akelleh/do_sampler_demo.ipynb

## do_sampler_demo.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import dowhy.api\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as pp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={'figure.figsize':(9, 7)})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, there is no causal relationship between title length and click-through rate. For each population, the click-through rate is independent of title length. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "n = 30000\n",
    "\n",
    "# normally, ctr is independent of title_length\n",
    "title_length = np.random.choice(range(25), size=N) + 1\n",
    "click_through_rate = np.random.beta(5, 100, size=N)\n",
    "\n",
    "# but one quirky person prefers titles in a much narrower range, and is especially talented\n",
    "title_length_2 = np.random.normal(13, 3 ,size=n).astype(int)\n",
    "click_through_rate_2 = np.random.beta(10, 100, size=n)\n",
    "\n",
    "# let's concatenate these together\n",
    "all_title_lengths = np.array(list(title_length) + list(title_length_2))\n",
    "all_click_rates = np.array(list(click_through_rate) + list(click_through_rate_2))\n",
    "\n",
    "df = pd.DataFrame({'click_through_rate': all_click_rates, \n",
    "                   'title_length': all_title_lengths, \n",
    "                   'author': [1]*N + [0]*n})\n",
    "\n",
    "# restrict down to where there's enough data\n",
    "df = df[df.title_length > 0][df.title_length < 25]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAE8AAAAUBAMAAAA+U32BAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAImYQu82Z3XZU70SJMqtHo0JKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABS0lEQVQoFY2SMUvDYBCGnyRNY2ii3RyNqwg6CJUiVsVZ6ugW/AV1km7iIE5Sf0EFRxHcHO2idFH6D9pZEASpKIh4l4SWrxbrO3x3997L3XfHgRUxEeeqWIHg+PG3tr3YGZDTRXF3YY+SegbchnOQEP4lOFWwY+wLAvFMLMC9MutPb/K2wSuS6+FrZGAHmsvK2JoqwyzM9Ch8GioJvuGuoWQidG84gUpM4WNE6LyL8GggDGK2YK5FOCq0vqA7FNpV6iKMCYU3YEnFbkuppLVf4+pfrcNXFcow/rhhmsNhRCitcxtpfaP3NiwN11Oo6TB/LjytIcM8SJ065Y7TNwriNZxrKlWY0oXnW9yKWTt9gcOkkUQpnPZZB69H/rm/D674xSzjm8KMzUxJPhhnfmCmzEiOgs2MWjVTRqRnpoebIMrsOJPT/1mRPBMwDz+CNk+5m8RT8QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( 0, \\quad 0.1\\right )$$"
      ],
      "text/plain": [
       "(0, 0.1)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGyCAYAAAAYveVYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmUnNV95//389S+dVfv3epWaxfaEGITQgbZmEUghMGAY5ZYxhPjQ44nJD6TOJkQjo/tOcEzIeNxnBkzmYyXEEjA83OQsYkgxOAFCYEEWCwSQrvU3ep9qb3qWX5/VHVJzSJaQi2VSp/XOTrdter2ud1Vn7r3e+81XNd1EREREakS5ulugIiIiMjJpHAjIiIiVUXhRkRERKqKwo2IiIhUFYUbERERqSoKNyIiIlJVFG5ERESkqijciIiISFWZ0nDz5JNPsmbNGq6++moeeeSRD7zfn/7pn/KTn/ykfLm7u5s777yTa6+9lt///d8nlUpNZTNFRESkikxZuOnt7eXb3/42jz76KOvXr+exxx5j165d77nPPffcw4YNGyZc//Wvf5077riDDRs2sGTJEv7X//pfU9VMERERqTJTFm42btzIihUriMfjhMNhVq9e/Z4Q8+STT3LllVdy3XXXla8rFAq8/PLLrF69GoCbb775PY8TERER+SDeqXrivr4+mpqaypebm5vZtm3bhPt88YtfBGDr1q3l64aHh4lGo3i9xaY1NTXR29s7Vc0UERGRKjNlIzfvdx6nYRhT9jgRERERmMKRm5aWFrZs2VK+3NfXR3Nz84c+rr6+nmQyiW3beDwe+vv7J/W4ow0Pp3Acl4aGKIODyeNuu0wN9UflUF9UDvVF5VBfVJbx/jBNg7q6yHE9dsrCzcqVK/nud7/L0NAQoVCIZ555hm9+85sf+jifz8dFF13EU089xQ033MATTzzBqlWrjuv/dhwXx3HL30vlUH9UDvVF5VBfVA71RWU50f6YsmmplpYWvvKVr7Bu3Tpuuukm1q5dy9KlS7n77rt5/fXXj/nYr33tazz++OOsWbOGLVu28Ed/9EdT1UwRERGpMob7fkUuZ7jBwSSO49LUFKO/P3G6myMl6o/Kob6oHOqLyqG+qCzj/WGaBg0N0eN6rHYoFhERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqorCjYiIiFQVhRsRERGpKgo3IiIiUlUUbkRERKSqKNyIiIhIVVG4ERERkaqicCMiIiJVReFGREREqsqUhpsnn3ySNWvWcPXVV/PII4+85/bt27dzyy23sHr1au677z4sywLg0KFD3Hnnndx444187nOfo6urayqbKSIiIlVkysJNb28v3/72t3n00UdZv349jz32GLt27Zpwnz/5kz/h/vvv5+mnn8Z1XR5//HEAvvOd73D99dezfv16rrnmGr797W9PVTNFRESkykxZuNm4cSMrVqwgHo8TDodZvXo1GzZsKN/e1dVFNptl2bJlANx8883l2x3HIZlMApDJZAgGg1PVTBEREaky3ql64r6+PpqamsqXm5ub2bZt2wfe3tTURG9vLwB/+Id/yG233cbDDz9MoVDgsccem6pmioiISJWZsnDjuu57rjMMY1K3/+mf/inf+MY3uOqqq3j66af5j//xP/LTn/50wuOPpaEhWv6+qSl2vE2XKaT+qBzqi8qhvqgc6ovKcqL9MWXhpqWlhS1btpQv9/X10dzcPOH2gYGB8uX+/n6am5sZGhpiz549XHXVVQCsXr2ar33tawwPD1NfXz+p/3twMInjuDQ1xejvT5ykn0g+KvVH5VBfVA71ReVQX1SW8f4wTWPCoMVkTFnNzcqVK9m0aRNDQ0NkMhmeeeYZVq1aVb69vb2dQCDA1q1bAXjiiSdYtWoVdXV1BAKBcjDaunUrkUhk0sFGREREzm5TOnLzla98hXXr1lEoFLj11ltZunQpd999N/feey/nnnsuDz74IH/xF39BKpVi0aJFrFu3DsMw+Nu//Vu++c1vks1miUQifPe7352qZoqIiEiVMdz3K345w2laqjKpPyqH+qJyqC8qh/qislTktJSIiIjI6aBwIyIiIlVF4UZERESqisKNiIiIVBWFGxEREakqCjciIiJSVRRuREREpKoo3IiIiEhVUbgRERGRqqJwIyIiIlVF4UZERESqisKNiIiIVBWFGxEREakqCjciIiJSVRRuREREpKoo3IiIiEhVUbgRERGRqqJwIyIiIlVF4UZERESqisKNiIiIVBWFGxEREakqCjciIiJSVRRuREREpKoo3IiIiEhVUbgRERGRqqJwIyIiIlVF4UZERESqisKNiIiIVBWFGxEREakqCjciIiJSVRRuREREpKoo3IiIiEhVUbgRERGRquI93Q0QEQFwXZds3iaVLZDL2xgGGIaBaRgYJpgYGGbpsgEGxds56vvx66H0vQEGBl6vgcfUZzmRs4XCjYicNpbtkM5ZjCXzjKVyWG5xONnjKUYU1wVwcce/d4uXDIwJtxnjCcd1i4mGd30PhANeaiN+wkEvQb8X0zQQkeqkcCMip8zRozOjyRyprFUeWTFNg57eJAd6Eziui2kUrxsfqRn//shX3nXZwBwf7TnqPn6vh7aGMLbj0DuUwcEFIBr0EYv4CAd8BP0ehR2RKqJwIyJTasLoTDqP5biYgN9nkCvY7Do0yq6uUQ70JrEdF49p4PEYOE4xDDmuWxqlOXEe06CjOcrsthiz2mpoawhTsG0ODxZwSyM80ZCPmoifcMBLwO/BNBR2RM5UCjciclIda3TGMKC7P8nurjF2dY0ymsoD0BQPsnxhM3M7aulsjuLxmO95Tsd1cRyKYceZeNkZv+xy5PvS13TWYv/hBHt7xnju1W6ee7WbgM/DjNYYs6cVw05DTYC8ZdMzkMJxXQzDIFYKOyGFHZEzjsKNiHxk46MziVSe0dSR0Rmf1yCbt9nVNcquQ6Mc6EviOC5+r8msaTVctrSNue011EYDE57PLQWMcYZh4DEMPCdYEzx/ehyAVLbAvp5i0Nnbk2DnwREAoiEfs6fVMKs0shMNeckVbLr6U7gGmAbEwj5qwqWw4/NMaJ+IVBaFGxE5Ydm8Rf9IhpFEcQTG6zUwKI7O7OoaZVfXGGOl0ZnmuhArFrUwt72G6e8anXFdl1zeJm855ZVPzng9cCnouMVq4glFwhMUK4snFBK7rosBmB6DgNdDJOhj8ax6Fs+qB2A4kSsHnV1do2zbPQhAQ02wPKozozVGwOchk7MYSxYACPg9NNYGqYn48Z5o4hKRKaNwIyLHLZMrhZpkDo9ZHBHZ3T3Grq4xDvYmcVyXgM/DrGkxVp3Xxtz2Wmoi/gnPYdkOubyN7RRHRmoiftqiAcIBDz6vByiGE7e0Qmq87ubIdYBbnJYqXl+8wjn6fqXrM7kCY6lCOWh5SmGnLhagLtbEBfObcF2XvuEMe0ph57Vdg7y8ox/DgLaGMLPaapjVVkNncxQX6OpP0tVvEI/5qa8JEg54NZojUiEUbkRk0tLZYqgZTeVwXdhxYJjNb/UykiyGhpa6EJcuaWFuey0dzZEJe8u4rkvecsjnHcDF5zVpqA0SC/sJfcDS7PG9a8Z3rzlRtRE/rfXFQJXN2yQz+fcNOy31YVrqw1y6uBXbdugaSLGnuxh2Nr3RywuvH8bvNVk8q57z5zcyrSFMMlNgOJHD7zVpjIeojfjL4UxETg+FGxH5UOlsgd7hDIl0Hsd1eWPPEJvf6iWVtehsiXL5edOY215DLDxxdMZxXHIFm4LlYBgG0aCPpqYQkaCPgP/UBwCvxyQaMomGfJMKO50tMTpbYnzifMgVbA70Jti+b5g39g7x6jsDNMdDLJvXyNI59Xi9Jj0DaXoGUtREAjTUBgkHvSpEFjkNFG5E5AOlsgV6h9IkMxa2bfPb3YO8tL2PbN5mzrQaLjuvjRktsQmPKZRGZ2xcvKZBPOovjs4EvBVXn/JhYSeRzuO6R8LOvI448zrirF7eyZv7hnh15wDPvHyQf996iAWdcc6f38TM1iiZXIE93Xm8HoPGeJDaSICAT6M5IqeKwo2ITOC6LqmsRe9QmlS2QN6yefWdAbbu6CdvOZzTGeeypW20N0bK988VbPKFYjFwMOClpT5EJFTcHO9MqkOZbNjxeU3On9fIBfOb6B1O89rOAbbtGeTNfcPEo36WzWtk2dxGAj4PfUMZDg+miYX8NMSDRIM+bRgoMsUUbkQEOBJqDg+mSGUt8gWbLW/38eo7A9iOy+JZ9Vx2bhvNdaHy/dNZG9txqI0EaGuYWAxcDd4v7GRyFsOJHKPJPOBSHwuy+pJOrrywg7cPjvDKzn6ef7WbX77WzZxpNZw/v4l57TXkLZt9PWN4TIOG2iDxaICgXy/BIlNBf1kiZznXdUlkitNP6ZxFOlvg5R39bNs1CAacN6eBj53bSn1NECjW0aSzFq4LjfEg9TXBs2bKxesxiYWL02yFBodEOs/AaJaxVAGPBxbOqGPxrHqGEzlee2eA13YN8OPndhMJelk6p4Hz5zVSXxNkcDRH33CGSNCHJ+AtHzchIieHwo3IWcpxXRKpPL3DGbI5i0Qmz+a3+nhz3xAe0+CiBU1cuqSV2tISbtt2SGctDMOguS5EXSyIz1tZNTSnks9rUl8TpC4WIJu3GU3mGBjN4rrFQzo/cf40Pr5sGru7R3l15wCb3+pj05u9TG+OcsH8RhbOqMN2HHZ3jWFlC7Q3RQgF9JIscjIYrvtRT22pPIODxV1Qm5pi9PcnTndzpET9URkc18UX9PHWOwPkChbDiSwvvtnH2wdH8HtNLlrQzIrFLURDPgAsyyGds/Ga0FwfJh4NVFxhcKWwHYdU1mJwJEsik8c0DIIBD16PSTJTYNvuQV7d2c/gWA6/z+TcWQ18cnknhlOsWWquC9EUD01YQi+njl6jKst4f5imQUND9Lgeq48JImeJ8emnnoE0wbCfQ/0JNr3Zy57uMYJ+Dx9fNo3lC5vLowf5gk02b+PzmkxvjlAT8etN90N4TJOasJ+asJ98wWYslad/NEM6a+Hzmly6uIVLF7dwoC/JqzsH+O3uQV55p59LFrXw8fPa6B/OMJLM09EULYdLETl+CjciZ4FkplAuFO4dTrF5ez97u8eIBL1ceWE7Fy1oLtfN5EqhJuT3MqMlRiziVz3ICfD7PDTGQzTUBknnLIbHcgwnc7iuS1t9mBmXz2L18ulsfLOPF7Z1s33fMGsuncGMlii7u0ZoqA3SUhc5q6f+RE6Uwo1IFUtnLQ4PpYq76CZz/Oq1bvb2JIhHA1x7SSfnz2ssv3lmcxZ5yyES9DFnWpRIUMcJnAyGYRAJ+ogEfbQ2hEmWi5DzeEyTT398DvPaY/xs437+6dl3WDyrnmsu7mA0WWA0OUJ7Y5jaaEB9IXIcFG5EqlA2b9E7nGE0kSOdK7DxzV7e2DNEKOBl9fLpXHXJTBKJTGk5t4VlO9RE/HS2hAgHNR0yVbwek3gsSDwWJJu3GEnmSOct6mNB7r5hIZve7OXXv+1hd9coV1/UwZLZ9RzoTRJL5pjWED0tuzqLnIkUbkSqSK5gMzCSYXAsS8Fy2Pp2X+nwR4PLlrayckkrwdI5TumMheU41MUCNNaGtFLnFAv6vbTWe6mrj/DG270MjGa4dHELi2bW87ON+3hy43627R5k7cqZZPM2Ow8N01Yfob42qGlCkQ+hVzORKlCwHAZHM/SPZnBseGPfIL/Zdpi8ZbNsbiMfXzaNmoi/uFFfxsL05YnHiucfnS171FQqr8ekvSlKbTTAwd4EQZ+Hdavn89quQZ7dcoiH1r/J5ee1cemiFroH0wwnc7Q3RgkH9fIt8kH01yFyBrNsh6FElr6hDI7jsqdnjOdf7WIsXWBeRy1XXthR3lE4k7MoWA6NtSEWzmlkdCR9mlsvR4uGfMybHqdvuBhSl8yqZ35HnKdfOsDzr3bz5t4h1q6cScAf4p2uEZrjxWXjWpYv8l4KNyJnINtxGEnkOTycxrYdeofS/PsrXfQNZ5jWGOamy2cxs60GKK1+ytnEowGaW4vTT36N1lQkr8dkWmNx2f3BvuJ+Kzd/fDbnzhnlX188wA+e2sFF5zRxxQXtDIxkGUnk6GiOvuc0dpGzncKNyBnEcV3Gkjm6B9NYtsNYOs9zr3SxtydBXSzALR+fzaKZdRiGQcEqnoMUDniZ21FLRIXCZ4xoyMe8jjj9Ixn6hzPMbI3x+zct5rlXu3hpe3HDxesu6WROey17e8aIRwO0NYSr6lwvkY9C4UbkDOC6LmOpPD1DaQoFh1zB4lfbenhjzxDh0gqoi85pwuMxi7vkpgv4vCYzWmLURPxaRnwG8npM2hrGR3GSFCyHay6ezpLZDfzshX08/txuzumMc+3y6STTBXamR2lrDFOnZeMiCjcila5g2ew/nCSdK+Di8uKbfWzZ0feeFVCO45JMFzAMg/bGKHWxAKapN7kzXSToY15HLf3DGXpHMjTVBvniDQvZ/FYfz7/azffWv8mVF3Zw/twGDvUlGR7LMb05qqlHOasp3IhUsFzeZu/hMfIFm9f3fPAKqHTWwnYcmuIhGmtVZFptPKZJa3kUJ0U6Y3Hp4hYWdMZ56sUD/OuLB3i9tGzc6zF559Bo6ciMwOluushpoXAjUqHSWYu9PWPs7h7hF1vffwXU+K7CdbEAzXVhLeuucuGgj7kdNfSPZOgbzhAJ+bjz6nm8vmeIZ146yN/99K3ysvG9PWO01IVprgtrBE/OOpMON4cPH+btt9/msssuo6+vj7a2tqlsl8hZLZkpsLdnjNd3D/Ds1i7aGiaugMoXbDI5m1jYx4zWGm3AdxbxmCat9RFqwgEO9iVJpAssmVXP3PYaNrx0kF++1s3+3gSfvnwWfSMZUlmLjuaogq+cVSY1dv38889z22238fWvf53BwUHWrFnDs88+O9VtEzkrjSZz7O4a4aXtvTy7tYsFnXG+cN0CZrbVFFdIpfKAwZz2Wma1KdicrcJBL/M6ammuC5FMF/CYJjevms2nPjaTQ30p/s+T2xkcy5K3bHYdGmUslTvdTRY5ZSYVbv72b/+Wxx9/nJqaGpqbm3n00Uf5m7/5m6lum8hZZ3Asy96eMX69rYcXXj/M+fMaufUTczBNg0QqT77g0NkcZW5HLdGQT6tiznKmadBaH2He9DimaTCWyrN0TgNfXLuQoN/DPz6zk5d39OHzGuztSdAzmMJx3NPdbJEpN6lw4zgOzc3N5csLFy7Ui6rISeS6Ln3DafYfHuPZLYd4ZecAK5e0snblDAq2Qypj0dYQ4ZzOOPGYzhaSiUIBL3Om1dJaHyaZKRCP+vni2oUsmVXP86928/hzu/GY0D+SYU/3GLmCfbqbLDKlJhVuQqEQ3d3d5UCzZcsWAgFV4YucDI7r0jOY5kBvgqdePMBb+4e56qIOrrqog1zexrZd5nbU0hgP4TG1Ckren2kaNNeFmdsex3ZcCgWXmy6fxdqVM9h/OMHfPbmdoUSOgm3zzsERRpOappLqNalXyv/0n/4T/+E//AcOHDjAZz/7Wb785S/zx3/8xx/6uCeffJI1a9Zw9dVX88gjj7zn9u3bt3PLLbewevVq7rvvPizLAqCvr48vfelL3HTTTdx2220cOnToOH8skTOD47h09SU50Jtg/W/2sbdnjE99bCYrl7SSzlpgGMyZVqu6Gpm0Yi1OnGjYx1i6wLK5jfze2oX4fSYPP/02W97uJ+A32Xc4QfeApqmkOhmu607qN3tsbIxXX30Vx3E477zzqK+vP+b9e3t7uf322/nJT36C3+/ntttu47//9//O3Llzy/dZu3Yt/+W//BeWLVvGn//5n7NkyRLuuOMO7rrrLlavXs3tt9/OP/3TP7F582b+x//4H5P+oQYHkziOS1NTjP7+xKQfJ1NL/TGRZTsc6kvSPZjiiV/vZTiR45aPz2bBjDpSGQu/z2Rma2xKttRXX1SOqeoL13UZHMvSPZAi6PfgAj/buJ839w4xZ1oNN14+E9eBYMBLZ3OMgF+rqfR3UVnG+8M0DRoaosf12EmN3Hzxi1+kpqaGj3/841xxxRXU19fzO7/zO8d8zMaNG1mxYgXxeJxwOMzq1avZsGFD+fauri6y2SzLli0D4Oabb2bDhg0MDQ2xY8cObrvtNgBuueUW/uiP/ui4fiiRSlewHPb1JDjYl+THz+1mLJ3nzqvnsWBGHcl0gaDfw6y2Gp0VJCfMMAwaa0PlaSqr4HLzqlmsWdHJvsMJ/s+T2xlO5rBth3cOjTCiaSqpIscc67733nvZu3cvBw8e5IYbbihfb1kW5ofM/ff19dHU1FS+3NzczLZt2z7w9qamJnp7ezl48CDTpk3jL//yL9m8eTPTpk3j/vvvP+4fTKRS5Qs2e3sSdA8k+f9+tQeAz197Dq31YcZSeWrCfqa3RFVfIyfF+DRVV3+KkVSOC+Y30d4U5f89v5sfbXibT17QziULm9l/OEGqtkBrQ1i/e3LGO2a4+epXv0pXVxf333//hIDh8XiYN2/eMZ/4/Wa7jl5h9UG3W5bFW2+9xR/8wR9w33338eMf/5g/+7M/4+GHH/7QH2bc0cNXTU2xST9Opt7Z3h+ZnMXOA8MMJHI8/vxuwgEv99y8lKZ4iJFkjjmdMaa31uA5BTvKnu19UUlORV+0ttTQN5zhYG+COdNDfPV36/jnZ3fy71u76BnKcPvV5+C4LkMpi1ntZ2+dl/4uKsuJ9scxf3s7Ojro6Ohgw4YN7xmpSafTx3zilpYWtmzZUr7c19c3YTl5S0sLAwMD5cv9/f00NzfT1NREJBLhiiuuAI7U5RwP1dxUprO9P9LZAnt7EuzpHuWnL+yjvibAnVfPx4vD/q4RmuIhQl6DocHklLflbO+LSnIq+8IEmqJ+9veO4djwqZWdTKsP8czLB/mrR7Zyy8dn0xQPcrgvQUdzhLpY8JS0q1Lo76KyTHnNzS9+8Qs+9alPcdVVV3HllVdyxRVX8LGPfeyYj1m5ciWbNm1iaGiITCbDM888w6pVq8q3t7e3EwgE2Lp1KwBPPPEEq1atorOzk5aWFn75y18C8Nxzz7F48eLj+qFEKk0inWd39yg7DgzxxG/20toQ5q7rFhAN+RhLFWhtCNHWENb+NTLl3r2a6sJzmvjCmgWYBvzoX9/m1Z0DhAIm+3uTHOpPYDvO6W6yyHGbVLj5b//tv3HPPffQ1tbG1772NS6//PJywe8HaWlp4Stf+Qrr1q3jpptuYu3atSxdupS7776b119/HYAHH3yQBx54gOuuu45MJsO6deuA4o7If//3f8/atWv5h3/4B/7yL//yI/6YIqfPcCLL3u4xfvvOID/fdIDZbTV87pr5+H0miXSe9sYILXURbYwpp4zXY9LZEqW9MUIyU6CxNsiXbljE/M5a/m3LIf7fL/fg9xoMj+XZdWisuC2ByBlkUkvBb7zxRtavX89f//Vfc9lll3HxxRdzyy238C//8i+noo3HTdNSlels7I+B0QyH+pO8vL2PTW/2snhmHTddPgsXSGUtpjdHqT8NQ/9nY19UqtPdF+msVZ6mCoc8vLS9j3/bcohYyMctn5hNUzxELmfT2hCmMR6q6tHF090XMtGUT0v5/X7y+TydnZ1s374d0zTJ5/Mn1FiRs4HruhweSnGwN8mvXutm05u9XLSgiU+vmo3jFt9QZrXGTkuwETnau6epLl7QzBeuOwfDgB/+69tsfbufcNBLz2CaPV2jZPMaxZHKN6lwc+WVV/KlL32Jyy+/nB/+8If8wR/8AfF4fKrbJnJGsmyHroEUXQMpntlykNd2DbLqvDauu6QTy3bI5W1mT6ulJqIjTKQyvHuaqike4u4bFjF/ei3PbjnEP/7bThzXxbId3jk4wsBoBmdy+7+KnBaTmpYaGBggn88zbdo03nrrLV5++WXWrl1LQ0PDqWjjcdO0VGU6G/pjLJXjUH+KTM7i55v2s+9wgmsvmc7yhS3k8jaW4zK7rea0L7M9G/riTFFpffHuaaptuwfZsPkgrutyzfLpnDengVTGIhr20d4UJeCrno0mK60vznYfZVpqUq+wv/u7v1veXXjRokUsWrTo+FspUsXyBZueoTQjiRwuLj/51R56hzJ8+vJZnDungUzOAhfmTKsh6D879w+RM8O7N/07d3YDM1tjrH9hHz/buJ+3D4xww8qZ5ArFAzinNUWoiwZUEC8VZVLTUu3t7bzyyis4WhIoMoHjugyNZdl5cIRkOs+B3gR//7Pt9I9k+eyVczh3TgPprIVpmMxpr1WwkTPCu6epgn4vn7tmPquXT2dvzxjfW/8m+3oShAIeDvUmOdCbpGDZp7vZImWTeqXdvXs3d9xxB16vF7/fj+u6GIbBK6+8MtXtE6lYmZxFV3+KdLZA3rbZsPkge7rHaG+MsPZjM2ipC5fOifIyozWGz6st7eXMMX42VSTo41B/kkS6wPKFzcyeVsP6X+/lx8/v5tzZ9Vx3SSepbIGdB0fpaIpQG1UtmZx+kwo3jzzyyAfe9tJLL7F8+fKT1iCRSmc7DgMjWXqH03hMg9f3DvLL13owTbjukk4uPKcJ0zSHtRiPAAAgAElEQVRIpPJEwz6mN8fwehRs5MwUCniZM62WvpE0vcMZaiN+vnD9An6z7TC/+m03+w8n+NRls+hsjrLvcIKGmjytDRH9zstpNalw097e/oG3PfDAAxW7343IyZbMFDjUn6RQsBlN5nnqxf30DmdY0Bnn2ks6qYn4cRyXsVSB2miA6U1RzFNwTpTIVDJNg9b6CDVhPwf7UqSzNqvOa2NeRy1P/Hov//jMTpYvbOaTF7Qzmiwwlh6hsyVGNOQ73U2Xs9RHLgCYxGIrkTNewXLoHU4xOJbDNOA3rx/m5R19xMI+PvvJOZzTWQdANmeRtxxa66t/wzM5+4SDPuZ2FA/g7BvO0BQPcvcNi/jFK4fY/FYfu7tGuenyWTTXhdjdPUpTbYiW+pBOGZdT7iOHG1XISzVzXZfRVJ6ugRQ40D2QYsPmA+X6gysuaCfg82A7DqmMRSToY2abVkRJ9fKYJm0NEWJhPwd7i2dPXXPxdOZPj7P+N/v4/lM7uGxpG5cvbWNoLMtYOk9nc5RwUKM4curoFVjkA+QKNt0DKRKpPJbj8G9bDvH2gRFa6kL8zhVzaG8q7ruQzlo4jktHU5S6mJbEytkhGvIxb3qcw0NpBkezTG+Kcs+Ni9iw+SC//m0Puw4VR3FqfX7e6RqlpS5EczysaVo5JRRuRN5lfHl3z2AKA4MdB0f4xSuHcBy48sJ2VixuwWOaFCyHdNYiHgswrSGMz1s9m5mJTIbXY9LRFKUm7ONQfwrXhRsvm8k5nXF+vnE/f/fTt7jywnaWL2ymbzhDIl2goyl62jexlOqnmhuRo6SzFl0DSTI5m1S2wFObDtA1kGL2tBquv3QGdbEAruuSTBcwTYNZbTEdoyBnvZpIgHkBHz2DKYYSOeZ11HLPTYv52cZ9PPNyccTzxstm4Tgu7xwaKdak1YY0iiNT5iOHm0suueRktEPktLJsh/6RDH0jGTyGwUvbe9n05mFCAS+fXjWLJbPqMQyDXMEmm7dpjodoioe03FWkxOc16WyJURPx09WfwjQMPvvJufx21yAbXjrAQz99k2uXd3Lu7HoOD6UZTuTpaI4QUS2OTIFJhZs9e/bw/e9/n8HBwQkjNQ899BD/+T//5ylrnMipkMwUONibwLJd+oczPPXifkaSeZbNa+TqizoIBbw4jksyUyDg8zCvvVbFkSIfIB4NEAl66epPMZbKc+6c+uLxDb/Zy09f2Mebe4dYc+kMDMNl96FRGuMhmuv0QUFOrkmFmz/+4z/mwgsv5Oqrr1axpFQNy3boHUozMJrBxeW5V7p5fc8QDTUB1q2ez8y2GqC4E7FlO7TVR6ivCWooXeRD+LweZrTGGE7m6O5PEfB7WHftOby0vY/nXunie0+8ySfOn8YlC5sZGssykszR0RTRFK+cNJMKN4VCgfvuu2+q2yJyyiTSeQ72JbFthz3dYzy79RC5gsOq89q47Nw2vF4T23ZIZS2iIR+z22oJ+FUwLDJZhmFQHwsSCfro6k+SSBW4eEEzC2bU8a8vHuDZLYd4Y88g16+cSXM8xN6eBHWxPG0qzpeTYFLhZtq0aRw8eJDp06dPdXtEppRlOxweSjE4miNXsHn6pQPs7UkwvTnK2pUzaIqHcF2XVMYCoLM5Sq1OPBY5YQGfh1ltNQyN5egeTOL3mnz2k3PYcWCEDZsP8P2fb+fiBc184vxpJNMFdqZHmdYYJq6/O/kIjhlu7rnnHgD6+/u59dZbOffcc/F6jzzkoYcemtrWiZxEY6kch/pT2JbLG6XzoDymwZoVxfOgDMOgYNmkszb1sQCtDREddilyEhiGQUNtkEioVIuTLjC/o5ZZbYv5xStdvLS9jx37h7luxQzmttdwsC/JcCJHe2NUI6ZyQo4ZblavXn2q2iEyZQqWw+HSEtVEKs9Tmw9weCjNOdPjXLeieB6U67ok0gW8HoM57bU6E0dkCgT9XmZPq2E4maNnIA0UD5s9d3YDP9+4n8d+sYuFM+q49pLp5Ao2Ow8N09ZQqnXTKI4ch2OGm09/+tMAdHd3T7jeMAyCweDUtUrkJBlN5jg0kKJQcHh5Ry8vvtVLJOjjM5+Yw4IZcQzDIJu3yOcdmutCNNXpHByRqTReixML+cvntbXUhbj7UwvZ9EYvv/ptN3u6x/jkhe1cMK+R7oE0I4kc7dr8T47DpH5Tbr/9dvr6+ohGoxiGQSKRwOPxUFdXx3e+8x0uuOCCqW6nyHEpWDbdg8UXxb7hNBteOshwIscF8xu58sLi8m7bdkhlCoSDPmZMr9ELp8gp5POadDTFiEeDdPWnSGUsVi5pZdHMOn6+aT//+uIBtu0eZO3KGQR8Ju8cHKG5vri/lD6AyIeZ1Kv5ypUrueSSS7jpppsAePrpp3nhhRe47bbb+NrXvsaPf/zjKW2kyGSVD7rsT5HNWfzm9R5+u3uQ+qOWd7uuSzpj4QIdzVHisYCGvEVOk2jIx7yOWgbGMvQOpQkFvPzuNfN5fc8Qz7x0kP/z0+1cuqSFy5e20T+cYSSZp6MpqqljOaZJhZsdO3bwwAMPlC+vXr2a//2//zeLFi2iUChMWeNEjkfBsukeSDOczLL/cIJntxwik7P52LmtrDpvGj6vSb5gk8mNFwxryalIJTBNg+Z4mNpwgO7B4rLxRTPrmNtey79tOcgLrx/mrX3DXH/pDDqaIuzuHqWxJkhLfVib/8n7mlS4sSyLnTt3Mn/+fAB27tyJ4zjkcjksy5rSBop8GNd1GU3m6BpIk0jnee7VLt45NMq0hjB3XjOT1vowjuMylsrj93lUMCxSoQJ+DzNbaxhN5ekeSOE4Lp/62EyWzmng55v284/P7GTpnAauuqid4USO0VRxX5xY2K+QIxNMeofiz33uc8ybNw/Hcdi/fz8PPvggf/M3f8NVV1011W0U+UD5gk33QIqRZI4dB4Z5/tVuXOCaiztYvrAF0zSKOwxbDq31YRpqg5qvF6lghmEQjwaIhnz0DWcYGMkwrTHCPZ9azK+39fDCG4d559Ao11zcwaIZdRzsS2IaBnWxAHWxAKGAV/vjCIY7yWO9R0ZG2LJlCx6PhwsuuIDa2lqSySTRaHSq23jcBgeTOI5LU1OM/v7E6W6OlJzM/nBdt7y1+1Aiy79tOURX/8TTuy3bIZ2xiIZ9TGuMEPSrYHic/jYqh/ri2NLZAl0DKTI5m2jIy8Bolp9v2s/BviQzW2N8bGkrs1pjZPMOlu0Q8HlpjAepCfuPe58q9UVlGe8P0zRoaDi+rDGpV/tnnnmm/L3jOGzevBmAa6655rj+M5ET5bgu+YJNLm8zli6QyOTJ5Cx+u2uAjW/0EvB5uOnyWZw7ux6AZLqAYRp0tsaojfj1SU7kDBUO+pjTXsvQWJaegTTRkJe7rjuHV3YO8ItXunjkmXeIR/0sm9fIsrmNeEzoHkjR7aaojfqprwkSCWo052wzqXDz8MMPl78vFAq8/fbbLF++XOFGpoztOOTyNumcRSJdIJUtMD7G6PMY9A2neWrTAQZGs5w7u55rlk8nEvSRK9jkcjYNtSo2FKkWpmHQWBsiFvZzeDDNSDLH0jkNnDengR0HRnjlnX6ef7WbX77Wzdz2Wi6Y38jc9lpS2QIjyTw+r0FTbYiaiB+/T4sIzgbHHW4ADhw4wF//9V9PSYPk7FSwbLJ5m1S2QCJdIJu3cV0wDfB6DVKZ4tD0wd4kB/uSDCVy1Eb83HHVPOZ21JYKhguEAh7mdtQSDqpgWKTaBHweOlui1MX8xaNUbJfFs+pYMrueobEsr+0a4LV3BnnsF7uJhnycN7eB8+c1Eg76OTyUpmcwRSzip7E2RDjo1RYQVeyEihA6OzvZs2fPyW6LnCVc1y2OsBQckuk8Y+k8BcvBMAxMEwwMBkYzHOpLcrA/xaG+JNm8DUA44GV6c5QLFzRx4fwm/D4P6WwB24ZpjWFt0y5S5QzDoCYSYH7Qx8BIhoHRLI7rEg35+OQFHXxiWTvvHBrllZ39bHzjMC+8fpiZbTEumNfIOZ1xsnmLPd2jeD0mjfEgtZEAAY3mVJ3jrrlxXZc33nhjwgGaIsfiOC7ZvE3/cJp9PWOksgUcF3BdvF6TbN6iZzBdHpU5PJTBKc1BNcWDLJxRx/TmKNNbotTHjpwUXLAcRpN54tEAbQ1hDTeLnEW8HpPWhgiN8RCJdJ6B0Sxj6Txew2ReRy3ndMYZS+VLozkD/ORXewkFPCyd3cD585toqPXSN5Th8FCGWNBHQzxIvTOp9TVyBjjuaSnDMKivr+db3/rWlDVKzmyW7ZAr2GSyFmOlehlcl9p4mEzOKu8gfLA/ycHeJKOpPFB8sWpvDLPy3BY6mqJMb37vWTKW5ZC3LGwbvB6DWW0xaiKB0/FjikgF8HpM6mJB6mLB0utLjqHRLJbj4veZXL60jcuXtrG3J8ErO/t5+e1+Nm/vo70pwgXzGlk8q568ZbOvZ4yxnI2dLxAKeAn6vfg8Jj6vqdq9M9Ckl4JDcTM/13Xx+Sq7nkFLwU8ty3bI5m3S2QJj6TyZnA0uGAa4uPQNZ4ojMsMZ9nWPkbccAGJhH9Obo8Ug0xKltX7imTGu61KwHPKWUywmdoubfNVEfESCPsJBr/asOUH626gc6ouTz3FcUtkCA6NZkuk8YBAKevB6TNLZAr/dPcirOwcYGM3i95osmV3P+fObWDCrkYHBJJbtYDsuhmHguuDzGoT8XsIhL0GfF59XoedUmPKl4IODg/zZn/0ZmzZtwrZtLr74Yv7qr/6KlpaWE2qwnNnKxb+ZAmPpAtm8hUFxC/VcwS5OMfWNTzGlcUtBp60hwtK5DcUppqYotdGJS7QdxyWXtylYDm7pQeGAl5a6IOGgl4DPoxcTEflQpmkQC/uJhf3kCzZj6TwDIxnSWQuvx2TFohZWLGrhUH+KV3b2s233EK/sHKC5LkRDTZC6WID6mkDpaxCf1yRn2aSGC8XQA4CBx2MQCngJB7wEA178Cj0VY1IjN3/4h3/IvHnzWLduHbZt8/DDD7N9+3a+973vnYo2HjeN3Jw8ruuSt4rLspOZPGOpPHnLwcDANIsrlLpLYeZAb4KRZHGKyec1aW+MML0lSmdpdKa1uYbhkVT5uW3HIV9wsCwHDAPTKB6iFw37Cfm9BP0eTFPFwVNBfxuVQ31xariuSyZnMZzIMZTIARDwmfh9HrJ5izf3DrOnJ0HfUIrhZB7nqPobj3lkB+Ri6AlSHwsQj/qJhHw4rotjl+5vgMc0CQW8hAIegn4vHtPA6zHxegw8HlOLHiZpykdu9u3bx3e+853y5XvvvZfrr7/++FopFct1XWzHLQ/FWrZLwSoeMJlI57EdF9d1MQwYHM3SNZDiQKn4d3wVUyToZXpLlOULm5neEnvPFBMUC4DT2QKWDRjFmpmakJ9Y2EfA7yHg82ijLRGZEoZhEA76CAd9tNSHSabzDI5lGUvl8Zgmy+Y1cNUlMxkeSZXPohtK5IphaCxbDkX7DicolKbWi88L8Wgp+JQCUDwaoCbiIxry4THN0h5d43P1Lh6PWQxWXg9+n4eAvzgqXQxBCkAnw6QPzszlcgQCxcLNTCajN6EzhOO62LaDZRcDjG075CybfKE4apK3itNA5bmj0kCeYRoUCjaHh9N09RdHZroHUtilTzONtUdWMRX3nQhM+J1wSzsK5wsO43/XoUiQumiQSNhH0OfR6iYROS28HpN4LEg8FiSbtxhN5hkczTCSzJHKFPB4DKIhH7VRP7On1Ux4rOu6JDOFYvAZy5UD0PBYljcGUuUPfOMiQS/xaIB4aaQnHg1QGyl+qIuGvGAYOOWHvE8A8hWnu44OQKZp4Cn903vx+5tUuFmzZg133XUXN998MwA/+clPWL169ZQ2TCbHsp3Sv1KYsGwKhSMBxradI6GlWOHL+IBKOmeRzhVIpi2SmeLmeYl0nrF0gbFUnuHS0K1pGkxrCLN8YTOdLVE6mqNE3rVJnu045PM2ln1kKDca8lEfCxIq1cu0tdZq+F1EKkrQ7yVY76WpLkQoEqSr2yCTs8nkLbKZ4iIaAJfia6HXNAj6vXQ2+5jREnvP82VyVin4FEd7RpJ5RpI5ugdSbN83XN7mYlxN2EdtdHzEpxh+4lE/NRE/HtNLwc4z5rhHApAx3priBY9p4PeZeE2z+NVbHBEab2sxCB0JRWeLSYWbL3/5y7S2tvLrX/8ax3G4+eabufXWW6e6bRWnYNk4DuXUfCp+UcanjAqWQ8F2yucrZfMW2YJTnOctNcMADMOlYLmkshbpbIFk1iKRzpNIFVcyjQeYVNZ6z//lMQ1iYR+xsJ+2hjDL5jXS2RxlWmNkwgF0xVVMNnnLwXGKl31ek9hRU0x+n0fDqiJyxjANg5qIn1w8VL7OdYvT9JZ95PU3k7PJFSzSWasYVFzAMDAAj6c4rdRWH6a9MfKe/8NxXBLpPMPJPKPJ4ojPaDLPcLI43TVW2hZjYpt85ZGfWNhHOOAtLVX3EA54Cfg9WLYHn9ckky/+H+PNmsAtrv7ylxZmFFd8Gfi9niP1QKaJx1MdI0KTCjef//zn+dGPfsQtt9wy1e2paIf6UyTSeUzDKP4+A16vic9j4PUWk3PxcrFa/ugQNP71/d7wi39ADgWr+DVfsMnki39A2Zxd/uQwPgKTHS/uTRcYLY2wjCSLxb6JdH7C6Mm4cMBbDC4RP9Maw0RDPmpKQ6PFVQXFP5r3+4UubsJnYVnuUauYPDTWhIiEfAR8Jj6vpphEpLoYhoHPa+DzmoTeddvRHzzH9/Ya/+CZydvYtnPk9bQ0Yu4xi6uromEfM1qi73m9tW2H0VS+PNpT/JdnJJFj16FRkpnCMdsb9HvKhcyhUggK+Y++XKxtDJTqfEJ+Dx6PWRrYH29r8TXe6zHwl0aBfKX6IK/nSPgZD0KV+iF2UuEmkUiQTqcJh8NT3Z6KNZbKs/Xtfrye4im1wVIBrN9nYnk9eO3iXiyO6xar7Etpvmy8lsUwyoHIY5rFlUgF+6i7ueQtm2S6OE00ls6Xk/1IovgLb79rF81YuJjs2xsjxMLxCYGlJuwnGvZ96NJE13WxbRfbdXAct1REDOBimibRkJdY3E/Q7yXgN7W/jIic1YxSABh/bX33BNXRJQPjtY6FgkOu4JArWMUPoRNqHYvHz4QCXqIhHzNao5iG8d7tMgo2mZxFOmcVp89y1lH/itNpmWzx8tBYjkzOek8d0NFCAQ+1ET81kWItUG3Uf1RNUPFDr5MZHxFyj7yvlb4fnxbzeU1qI8Vi6kowqXATCoW44oorOOeccyYEnIceemjKGlZptu0e5KkX93/g7YZRPNQt6C/+C/i9BCdc9kwIRH6vh4LtMJYqjriMJ/XhRH5C2IFiGq+LBWiuCzF/em2pMK00Rxvx4/V+cNBwS2ErX7CxneL3jkNx2OmoX1TTNPB5zGKiH0/qpYTv95pn/BCliMip5PUce7+b8Skv2zkSgPKlWslcebGH/Z4AZBjFOqHiJqbFkRPzQ6aRxkffMzm7FIqK/xKl+sriaFGO/YcT73n/MU2jWBcU8VMbDVATKYaf8SAUC/lwXUikC7guZ1a4ORvra97tsqVt1ET8pLN5rNJZScUhSJtcoTgUOX55/N9wMle+37t/YY7m9RjlOdXpzdGjKusD1EX9BAPv303joy25vI3tlOpfoDRq5IJb/EPweYtBJeI1y8HKc/SyQ9M8qwrNREROt/KUF8cOQOPbdIyveC2UVrhalkPBLo70Wzn7vTU2jNdhlupEPSa1keIH5WO93mfzFmOpAqOpYj3QaOkD+Ggqz/7DCcbS+fGJiLKgvzj6c8PHZjKzteb9n/gUm1S4+fSnP41t2ySTSY7jtIaqU9yzoBgWjNKmc4YxucLi8dGTbN4mW/rqMQ3iUT/RkO8Da11sxyGbs4rLuMd3xiytejIM8Ps8BAMeAj5/cW7Ua04oCtNOmSIiZ6Z3T30dy/j7xfgIvVX6Or6CtmA5WE7xqJxiaUPpvdw1cHExDaP8vtFQE6ApHvzA96VEpsBoMlcOPaOp4mKVSvqQPKlw8+ijj/Ktb32LQqFYzOSWqq63b98+pY2rNHUxPz6vWdw3ximm5uKISamehvGRkyPLrosppBSCDKNY8BX0FoMRBrbrkskVR17Ky55KzzG+z0GotC+Mz+cp/qKXwouCi4iIQHF0xjQ9TObkx/ERofG9zyx7PATZpZogm0LOxh1fOnNUTZDHLC4oiYWimC1H3rNyBZtABe1dNqlw83//7//lscceY+HChVPdnopWPHn2vde7rlsqJD5SUFy8jnKKHg9Etu1SsIuByLJdAl4PgVBph8rxFVel4FJJKVhERKrDkREh4AMCyXumxCZsAFsstSivzC29/wX9Z1i4qa2tPeuDzbEYhoHHMNBAioiIVIPJTIkdvQO+ZTsVNZtwzHAzMjICwLJly/jhD3/I2rVr8XqPPCQej09t60RERKQimYaB6fXgm9Qwyal1zCatWLECwzDKRcTf+ta3yredjTU3IiIiUvmOGW527NgBgOM4mO/atG18VEdERESkkkxqguz9jl244447TnpjRERERD6qY47cfP7zn+f1118nm81ywQUXlK93HEcFxiIiIlKRjhlu/uf//J+MjIzw53/+5zzwwANHHuT10tTUBHDWnzklIiIileWY4SYajRKNRvmHf/iHD7zPnXfeyb/8y7+c9IaJiIiInIiPvCj9bD6OQURERCrPRw43Oi1aREREKknlbCcoIiIichIo3IiIiEhVUbgRERGRqnLC4cayLAAtAxcREZGKMqlw88gjj0y4vGvXLj7zmc+8720iIiIip9Okws0///M/8/Of/xyAH/zgB9xxxx3ceOONU9owERERkRMxqYPKv//97/OFL3yBH/3oR5imyWOPPcasWbOmum0iIiIix+2YIzcjIyOMjIzg8/l48MEHOXz4MPfccw91dXU6FVxEREQq0jFHblasWIFhGBN2Ib7nnnuA4uZ927dvn9rWiYiIiBynY4abHTt2lL93XRfDMLBtG8dx8Pl8U944ERERkeM1qYLizZs3lwuI9+zZwyc+8QleffXVKW2YiIiIyImYVLj5r//1v/LAAw8AMG/ePP7u7/6ufPlYnnzySdasWcPVV1/9vkvGt2/fzi233MLq1au57777ynvnjHvrrbdYsmTJZJooIiIiAkwy3BQKBRYvXly+vHjxYvL5/DEf09vby7e//W0effRR1q9fz2OPPcauXbsm3OdP/uRPuP/++3n66adxXZfHH3+8fFsmk+Eb3/gGhULheH4eEREROctNKtyEQiF+9atflS9v2rTpQ3cm3rhxIytWrCAejxMOh1m9ejUbNmwo397V1UU2m2XZsmUA3HzzzRNu/9a3vsVdd911PD+LiIiIyOT2ubnvvvv48pe/jNdbvLtpmnz3u9895mP6+vpoamoqX25ubmbbtm0feHtTUxO9vb0A/Pu//zvZbJZrr7128j/JURoaokc9b+yEnkOmhvqjcqgvKof6onKoLyrLifbHpMLNeeedx/PPP8/OnTvxeDzMmjULv99/zMccvXx8nGEYH3p7f38/3/ve9/jhD384maa9r8HBJI7j0tQUo78/ccLPIyeX+qNyqC8qh/qicqgvKst4f5imMWHQYjKOGW7Wr1/PjTfeyA9+8IMJ12/cuBGAL3zhCx/42JaWFrZs2VK+3NfXR3Nz84TbBwYGypf7+/tpbm7m+eefZ2RkhDvvvLN824033sgjjzxCNHp8P5yIiIicfY4Zbvbv3w/Azp07j/uJV65cyXe/+12GhoYIhUI888wzfPOb3yzf3t7eTiAQYOvWrVx44YU88cQTrFq1is985jPlQzkBzjnnHNavX3/c/7+IiIicnY4Zbu69916ASS37freWlha+8pWvsG7dOgqFArfeeitLly7l7rvv5t577+Xcc8/lwQcf5C/+4i9IpVIsWrSIdevWndhPISIiIlJiuO9X/FJyww03HPPBTz755Elv0MmgmpvKpP6oHOqLyqG+qBzqi8oyZTU3999/PwA9PT089NBDfP3rX2f37t08+OCDfOMb3zjxFouIiIhMkWOGm+XLlwNw1113cfvtt7N8+XLOO+88crkcTzzxBNdff/0paaSIiIjIZE1qE7/h4eFyPUwgEOCuu+6iv79/ShsmIiIiciImFW5s2y5vsAcwMDDwvvvUiIiIiJxuk9rE76677uKmm27i8ssvxzAMNm7cyFe/+tWpbpuIiIjIcZtUuLn11ltZsmQJL774Ih6Ph9/7vd9j/vz5U902ERERkeM2qXADsGDBAhYsWDCVbRERERH5yCZVcyMiIiJyplC4ERERkaqicCMiIiJVReHm/2/v/mOqrvc4jr+OHtG43jL1wBxTW2YyU9PNNdSFs0KcgJjQKit0JLlW2awxceosp2bOSf5k5dz6I12wtjSaQ62mbUhbsiVtWUhpFhAgOJOT4kE+94+7zlIxf8DhHN/n+fjL8/1yzvnge5/t6TnH8wUAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKaENG5KS0s1c+ZMpaSkaNeuXecA5+UAAAxHSURBVNecP378uLKyspSamqply5apvb1dklRZWamsrCxlZmZq3rx5qq2tDeUyAQCAISGLm4aGBhUWFmr37t3au3eviouLVVNTc8XP5Ofna8WKFdq/f7+ccyopKQkeX7Nmjfbu3auMjAytXr06VMsEAADGhCxujhw5oqSkJA0YMECxsbFKTU1VWVlZ8Hxtba0uXryo8ePHS5LmzJmjsrIyXbp0Sa+//roSExMlSaNGjVJ9fX2olgkAAIzxhuqBGxsb5fP5grfj4uJUVVV13fM+n08NDQ2KiYlRZmamJKmjo0Nbt27VE088cUvPPWhQ/3887n9v91dACDCPyMEsIgeziBzMIrLc7jxCFjfOuWuOeTyemz5/6dIlFRQUqL29XQsXLryl525ublVHh5PP9181NZ2/pfsidJhH5GAWkYNZRA5mEVn+nkevXp4rXrS4GSF7Wyo+Pl5nzpwJ3m5sbFRcXNx1zzc1NQXP+/1+LViwQO3t7SoqKlKfPn1CtUwAAGBMyOJm8uTJqqioUEtLiy5cuKADBw4oOTk5eD4hIUF9+/ZVZWWlJGnPnj3B8/n5+Ro+fLg2bdqkmJiYUC0RAAAYFLK3peLj47V48WLl5OQoEAgoOztb48aNU15enhYtWqSxY8dqw4YNWr58ufx+v0aPHq2cnBz98MMP+vLLL/XAAw9o9uzZkv7/eZ0dO3aEaqkAAMAQj+vswy93OD5zE5mYR+RgFpGDWUQOZhFZIvIzNwAAAOFA3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYEpI46a0tFQzZ85USkqKdu3adc3548ePKysrS6mpqVq2bJna29slSXV1dXruuec0Y8YMvfzyy/L7/aFcJgAAMCRkcdPQ0KDCwkLt3r1be/fuVXFxsWpqaq74mfz8fK1YsUL79++Xc04lJSWSpLfffltz585VWVmZxowZo+3bt4dqmQAAwBhvqB74yJEjSkpK0oABAyRJqampKisr06uvvipJqq2t1cWLFzV+/HhJ0pw5c7R582Y99dRT+vbbb7Vt27bg8eeff175+fk3/dy9enk6/TPCj3lEDmYROZhF5GAWkaVXL89tzSRkcdPY2Cifzxe8HRcXp6qqquue9/l8amho0NmzZ9W/f395vd4rjt+Ke+/9T/DPgwb1v91fASHAPCIHs4gczCJyMIvIcrvzCNnbUs65a455PJ4bnr/R/QAAAP5NyOImPj5eZ86cCd5ubGxUXFzcdc83NTUpLi5OAwcOVGtrqy5fvnzFcQAAgJsRsriZPHmyKioq1NLSogsXLujAgQNKTk4Onk9ISFDfvn1VWVkpSdqzZ4+Sk5PVp08fTZw4Ufv27bviOAAAwM3wuM7eB+ompaWlev/99xUIBJSdna28vDzl5eVp0aJFGjt2rH788UctX75cfr9fo0eP1jvvvKOYmBjV1taqoKBAzc3NGjJkiDZu3Kh77rknVMsEAACGhDRuAAAAehrfUAwAAEwhbgAAgCnEDQAAMIW4AQAAppiNmxtdtBM9JycnR2lpacrMzFRmZqaOHTsW7iVFndbWVqWnp+v333+X9P/Lo2RkZGj69OkqLCwM8+qiy9WzWLp0qaZPnx7cHwcPHgzzCqPH1q1blZaWprS0NK1fv14SeyNcOptFl/aGM+iPP/5w06ZNc2fPnnV+v99lZGS4EydOhHtZUamjo8NNmTLFBQKBcC8lan333XcuPT3dPfTQQ+63335zFy5ccFOnTnWnT592gUDA5ebmukOHDoV7mVHh6lk451x6erpraGgI88qiT3l5uXv66addW1ubu3TpksvJyXGlpaXsjTDobBYHDhzo0t4w+crNPy/aGRsbG7xoJ3reL7/8Io/Ho7y8PM2aNUsfffRRuJcUdUpKSrRy5crgN31XVVVp+PDhGjp0qLxerzIyMtgfPeTqWfz111+qq6vTihUrlJGRoc2bN6ujoyPMq4wOPp9PBQUFiomJUZ8+fTRixAidOnWKvREGnc2irq6uS3vDZNx0dtHOW734JrrHn3/+qUmTJmnbtm368MMP9fHHH6u8vDzcy4oqa9as0cSJE4O32R/hc/UsmpublZSUpLVr16qkpERHjx7VJ598EsYVRo+RI0dq/PjxkqRTp05p37598ng87I0w6GwWjz76aJf2hsm4cVx8M2JMmDBB69evV2xsrAYOHKjs7GwdPnw43MuKauyPyDF06FBt27ZNgwYN0l133aUXXniB/dHDTpw4odzcXC1ZskTDhg275jx7o+f8cxb3339/l/aGybi50UU70XOOHj2qioqK4G3nnLxebxhXBPZH5Pjpp5+0f//+4G32R8+qrKzU/Pnz9eabb+rJJ59kb4TR1bPo6t4wGTc3umgnes758+e1fv16tbW1qbW1VZ9++qlSUlLCvayo9vDDD+vkyZP69ddfdfnyZX3++efsjzBxzmnt2rU6d+6cAoGAiouL2R89pL6+Xq+88oo2bNigtLQ0SeyNcOlsFl3dGyb/iRAfH6/FixcrJycneNHOcePGhXtZUWnatGk6duyYZs+erY6ODs2dO1cTJkwI97KiWt++fbVu3Tq99tpramtr09SpUzVjxoxwLysqJSYm6qWXXtKzzz6r9vZ2TZ8+Xenp6eFeVlTYuXOn2tratG7duuCxZ555hr0RBtebRVf2BhfOBAAApph8WwoAAEQv4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAHSb3NxctbS0KC8vTzU1NVcck6THHntM33///W0/flfv/2+6c50Awsvk99wACI+/rxu2Y8eOa45FujtlnQBujLgB0C2WLl0qSZo3b55qampUUlKi3bt3B4998MEHV/z8V199paKiIgUCAfXr109Lliy5pS94vN79t2zZotraWjU1Nam2tlYDBw5UYWGh4uPjVVVVpbfeekuBQEDDhg1TXV2dCgoKtGfPnmvWWVxcrJUrV6qlpUWZmZlavHhxd/w1AegJDgC6yYMPPuiam5vdtGnTXFVV1RXHnHPB4ydPnnTp6emupaXFOedcdXW1mzJlivP7/f/6+Ddz/82bN7vHH3/cnT9/3jnn3MKFC92mTZtcIBBwycnJ7tChQ8455yoqKtyoUaPcN9980+k6V61a5ZxzrrGx0Y0ZM8bV1dV1518VgBDilRsAPa68vFyNjY2aP39+8JjH49Hp06eVmJjYpftL0iOPPKL+/ftLkkaPHq1z586purpakjR16lRJUlJSkkaOHHnd5/j7q959Pp8GDx6s5uZmDRky5JZ+TwDhQdwA6HEdHR2aNGmS3nvvveCx+vr6m74C87/d/+DBg+rXr1/wuMfjkXNOvXv3lrvqajO9e/e+7nP88wrEfz8GgDsD/1sKQLfp3bu32tvbb3gsKSlJ5eXl+vnnnyVJhw8f1qxZs9TW1nZTz3M79x8xYoRiYmL09ddfS5KqqqpUXV0tj8dz3XUCuDPxyg2AbpOSkqK5c+fK7/dfc2z79u3BYyNHjtSqVav0xhtvyDknr9eroqIixcbG3tTz3M79vV6vtmzZopUrV2rjxo267777NHjw4OCrPJ2tE8CdiauCA4ga7777rl588UUNHjxY9fX1yszM1BdffKG777473EsD0I145QZAxPjss8+0c+fOTs9lZGRowYIFXXr8hIQEzZ8/X16vV845rV69mrABDOKVGwAAYAofKAYAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABM+R9B2RBy/jK5jAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lineplot(data=df, x='title_length', y='click_through_rate'); pp.ylim(0, 0.1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see now the click-through rate depends on the title length! From the data-generating process, we know the dependence isn't causal: it was introduced because one quirky author happened to like 12-16 word titles, and happened to be really good at writing titles. If we were to control for author, this effect would go away, and we'd recover the underlying causal relationship! Let's use the do-sampler to do the adjustment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:dowhy.do_why:Causal Graph not provided. DoWhy will construct a graph based on data inputs.\n",
      "INFO:dowhy.do_why:Model to find the causal effect of treatment ['title_length'] on outcome ['click_through_rate']\n",
      "INFO:dowhy.causal_identifier:Common causes of treatment and outcome:['author', 'U']\n",
      "WARNING:dowhy.causal_identifier:There are unobserved common causes. Causal effect cannot be identified.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARN: Do you want to continue by ignoring these unobserved confounders? [y/n] y\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:dowhy.causal_identifier:Instrumental variables for treatment and outcome:[]\n",
      "INFO:dowhy.do_sampler:Using WeightingSampler for do sampling.\n",
      "INFO:dowhy.do_sampler:Caution: do samplers assume iid data.\n"
     ]
    }
   ],
   "source": [
    "causal_df = df.causal.do('title_length', \n",
    "                          method='weighting', \n",
    "                          variable_types={'title_length': 'd', \n",
    "                                          'click_through_rate': 'c', \n",
    "                                          'author': 'd'},\n",
    "                          outcome='click_through_rate',\n",
    "                          common_causes=['author'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAE8AAAAUBAMAAAA+U32BAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAImYQu82Z3XZU70SJMqtHo0JKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABS0lEQVQoFY2SMUvDYBCGnyRNY2ii3RyNqwg6CJUiVsVZ6ugW/AV1km7iIE5Sf0EFRxHcHO2idFH6D9pZEASpKIh4l4SWrxbrO3x3997L3XfHgRUxEeeqWIHg+PG3tr3YGZDTRXF3YY+SegbchnOQEP4lOFWwY+wLAvFMLMC9MutPb/K2wSuS6+FrZGAHmsvK2JoqwyzM9Ch8GioJvuGuoWQidG84gUpM4WNE6LyL8GggDGK2YK5FOCq0vqA7FNpV6iKMCYU3YEnFbkuppLVf4+pfrcNXFcow/rhhmsNhRCitcxtpfaP3NiwN11Oo6TB/LjytIcM8SJ065Y7TNwriNZxrKlWY0oXnW9yKWTt9gcOkkUQpnPZZB69H/rm/D674xSzjm8KMzUxJPhhnfmCmzEiOgs2MWjVTRqRnpoebIMrsOJPT/1mRPBMwDz+CNk+5m8RT8QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( 0, \\quad 0.1\\right )$$"
      ],
      "text/plain": [
       "(0, 0.1)"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGyCAYAAAAYveVYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYFFeh//9P9d49PSuzEIZANhLMQjDGiFGJaMhEApIAapYr4oLGJxrN1Zj4TXhiEpfoRblJ3P25XW+4gvdREBfALfEaMAaMEg1ICPsAsw8z09N71e+P6q6ZgWEYlmaamvfrMdJV1cvpOdVVnzrnVJVhWZYlAAAAl/CMdAEAAABOJ8INAABwFcINAABwFcINAABwFcINAABwFcINAABwFcINAABwFcINAABwlYKGmzVr1mjWrFmaOXOmnnrqqWM+77777tNPf/pTZ/rAgQO64447dOONN+rDH/6wYrFYIYsJAABcpGDhpqmpScuWLdPy5cu1evVqrVixQjt27DjqOXfeeafWrl07YP7DDz+s22+/XWvXrtXll1+ur3/964UqJgAAcJmChZsNGzZo2rRpqqioUCQSUUNDw1EhZs2aNXrrW9+qt73tbc68dDqt559/Xg0NDZKkefPmHfU6AACAY/EV6o2bm5tVU1PjTNfW1mrLli0DnvOBD3xAkrR582ZnXkdHh6LRqHw+u2g1NTVqamoqVDEBAIDLFKzlZrD7cRqGUbDXAQAASAVsuamrq9OmTZuc6ebmZtXW1h73dVVVVerp6VE2m5XX61VLS8uwXtdfR0dMpmlpzJio2tp6TrjsKAzqo3hQF8WDuige1EVxydeHx2OosrLkhF5bsHBz7bXX6sknn1R7e7vC4bDWr1+vRx999Liv8/v9uvrqq/WrX/1Kc+bM0apVqzR9+vQT+mzTtGSalvMYxYP6KB7URfGgLooHdVFcTrY+CtYtVVdXp3vuuUcLFy7UzTffrNmzZ2vKlClavHixXnzxxSFf+9BDD2nlypWaNWuWNm3apI9//OOFKiYAAHAZwxpskMtZrq2tR6ZpqaamVC0t3SNdHORQH8WDuige1EXxoC6KS74+PB5DY8ZET+i1XKEYAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4SkHDzZo1azRr1izNnDlTTz311FHLt27dqvnz56uhoUEPPPCAMpmMJGn//v264447NHfuXL373e9WY2NjIYsJAABcpGDhpqmpScuWLdPy5cu1evVqrVixQjt27BjwnHvvvVdLlizRunXrZFmWVq5cKUl6/PHHddNNN2n16tW64YYbtGzZskIVEwAAuEzBws2GDRs0bdo0VVRUKBKJqKGhQWvXrnWWNzY2KpFIaOrUqZKkefPmOctN01RPT48kKR6PKxQKFaqYAADAZXyFeuPm5mbV1NQ407W1tdqyZcsxl9fU1KipqUmS9LGPfUy33nqrfvSjHymdTmvFihWFKiYAAHCZgoUby7KOmmcYxrCW33fffXrkkUd0/fXXa926dfrIRz6in//85wNeP5QxY6LO45qa0hMtOgqI+ige1EXxoC6KB3VRXE62PgoWburq6rRp0yZnurm5WbW1tQOWt7a2OtMtLS2qra1Ve3u7du7cqeuvv16S1NDQoIceekgdHR2qqqoa1me3tfXINC3V1JSqpaX7NH0jnCrqo3hQF8WDuige1EVxydeHx2MMaLQYjoKNubn22mu1ceNGtbe3Kx6Pa/369Zo+fbqzvL6+XsFgUJs3b5YkrVq1StOnT1dlZaWCwaATjDZv3qySkpJhBxsAADC6FbTl5p577tHChQuVTqe1YMECTZkyRYsXL9bdd9+tK664QkuXLtWDDz6oWCymSy+9VAsXLpRhGPrqV7+qRx99VIlEQiUlJXryyScLVUwAAOAyhjXY4JezHN1SxYn6KB7URfGgLooHdVFcirJbCgAAYCQQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsQbgAAgKsMO9wcOnRIzzzzjLLZrA4ePFjIMgEAAJy0YYWbp59+WrfeeqsefvhhtbW1adasWfrtb39b6LIBAACcsGGFm69+9atauXKlysrKVFtbq+XLl+uJJ54odNkAAABO2LDCjWmaqq2tdaZf9apXyTCMghUKAADgZA0r3ITDYR04cMAJNJs2bVIwGCxowQAAAE7GsMLNJz7xCb3vfe/T3r179a53vUt33XWXPvnJTx73dWvWrNGsWbM0c+ZMPfXUU0ct37p1q+bPn6+GhgY98MADymQykqTm5mZ98IMf1M0336xbb71V+/fvP8GvBQAARqthhZurrrpKK1eu1Fe+8hXdeeed+vWvf61rrrlmyNc0NTVp2bJlWr58uVavXq0VK1Zox44dA55z7733asmSJVq3bp0sy9LKlSslSZ/61Kc0Y8YMrVq1SnPnztXSpUtP8usBAIDRZljh5gMf+IDKysp03XXXacaMGaqqqtI73/nOIV+zYcMGTZs2TRUVFYpEImpoaNDatWud5Y2NjUokEpo6daokad68eVq7dq3a29u1bds23XrrrZKk+fPn6+Mf//jJfj8AADDK+IZaePfdd2vXrl3at2+f5syZ48zPZDLyeIbORc3NzaqpqXGma2trtWXLlmMur6mpUVNTk/bt26dx48bp85//vJ577jmNGzdOS5YsOeEvBgAARqchw82nPvUpNTY2asmSJQMChtfr1aRJk4Z8Y8uyjprX/wyrYy3PZDJ66aWX9NGPflQPPPCAfvKTn+j+++/Xj370o+N+mbwxY6LO45qa0mG/DoVHfRQP6qJ4UBfFg7ooLidbH0OGm/Hjx2v8+PFau3btUS01vb29Q75xXV2dNm3a5Ew3NzcPOJ28rq5Ora2tznRLS4tqa2tVU1OjkpISzZgxQ5I0e/Zsffaznx3+N5LU1tYj07RUU1OqlpbuE3otCof6KB7URfGgLooHdVFc8vXh8RgDGi2GY8hwk/f73/9eTzzxhHp7e2VZlkzTVGdnp1544YVjvubaa6/Vk08+qfb2doXDYa1fv16PPvqos7y+vl7BYFCbN2/Wa17zGq1atUrTp0/XhAkTVFdXp2eeeUbXXXed/vCHP+iyyy47oS81GqUzWTV3xGVakseQPIYhwyN5PIY8MmR4DHkMQx5DkmHIMCT74cDHkv1aGfb7GIYhr8fgukY4JtOyZJqWva65eD1JZ7JKpk0lUhkFfF6FAl4F/N6RLtZZxTQtZU3TWVfYrqBQhhVuvvSlL+njH/+4/ud//keLFy/Wb3/7W5WUlAz5mrq6Ot1zzz1auHCh0um0FixYoClTpmjx4sW6++67dcUVV2jp0qV68MEHFYvFdOmll2rhwoWS7CsiP/TQQ/qP//gPRaNRPfbYY6f+TV0smc5q18EuZbKmfF5D+R4/y7JkSfa0JVmyJMtQ7pE8hiFL+SfkEo4jNy/3j8djyOcx5PMa8no98ns98noN+bwe+bweZ2PlzYcojyGPR2zAXMSyLGWyplIZU8lUVolURr3JjOLJrPIrisdjyJ9bJ/y59cPv88jns9cRb//1xFO8wTmdMZXKZBVPZhRLpBWLZ5TJmvZPpO//5PcaKisJKhrxKxzwyu8j7OSZlqV02lQynVVvIq2eRFrxREbK1bchyesx5PN5nO2I35dbZ7z91pd+/57s9sS0LHt7aNnrsZn798hpf3dSPfG0fP22bTg7GdZgg1+OMHfuXK1evVpf/vKX9cY3vlGvfe1rNX/+fP3sZz87E2U8YaOpWyqezGjXwS51dCXUk0jL57M3DP7cvz5f32O/z3PSOxLTtHIbAXtD0DctWaYdl+xt/tHv7/XYG4q6mqjSibRCQZ+CfnvDUYw7ttHgeL+NdMbM/ZdVbyKj3mRa8VTWCc6G1C/c9u1wTNNyWnLsf+0diyz7X6e2+z3Irx8DApHfI7/Pm1tmyOuxw3QhWoYyWVOpdFaJVFaxuL0TzmRMu4VTks9rKJHK6lB7rxpbY0oks5owNqrzzylTOOBVMm0qm7VkyVLA71V5JKBoxK9QwCe/7/g7RzdspyzLUjpjB5l4KqOeXnudyYcHn9fI1aecIQ75cDHYOuPslfpXd+5AK7+++HOhyGPkW4Ry2ybLlGlKWdN0tlX5t7Lyb3IMFeVhdR6OOx9oyFDA71U44FUw6FXQ57XDmMcjn89eL1E4Be+WCgQCSqVSmjBhgrZu3arXve51SqVSJ1VYnD69CTvYvLy/U7/YuMf5EQ/F6zHssNMv8OQ3FEc+Dvg9Cvq9ff8F8o89A6aPF5ryG61EKquWdrtrU7nWnZKgT9Gwf0QDj2lZymYt5+iwWORbSpJpO2gYRl83Yr6L0elWHHT6+N8lkzVzrRSmEkm7JaY3mZGZNXObdjvE+H2GSkI+GYahVDqr1sMJtR5OqKUzrvaupAxDCgW8CgV8CgXsdSPk/OdT0O915h+5vuTXj0w2q1SmLxBZ+fXZMOy9Xa7FJ+D3KOjzyu/3KODsbPoC0FCtQVnTVDJlKpnOqCdut8qk0lkZsv9uPp+hdDarpo6EDrTGdKA1poNtveqJpyXZf3e/z6PN21skSXWVYZ0/rkznn1OmiXVReTyGOnrsv41kKRjwqSziV0l4+GHnbDCgZSueVk8iY29/LMnjUW5bYqi5I65D7b061Narg229ajmckNdjKBL0KRLyDfg3fMR0JOTPzfc6IaJ/IEpnskqmMk4rtL3eS4YMeT2Sz+t1fgtHsnK/+XRu/c//6/X55PMaCgd9A54XS6bVFU/JzGpANvJ4DIX8HmcdD/i9fa0+Po+ru2mL3bDCzVvf+lZ98IMf1GOPPaZbb71VmzdvVkVFRaHLhiH0xNPadfCwtrzSpt9tbtSEuqhmXj1eWdP+wWYyptKZfj/eTFaZrOUckWec+X0/7GQ83e+I3d54ZbLHD0yGoQEhKBA4OhSFAl5NPKdC0ZBHVaVBGYYh07SUzGTV054uaODJZO0j64xp2n+DdFbJtD1+IpnO2N8xd6gY8HudDW1+Z+Q/AxupfBN+KmMf+cbiuZBht5fLMPLdjdagrWOOfBdjri3Nbu2QvB5PrpvQPtI9nMio8VCXstlcn6X6WmLCAa88Hp96E2m1HE6otTOh1sNxtXYm1HI4oa5Y34GNxzBUWRqQJCVSdutH9jgh22MYR4Sf/Dric6ZLQn5FI36Vhv2Khv0qCfvk9XicI/RYMi0zbj+2LHuHZu907O/v99nBJ+C3Q3o6bSqWyCiZzub+TnaQMU1LbV35INOrA60xHe73/arLQ7pgXJnGVZdoXHVEdZUReT2GDrX3aueBLu082KXntzbrz/9sksdj6NyaEl0wrkznjyvTuDElMi1L7d32302yFPT7VFbiVzQcUCjgLbpuj8Ea8rOmpVRurFEskVFPPOW0bHlk5P6OphNkDrb16lB7r9q7ks57REI+nVMV0QXjymRadotzbyKt3mRGbV0JxZPZvroZRNDvVSTkU/iIUBQKeO1t3hHbtEy/0JJxtnMDnzeUoN+rytKgKkoDqowGVVEaVEU0aM8rCciXC6n5MUTdvSl1mPlAnlsXLUs+n0flJQGVRgIKB31FV99uNqxuqdbWVqVSKY0bN04vvfSSnn/+ec2ePVtjxow5E2U8YW7vljrck9TuQ136y9ZmbfxnkyZPqNAt0y8oyFFhNtdykHQCQVbJVFYpJxzYO7RUbv7A5/W9rv/GJODzqK4qorFVYY2tiqiuKqLairB8Pk/uiMze+Fj9jtZLQj6VhPoCT35sg93qYoeWvlaIvs9PpbN9+3vJDgq5PnxvblxQIve8imjQbl7PmspkcuOVcu0XQb9XJbmNa8DvVcB38t18WdNUKm23lvTG0+pNptWbzNqfZdllyncvGoYdZOPJbK5pv298Qv/uoGPpP67Ayo27siypvDyiWE9ChiF196bV0hm3W2M6E2rJBZneZMZ5H5/Xo+rykGoqQqouD6m6Iqzq8pCqyoJHNc1nMqYddNL2uJxEKpsbo3PEdLpvXt/y7DF3PJGQT6VhuxWkNGyHn2gu/PQPQn5fXwjK/5tvkWs93Ncic6C1V21dCef9K0uDOmdMJBdkSnTOmIiCwxgwnM5ktbepRzsPdmnXgS4dare7NUIBr84bW+qEnarSoLJZO9BncwcNoYBP542vVOfhmNO96/xrSWZu0JzpzFeuPk1Zkszcn2rAmJKBK4AG28IfOevI5/Rfrax8nrYkn1eKpzJq7kjoYK5F5lB7r7p7087zK6IBja2KaOyYiP1vVUSlEf9x19Vs1nRaDnsTmVwA6psebH5+XRnQDe+M8zIGdNHnx34d+bz+ry2JBNXY3KXO7pQ6epLq7E6qsyd51EFeacTfF3aigVwQCqoyGhzwXfMthZmsKcMwFA35VR4NqCTkV8BPt/zxnEq31LDCzY033jjg6sLFzs3hpr07od2HuvXMC416cWe7XnNJjd72ugkyDOUGduabZpU7gDCcDVW+yTZ/YJH/YR2r6fZ0ymRMJbLSy3vadKjdPsJrau9VKrdx8hiGqitCubATdjaK4aDvmIHHMNS30cklGMOwm8W9Ho8TXnriGXX2JHU4llJnT0qHe5I63JN7HEspk+3bQNZW9n12XVVYdZVh+bwepwk7m7WcMhiGFA54FQ7aocfuovMOCB35VrNEKutskJPprKzc+BNvv0GUliW1dSXsZvz2XjXl/k69icwgf1FbvpuxbwyCMSD89A9D9gbcXp7MWDrQ0q3Wwwml0n1hIhTwqroipJrysKorQqouD6umPKTyaOCMbYizWVM9iYx6elPqiWfUE0+ruzelWDyj7nhaPfGUenrT6oln7PE8Rwj4PAOCj9/r0aH2XjV3xp2deGnEr3Fj7NaYc6pLNG5MiSKhwRuy892D9lF/vhXNktdrOGMw+osl0tp9sNtu2TnQ5bQElZcE7KBzTpnOH1eqkpBf6YypSElInYdjyg+zdX6vuf/r/5uV8sHD6LdczgsMGf2eM7R8GMrmAmD/MJjNtUhkspbaD+eCTC7MJFJ925nq8pDze8mHmXyXzmCfl83aXZD5blTPKQ4oz5rmaT1pobKiRB2dsQHzLMtSTzytju6kOntSuX+TznT/lkzJ/k1WRAOqKA2qujysC8eVaeLYUvm8ht1CnjJlGZb8Xo8qo0FFIwGFg17G7wyi4OHm/e9/v+666y5NnTr1uFcmLgZuDTetnXHtOnRY6/6yXzsPdOnNrx6nN005R5YldfWmNaY0KJ/P069fWjJlyjKVOxug72iv/5GeaVp9R3LOYVr+sb3BzJ9Vld+w64jxH4ON+zhyo3PkhsOyLHV0J52d+aH2uJqOOAosLwkMaOUZWxVReTTgtEIYMtQVS6kzZoeWfHixp+3wcuRYpEjQp4poQOXRoMqjAVVEA/J7PWrujOeORONOE7lhSGPKQk4Z6nJliIb9uZ2evePLZOzPsGSP3Qk5A01N5XdY3n5hI5XOqqmjL+Qdau9Vc0fcCWtej+EErbqqiEpCPvtz+jWtZ7KmE7j65uVannLN8Ue/xl5eGvarqiyomlwLTL41Jj+u5liyWVMZ0/7OZtZeV6zceuGM+FTfQzvE9YXO/FggZ/04hTPqLMtSPJkPPOlc4Ol73B1PKxZPK5UxVVsRdlpkxlVHVBoJHPM9890ZZr8GpEjQq5KwX5GgfcSdyVrqiafUFUsrmcpIMuT1KjfmwjPg/dq7k9qV68LafbDbCQdjq8I6/5wynVdfoZ6exKAB4+jH5nGWD3zuUYElFy7y/w6X12OorjLc1xozJqK6yvAxzw4z+3WPm/mez1z3dcDncbqI8+U7apxvfjuU+ye/nvRfdzyevu1O/5bJ3P+cblzLmT72svzZpGXlEXUfjjtjt7ze4w8azmRNHY7lQk93Mtfik1JnT1LNHXFlTXtA9cSxpbqovlwX1Zerqiwo07SUSGeVzdqX3SiN+FVeElAk5OcSAzkFDzdvfvObdejQIfl8PgUCAWdj9te//vWkC11Ibgs3lmWpubNXOw9065cbdutge69uev1EXXVxjTJZU73xjMbVlKi6PHzKn5PfCFhWvvm73xkMyp/R0Hc0m83a1z3KmPaOz8zt8PMbYSdXWJYqKiLq6Oy1d/JDnL0Vi6ePar2wB2jaQgGvqsqCisUz6upNHdWkHg37c6ElqPKSgBNkKqIBlZcEjrvhsCxLh3tSR5Wh/1iMaNg/oIWpriqiqtKgPB7D2bF4PX3dSvnglg8zbf3GI4SDXic05VuMqstDBT2SG+wItf/3zwch+ywgZ4H8PnvwZDjoUyjoy9WjvQMY7CypvsHC9rqSNXNhLBcK89NHboRyQxYk2aGo/2nkntM08PvYQcankrA9oDU/cHmoMVfpjD0epac3rcO9KaXS9qhTn9dQ0O+Rt1/YMU1LB9tiufE63drX3HPcEwHyXaj5s8Xsx55+j/v+LgOemx9kPaAbduC/A97Le/TzvB6PKkoDx1wfBwv4kiWv12OPiwnbY9eCPnsg+WD11n8b46w/Vm48lWk5YTqb7Vt37FalvnXHkJywkw/L+TFm+YHGXiN/vS/7gM3jMQaEJo/HUPWYqA41dznd7olUZkDrsGEYuQMYz4C/17HCeTqT1Z5DPdrReFivNB52fvcV0YAuGl+uC+vLdf7YUvl9HrurOp2VZdgtgpWlAWfs4WgdmFzwcNPY2Djo/Pr6ev3lL3857h3CzzQ3hRvTsnSorVc79nfq5xt2qyuW0vzrLtQlEyqccS4Tx5aqrCQ40kUdlNM6ZFmqqCzRgYOH7TEmiax9qqjZ18WT7/sebAOYSmfV3Bm3g0Zbr9q7E4qG7X5vJ8hEAyqPBI7qJhiqbNlsXxN8/kyY/juj/uLJTC6g9AWVls6EcwTcv1sr4PeoKfe8WL9upYpowDn6zQeasmOMR+jfHZLNj7U0+lpC+nc7Ot2MzmMN/rzc51RWlKi1rdvZMWXNvtO0PbmWp3DQ3jEF/F4njBbqbLJ8ELJyO7lsv5aHdMYeh5Mfp5TOmjJzYxj6srNl78A8A3f0+W4PJ8hkTGVNyTDsHVUk6LfHUoV8wwoyw5HOZBVPZdXTa7ccpnMD1vNnIPYPCal0Vl6/Xz2xhLz5cnsN53Exnb2XNe0Ak87arcFWbl0MB3OD8IP9zxgq/hb+wQy2zxhsXF9+fGEqkx/X17/FW/J6+7rG+2+POrqTTtDZdbDbPkPLY2hCXVQX5lp1aipCymbtVh3LtMcIlkfsA7Rw0D1n3A1HwcPNUG655Zaiu96NW8KNaVo60NqjbXs6tfrZXcqalm67fpLOrY0qkcwoa0nnjy075liBYnNkfViWZe+scqch98Ttsyey2YHXGBkqcAwl3yqQb5LPByn7w/OnrHoVzJ3dlU6biiUzSqXNfEdcbkzMsQfvZrKmWjsTA7qXDrXHlc7aXSH9u7LqKsMKDTEe4ZjdISF/36mpkvNdsvmuRdNS1lK/eZKZu8ZHPliaZl+rm2FI5eVhJXpTCgY8dojx950ZdjbsmEwnmJpO/ebPismvU/lQZMn+zoUIMsORv4ZOdy7sZLN2mfw+e70bUxU9Ziva8eQ33/mt+Il0yzjPdR4bzns463r+/XPljYR8Kgn6FcxdnflMnEl4Jp3MPsPKtTLZwTl3+YZUVol0xrnsgGTI77XPEsz/bTNZU3ubevRK42G90til5k57IHpZxG8HnfHlOv+cUgV8XqXS2VxIls6pjmhMWWhUDEYu+HVuhnKK2QjHkDVN7W+O6cWdrfrFxj0KBXxaeOMk1VSE1RvPyOv16KJxpcM6m6NYGYbhnC4eDftVXWF3qzlnPKWyiiXsa5H0JtLKd8znB8UahpxxBFnTtMcG9TuCsrsFfCoJ2TsRf/4aFJ6+qysPxjnrKm0PAo4l7cGr+bOZDE+/szO8HrslZkzEeX2+tepYR9358Qj2AOn838Le+ZaXBAu+8zUtSzXVpWpr6znt732meAxDHp8hv44fxE73oNMTFci1ZpSVBDSu2g70iVRWXbFULvAk+40z6z9uyeobVZyfbeXHOfWdOGB47AMBZ4Cx080ieYz8tV6O6Jbp95y+x/bfyGn56ze+5WwJvSPBMAxnIP9g0hlT8WTfSQ1WrmUn5PfqgnFlumBcmWa+VjocSzlB56XdHXrh5VYZhnRubVQX1dtdWLUwYG3HAAAgAElEQVQVIR1oiam7N6X66ihjc4ZwyuFmNKTHM81O9N366/YWrf3LPlWXhXT7zEkqjfjVHUupJOzXhLpS125s8i0IJSG/qsrsef0DRyyRViyRkWna42/sQYrefmcG9Y0jOBn5s4vCQZ/Ko3Z335HXoenNXYcmm7vYnceww1T+UvL5j86aubNsMpYzpsnr9agk6NOY0pBCQd8pnVZ+MvI7stGimM5C6R/oy0vs8YtlFRG1tvQM6FLMBxznzKlBz5oaPXV4NrO3ZwGVlQRkmpbiqYx9XZzupNLJjAzLUDDgUVnEr6surtFVF9fYB7ctMb3SeFg79nfp939t1O//2qjairDe/sbz5PMZenn/YY2vKXG2URjo7OjPGEXSGVN7DnVr4z8O6Q9/sy/O9663XKRQwKvuWFpVZUGNq46Oqp2TNHjgOJPyF54LBrz2mTa5a1j2b2XKn+rdnew7Cvd5PYqG/AqX5wZW9rtGD2AYhr1eBFgnRgOPx1BJyK+SkF91lRGl0qZ6Eml1dudb7+yrgQf9Xk2sK9XEulK95Sqppzetl/d36um/HdB3f7lVb5pyjl5/WZ12H+rWmPKUxlaVuPZg92QRbopIKp3VzgNd+sMLjXp+W7MmT6zQvDddIMNjX2ittiqsusoIR2xFZLBWpnxrDTfeA3AsRr8DpjFlIfvM12TGvg5XLCUrd/HJYMCraMSvV19co8kTK7XuL/v0x78f1PZ9h/X2N56nzm5DPfGMJtRGFQn5R/prFQ3G3BSJZCqrlxs7te65vfrn7g5dfUmNbnzdBJmmpVhvRuProqoqDY10MTEMXo9H3gChBsDw+bwelUUCKosEVG9ZSiQz6u5Nq6MnqXgiJcnuvrr5Tedr8oQK/XLjHv1/v9iqN08dp6svqdGOxsMaWxVRdUXYVYO8T9Yph5vXve51p6Mco1o8mdG2vR36xYbd2nWw27k4XypjKp02dUF9uaJhEjkAjAYew7BvHBryq64qkju5Iq32LvvebpPOLdeddZfp13/eq9//tVHb9nbq7W84T4fa4+ruTWt8bfSsPtnkdBjWqeA7d+7U9773PbW1tQ1oqfnmN79Z0MKdrLPpVPBYIq1/7mzX6md36VC/i/PFkxnJks47p+yYlzM/25wN9TFaUBfFg7ooHsVeF/mruh9ojcnjMRQJ+fTPXe361Z/3KJ0x9Zar6jXlwjEyTam+pkQV0eBZPYyh4KeCf/KTn9RrXvMazZw586z+QxWbrlhSf9/RqlV/2q3u3pTe+ZaLdMm5FerpTSsYsAeUcaofAECyx+lUlYUUDft1sK1XnT1JXXJuhSbWleoXG3dr/fP7tW1vp2a/fqL2NveouzetcdWjc7DxsMJNOp3WAw88UOiyjCrt3Qn9dVuLVj+7S6Zl6d0Nl2h8TYm6YimVRQIaXxsdlSskAGBoAb9XE8eWqjIW0P4W+wKQ75xxoV7c2a61z+3Td36xVddfPV6vmlCpnnhaE+pKR93QhmHtPceNG6d9+/YVuiyjQv4+Uc9uOaj//eMr8no9WvS2yaqvLlF3LK0xZSFNGOvea9gAAE6PspKgLj63QpWlQXXHM3rVxErdefNlOrc2ql//ea9+9n871ZtMa2fjYR1qix33PmZuMmTLzZ133ilJamlp0YIFC3TFFVfI5+t7SbGOuSlWpmXpYGtMf9pyUL/ZvN+5OF9JyKfu3pTGVUdVXT46LqsNADh1Pq9H9TVRVZQGta+5Rx5Jt731Iv1tR5vWP79P3/3FNt3w2vGyZKk7nta5tVGFAu4YxzmUIb9hQ0PDmSqH62VNU/uaevTbzfu14R+HdN7YUr1zxoXyeg3FEhlNHFumCq40CQA4CSUhvyaNL1dLZ1zNHXFddn6lLhhXpp//aZfWbNijSePLNfO145XYl9G4mqiqSs/uwcbHM2S4ueWWWyRJBw4cGDDfMAyFQlxzZbjSGVO7DnTpV8/t0ZZX2nT5BVV6+xvOUyZrKpW2dGF9uUq4+BIA4BR4PR6NrSpReUlQjS098nos/dsNF2vTv1r0u82N+t4vt+nGayYoa1nqjqVUX1Pi2iumD6tt6rbbblNzc7Oi0agMw1B3d7e8Xq8qKyv1+OOP66qrrip0Oc9ayVRW/9rfoTV/2q1XDnTpDVeM1Vuuqlc8mZXHY+ii+jIuvQ4AOG3CQZ8uqC9Xe1dCB1t7dcUFY3RRfblW/WmXVv1plyZPrNCMV9erO55WWcSv8mhQJSGfq4LOsMLNtddeq9e97nW6+eabJUnr1q3Ts88+q1tvvVUPPfSQfvKTnxS0kGer3kRa/9zVrlV/sq9h87ZpE/TaybXq7k0rEvRpQl2p/D4GDgMATi+PYai6PKzSSEAHWmPKZEy9+4aL9fy2Zj39wgHtberRm6acowvGlepwLCVDUijoU2VpUCUhv0IB71ndbTWscLNt2zZ94QtfcKYbGhr0rW99S5deeqnS6fQQrxy98tew+dn/7VJPPKN3zrhQF59boa6elCpKg6qvKSmquxUDANwn6PfqvLGlOtyTVGNrr6ZeVKOL6sv1iw17tO4v9lnQ42tKNHlipS6qL1MylZVlST6focpoUKUlAYUDvrPuZs3DCjeZTEbbt2/XxRdfLEnavn27TNNUMplUJpMpaAHPRm1dCW3a1qyfP7tbkrSw4WKNqy5RVyytusqwaqsi3PsDAHBGGIahitKQSsIBHWqLKZM1tfDGi9XZk9K2PR3auqdTv920X7/dJI2tCmvyxEpdPL5cpmWppTMhw5DKowGVlwQVCfnOikuVDPsKxe9+97s1adIkmaapPXv2aOnSpXriiSd0/fXXF7qMZw3LstTU0auN/zikXz+3V6WRgG6/fpLKowH7fh81UY0pZyA2AODM8/s8OreuVBWlQe1v7lHI79PrLxurN105Th3dSTvo7O3Q0y8c0NMvHFB1eUiTJ1Zo8oQKeT1SR3dKkqWSkN/pvirWMaPDureUJHV2dmrTpk3yer266qqrVF5erp6eHkWjJ3a/hzNhJO4tZZqWDrT16A8vHNDTLzTqnDER3fbWSfL7PEqms5pYV6qyktF9qnex37dlNKEuigd1UTxGU11ksqY6u5Nq60oqmc7IYxgKBrzy+zzq7k1p295ObdvTod2HumVZUnlJQJMnVupVEytUVxlWKm1JsuT3e1VVGlQ0HFAo6D2tvRIFv7fU+vXrncemaeq5556TJN1www0n9GFulcma2tvUpbV/2adN21o0aXy55l93gbKmpUzW0kX1Fa65+SUA4Ozn83pUXRFWdUVYiVRGPb1ptXcn1R1LyTAMTb2oWq+dXKveRFrb9x3W1j0d2rStWc+91KRo2K/JEyo0eWKl6mtK1NwZ16H2XlWXhzWuumSkv5qkYYabH/3oR87jdDqtf/3rX7rmmmsIN5LSmaxeaTysNRv2aOueDl11cbVmTZuoeCKrQMDDzS8BAEUtFPApFPCpuiKsZCqr7nhK7V19Qeey86s0dVK1kqmsXt7fqa17OvX3V9q06V8tCge9uvjcCl1UX67SSPFcr+2Ew40k7d27V1/+8pcLUqCzSSKV0bbdHVr17C7tberRjKvq9YbL69QTT6s04te5tdwjCgBw9ggGvAoGwqouHxh0umIpGTJ0yYQKXX7BmNyBfZe27unQtj2d+vuONlnWBTr/nPKR/gqShhlujjRhwgTt3LnzdJflrBJLpLVlR6tW/WmX2g4nNfeN5+mKC8aoO27f/HLcmJKz7tQ5AADyBgs6Hd1JdfWmZVjShfVlmjyxUtmsqT3NPRpfJF1S0kmMubEsS//4xz8G3EBztOnsSeqv21u0+k+7lExnddvMi3Te2FJ196Z1zpiIairCZ/XFjwAA6G9A0Eln1RNPq70r4QSduopwUQ3BOOFuKcMwVFVVpccee6xghSpmrZ1x/fmlJv1y4x4F/B4tettkjSkLKRbPaEJdVJWlnOoNAHCvoN+roN+rMWWhAUHHX0TDME4o3GQyGVmWJb+/eAYNnUlNHb16+oVG/WbTfo0pC+r26ycpFPApmcpy80sAwKjTP+gM88oyZ8SwYlZbW5sWL16sqVOnasqUKVq4cKGampoKXbaiYlmWfvP8Pq37yz6dW1OiRW+brIDPK0vSReMJNgCA0a2YhmMMK9w88sgjuvLKK7VhwwZt2LBBV199tT7zmc8UuGjF5c//bNLv/9qoV02s0B03XKxs1pLf59WF48oUCoze8UcAABSbYe2Vd+/erccff9yZvvvuu3XTTTcVrFDF6IL6Ms2+9jxddl6FehNZlZcENL6Wm18CAFBshrVnzmQySiaTznQ8Hi+q5qczoa4yoqkXVSuWyKqmIqRz66IEGwAAitCwWm5mzZqlRYsWad68eZKkn/70p2poaChowYpRwOfRhLqoqsvDI10UAABwDMMKN3fddZfGjh2r//u//5Npmpo3b54WLFhQ6LIVnfG1xXeTUAAAMNCwws173vMe/fCHP9T8+fMLXR4AAIBTMqxBI93d3ert7S10WQAAAE7ZsFpuwuGwZsyYoUsuuUSRSMSZ/81vfrNgBQMAADgZwwo3o3F8DQAAODsNK9zccsstymaz6unpKarLKwMAABxpWOFm+fLleuyxx5ROpyXZtyIwDENbt24taOEAAABO1LDCzXe/+12tWLFCr3rVqwpdHgAAgFMyrLOlysvLCTYAAOCsMGS46ezsVGdnp6ZOnaof/OAHam1tdeZ1dnaeqTICAAAM25DdUtOmTZNhGM4g4scee8xZxpgbAABQjIYMN9u2bZMkmaYpzxE3iaTlBgAAFKNhjbkZ7LYLt99++2kvDAAAwKkasuXmPe95j1588UUlEgldddVVznzTNBlgDAAAitKQ4eZrX/uaOjs79f/+3//TF77whb4X+XyqqamRJPX29g64JQMAAMBIGjLcRKNRRaNR/dd//dcxn3PHHXfoZz/72WkvGAAAwMkY1piboXA7BgAAUExOOdwYhnE6ygEAAHBanHK4AQAAKCaEGwAA4CqEGwAA4ConHW4ymYwkcRo4AAAoKsMKN0899dSA6R07dugd73jHoMsAAABG0rDCzY9//GP98pe/lCR9//vf1+233665c+cWtGAAAAAnY8iL+OV973vf03vf+1798Ic/lMfj0YoVK3T++ecXumwAAAAnbMiWm87OTnV2dsrv92vp0qU6dOiQ7rzzTlVWVnJXcAAAUJSGbLmZNm2aDMMYcBXiO++8U5J98b6tW7cWtnQAAAAnaMhws23bNuexZVkyDEPZbFamacrv9xe8cAAAACdqWAOKn3vuOWcA8c6dO/XmN79ZL7zwQkELBgAAcDKGFW6++MUv6gtf+IIkadKkSfr2t7/tTA9lzZo1mjVrlmbOnDnoKeNbt27V/Pnz1dDQoAceeMC5dk7eSy+9pMsvv3w4RQQAAJA0zHCTTqd12WWXOdOXXXaZUqnUkK9pamrSsmXLtHz5cq1evVorVqzQjh07Bjzn3nvv1ZIlS7Ru3TpZlqWVK1c6y+LxuB555BGl0+kT+T4AAGCUG1a4CYfD+uMf/+hMb9y48bhXJt6wYYOmTZumiooKRSIRNTQ0aO3atc7yxsZGJRIJTZ06VZI0b968Acsfe+wxLVq06ES+CwAAwPCuc/PAAw/orrvuks9nP93j8ejJJ58c8jXNzc2qqalxpmtra7Vly5ZjLq+pqVFTU5Mk6Xe/+50SiYRuvPHG4X+TfsaMifZ739KTeg8UBvVRPKiL4kFdFA/qoricbH0MK9xceeWVevrpp7V9+3Z5vV6df/75CgQCQ76m/+njeYZhHHd5S0uLvvGNb+gHP/jBcIo2qLa2HpmmpZqaUrW0dJ/0++D0oj6KB3VRPKiL4kFdFJd8fXg8xoBGi+EYMtysXr1ac+fO1fe///0B8zds2CBJeu9733vM19bV1WnTpk3OdHNzs2prawcsb21tdaZbWlpUW1urp59+Wp2dnbrjjjucZXPnztVTTz2laPTEvhwAABh9hgw3e/bskSRt3779hN/42muv1ZNPPqn29naFw2GtX79ejz76qLO8vr5ewWBQmzdv1mte8xqtWrVK06dP1zve8Q7nppySdMkll2j16tUn/PkAAGB0GjLc3H333ZI0rNO+j1RXV6d77rlHCxcuVDqd1oIFCzRlyhQtXrxYd999t6644gotXbpUDz74oGKxmC699FItXLjw5L4FAABAjmENNvglZ86cOUO+eM2aNae9QKcDY26KE/VRPKiL4kFdFA/qorgUbMzNkiVLJEkHDx7UN7/5TT388MN65ZVXtHTpUj3yyCMnX2IAAIACGTLcXHPNNZKkRYsW6bbbbtM111yjK6+8UslkUqtWrdJNN910RgoJAAAwXMO6iF9HR4czHiYYDGrRokVqaWkpaMEAAABOxrDCTTabdS6wJ0mtra2DXqcGAABgpA3rIn6LFi3SzTffrDe96U0yDEMbNmzQpz71qUKXDQAA4IQNK9wsWLBAl19+uf785z/L6/Xq/e9/vy6++OJClw0AAOCEDSvcSNLkyZM1efLkQpYFAADglA1rzA0AAMDZgnADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABcpaDhZs2aNZo1a5Zmzpypp5566qjlW7du1fz589XQ0KAHHnhAmUxGkrR582bNnz9fc+fO1Xve8x41NjYWspgAAMBFChZumpqatGzZMi1fvlyrV6/WihUrtGPHjgHPuffee7VkyRKtW7dOlmVp5cqVzvzPfe5zWr16tebMmaPPfvazhSomAABwmYKFmw0bNmjatGmqqKhQJBJRQ0OD1q5d6yxvbGxUIpHQ1KlTJUnz5s3T2rVrlUql9LGPfUyTJ0+WJF1yySU6ePBgoYoJAABcxleoN25ublZNTY0zXVtbqy1bthxzeU1NjZqamhQIBDR37lxJkmma+upXv6rrr7/+hD57zJhov/ctPdmvgAKgPooHdVE8qIviQV0Ul5Otj4KFG8uyjppnGMawl6dSKd1///3KZDL60Ic+dEKf3dbWI9O0VFNTqpaW7hN6LQqH+ige1EXxoC6KB3VRXPL14fEYAxothqNg3VJ1dXVqbW11ppubm1VbW3vM5S0tLc7yWCymD3zgA8pkMvrGN74hv99fqGICAACXKVi4ufbaa7Vx40a1t7crHo9r/fr1mj59urO8vr5ewWBQmzdvliStWrXKWX7vvfdq4sSJevzxxxUIBApVRAAA4EIF65aqq6vTPffco4ULFyqdTmvBggWaMmWKFi9erLvvvltXXHGFli5dqgcffFCxWEyXXnqpFi5cqJdeekm/+93vdNFFF+nmm2+WZI/X+c53vlOoogIAABcxrMEGv5zlGHNTnKiP4kFdFA/qonhQF8WlKMfcAAAAjATCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcJWChps1a9Zo1qxZmjlzpp566qmjlm/dulXz589XQ0ODHnjgAWUyGUnSgQMHdMcdd+jGG2/Uhz/8YcVisUIWEwAAuEjBwk1TU5OWLVum5cuXa/Xq1VqxYoV27Ngx4Dn33nuvlixZonXr1smyLK1cuVKS9PDDD+v222/X2rVrdfnll+vrX/96oYoJAABcxleoN96wYYOmTZumiooKSVJDQ4PWrl2rj3zkI5KkxsZGJRIJTZ06VZI0b948PfHEE3rHO96h559/Xl/72tec+f/2b/+me++9d9if7fEYgz7GyKM+igd1UTyoi+JBXRQXj8c4qTopWLhpbm5WTU2NM11bW6stW7Ycc3lNTY2amprU0dGhaDQqn883YP6JqKwscR6PGRM92a+AAqA+igd1UTyoi+JBXRSXk62PgnVLWZZ11DzDMI67/HivAwAAGErBwk1dXZ1aW1ud6ebmZtXW1h5zeUtLi2pra1VVVaWenh5ls9kB8wEAAIajYOHm2muv1caNG9Xe3q54PK7169dr+vTpzvL6+noFg0Ft3rxZkrRq1SpNnz5dfr9fV199tX71q18NmA8AADAchjVYP9BpsmbNGn3rW99SOp3WggULtHjxYi1evFh33323rrjiCm3btk0PPvigYrGYLr30Un3hC19QIBBQY2Oj7r//frW1temcc87RV77yFZWXlxeqmAAAwEUKGm4AAADONK5QDAAAXIVwAwAAXIVwAwAAXIVwAwAAXMW14eZ4N+3EmbNw4ULddNNNmjt3rubOnau///3vI12kUaenp0ezZ8/W/v37Jdm3R5kzZ45uuOEGLVu2bIRLN7ocWRef/vSndcMNNzi/j9/85jcjXMLR46tf/apuuukm3XTTTfrSl74kid/GSBmsLk7pt2G50KFDh6wZM2ZYHR0dViwWs+bMmWO9/PLLI12sUck0TesNb3iDlU6nR7ooo9bf/vY3a/bs2dZll11m7du3z4rH49Z1111n7d2710qn09b73vc+6+mnnx7pYo4KR9aFZVnW7NmzraamphEu2ejz7LPPWu9617usZDJppVIpa+HChdaaNWv4bYyAwepi/fr1p/TbcGXLTf+bdkYiEeemnTjzdu7cKcMwtHjxYr397W/Xf//3f490kUadlStX6qGHHnKu9L1lyxZNnDhR5557rnw+n+bMmcPv4ww5si56e3t14MABLVmyRHPmzNETTzwh0zRHuJSjQ01Nje6//34FAgH5/X5deOGF2r17N7+NETBYXRw4cOCUfhuuDDeD3bTzRG++idOjq6tLr3/96/W1r31NP/jBD/TjH/9Yzz777EgXa1T53Oc+p6uvvtqZ5vcxco6si7a2Nk2bNk2f//zntXLlSm3atEn/+7//O4IlHD0mTZqkqVOnSpJ2796tX/3qVzIMg9/GCBisLt70pjed0m/DleHG4uabRePVr361vvSlLykSiaiqqkoLFizQM888M9LFGtX4fRSPc889V1/72tc0ZswYhcNhvfvd7+b3cYa9/PLLet/73qf77rtPEyZMOGo5v40zp39dXHDBBaf023BluDneTTtx5mzatEkbN250pi3Lks/nG8ESgd9H8fjXv/6ldevWOdP8Ps6szZs3a9GiRfrEJz6hW265hd/GCDqyLk71t+HKcHO8m3bizOnu7taXvvQlJZNJ9fT06Gc/+5lmzpw50sUa1a688krt2rVLe/bsUTab1S9+8Qt+HyPEsix9/vOf1+HDh5VOp7VixQp+H2fIwYMHddddd2np0qW66aabJPHbGCmD1cWp/jZceYhQV1ene+65RwsXLnRu2jllypSRLtaoNGPGDP3973/XzTffLNM0dfvtt+vVr371SBdrVAsGg3rsscf00Y9+VMlkUtddd51uvPHGkS7WqDR58mR98IMf1G233aZMJqMbbrhBs2fPHulijQrf/e53lUwm9dhjjznzbr31Vn4bI+BYdXEqvw1unAkAAFzFld1SAABg9CLcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcADht3ve+96m9vV2LFy/Wjh07BsyTpLe85S168cUXT/r9T/X1Qzmd5QQwslx5nRsAIyN/37DvfOc7R80rdmdLOQEcH+EGwGnx6U9/WpL0nve8Rzt27NDKlSu1fPlyZ963v/3tAc///e9/r2984xtKp9MKhUK67777TugCj8d6/ZNPPqnGxka1tLSosbFRVVVVWrZsmerq6rRlyxZ95jOfUTqd1oQJE3TgwAHdf//9WrVq1VHlXLFihR566CG1t7dr7ty5uueee07HnwnAmWABwGly8cUXW21tbdaMGTOsLVu2DJhnWZYzf9euXdbs2bOt9vZ2y7Isa/v27dYb3vAGKxaLDfn+w3n9E088Yb31rW+1uru7LcuyrA996EPW448/bqXTaWv69OnW008/bVmWZW3cuNG65JJLrD//+c+DlvORRx6xLMuympubrcsvv9w6cODA6fxTASggWm4AnHHPPvusmpubtWjRImeeYRjau3evJk+efEqvl6RrrrlG0WhUknTppZfq8OHD2r59uyTpuuuukyRNmzZNkyZNOuZn5C/1XlNTo+rqarW1temcc845ocqcVcQAAAFrSURBVO8JYGQQbgCccaZp6vWvf73+8z//05l38ODBYd+BeajX/+Y3v1EoFHLmG4Yhy7Lk9XplHXG3Ga/Xe8zP6H8H4vx7ADg7cLYUgNPG6/Uqk8kcd960adP07LPP6pVXXpEkPfPMM3r729+uZDI5rM85mddfeOGFCgQC+uMf/yhJ2rJli7Zv3y7DMI5ZTgBnJ1puAJw2M2fO1O23365YLHbUvK9//evOvEmTJumRRx7Rv//7v8uyLPl8Pn3jG99QJBIZ1ueczOt9Pp+efPJJPfTQQ/rKV76i8847T9XV1U4rz2DlBHB24q7gAEaNL37xi3r/+9+v6upqHTx4UHPnztVvf/tblZWVjXTRAJxGtNwAKBo///nP9d3vfnfQZXPmzNEHPvCBU3r/+vp6LVq0SD6fT5Zl6bOf/SzBBnAhWm4AAICrMKAYAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4yv8Po6oMjkDDWewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lineplot(data=causal_df, x='title_length', y='click_through_rate'); pp.ylim(0, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we see the bump in the middle goes away! We've recovered the true, underlying causal effect after removing the confounding by our author's title length preference. Go back to the article for more explanation!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}