From 513946e463979d368480eb554ee957043e2a274b Mon Sep 17 00:00:00 2001
From: Riku-Laine <28960190+Riku-Laine@users.noreply.github.com>
Date: Mon, 27 May 2019 17:14:23 +0300
Subject: [PATCH] Analysis with unobservables Z added

---
 .../derivation_validation.ipynb               | 268 ++++++++++++------
 1 file changed, 182 insertions(+), 86 deletions(-)

diff --git a/analysis_and_scripts/derivation_validation.ipynb b/analysis_and_scripts/derivation_validation.ipynb
index 46683e3..3dc1e21 100644
--- a/analysis_and_scripts/derivation_validation.ipynb
+++ b/analysis_and_scripts/derivation_validation.ipynb
@@ -2,12 +2,13 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 317,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Imports\n",
     "\n",
+    "import warnings\n",
     "import numpy as np\n",
     "import pandas as pd\n",
     "from datetime import datetime\n",
@@ -32,8 +33,6 @@
     "\n",
     "# Suppress deprecation warnings.\n",
     "\n",
-    "import warnings\n",
-    "\n",
     "\n",
     "def fxn():\n",
     "    warnings.warn(\"deprecated\", DeprecationWarning)\n",
@@ -42,27 +41,32 @@
     "with warnings.catch_warnings():\n",
     "    warnings.simplefilter(\"ignore\")\n",
     "    fxn()\n",
-    "    \n",
+    "\n",
+    "\n",
     "def sigmoid(x):\n",
     "    '''Return value of sigmoid function (inverse of logit) at x.'''\n",
     "\n",
-    "    return 1 / (1 + np.exp(-x))\n"
+    "    return 1 / (1 + np.exp(-x))"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the code below\n",
+    "In this notebook we try to validate the derivation of the following equation:\n",
+    "\n",
+    "\\begin{equation}\\label{int}\n",
+    "P(Y=0|do(R=r))=\\int_x P(Y=0|T=1, X=x)P(T=1|R=r, X=x)f_x(x)~dx\n",
+    "\\end{equation}\n",
     "\n",
-    "$$P(Y=0|do(R=r))=\\int_x P(Y=0|T=1, X=x)P(T=1|R=r, X=x)f_x(x)~dx$$\n",
+    "where $f_x(x)$ is the pdf of x. (Now $X \\sim U(0, 1)$, so $f_x(x)=1$.)\n",
     "\n",
-    "where $f_x(x)$ is the pdf of x. (Now $X \\sim U(0, 1)$, so $f_x(x)=1$.)"
+    "First we generate data and estimate the causal effect from a causal model without unobservables Z."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 320,
    "metadata": {
     "scrolled": false
    },
@@ -71,29 +75,52 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Bias: (0.07131297477417739, 7.917330654880471e-16)\n",
-      "Bias: (0.035957051207116175, 3.9920346147766154e-16)\n",
+      "0.0 % 10.0 % 20.0 % 30.0 % 40.0 % 50.0 % 60.0 % 70.0 % 80.0 % 90.0 % \n",
       "Analytical: 0.03709585053394618\n",
-      "Estimated:  0.035378945320598335\n",
-      "Difference: 0.001716905213347844\n",
+      "Estimated:  0.035394276460008964\n",
+      "Difference: 0.0017015740739372148\n",
       "\n",
       "Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\n",
       "\n",
       "Analytical: 0.14973849934787756 0.11264264881393138\n",
-      "Estimated:  0.1432745402973986 0.10789559497680028\n"
+      "Estimated:  0.14676199070433749 0.11136771424432852\n",
+      "\n",
+      "Average difference: 0.00045466276498828695\n"
      ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.lines.Line2D at 0x1a2103aa90>"
+      ]
+     },
+     "execution_count": 320,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0QAAAG1CAYAAADZSDS7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcFMX5x/HPwy4sICCgHCIGEFRAkBtUUPHCI8YzKh5RjHj9vBJjjEajJlFjTk003gcqGo3ReARvBTyRGwQEFUQ8QFaQ+2br90f1yDjM7M7szlJzfN+v17xmt7um55np7pp+uqqrzTmHiIiIiIhIMaoTOgAREREREZFQlBCJiIiIiEjRUkIkIiIiIiJFSwmRiIiIiIgULSVEIiIiIiJStJQQiYiIiIhI0VJCVETMbL6ZzQ8dh2zNzMrM7CYzm2tmG83MmVnPai5rWPT6YXHT2kfTRiQp/0Mze9/MVkRlbo2bd6aZTTez1dG8n1UnJpFCZWYjon2jfehYoPbjqawuyQdmdq+ZfWFmZaFjSaW2f6vNbIyZFe09V3Jtn80mM3vZzCaYmYWOJd8oIcozcT9Gz1RS5tiozPVZfN+CrUByxC+Bq4DPgD8CvwUW1fabmtmuwFPALsA90fu+FM0bBIwAyoDbonnjajsmyR9mtrOZ/dzMXjOzz81sg5l9aWaPmVm3Sl7XxcyeMrMlZrbGzCaZ2RlVvM8IM1tkZuvMbGb0vkl/w8ysqZn9I4ppfXSi4XdmVj8bn7vQFWp9b2Z7AmcBNznn1oeOp7aY2eBsHwPkk3xL2qNYx2RpcdcDfYGTs7S8olEaOgDZpg4OHYCkdASwCjjMObexFpb/JdAFWJ4w/WB8wnOZc+7xJDEBnOmcUyIkyVwM/Ar4GHgRWAp0A4YCx5vZ4c65MfEvMLPuwDtAPeAJ4GvgGOAhM+vgnPttQvmdgfFAa3zyPg84BPgbsBf+ADe+fGPgTaA7PrmfBvQHfgPsE8W0OUufH/yJjJvx+1gxSFWX5IPrgBXA/aEDCewMoGHoIAIq2H3WOfeemb2LT4wSf9OlEkqIiohzbm7oGCSlnYAltZQMES13dor3heStUZXNEwGfqOznnHs7fqKZnQj8G7gD6JrwmjuBRvjk/9Wo/HXA28BvzOxx59ycuPJ/BNoAw51z90flS4BngWFmNtI593pc+V/hk6EbnHO/iYvpTuB8YBhZPCB2zi0EFmZrebmukrokp5lZK+BYYEQhtw6lwzm3IHQMIRXBPvso8E8z298592boYPKGc06PPHoA7QEHPFNJmWOjMtcnTJ8PzE+Y1gy4Cf8Dtwb4FpiBP5BpFPc6l+QxImFZZwMTgNX4s3BvAsekiLEj/mzv8qjsK0APfBctB7SPKzssmjYMOA7fbWs1MCaavz1wJfAW/uB9A7AAuBtoneS9Y++xK/7gaS6wFpgKHB6VaQL8E19prgXeAPbIcF3tjz9D/W20jA+Ay4HSuDLXp/hux6Sx/O2AvwJfRcufDJwQ/30l2W5GJPyf7DEs1byE9z8QeAFYAqwDZkXroTShXKXrLypTBzgHeB/fUrYKeBc4vpL11wG4BPgIWB+tx0tTfFf18d0SJ0fLXoFvObgBqJtQthfwJL7lIrbcPxDtDwllT8IfyH8TrYPPgGeAQRlsJ1XuN/gz2w74cYpl/Dqaf0ptraNq1FVzomXuGDetczTt1STlj4nm/SFuWpNoHXyUpHyPqPy/4qYZfn9YDjRMKL8jvm54L834q6wbE7bHVHXWMcDEaBmfAb+Ii/WyaPtdB8wEjkoSx3wS6u3K5qWIJ+06kirqexLqkoTXVlnnReUGR8u4HugHvIbfL5fiD+ZaJFn2IcCrUfzrgC+Al4Gj01yfP4/ec0iSeW2A3+ET/PJom/sE+AvQOEn5MdGy6uK7En8WvWYmCftgVH534M/435hv49b3NSTUP8nWK/BA9H69Uny2h2PzSf2b4hLjT7KctOvJbD7wxwMPRut0Q/R8R3W2Ayr5/SLaJ6j9fbYP/vhhZvQdro6+0wsBS7IfJHsMjitXFq2XaVFMy/H7zAEpvs9W0TLur611VogPtRAVseiiu5fx/U1fAZ7Dd2PpiK8UbsJXirdG//cA/g4sixYxNW5ZtwEX4SuPu6PlnAg8Y2a/cM79La7sLvguMy2j95wF9MQfCH63zCROxnfxehZfoW+KpnfB/yi9AfwH/0PcAzgXGGJmvZ1z3yZZ3i34iut5fGvpacBzZjYQXxmX4puc2+OTzP+ZWWeXRncbMzsJeAxfeT2B/xE8Ev+jOMjMjnO+5hoTvSQ2WEFsQIP5VSy/DvA/fIU6Gf+D2CZ6z1erig+/Dn8bvf4A4KG495wazTuWrdd57P0vAv6BTwSejT7fIHzi0B84Psl7Jl1/0Xb4L3xyMSuKBeCHwFNm9jPn3N+TLO8v0Xv+D7/9ngjcambrnXN3xcXaEHgd2Dta/r34H4vOwBXRcpZFZY/Dr/MN+MRmEdAbfzB5YHTGbUNU9kLgdnzC9Dh+X2kTfZ+D8YlSpTLYb0biD3ROw2/jiU6N3v/ZuGVnbR1VU6y1M34ZB0TPybbR14CKuDIA++C/k9cSCzvnppnZ1wnld8e3bL7onFuTUP4bM5sM9Dez+s65dakCz6BurMrx+O/zv/jt4XjgL2a2Bl9vHY/ffusApwNPm1lX59wnaSw7E5nUkVXW98lkUOfF64ffB18D7gIG4rflXc1s31h5MzsKvw4W4bfNpfgulAOAo6N5VTkQv9+/n2Te/viE6XX8iRgXLfsXwP5mNtAlb71/HF8/jML/XpwCPGZmy5xzL8aVOx7ftfON6LOW4bfb3+O3sWOriP2e6PU/xXdT/Y6ZNcGfCJvinJtiZtvjf7POBMay5TemUpnUk9lkZvvgk+j6+PX4afSe5wOHmVk/59zSqGw628FU/HZ7KT6BiL/eOp34s7HPngMchT+mGQU0Bobgfy92Y8vv/Xz8fnkd/jdgRNwy5kefuT6+DtoPf+LsHvzJ0GOA183sJOfc0/EfwDn3dTQox0FpfF6JCZ2R6ZHZgy1n52bjD5CSPR4njRYifP97B/wtyfs0AerF/T+ChDMqcfMOiOZN4/tnTlvj++huBHaNm/5oVP7ChOVcQ8KZnGj6sGjaJnz3nMT33x5onmT66dHrrkmYHvssHwI7xE0/IZr+Lf4AvSRu3m3RvK1aLFJ8d8vwB0xd4qaX4g+yHHBGZesmjff4abScZ4A6cdMPjvsOhyXZbkYkLOd6Es5GVbXOgT2jdToO2D5uuuErfEdcS0Ya6++8aP4/E77z7fAHL+uBNkni+hhoFTe9UxTXnITl/y0qfydxZ+eiea2Izl7jWxBW4K9RaZNQ7pfRMi6PmzYZf3YysSXCkm2PWdhvxkXfRdOE5cRaSh6urXWU6QN/osEBExKm/4VK9iPgc2Bx3P8XReUvS1H+rWj+dtH/R0X//yNF+Uei+V2riL9GdWPc97meuLP6wM74M8vL8AedyeqffyS833yy00JUnTqyfZLXtCehLiHDOo/vnxk/IW56HfxBuQP2iZv+dPRdJmsx2CFxWorvqhz4JMW8lrFtKGF67Dfp9ITpY6Lp7/L9fTe2T7+cUH7n+G0mbl+MJR2DEuYlW68f4BOAsoTpsfrzwrhpse/3+hSfdwxbt/inVU9m84E/yfBZ9Lm6JMw7MYrn9ky3g2TbaBr7yDCyt8/+gLjf5bh94SVgM9AuYZ4jRWs8/gSWA65MmN4i2k7KgQZJXvef6HVtki1Xj60fGmUuf+2BP6uQ7JHp6CJrEyc451a46Ex4Gs6Mnn/jnPvurKlzbhH+7GCs9YVoqNMT8Ac+dycs56/4ijGV/zrn3koS63IXnUFK8Ci+afmQFMu7yTm3JH75+JaBpsAV7vstQU9Ezz0qiS/mWPwByD3OuQ/j4tyE76IHW76z6jo9ev61c64i7j1eJ70Wopo4D79OL3LOfXdhtfO1cKzr1tAkr0u6/vDdCL4Ffh7/nTvnVuPPoNYjeWvGjc65r+PKf4I/o7d7dGE9ZlYKDMe3klwRxUjca76O1gv4C40bR+W+Snivv+J/eBI/1wYSWlGcV9l2HJP2fhN5FP9d/DhhOafFzY/J9jpKm5k1YssBx68SZjeJnlekePkK/L6TSfn4cumW3z7F/EQ1rRtHOuemxL32S/w2uj2+a2Cy+iedOiYjNagj01XdOm+sc+6puPIV+NZu8C0n8TaypdXxOwnfYVLR786O+G6wW3HOLY7qm0R3RM+pvp9fJ+y7Y/EHqd+L3Tn3ZeI2E+2Ld1ax/Hj34rtxHpcw/af4A/bH0lhGUhnWk9l0FD55uCl+u4ne80lgElvXU9XeDtJU433WObcg/nc5mrYJ37pTB99aWaWoJ8j5wEzn3M0JyyvHn2DakeQDZi2OnndO571Egyrks2edc0mb2c3sWPyOWpVZ+D7xV5m/580o/BnXGYkVYhVi98sZk2TemIQynfHdBcYnVrDOubVmNpXUzbwTUwVgZgfjm6H7AzsAJXGzd0r6In9mPv79K8ysHH/G//OEsrGBBdqkiiFOyu/DOTfVzJbHlamuHsBS59ysJPPeBg6t4fIrMwB/sHt01IUh0Vr8ek601fqLuml0w1/P8Gvb+tYJLaLnZMubkmRabNSgpsDK6HWN8WffViYpH29A9DwoxZDRGxPieAI/UtEMM3sC3z3lvRQHVslkst+Ab/n9Gz4Bug++69p1Cn77jO9WlrV1lAkzq4sfTKEbcJ1z7o3EItFzuvVLbZdPJVt147Qk0xYlmxdX/6RTx2SsmnVkuqpb51W1D8c8gU8EZpjZ49H7vO2cS7cL1w7Rc7Ku08B3A4GcF8XZjO/fliTV95Mq/h8kLLsO/lrBYfjW2yZs2VYrW368R/ADjJxNNHpYVE/1Bx5zybuFpyuTejIp8/e9a58weYRzbn4lL4vVud1SDBHeANjBzHZ0zn1DzbeDdNR4n40S8EvwJ6f3wA8iEy/d/W0P/H7wWYrvZ7fouTO+G1+82PawY5rvVfSUEBUx59wmMzsIfzHp8fj+3gCfm9kNzrl70lxUE2Cdcy7ZWdlFcWXAV7rgz7YnszjF9JTzor7rj+MPgF/Cn6GLndn9GT4BSyZZvJsqmQ7+ItqqxD5r0rOR+O+kYxrLqcz2+C5jyaR632xpjv8x/00lZbZLMi3Z+msWLasdvnUzk+UlG/Y3tp5iB3uxloDEFp9kmkfPl6ZRFuBP+B+dC/Bda64B1kU/1JelcYCSyX6Dc67czF7F96vfOTpzeQDQFrg1oUUzm+soLdFZ5n/hh2v/q3Pud0mKxdZZqhaaJnx/vaZTHrbss+mWr3TI6CzWjZXVJanmpVPHZKQGdWS6qlvnpbMP45x7wsw24S9o/zn+2p5NZvY88DNX9ahpsc/aINlMM/slfn9ejB8+/kt8qwv4einp9xPf+poQf2Lvm9uA/8N3D3uaLQNbNMXXN1V+/865b83sKeBUM2vnnPsMnxxBdIKkBjKpJ1MZxvev5wOfsMyv5DWxOreqHhPbAd9kYTtIRzb22afw18DOxrfclUfl2uM/a7r7W+z76UHlLcfJ6vLYtr4myTxJQglRkYuaXS+ILhDvhm9Z+Blwt5mVO+fSaWlaAXQ0syZJDu5axZUB/4MMW878J2pZWbgppl+H/8Hr7eKGFo/Onl9RyfJqS+yztkoxvxWpu/Skazmpv6tU75stK/D9oLdzmQ1fm2z9xb6Hd5xzg2oc2dZiZw7TOesei2U3l8ZF7VFLwT3APdGQvgfgu68Mw/+QHZPG+6W738Q8ik84TsF3l0jWXS72umytoypFw2CPxHeHvc05d3mKorEkvlOSZWyHX0/vp1M+shuwMK5VLp3yFfjrxCqVpboxGypInSQ1Ib26pLbryFqv86KudU+ZWTP8Bean4LtT7WpmvapouVuGb+FtnjgjSuSvwScDPaKWiNi8VlR+oqZK0TIuwLcu7OOcWxs3bwDpn4AB323uNOAsM7sJ33V6HmkOnFCJTOrJpJxzg6vxstg2cahzbquBU1K8T022g1pnZv3wydBLwA/ju86Z2clk1l0+9v086pw7vdKSW2sWPac6+SwJdA2RAL7p1zk33Tn3V3wFA37UlpjY2ecSthYbfWj/JPMOSCgzG3/RYv/oh+g7ZtaA6vWf7wh86La+z1IvUpwRrGUpvw8z2wt/VrDSEZvSMA1obmaJ93gBP5JYbRqP3w761XRBUfeM2fguE4ndCrJhDj4J3yeN5Y+PnvfO9E2iPvb/xrckfAwckbh9J5HJfhPzDH4I11Ojbhk/xg8ikdjVLWvrqCpRMvQIvnvIXc65SyopPjZ6Ttal8xD8b1L8fTPG4c+kb3WNhZn1wB9ox5f/CD9U/sCoPokvvyN+RLAJrpIR5hKlUTfWtmVAq+h7/o6ZtWPLQU9VMq0jK6vvk9kWdR7gW0qcc885507BD8DQg4Quakle4/BDIHeMuq/F2xGfWL4XnwxFBmYh5A741trX4pOh6iw/ukZpDn7EuWPxsT+QJAnIdP1lUk9mU03q3Mq2g0w/fzbFWkJHJV5HROr1XUHyWD/Er5d+ift/GnbHnwTJ9oiVBUsJUREzsw5mluw6gthZvvjKO3ZBbrIL9GIXwf42uiYktvyW+LOPm4gu+IzOVv8X2AXfXzveZWzp652JBcBu0fvF3rsJfujNEJ7Fn9k518y+O1MdVWh/jP59ONkLMzAyer4x/gc+uk6gNq8fAn+h8Wb8jd+26gttZq3MrEsGy7sN32XjjmiI0cTl7Rm/bjMRXad2L/7A4c+WcJGSmbWMS1wexI+S9Ucz240EZtbUzHrF/T8kyY9UQ3x/8Q34H7nKpL3fxH2e1fikqBe+u0hTtmwL8bK9jpKKtr0H8YnC/fhuQSk552bjR+Y62My+206j5OXaKOYH48ovx4+WtJuZnR1XvgR/bxSI6y4UHRg+iD/A/XXC2/8O39JSZfeiDOvG2jYJH3csGYtdq/WXDJaRaR1ZWX2fTK3WeWZ2UHQCIH5aKVtafNJJcN/Cd9lOXK+LiVrP4pPoaL+5qboxx4l149onvv4xs92Bq6qxvPvwB/634PeXEUnKZLT+Mqwns+kZ/ABLV5pZ/8SZZtYgakWL/Z/udhDrrhxiQIHY+v5e8mNme+OHuU9mKUlijdbLXfjk5oZkSZGZDYj//YimleJ/I8YlDuYhqanLXHHrAfzXzMbhz54txp/NOhbf7/SuuLKj8TfYu8vM/hPN/8A5N8o5N8b8XeAvwF/s+F/8aFgn4bt1XZFwZvJK/Kgot5nZIfizID3xFchb+Gbwqg4m492O/2GfHPWxLsN3K/qKmvWJrhbn3HIzOx9/oDrJ/DUly/CtB93wFz/WNCEaAfwEv64mmL+2pA3+LP0ofJN9rXDOfWBmF+O/94/M7AV8P/Fm+C5Jg/DXrnyYciHfdyewL/7zDDazN/B97HfCD3/cE38/mupe3/Ib/LZ1PrCfmb2C3752Bw7DH+Quc84tNrPT8Bfuzow+18f4/tm7suV+TedHy/03sMrM3sZfG9AQ/73vBNyQ5Ozg91Rjv4l5FN9t5vro/61Gl6qFdZTKtfj1tgy/r11nWw+MkXhh9fn4+5A9H+0bi/HdC3cHfuucm5Pw+ivwwwjfY2aH4bsIHYIf2vuhJF1tbsa34FxjZn3wrakD8CM7vU5cwlWJTOrG2vZPfDfMB6Ikchm+/lyJbw1LR6Z1ZMr6PtnCt0Gd9zegrZmNwW/HJfhtoBt+VLB0rpt8Fn8Pn0Pwg2bEYq8ws7vw16RMMbNR+APs2H1k9qhB3Djnvor27ePwdfVofF19NH5I8hMyXORDwI3RMkZF1xImmoPfNoaa2WqigSpcwkhlCdKqJzOMtVLOufXmB7N4ERgXvecs/LFpe3yd+x5wePSStLYD59wqM5sAHGBm9+HvFeeAO1Nc95VN7+MHpxlqZq3x9w7alS33SUq2vkcDJ5ofnGc6PtF9LLom6lr8qIVXAseZ2Vv4BKptND1277X4a4X2xbf8Poukz+XA2N96pP9gy/j6z1RS5liS3IOAre9D1BZ/8PA+/gd/Hf5g4yES7gkQlb8KX7FsZOv7UBj+ZmST8DvmKnxyk+p+I53wLUUroscr+IOQ56NlN4srOyyaNizFsgw/dPMs/Jm+z/E//o0TP3NUfgSp77GxVfmE731EshhSxDUY/4O3jC13tf4VadyZPM3lN8L/QCyMPvcUfGW71feVKn6qcR+iuPn7AE9G778Bn8SMw1fgP0h3/cWVOw3/w/Atvlvlguj7u4C4e4RUsf6SzsP/OFyFv5fHWvw1WFOJWg0SynaNlvN59Lm+wd9z6Gagc1y5C6Lt9bNo/S7GH0CdnME6zGi/iV5Tir943eGvvaps+VldR5V835U9km1bXfEXly+N1sdk4MxK3qctvl5aHG0bH+JblEtSlG+Gb3n8Iio/Dz+Ee/00P1fadWOyba6y77OK7Xc+yeufIfiDrPXRur+dDOo3Mqwjo9ckre+ppC4kzTqPSu6Tk2we/kTPv6N41gBL8N2tziPN++NE38HH+BFOE+fVw+8Tn0Rxz8XXjfVIco8YktzHp7J50fd8K1vqipn4Ucg6JPsuU62TuPnPRa87rpIy++LrklVRWVdV/GRQT2bzgW/xit3kej3+N+AD/D7crzrbAb4l8OXoM8TqovbbYp/FJ48j8Ccb1uDr99OSbdtR+Tb4lvAl+CT0e/Umvs6/EF93r4jWzTx8C9sZST773dF2VuX98PTY8rDoyxMJLup+Mw8/7HW1ukiJiIgkY36AjNvxgydMDx1PdUS/k3PxLdJtnXNb3ZNHipf5e/B9BjztnBseOp58omuIJIhk1zXgu8a0Q828IiKSfffgW4GuDR1IDRyHb6W7X8mQJHERUJ8tXaolTbqGSEL50MzG47tw1MHfXG4AvjvI9QHjEhGRAuSc22j+BqIHm1mZy2xI+qDM7AJ817Jz8d2mQg0aJLltJb778RehA8k36jInQZjZH/EXoO+Cb/pfhO/v+zuXnZuriYiIFAQzm48fiexD/I2f07pvj4ikRwmRiIiIiIgUrbzsMrfjjju69u3bhw5DRKRwzYlGv96jRqMOi4iIBDNp0qRvnHMtqiqXlwlR+/btmTgx8cbsIiKSNYMH++cxY0JGISIiUm1m9lk65TTKnIiIiIiIFC0lRCIiIiIiUrSUEImIiIiISNFSQiQiIiIiIkVLCZGIiIiIiBQtJUQiIiIiIlK0lBCJiIiIiEjRUkIkIiIiIiJFSwmRiIiIiIgULSVEIiIiIiJStJQQiYiIiIhI0VJCJCIiIiIiRUsJkYiIiIiIFC0lRCIiIiIiUrSUEImIiIiISNEqDR2AiIjkkdVLYOFUMIMWXaBxa/+3iIhInlJCJCIiyTkHn7wOX03xSdBXU2H5598vU78ptOwKLTtDqz2h+4lQf/sw8YqIiFSDOedCx5Cxvn37uokTJ4YOQ0SkMDkH/bvDsvlwWomf1nxX2KkntOnpn81g8YdbHuUfwrrl0HAHOOBX0PenUFI36McQEZHiZmaTnHN9qyqnFiIREdli3lh47Toonw31GsJJD0OHA6BB063Ldth/y9/O+ZakV6+FF6+A9++GQ38LnY9SlzoREclpGlRBRERg0QfwyPHw8NGwqhx23B3a9IKuxyRPhhKZwc694czn4dR/Q51SeOJ0ePAI+HJy7ccvIiJSTUqIRESK3aSH4J7B8NVkGHIjXDwJGrUEqtGyYwa7HwYXvAtH3QJL5sIDh8HMZ7IdtYiISFYoIRIRKVYVm+Glq+D5S3y3uEumwL4XQd36NV92Sam/juii8dCmNzw5DCY+UPPlioiIZFlaCZGZtTWz28zsPTNbY2bOzNpX8ZqronJvJ5lXJ5o/38zWmdk0Mzuheh9BREQytm4FPHYyjLsDBlzgu7k1aJb992nQDH7yX9htCPzv5zD2z/56IxERkRyRbgtRJ+Ak4FvgraoKm9muwNXA4hRFfg9cD9wOHAGMA540syPTjEdERKpr6adw/6Ewb7Tv1nbEzb5Fp7bUawhDH4W9hsLoG+ClK6GiovbeT0REJAPp/gK+6ZxrBWBmw4EhVZS/E3gU2CPxPcysJXA5cLNz7i/R5NFm1gm4GXghzZhERCRTn4/3LUOuwrfcxI8UV5tK6sKxd0LD5r5Vas1SOPYODc0tIiLBpdVC5JxL+1SemZ0K9AauSlHkMKAeMDJh+kigu5l1SPe9REQkA4tmwMgf+25s57yx7ZKhmDp14LCb4KDfwAf/hucuVvc5EREJLqt9JMysGXALcIVzbqklv/fEnsB64JOE6TOj567Ap9mMS0Sk6H07H0YeD2WN4IxnoekuYeIwg/0vh4pNMOYP0Lo77HNhmFhERETI/ihzfwY+AkZUUqY5sMy5rU4LLo2bvxUzO9fMJprZxPLy8hoHKiJSNFaVwyPHwab1cPrT4ZKhePtf4W/a+so1MPeN0NGIiEgRy1pCZGb7AWcAFyRJdr5XFEg2v9IbXjjn7nHO9XXO9W3RokUNIhURKSLrVsCjJ8CKhXDak9Cyc+iIvDp14Li7oUVnePIsf78iERGRALLZQnQ3cD/whZk1NbOm+C55JdH/ZVG5pUAz27o/XbO4+SIiUlOb1sMTp/lrh056GHbpHzqi7ytrBEMf893oHj8V1q8MHZGIiBShbCZEXYDz8UNzxx4Dgb2jvy+Iys0EyoCOCa/vGj3PymJMIiLFqaICnj4HPn3Tj+a2e1WDgwbSvAOcOAK++RiePk/DcYuIyDaXzYTowCSPacCM6O//ROVeAjYApyW8/nRghnNOAyqIiNTUW3+FWc/CkBugx9DQ0VRu18F+9Lk5o/xACyIiIttQ2qPMmdmPoz/7RM9HmFk5UO6cG+ucG5PkNcuA0vh5zrnFZnYLcJWZrQQmAycDBwHHVOtTiIjIFp++CWNugu4nwj4XhY4mPQPOg68/gDf/BO32hY4Hho5IRESKRCbDbj+Z8P8d0fNYYHCG73s1sAq4FGgNzAFOcs49n+FyREQk3qrF8NRwaN4RjrrVX5+TD8zgyL/Agvfh2Yvg/96F+tuHjkpERIpA2l3mnHOW4jG4ktcMds7xjdJ9AAAgAElEQVQNSjJ9s3PuBudcO+dcmXNuL+fcf5ItQ0RE0lSx2V83tG65vy6nrFHoiDJTtwEcdxes/ApeSnVvbxERkezK9n2IREQklDf/AvPGwBF/gtbdQkdTPW37wqDLYOqjMPuF0NGIiEgRUEIkIlII5o31AxLsdTL0PiN0NDVzwK+gVXd4/hJYvSR0NCIiUuCUEImI5LuVX/vrhnboBD/8W/5cN5RKaT3fdW7tMhj1c6j0Xt8iIiI1o4RIRCSfOQfPnA/rV8BJD+XfdUOptO4GB/7aDx3+gS4xFRGR2qOESEQkn019DOa+AYf+HlrtGTqa7Nr3EmjbD174Baz4KnQ0IiJSoJQQiYjkq1WL4eVfwy57Q7/hoaPJvpJSOO5u2LQB/qeucyIiUjuUEImI5KsXfgkb18DRt0GdAq3Od+jou8599BLMHhU6GhERKUAF+gsqIlLgZo+CWc/AAVdAi91DR1O79r4AWu4JL14B61eFjkZERAqMEiIRkXyzbjmM+gW06gYDfxY6mtpXUheOugVWfOmHFhcREckiJUQiIvnm1Wth1de+q1xJ3dDRbBs/GAC9z4Rxd8KiGaGjERGRAqKESEQkn8x/GyaNgL3/D3buHTqabeuQ66FBUz/AQkVF6GhERKRAKCESEckXG9fCcxdDs/Zw4NWho9n2GjaHITfAF+NhysOhoxERkQKhhEhEJF+883dYOg9+9Heo1zB0NGH0OAXaDYJXr4NV5aGjERGRAqCESEQkHyz/Et6+FfY8DnYdHDqacMzgqL/BhtXw6m9CRyMiIgVACZGISD54/bfgKuCQ34aOJLwWe8DAS2Dav2D+O6GjERGRPKeESEQk130xCaY/AftcCM3ahY4mN+x3OWy/C7z4K6jYHDoaERHJY0qIRERymXPw0pWwXUvY77LQ0eSOeg3h0N/B1x/4UfdERESqSQmRiEgum/GUH1Xt4GuhrHHoaHLLnsf5ARbeuAHWfhs6GhERyVNKiEREctXGtX40tdbdoeepoaPJPWZwxM2wbhmM/kPoaEREJE8pIRIRyVXv3Q4rvoDD/gB1SkJHk5tad4c+Z8GE++DrWaGjERGRPKSESEQkF61YCG/dAl1+BB32Cx1NbjvoGt+d8KUr/TVXIiIiGVBCJCKSi974PVRs9AMHSOUaNocDr4ZPx8KHz4eORkRE8owSIhGRXLP4Q5j6GPQ/F5rvGjqa/ND3p9CyK7xytb/2SkREJE1KiEREcs2YP0C9RrDfL0JHkj9KSuHwm2HZAnj39tDRiIhIHlFCJCKSSxZOh1nPwt4X+K5gkr5dD4AuR8Pbf/PXYImIiKRBCZGISC4ZfRPU3x72uTB0JPnp0N9BxSZ/byIREZE0KCESEckVX0yCj16EfS6GBk1DR5OfmneAAefD1Efhq6mhoxERkTyghEhEJFeMvhEaNIe9zw8dSX7b7xe+u+Er12gYbhERqZISIhGRXLBgHMx9HQb9zN9TR6qvQVMYfBXMfwvmvBA6GhERyXFKiEREcsEbN8B2LaHfOaEjKQx9zoId9/CtRJs2hI5GRERymBIiEZHQ5o31rRn7XQb1GoaOpjCUlMJhN8LSeTDhvtDRiIhIDlNCJCISknP+2qHGbXyrhmRPp0Og40Ew9o+wZmnoaEREJEcpIRIRCemT1+Hz92H/X0Dd+qGjKSxmMORGWL/CJ0UiIiJJKCESEQnpzT/D9rtArzNCR1KYWnWF3mf6bnPffBw6GhERyUFKiEREQvnsPfh8HOx7MZTWCx1N4TrwaihtAK9dHzoSERHJQUqIRERCefsWaLgD9PpJ6EgKW6MWMPBSmP0/P7y5iIhIHCVEIiIhLJoBH78MAy7QyHLbwj7/B41awyu/0c1aRUTke5QQiYiE8M6tUK8R9B8eOpLiUG87OPAq+GK8bykSERGJpJUQmVlbM7vNzN4zszVm5sysfUKZvmZ2j5nNjsosMLNHzaxDkuXVMbOrzGy+ma0zs2lmdkJ2PpKISI5b+inMeAr6DIMGzUJHUzx6nu5v1vra9bB5Y+hoREQkR6TbQtQJOAn4FngrRZmhwJ7AP4AjgCuB3sBEM9sloezvgeuB26Oy44AnzezITIIXEclL790OdUphn4tCR1JcSkrhkOthyScw+eHQ0YiISI4oTbPcm865VgBmNhwYkqTMH51z5fETzOwd4FPgHODaaFpL4HLgZufcX6Kio82sE3Az8ELGn0JEJF+sWgxTRkKPodBkp9DRFJ89joAf7Atj/gB7nQRljUNHJCIigaXVQuScq0ijTHmSaZ8B5cDOcZMPA+oBIxOKjwS6J+tiJyJSMMbdCZvWw76Xho6kOJnBkN/D6nJ49/bQ0YiISA6o1UEVzKwL0BL4MG7ynsB64JOE4jOj5661GZOISDDrlvsbhHY9BnbsFDqa4tW2r18H794GK78OHY2IiARWawmRmZUCd+FbiO6Pm9UcWObcVuOeLo2bn2x555rZRDObWF6+VWOUiEjum/ggrF8Bg34WOhI5+DrYvB7G3hw6EhERCaw2W4huB/YFTnfOfRs33YBkN4GwyhbmnLvHOdfXOde3RYsWWQxTRGQb2LgOxt0Bux4IbXqFjkZ26Ah9zoJJD8GSuaGjERGRgGolITKzPwDnAj91zr2SMHsp0MzMEhOgZnHzRUQKy4z/wKqv1TqUS/b/JZSWweibQkciIiIBZT0hMrOr8UNuX+qceyRJkZlAGdAxYXrs2qFZ2Y5JRCQo5/xgCi33hA4HhI5GYhq3ggHn+XtCLZoROhoREQkkqwmRmV0C3ABc7Zy7LUWxl4ANwGkJ008HZjjnPs1mTCIiwc1/C76eAXtf4Ec5k9yx7yVQ1gRG3xg6EhERCSTd+xBhZj+O/uwTPR9hZuVAuXNurJkNBW7FJzxvmNnecS9f4ZybBeCcW2xmtwBXmdlKYDJwMnAQcEzNPo6ISA4adyc03AG6nxg6EknUsDnsezGMvgG+mOhHoBMRkaKSdkIEPJnw/x3R81hgMHA4fmCEw6NHvFiZmKuBVcClQGtgDnCSc+75DOIREcl9S+bCnBf99Sp164eORpLZ+3x4/y54/Xdw5nOhoxERkW0s7S5zzjlL8RgczR9WVZm4ZW12zt3gnGvnnCtzzu3lnPtPdj+aiEgOeP9uqFMK/c4OHYmkUtYY9rsMPh0L88aGjkZERLaxWr0xq4hIUVu7DKaMhG4nQOPWoaORyvQ9G5rsDG/83g+CISIiRUMJkYhIbZkyEjau9oMpSG6rW993a/xiAnz0cuhoRERkG1JCJCJSGzZv8t3l2g2ENj1DRyPp6HU6NOvgW4kqKkJHIyIi24gSIhGR2jDnBVi+QK1D+aSkLhz4az9E+synQ0cjIiLbiBIiEZHaMO5OaPoD2OPI0JFIJrqdAC26wNg/qZVIRKRIKCESEcm2r6bAgndhwPlQpyR0NJKJOiVwwC/hmznw4bOhoxERkW1ACZGISLaNuwvqNfLXpEj+6Xos7Lg7jP2zWolERIqAEiIRkWxa/Y2//qTHKVB/+9DRSHXUKYH9LofFM2HOqNDRiIhILVNCJCKSTVMegc0boN/w0JFITXQ7AZrv6q8l0n2JREQKmhIiEZFsqdgMEx+A9vtBy86ho5GaKCn1rUSLpuu+RCIiBU4JkYhItnzyGixboNahQrHXSX6kwLF/VCuRiEgBU0IkIpIt4++FRq2h8w9DRyLZUFIX9vsFfDUZ5r4eOhoREaklSohERLJh6TzfQtRnmD+QlsLQ41Ro0hbGqJVIRKRQKSESEcmGiQ+C1YE+Z4aORLKptB7s93P4Yjx8OjZ0NCIiUguUEImI1NTGtX50uc4/hCZtQkcj2dbzdGi8k78vkYiIFBwlRCIiNTXzGVj7rQZTKFR168PAn8Fnb8Nn74aORkREskwJkYhITU24D3bcHTrsHzoSqS29z4AGzeGdf4SOREREskwJkYhITXw1Bb6c6FuHzEJHI7WlXkMYcB589CIsnh06GhERySIlRCIiNTHhPqjbEHoMDR2J1LZ+50BpA3j3ttCRiIhIFikhEhGprrXfwgf/8TfwrL996Giktm23A/Q6HaY/ASsWho5GRESyRAmRiEh1TXscNq2DvmeHjkS2lX0uBLcZ3r8zdCQiIpIlSohERKrDOZg0AnbuAzvtFToa2Vaad4Cux/r7Tq1bEToaERHJAiVEIiLV8fn7UD4b+gwLHYlsawMvgfUrfEIsIiJ5TwmRiEh1TBoB9RrDnseHjkS2tTa9/BDr4+6ATRtCRyMiIjWkhEhEJFNrv4WZ/4W9ToSyRqGjkRAGXgorF8IHT4aOREREakgJkYhIpqY/6QdTUHe54tXxYGjVzQ/BXVEROhoREakBJUQiIpmIDabQphfs1CN0NBKKGex7CZR/CJ+8GjoaERGpASVEIiKZ+GIiLJ6p1iGBbsdDk7bwzt9DRyIiIjWghEhEJBOTRkDd7aDbCaEjkdBK6sKA8+Czd2Dh9NDRiIhINSkhEhFJ17rlMOMp6P5jKGscOhrJBb1/AnUbwvt3hY5ERESqSQmRiEi6pv8bNq1VdznZokEz6HGKH21uVXnoaEREpBqUEImIpMM5mPQQtN7LD6ggEjPgPNi8QTdqFRHJU0qIRETS8dVk+PoD3zpkFjoaySUt9oCOB8GE+3SjVhGRPKSESEQkHZNG+GtFup8YOhLJRQMugFWLYNazoSMREZEMKSESEanK+lXwwVN+mOX6TUJHI7mo0yGwQyd4/87QkYiISIaUEImIVGXWM7BxNfQ6I3Qkkqvq1IH+58GXk/y9qkREJG8oIRIRqcqUkbDDbrBL/9CRSC7reQqUNYFxaiUSEcknaSVEZtbWzG4zs/fMbI2ZOTNrn6RcfTP7s5ktNLO1Ufn9k5SrY2ZXmdl8M1tnZtPMTHc5FJHc880nsOA96HW6BlOQypU19tvJrGdgxVehoxERkTSl20LUCTgJ+BZ4q5Jy9wPnANcCRwELgZfNrGdCud8D1wO3A0cA44AnzezItCMXEdkWpo4EK4EeQ0NHIvmg/7lQsRkm3B86EhERSVO6CdGbzrlWzrkjgSeTFTCzHsCpwM+dc/c6517HJ1ELgN/FlWsJXA7c7Jz7i3NutHPuPGA0cHMNPouISHZt3gTTHofdDoXGrUNHI/mgeQfY4wiY9CBsXBc6GhERSUNaCZFzriKNYkcDG4En4l63CXgcOMzMyqLJhwH1gJEJrx8JdDezDunEJCJS6+a+ASsX+m5QIukacD6sWQIz/hM6EhERSUM2B1XYE/jUObcmYfpMfALUKa7ceuCTJOUAumYxJhGR6pvyCDTcEXY7LHQkkk867A8tOsP4e8G50NGIiEgVspkQNcdfY5Roadz82PMy57b6lUgs9z1mdq6ZTTSzieXl5TUOVkSkUquXwJwXYa+TobRe6Ggkn5hBv+GwcCp8OTl0NCIiUoVsJkQGJDsVljgsU7rlvsc5d49zrq9zrm+LFi2qGaKISJo++DdUbFR3OamevU6Geo1gwr2hIxERkSpkMyFaSvLWnWZx82PPzcy2Gr82sZyISBjOweRHoE1vaKVevFIN9Zv4pGjG0761UUREclY2E6KZQAcza5gwvSuwgS3XDM0EyoCOScoBzMpiTCIimVs4FRbPVOuQ1Ez/c2Dzen8tmoiI5KxsJkTPAXWBE2MTzKwUOBl4xTm3Ppr8Ej5BOi3h9acDM5xzn2YxJhGRzE0ZCaX1oZvuFy010LILtBsEE+/39yYSEZGcVJpuQTP7cfRnn+j5CDMrB8qdc2Odc1PN7AngVjOrC3wKXAB0IC75cc4tNrNbgKvMbCUwGZ80HQQcU+NPJCJSExvXwQdPQpcfQYOmoaORfNd/ODw5DD5+FfY4PHQ0IiKSRNoJEVvfkPWO6HksMDj6+yzgRuAGoCkwDTjcOZc4zM7VwCrgUqA1MAc4yTn3fAbxiIhk3+z/wbrl6i4n2dH5KGjUGibcp4RIRCRHpZ0QOecqHQUuKrMWuCx6VFZuMz5puiHd9xcR2SamjITtfwDt9w8diRSCkrrQZxiM/SMsnQfNdw0dkYiIJMjmNUQiIvlt2QKYNwZ6ngp1VD1KlvQZBlYHJtwfOhIREUlCv/giIjFT/wU4nxCJZEuTnaDLUb71ccOa0NGIiEgCJUQiIgAVFTD1UehwADRrFzoaKTT9zoF1y2Dm06EjERGRBEqIREQAPnsbln0GvX4SOhIpRO0HQYvOMP5ef+NfERHJGUqIRETAd2cq2953bRLJNjPoN9zf9PfLSaGjERGROEqIRETWLYdZz0L3E6Bug9DRSKHa62So18gPwS0iIjlDCZGIyIynYdM63XtIalf9JtBjqN/eVi8JHY2IiESUEImITBkJLbtCm96hI5FC1284bF4PUx4OHYmIiESUEIlIcVv8IXw5EXqe5q/zEKlNLbtAu0Ew8QGo2Bw6GhERQQmRiBS7KSOhTqm/vkNkW+g/3N8E+ONXQ0ciIiIoIRKRYrZ5I0x/AnY/HBq1CB2NFIvOR0Gj1jDh3tCRiIgISohEpJh9/AqsLte9h2TbKqkLfYbBJ6/BkrmhoxERKXpKiESkeE0ZCY1aQadDQkcixabPMN9Vc+IDoSMRESl6SohEpDit/Bo+etkPg1xSGjoaKTZNdvJd56aMhA1rQkcjIlLUlBCJSHGa/gS4zdBT9x6SQPqfA+uWwYynQkciIlLUlBCJSPFxzp+Z32UAtNg9dDRSrNoNhBZd/OAKzoWORkSkaCkhEpHi88VE+GYO9FLrkARkBv3OhoXT4MtJoaMRESlaSohEpPhMeQTqNoQ9jwsdiRS7HkOhXiMYryG4RURCUUIkIsVlwxqY8TR0PRbKGoeORopdWWOfFM18GlYvCR2NiEhRUkIkIsXlw+dgw0p1l5Pc0W84bN4AUx4OHYmISFFSQiQixWXKSGi+K7TbN3QkIl7LLtB+P5jwAFRsDh2NiEjRUUIkIsVj6acw/y3oeaq/oF0kV/Q7G5YvgI9fDR2JiEjRUUIkIsVj6mOAQY9TQ0ci8n2dj4LGO/khuEVEZJtSQiQixaFis0+IOh0M2+8cOhqR7yupC32GwSevwZK5oaMRESkqSohEpDjMGwMrvtBgCpK7ep8JdUph4gOhIxERKSpKiESkOEwZCQ2awR5Hho5EJLkmO/muc1NG+uHhRURkm1BCJCKFb81SmD0Kup8EpWWhoxFJrf85sG4ZzHgqdCQiIkVDCZGIFL4ZT8Hm9eouJ7mv3UBo0cUPruBc6GhERIqCEiIRKXxTHoHWe8FOe4WORKRyZn4I7oXT4MtJoaMRESkKSohEpLAt+sAfXKp1SPJFj6FQrzGM1xDcIiLbghIiESlsUx6FknrQ/cTQkYikp6yxT4pmPg2rvwkdjYhIwVNCJCKFa9N6mP4EdP4hNGweOhqR9PU7GzZv8N09RUSkVikhEpHCNedFWLtU3eUk/7TsAu33gwkP+JsKi4hIrVFCJCKFa8pIaLIz7Hpg6EhEMtdvOCxfAB+/EjoSEZGCpoRIRArTiq9g7uvQ4xSoUxI6GpHMdf4hNN5JgyuIiNQyJUQiUpim/QtcBfQ8NXQkItVTUhf6DPOJ/ZK5oaMRESlYSohEpPA457vLtRsEO3QMHY1I9fUZBnVKYeIDoSMRESlYWU2IzGygmb1iZovNbIWZTTaznyaUqW9mfzazhWa21szeM7P9sxmHiBS5Be/B0nnQ67TQkYjUTOPW0OVHfrS5DWtCRyMiUpCylhCZ2V7Aa0Bd4BzgBGACcL+ZXRBX9P5o/rXAUcBC4GUz65mtWESkyE0ZCfUaQddjQkciUnP9hsO65TDjqdCRiIgUpGy2EA0FSoAfOeeedc696pw7D3gfOAPAzHoApwI/d87d65x7HTgJWAD8LouxiEixWr8SZj4D3Y6HetuFjkak5toNhBZdYMK9vjuoiIhkVTYTonrARmBtwvRlce9zdFTmidhM59wm4HHgMDMry2I8IlKMZj4DG1dDr5+EjkQkO8yg/3BYOA2+mBg6GhGRgpPNhGhE9PwPM2tjZk3N7BzgYOCWaN6ewKfOucSO0DPxCVWnLMYjIsVoykjYYTdo2y90JCLZs9fJUK8xTLgvdCQiIgUnawmRc24GMBg4BvgS+Bb4J3C+c+7xqFjzaHqipXHzkzKzc81soplNLC8vz1bYIlJIvvkYPh8HvU73Z9VFCkVZY+gxFGY+Dau/CR2NiEhByeagCrsBT+Fbe34EHALcBdxlZrGhngxI1gG6yiMX59w9zrm+zrm+LVq0yFLUIlJQpjwCVuIPHEUKTb/hsHkDTH44dCQiIgUlm13mbsJfH3SUc+5/zrnXnXOXAP8G/m5mdfAtQclagZpFz0uTzBMRqdrmjTD1MdjjCD9UsUihadkZ2u/nu81t3hg6GhGRgpHNhKg7MM05l1hLjwd2AFriW486mFnDhDJdgQ3AJ1mMR0SKyZwXYXU59D4jdCQitWffS2DFlxqCW0Qki7KZEC0CeppZvYTpA4B1+Naf5/D3KToxNtPMSoGTgVecc+uzGI+IFJPJD0PjNtDx4NCRiNSe3Q6Fll3hnb9rCG4RkSzJZkJ0O9ABeN7MjjGzIWZ2O3AKcKdzboNzbip+yO1bzWy4mR2MH3K7A3BdFmMRkWKy/Av45DXodRqUlIaORqT2mMHAS2HxLPj41dDRiIgUhGyOMvcf4EigDLgPP8DCIOBC4JdxRc8CHgRuAEYBuwCHO+cmZysWESkyUx4FnB9dTqTQdTsBmrT1rUQiIlJjWT2V6px7EXixijJrgcuih4hIzVRs9qPL7XogNGsfOhqR2ldSF/b5P3j51/5GrW37ho5IRCSvZbPLnIjItjdvDCz/XIMpSHHpfQbU316tRCIiWaCESETy2+SHoEFz6PzD0JGIbDtljf19iT58HpbMDR2NiEheU0IkIvlrVTnMfgF6nAKlZaGjEdm2+p8HJfXg3X+EjkREJK8pIRKR/DX9cajYCL1/EjoSkW2vcSvoeQpM/Res/Dp0NCIieUsJkYjkJ+f8vYfa9oeWXUJHIxLGvpfA5g0w/u7QkYiI5C0lRCKSnz5/H775CPqcGToSkXB26AhdfgQT7oP1K0NHIyKSl5QQiUh+mvww1GsMXY8NHYlIWAN/BuuWw4T7Q0ciIpKXlBCJSP5ZuwxmPA3dT4CyRqGjEQmrbR/oeJAfXGHD6tDRiIjkHSVEIpJ/pj8Bm9ZCn7NCRyKSGwZfBWuWwPh7Q0ciIpJ3lBCJSH5xDiY+AG16Q5ueoaMRyQ279IeOB/tWovWrQkcjIpJXlBCJSH5ZMA7KZ0Pfn4aORCS3xFqJJtwXOhIRkbyihEhE8svEB6CsCXQ7PnQkIrlll37Q6RC1EomIZEgJkYjkj9VLYNYz0GMo1NsudDQiueeAK6NWIl1LJCKSLiVEIpI/pj3mb0KpwRREkou1Er2jViIRkXQpIRKR/FBRARMfhF32hlZdQ0cjkrsGXwVrl8L4e0JHIiKSF5QQiUh+mP8mLJ2rwRREqtK2L3Q6FN69DdavDB2NiEjOU0IkIvlh4oPQoBl0PSZ0JCK577tWIl1LJCJSFSVEIpL7Vn4Ns/8HPU+DuvVDRyOS+9r2gd2G+BHn1i4LHY2ISE5TQiQiuW/qSKjYBH2GhY5EJH8cdI1Pht6+JXQkIiI5TQmRiOS2is0waQS03w923C10NCL5Y6cesNfJMO5OWPZ56GhERHKWEiIRyW1z34BlCzSYgkh1HHS1fx59Y9g4RERymBIiEcltE+6D7VpA56NCRyKSf5r+APY+H6Y9Dgunh45GRCQnKSESkdy19FP46GV/I9bSeqGjEclPgy6DBk3h1WtDRyIikpOUEIlI7pp4P1gd6HtW6EhE8leDprD/FTBvNHzyeuhoRERyjhIiEclNG9bA5Eegy1HQpE3oaETyW7+zoWk730pUsTl0NCIiOUUJkYjkphlPwbpl0P/c0JGI5L/SMjjkOvh6Bkx/InQ0IiI5RQmRiOQe52D8PdCyK7QbGDoakcKw5/HQpje8cQNsXBs6GhGRnKGESERyzxcTYNF06DcczEJHI1IYzGDIDbDiS3jv9tDRiIjkDCVEIpJ7xt8DZU38TSVFJHvaD4QuP4K3/qabtYqIRJQQiUhuWfk1zHwGep4GZY1CRyNSeA67yXdLfeXq0JGIiOQEJUQiklsmPwQVG313ORHJvqY/gP1+AbOehbmjQ0cjIhKcEiIRyR2bN8LEB6HjQbBjp9DRiBSufS+GZh3ghV/Cpg2hoxERCUoJkYjkjtmjYOVXGmpbpLbVrQ9H/AmWfAzj7ggdjYhIUEqIRCR3jL/Xd+fZbUjoSEQK3+5DYI8jYeyfYPmXoaMREQlGCZGI5IZFM+Czt6Hv2VCnJHQ0IsXh8D+A2wyvXBM6EhGRYEpDByAiAsD7d0LdhtD7jNCRSAbaXzkq7bLzb/5hLUYi1dKsPQz6OYz5A/Q9CzrsHzoiEZFtTi1EIhLeqnKY/iT0OAUaNg8djUhxGXgpNG0Hoy7XAAsiUpRqJSEysyPN7E0zW2VmK8xsopkdFDe/mZndZ2bfmNlqM3vNzLrXRiwikgcmPQib18OA80NHIlJ86jbwAyx8Mwfe/XvoaEREtrmsJ0Rmdh7wLDAJOA44EXgSaBjNN+A54HDgYuAEoC4w2szaZjseEclxm9bDhPug06HQYvfQ0YgUpz0Oh67H+gEWyj8KHY2IyDaV1WuIzKw9cCvwS+fcrXGzXo77+2hgEHCQc2509Lr3gE+BK4BLshmTiOS4mf+FVV/D3heEjkSkuB3xJ5g3Gp6/FIaNgjrqVS8ixSHbtd1PgQrgrkrKHA18FUuGAJxzy4HngWOyHI+I5DLn4L1/wo57+Juxikg4jVvBkBthwbsweUToaEREtplsJ0SDgNnAUDOba2abzOwTM7swrsyewIwkr50J/MDMGmU5JhHJVQveg0XTfeuQWehoRDDXqzsAACAASURBVKTX6X6kuVevgxULQ0cjIrJNZDshagPsBvwZuBkYArwK3G5ml0ZlmgPfJnnt0ui5WbIFm9m50eAME8vLy7MbtYiEMe4OaNAM9jo5dCQiAv7ExFG3wuYN8MLloaMREdkmsp0Q1QEaA+c55+51zr3hnLsAeAm4KhpQwQCX5LWVnh52zt3jnOvrnOvbokWLLIctItvct/Nh9ijocxbUaxg6GhGJ2aEjDL4KZv8PZj0XOhoRkVqX7YRoSfT8asL0V4BWwE74lqBkNxqJtQwlaz0SkUIz/l6wOtBveOhIRCTRPhdB6+6+lWjtstDRiIjUqmwnRDNTTI+1/lREZfZMUqYrsMA5tyrLMYlIrlm/EiY/7If53X7n0NGISKKSUjj6NlhdDq9cEzoaEZFale2E6L/R82EJ0w8DvnDOLcLfg2hnMzsgNtPMmgA/iuaJSKGb8iisXwF7/1/oSEQklTa9YN+LYcoj8HFixw8RkcKR7YToBWA0cLeZnW9mQ8zsHvzgCr+JyjwHvAeMNLOhZnZYNM2AP2U5HhHJNRWb/WAKbftD2z6hoxGRygz+NbT4//buOzyqYv/j+HvSIAGkFwsICNIUEFGqiIIiRcDGz4oFURF7Ra/Xa7/Ye8Ny7b1QFVSKSBHFCkixUBSVYpASStr8/piNrDGBTbK7s+Xzep7zbPbsOZtPwmQ53zPnzLSE8RfDNl3RLiKJKawFkbXWAoOB14CbgYlAZ+A0a+1zgW0KgQG4+4wew/UqFQBHWGt/DmceEYlBiyfAnyvdmWcRiW3pleG4J2DLWnj/Wt9pREQiIuzTUFtrN1lrR1pr61trM6y1ba21rxTbJttae461tpa1Nsta28ta+024s4hIjLEW5jwMNZtAy/6+04hIKPY6CHpcBd++Dosn+k4jIhJ2ab4DiEgS+XkerJ4P/e6BlFTfaSTKGo+aVKbtV4xW0RwzDrsKlr4HEy+DRp2hSh3fiUREwibsPUQiIqWa8zBUrgHtT/WdRETKIi0DjnvSDcE96QrX2ysikiBUEIlIdPzxo5uI9ZBzIaOK7zQiUlb128AR18F342Dh277TiIiEjQoiEYmOTx+D1HQ49DzfSUSkvLpeCnsf7CZs3bzGdxoRkbBQQSQikbc128091HYIVKvvO42IlFdqGgx+AvK2wYRLdOmciCQEFUQiEnmfPwP526DLRb6TiEhF1d0fet0IyybD16/sfnsRkRingkhEIitvO3z2JDQ7Cuq18p1GRMKh0who1BUmj4KNv/hOIyJSIRp2W0Qia8EbkLMOuqp3KB4UDY392k9/AHByGYfKliSRkgKDH4XHu8G4i+CMd8EY36lERMpFPUQiEjmFhTDnEWhwIDQ53HcaEQmnWk3h6Fvhp+kw/1nfaUREyk0FkYhEzg8fwvql0OVinT0WSUQdh0HTnvDBvyF7ue80IiLlooJIRCJn1v1QvSEccLzvJCISCcbAwEcgJRXGjXS9wiIicUYFkYhExqpPYdVcN7JcarrvNCISKTUawjH/hZWzYd4TvtOIiJSZBlUQkciY9QBk1oIOZ/hO4kXjMg5GsGJ0/wglEYmC9qfB4gkw9WZofhTUae47kYhIyNRDJCLht+Y7WPY+dDofMqr4TiMikWYMHPsgpFWGdy+AgnzfiUREQqaCSETCb85DkJ4Fh57nO4mIREu1BtD/Xlg9330GiIjECV0yJyLh9efPsOBNVwxl1fKdRij75Xsi5XbACbB4PEy/A/bvA/Xb+E4kIrJb6iESkfCa+6h77DLSbw4RiT5joP99kFnDXTqXn+s7kYjIbqmHSETCJ+cP+PJ5OHAIVN/Hd5q4ol6c6CrL71sDXpRRlTow4AF4/TT45B444nrfiUREdkk9RCISPp+Ngbyt0O1S30lExKdWA6DtyTDzHlj9pe80IiK7pIJIRMIjNwc+exJa9IN6LX2nERHf+o6GqvVh7AjI2+47jYhIqVQQiUh4fPkCbNsA3S7znUREYkFmTRj4MKxbAjPu8J1GRKRUKohEpOIK8mDOI9CoKzTq5DuNiMSK5r2hw1CY8zD88oXvNCIiJVJBJCIVt+At2PQLdL/cdxIRiTVH3wbV9oRxF0L+Dt9pRET+QQWRiFRMYSHMfhDqtYHmR/lOIyKxpnJ1N+rcuiXw8V2+04iI/IMKIhGpmO+nwLrF0P0yNweJiEhx+x8N7U6BWffDb9/4TiMi8jcqiESkYmbdD9UbQZvjfScRkVjW5w43R9HYkZqwVURiigoiESm/lXPh53nQ9WJI1TzPIrILWbVgwP2wZoE7kSIiEiNUEIlI+c26H7Jqw0Gn+04iIvGgZX844ESYeTesWeQ7jYgIoIJIRMprzSJ3/1CnCyAjy3caEYkXfe+CzBow9kIoyPedRkREBZGIlNPsByG9Chxyru8kIhJPqtSGfvfAb1/Dp4/5TiMiooJIRMphw0o399DBZ7n7AkREyqL1IGjRD6bfAdnLfacRkSSngkhEym7uo2BSoMtI30lEJB4Z43qJUtJg4uVgre9EIpLENCyUiJRNznr48gVoOwSq7+07jUjMaTxqUsjbrhjdP4JJYlz1vaH3f+C9q+Db16Hdyb4TiUiSUg+RiJTNvCcgfxt0vcR3EhGJdx2HwT6HwuTr3MkWEREPVBCJSOi2b4R5Y6DVQKjX0ncaEYl3KSkw8CHYsRmmXO87jYgkKV0yJyKh++wp2LERelzlO4kkAV16liTqtYLul8PMu9yluM16+04kIklGPUQiEprcHDeYQvOjYc92vtOISCI57Eqo3dwNsJCb4zuNiCSZiBZExpjJxhhrjLmt2PqaxpinjTHrjTE5xpiPjDEHRjKLiFTQ/P/BtmzocbXvJCKSaNIru0vn/lzlhuIWEYmiiBVExphTgH+cRjbGGGA8cAxwMXACkA5MN8bsE6k8IlIBedthzkPQpAc0PNR3GhFJRPt2hQ5nwqePw+8LfacRkSQSkYLIGFMDuB+4ooSXBwLdgTOsta9aaycH1qUA10Qij4hU0FcvwpY16h0SkcjqfRNk1nCXzhUW+k4jIkkiUoMq3AUssta+aox5pdhrA4FfrbXTi1ZYazcaYyYAgwCN5SsSS/JzYfaD0LATND7MdxqREpVlAAaJYVm14OjbYOwIdyLm4DN9JxKRJBD2HiJjTHdgKHBhKZu0AUrqC18ENDLGVA13JhGpgG9fh40/u94hY3ynEZFE1+4U2LcbfPQfyPnDdxoRSQJhLYiMMenAk8A91tqlpWxWC9hQwvrswGPNUt77PGPMfGPM/HXr1lU8rIjsXkE+zLrPjSqnoXBFJBqMgf73urmJPrzRdxoRSQLh7iG6FsgEbt/FNgawpawvlbV2jLW2o7W2Y926dSsQUURCtuhdyP5JvUMiEl31WkGXi+Drl2DlXN9pRCTBha0gMsY0Av4F/BuoZIypERhcgaDnqbieoFolvEVRz1BJvUciEm2FhfDJPVC3FbTQpJciEmWHXwPVG8GkK6Agz3caEUlg4ewhagpUBl7CFTVFC8BVga8PxN0r1KaE/VsDq6y1W8KYSUTK67uxsG4J9LgKUjSHs4hEWUYV6HsnrP0OPn3MdxoRSWDhPMr5GjiihAVckXQE8ANuDqK9jTGHF+1ojNkDODbwmoj4VlgAH98JdVpAm+N8pxGRZNWyH7ToBzNGw58/+04jIgkqbAWRtfZPa+2M4kvg5ZWB51twRc9c4CVjzMnGmD6BdQY3XLeI+LboXdc71PNaSEn1nUZEklnfO93j5FF+c4hIwor6dTDW2kJgAPAh8BjwLlAAHGGt1ekfEd+KeofqtoLW6h0SEc9qNHL3Ey2ZCEsn+04jIgkoUhOz/sVa+4+hqay12cA5gUVEYsnCd2D9MjjpOd07JCKxofNI+OY1eP9qaNIDMrJ8JxKRBBLxgkhE4khR71C9NtBqkO80Edd41KQybb9itEbbS0ZlbScSAWkZ0P8+eK4fzLwbev/HdyIRSSA6/SsiOy14C/74PnDvkD4eRCSGNO4G7U6FOQ/D2iW+04hIAtERj4g4Bfkw8y6ofwC0PNZ3GhGRfzr6Vjcc96QrwZY0x7uISNmpIBIRZ+Fb8McP0HOUeodEJDZVqQNH3QwrZ7l7ikREwkBHPSLieoc+vhMaHAgtB/hOIyJSuoOGwj6Hwgc3wNZs32lEJAFoUAURgQVvQPZPcPIrYP4xMGRc0Q3wIgkuJQUG3AdPHg5Tb4FjH/CdSETinHqIRJJd/g6Y/l/Ys72bEV5EJNY1OBA6XQBfPAe/zPedRkTinHqIRJLd58/AxlUw8KG47x2KNPU+icSQI66DRe/CxMtg+AxI1SGNiJSPeohEktn2jW5Oj6Y9Yb8jfKcREQldpWpwzH/h9wXw+VO+04hIHFNBJJLM5jwM27Kh902+k4iIlF3rQdCsN0y7HTb95juNiMQpFUQiyWrzGpj7KLQ5HvY6yHcaEZGyMwb63Q2FeTDlOt9pRCROqSASSVYf3wkFuXDkDb6TiIiUX62mcNiV7n6iH6b6TiMicUgFkUgy+uNHNzrTwWdB7f18pxERqZhul0LtZvDeVZC33XcaEYkzKohEktG0WyGtMvS4xncSEZGKS6sE/e9186nNut93GhGJMyqIRJLN6i/dpSVdRkK1+r7TiIiER9OecMCJMOs+WP+D7zQiEkdUEIkkE2vho/9AVm3oerHvNCIi4dXnDkjLdHMTWes7jYjECRVEIsnk+w9g+UzocTVU3sN3GhGR8KpWH46+BVZ8Al+96DuNiMQJFUQiySJ/B0weBbWbQ8dhvtOIiETGQUNh3+4w5QbY/LvvNCISB1QQiSSLuY+4G4773glpGb7TiIhERkoKDHwI8re7UedERHZDBZFIMti4GmbeAy0HQLNevtOIiERW7f2g5yhYPAG+G+87jYjEOBVEIsngw3+DLYQ+t/tOIiISHV0vhgYHwntXw7Y/facRkRimgkgk0a2YBQvfdhMX1mzsO42ISHSkpsPAhyFnLXx4o+80IhLDVBCJJLKCfHjvGqjeELpd5juNiEh07XUQdLkIvnweln/iO42IxCgVRCKJbP6zsHaRu1QuI8t3GhGR6Ot5HdRsAhMugdytvtOISAxSQSSSqHLWw/TboMnh0Gqg7zQiIn5kZLlR57J/0qVzIlIiFUQiiWrqzZCbA33vAmN8pxER8adJD+h8IXz+FPzwke80IhJjVBCJJKIfp8OXL0DnEVCvpe80IiL+9boR6raEsSNha7bvNCISQ1QQiSSaHZth/MVQuxkc8S/faUREYkN6Jhw/Brauh0lXgLW+E4lIjEjzHUBEwuyDf8PGX+CcKe4AQERiVuNRk8q0/YrR/WPivePWnu3cIAvTboUW/aHtSb4TiUgMUA+RSCL5cTp88T/oMhIadfKdRkQk9nS7DBp2gveuhI2rfacRkRiggkgkUfx1qVxzOPIG32lERGJTahoc94Sbp23sCCgs9J1IRDxTQSSSKIoulRv8mC6VExHZlVpN4Zg7YPnHMO8J32lExDMVRCKJIPhSuYaH+k4jIhL7OpwJLfq5uYl+me87jYh4pEEVROJdHF4qV5abvZPiRm+RGJB0f5fGuB71J3vAG2fC+TOhSm3fqUTEA/UQicQza2HCZbBptS6VExEpq8yaMOQFyFkL7wyHwgLfiUTEAxVEIvHss6dg4VuuZ0iXyomIlN1eB0Hfu+DHqTDzHt9pRMQDFUQi8ernz2HK9bD/MdDtct9pRETi18FnQduTYcZ/4YepvtOISJSpIBKJRznr4c0zYY+93PCxKfpTFhEpN2NgwP1QrxW8fa4bsVNEkkZYB1UwxpwInAJ0BOoBq4B3gDustZuDtqsJ3A0MBjKBucDl1toF4cwjkpAKC9x/2DnrYdgH7hr4BFaWG71FRMotIwuGvAhjesKbZ8FZ70Fahu9UIhIF4T6tfBVQAFwPHAM8DowAPjTGpAAYYwwwPvD6xcAJQDow3RizT5jziCSeGaPhp+nQ727Yq73vNCIiiaNOMxj0CPzyOUy41A1cIyIJL9zDbh9rrV0X9PxjY0w28DzQE5gGDAS6A0daa6cDGGPmAsuBa4BLwpxJJHEs+wBm3gXtT4cOQ32nEZEoi9ce07ga0rvNYFh3Pcy4A2o2hp7X+s0jIhEX1h6iYsVQkc8Dj3sHHgcCvxYVQ4H9NgITgEHhzCOSUNYugXfOhfoHQv973DXvIiISfodfA+1OdUXRN6/7TiMiERaNO7EPDzwuDjy2ARaWsN0ioJExpmoUMonEl42r4aXjIa0ynPyS5hsSEYkkY+DYB6HxYTBuJKyY5TuRiERQRAsiY8zewC3AR9ba+YHVtYANJWyeHXgs8Q5xY8x5xpj5xpj569aV1BElkqC2bYCXToDtm+C0t9wlHCIiEllpGfB/L0KtJvDaabD+e9+JRCRCIlYQBXp6xgH5wNnBLwEl3aW4y+t/rLVjrLUdrbUd69atG76gIrEsbxu8egpk/winvAJ7tvWdSEQkeWTWhNPehNR0ePlE2KITsiKJKNyDKgBgjKmMG0muKXC4tTZ4QP9sXC9RcUU9QyX1Hokkn4J8eGsYrPoUTnwWmvTwnahU8Xqjt4jEnrJ+nkR8EIaajeGU1+C5/vDKEBg6DirvEdnvKSJRFfYeImNMOvA2cCjQr4S5hRbh7iMqrjWwylq7JdyZROKOtfDelbB0EvS9Ew443nciEZHktU9HOPF/8Pu38PJJsEOHKiKJJKwFUWCuoZeBXsAga+2nJWw2HtjbGHN40H57AMcGXhNJbtbCRzfBF89B9yug0/m+E4mISMt+cMLT8Mtn8OrJkLvVdyIRCZNwXzL3KHAScDuQY4zpHPTaL4FL58YDc4GXjDFX4y6Ruw53D9FdYc4jEl+shcmjYN4TcPDZ0OtGb1F0GZyISDFtjnOXM78zHF471V1Kl17ZdyoRqaBwXzLXN/D4L1zRE7ycC2CtLQQGAB8CjwHvAgXAEdban8OcRyR+FBbAhEtcMdT5Qhhwv+YaEhGJNW1PgkGPwk/T4Y2hkJ/rO5GIVFBYe4istY1D3C4bOCewiEhBPoy9ABa8CYddBUfeoGJIRCRWHXQaFOTCxMvgrbPhpOfcSHQiEpeiMTGriOxK/g5480xXDPW6EXr9W8WQiEis63g29L0blkyE189w0ySISFxSQSTiU26Om/BvyUQ4ZjQcdqXvRCIiEqpO50H/e2HZ5MAE2ht9JxKRclBBJOLLxtXw7DHw41Q49iHoPMJ3IhERKatDznWjz/08D54boMlbReKQCiIRH379Cp7uBdk/uVGKDj7TdyIRESmvA0+EU16H9d/Ds31gw0rfiUSkDFQQiUTb4gnwv36QkgbDPoD9+/hOJCIiFdW8NwwdB1vXu97/tUt8JxKREKkgEokWa2HW/fD66VCvNQyfBvXb+E4lIiLh0qgTnPUe2ALXU7Rilu9EIhICFUQi0ZC7FcaOgI9uggNOgLMmQtV6vlOJiEi4NTgAzpniPuNfGAxfv+o7kYjshgoikUj740d45ij45jXoeR2c8AykZ/pOJSIikVKribsket8ubo65abdBYaHvVCJSChVEIpG0eCKM6QmbVsNpb0HPUZpjSEQkGWTWhNPfgYPOgJl3w9vDIG+771QiUoI03wFEElJBPky9GeY8BHt1gCHPQ41GvlOJiIRF41GTfEeID6npMPBhqNMcPrwRNv4CJ78CVev6TiYiQdRDJBJum36DFwa6YuiQc+GcySqGRESSlTHQ7VIY8iL8vgCeOgJ+/dp3KhEJooJIJJyWvg+Pd3XzDB3/lJvBPK2S71QiIuJb64HuBJm1bgS6b9/wnUhEAlQQiYRD3jaYdBW8ejJU3xvOnwlth/hOJSIisWSv9nDeDNi7I7wzHKb8y11iLSJeqSASqag138FTR8LnT0GXi+Dcqe56cRERkeKq1oWhY6HTBTD3EXjpeNia7TuVSFJTQSRSXtbCZ0+568Fz1sFpb0Of23WJnIiI7FpqOvS9EwY9Bqs+hTGH674iEY9UEImUx8bV7qzee1dB4+4wYg407+07lYiIxJODToOz34fCAnjmaPjiOXeyTUSiSgWRSFlY626EfbyLO6vX/143v1DVer6TiYhIPNrnYDj/E2jcDSZcCmNHQO5W36lEkormIRIJVc56mHg5LB4PDTvB4Meh9n6+U4VM84aIiMSoKrXdybWZd8OM0fDbN26Y7jrNfCcTSQrqIRIJxeIJ8FhnWDYZet/kLnGIo2JIRERiXEoq9BwFp78Nm3+HMT1h4Tu+U4kkBRVEIruyeQ28fga8fjpUbQDDp0P3y91/XCIiIuHWrBdc8AnUawVvnQ3jL9EldCIRpoJIpCTWwpcvwqOHwLIp0OtGOG86NDjAdzIREUl01feBs9+D7lfAly+43qLfF/pOJZKwVBCJFJe9HF4cDOMvgnptYMRsOOxKN0yqiIhINKSmQ+//uDmLtm90893NG6NR6EQiQIMqiBTJ2w5zH4aZ90JKGvS/Dw4+G1J03kBEJFaVZcCYFaP7RzBJhDTt6U7MjR0B718NP06DQY9AlTq+k4kkDB3piYC7LO6xzjDtNtj/aBg5Dw4ZpmJIRET8q1IHTn0DjhkNP06FRzu5wX5EJCx0tCfJLXs5vHIyvDLEXZ5wxrsw5AWovrfvZCIiIjsZA51HwHkfu/+jXj8d3h4OW7N9JxOJe7pkTpLT9o0w+0GY84grhI66FTpdAGkZvpOJiIiUrn5rOHcqfHIfzLwLln8Mxz4ELY7xnUwkbqmHSJJL3jaY/RA82A4+uRdaD4KLPodul6gYEhGR+JCaDj2vheHTIKsOvPp/8O4IyPnDdzKRuKQeIkkOBfnwzStuBvBNq6FZbzeU9p7tIv6ty3LDb1nF5Q3CIiKeJNzn8Z7t4LwZ8PGdMPsBWPqeG5muw1m6B1akDPTXIoktPxe+fhUe7wLjL4Y99oIzJ7qZwKNQDImIiERUWgb0+jdcMAvqHwATL4ene8GvX/lOJhI31EMkiWn7JvjyeZj7GGz+1c0n9H8vQ8v+7sbUCojkGUYREYlPZf2/Iew9SvVawVkTYcGbMOVfMOYI6HgOHHkDZNUK7/cSSTAqiCSxbPoV5j0B8/8HOzZB48Ng4MPQrFeFCyEREZGYZgy0HQL794Hpd8BnY2DBW9D1Yuh8AVSq5juhSExSQSTxLz8Xvp8CX70E338IWDdYQtdLYO8OvtOJiIhEV+Xq0PdO6DAUpt0O02+DeY/DYVdCx2GQXtl3QpGYooJI4teaRfDVy/Dta7D1D6jawJ0FO/hMqNXUd7qo0OV7IiJSqvpt4JRX4Jf5MO1WmHK9m26ix1XQ/lRIz/SdUCQmqCCS+FFYAKu/cKPoLH0f1i2BlHRo0RcOOgP2OxJS1aRFRET+Zp+OMHQcLJ8JU2+FSVfAtNvcCcSOw6BGQ98JRbzS0aPEtq3ZsHI2LJsMy6ZAzjpISYN9u7qbRQ84EarU9p1SRESkTLwMAd6kBwz7AFZ8AvOedBOUz34QWg6ATufDvt10v60kJRVEEjushT9Xwaq5blk5F9Yvda9Vqg7Nj3K9Qc16Q2YNv1lFRETikTGuMGrSw/2f+/nT8OULsHg81G4GrQdDm8FuCG8VR5IkVBCJH/k7YN1SWLPQ3QtU9Jizzr1eqTo06gTtToZGXVx3f2q638wiIiKJpEYjOOoW6HmdG41uwZsw6z745B53L27rQW5p0E4TvUpC81YQGWMaAvcDRwEG+Ai4zFq7ylcmCaPCQtiWDZt/c2egsn+C7OXuccNy+PNnsAVu27TKbv6E/fvAnu1dAVSvFaSkRi2uBicQEZGklZ4JHc5wS856WDIRvhsHsx+CWfdDZk1o2Bn27QKNurqJzdMyfKcWCRsvBZExJguYBuwAzgQscBsw3RjT1lqb4yOXhKAgH7ZtgC1rYPPvruDZ8nvg66BlyxoozPv7vpVruDNOe3eEA0+Ceq1dl3zt/aJa/IiIiEgpqtSBg89yy9Zsd//uylnuMvZl77tt0jJhr/ZQZ3+o28I91mkO1RvGzv/nhYWQt9UtuVsgN8ddnVKQBwW57hilID/wmBv0dWDBunuWU9PdY0q6G7gpJT1oXeD11EqQUQUqVYWMqm6+p1j5PUhIfPUQDQeaAi2stT8AGGO+Bb4Hzgfu85QrORTkQ+5m2BG0bN/kCp1tG1zPztbsnV9v27Dz+Y5NJb9nZk2oticzf0tlLU1ZazuwxtZkra3BaluHlbY+G7dXhT+Bn9wuZZ2lW704IiIiUZRVC9qf4haALWt33uP761euJ+nL53dun1YZqu0JVepC1XrusUpdV2SlZ7rX0yrv/Dol1d0/jA16LHSFS95WyNsWeNwe9HXxx62u2Cm+5Hk+t56WubNAyqgaVCxVhUp7uKKp0h5QOejrStUCz4PWpWfqXq4o8FUQDQQ+LSqGAKy1y40xs4FBxHtBZAN/0EWPBD2nhNcKC3eelSjMc8NLF31dkA+Fxc5a5G8LfDhshfzAh0RJj7lb/1n47Njs9tsl4yZ1y6oFmbUgq447+5NZyxU+WbXcB13VBlCtAVSt/9ckb0NVtIiIiCSmqvV23ldUJOcP+ON7WL/MLZt/d4VT9k+w6lM3TyA2PN8/Jd0VCEVLWqbrmcmo4o5VMrICz6vuXJ9RBdKrcMEbS9hOOvmkkWfTyCOVfFLJY+fX+aSRa9PIJxULpFFAOgWkGvdY9DyNfNIoII1C0kw+L5/ZPtALFeiJ2rElcPwVvG4zbF0PG1YEjsc2hXA8huuFqlQtsFQPPFYN/A6yAgVmVuB50NdpQb+nv3q6gnq3/noe3OuVHrhXzIBJCRRiQV+blH8+T5BizVdB1AYYV8L6RcBJUc5Sfh/d5CY4Cy56fCp+5iU9y/3hZNWBmk2C/qD2CPo6aF1RsVO5urp6RUREZPeq1HZLo84lv15YANv+3HkyN/ixsCDooDvo4Du9lDM2xAAACMdJREFUctABfdbOA/sKDK40+bUKnLDdVT1ngRbHlO99C/JdYVRUIBVdsbNjM+zYWOx50DZb1uzsNcvftrO3zNtxqPl7wTRynrsdIo4Ya8NUtZflmxqTC9xnrR1VbP1twChr7T8KNWPMecB5gactgKURDxqaOsB63yEkaai9STSpvUk0qb1JNKm9JYd9rbV1d7eRz2G3S6rESu13s9aOAcZELk75GGPmW2s7+s4hyUHtTaJJ7U2iSe1NokntTYL5GlR+A1CrhPU1A6+JiIiIiIhEnK+CaBHuPqLiWgPfRTmLiIiIiIgkKV8F0XigszGmadEKY0xjoFvgtXgSc5fxSUJTe5NoUnuTaFJ7k2hSe5O/+BpUoQrwDbANuAF3P9GtQDWgrbV2S9RDiYiIiIhI0vHSQ2StzQGOBJYBLwIvA8uBI1UMiYiIiIhItHjpIRIREREREYkFvu4hijnGmBRjzHXGmBXGmO3GmG+MMSeUYf/BxpivAvuuNMbcYIxJLbZNd2PMc8aYhcaYfGPMil28X0NjzFvGmI3GmE3GmHeMMY0q8CNKDIlGewts190YM8cYs80Y87sx5j5jTGaxbc4yxtgSlq/D8bNK9FTkc8MYU9kYc7cx5rdAe5lrjOlRwnYht11jzHBjzBJjzA5jzFJjzAUV/RkldsRSezPGzCjlc+yycPys4l+U2tsVxpgJge2sMeamXbxnSP8PS3xQQbTTrcBNwCNAX+BT4E1jTL/d7WiM6QO8DXwe2PdB3L1RdxTbtBdwGG6UvcW7eL8sYBrQEjgTOANoDkwP3H8l8S/i7c0Y0xb4EFgLDAhsczbwXClvfRLQJWg5o2w/kvgUhs+NZ4DhwI249vIbMMUY077YdiG1XWPMcOBJXFs9BngTeMwYM6I8P5/EllhrbwHf8vfPsC7Aa2X6wSQmRbG9DQfqAWN3kyfU4z6JF9bapF9wjX8HcHOx9VOBb0PY/yvg42LrbgRygQZB61KCvn4JWFHK+10KFADNgtY1AfKBK3z/vrRUbIlie3sX+B5ID1o3FDeISYegdWcF1jUr68+iJXaWinxuAO0CbeDsoHVpwFJgfNC6kNpuYN+1wPPFtnsWNzN8ell+Ni2xt8RSewusmwHM8v170RKZJRrtLbA+Jeh1C9xUynuG9P+wlvhZ1EPk9AEycEVKsJeAA40xTUrb0RjTEGhfwr4vAum4MwcAWGsLQ8wzEPjUWvtD0L7LgdnAoBDfQ2JXxNubMSYdd1b+DWttXtB2b+A+sNWOEk9FPjcGAnnA60H75uPOrvcxxlQKrA617XYB6paw3YtAbaB7iD+TxK5Yam+S+KLR3kI6TivLcZ/EDxVEThvcWagfiq1fFHhsvZt9ARYGrwz8oW7dzb67es+FJaxfVM73k9gSjfa2H1C5hO22Az+W8j1mGWMKAtdOP2GMqbW7H0RiSkU+N9oAy621W0vYNwNoFrRdKG23xHZawnYSv2KpvRU5KHB/SZ4x5ltjzLDd5JD4EY32VpYsFM9TweM+8SzNd4AYUQv401pbfMi97KDXd7UvwIYSXtuwm3139Z4lvV82ULMc7yexJRrtbVfbZRf7Hr8BtwDzcHODdQOuBboZYw4JFFES+yryubGrfYteL3oMpe2W1v5CaeMSH2KpvQHMxE3hsQyogbs8+GljzJ7W2tt2k0diXzTaW1myUMp7lve4TzxLyILIGNMbdzP57nxsre0JGNy1ov94q1C+XeCxvPuXJtzvJxESo+0t5HZprZ0CTAlaNd0YswB3U+npwNMh5JLYUJF2Fcq+ZdmutDySOGKlvWGtvbHYqnHGmHeBfxljHrCa4zARRLq9hSpSx33iUUIWRMAcoFUI2xV1n2YDNY0xptiZqJpBr5dmV2cYauxm39KUdoahJiWfkRC/YrG97Wq7muy85KQ044Ec4BBUEMWLinxuZAMlDV9bvE2G2naD299vQdvVKva6xK9Yam+leRUYDBwIzN3NthLbotHeQhWJ4z7xLCELosB1okvKsMsioBLuvovga5WLrgP9bjf7grum9K8PXGNMYyBrN/vu6j3blLC+dTnfTyIoRtvbj7hr7//WjowxlYGmuCGQQ6Ez/PGjIp8bi4DjjDFZxa6zb40bhOOHoO1CabvB7fS3XWwn8SuW2ltp1FOZOKLR3sqSBcJ73CeeaVAFZzLuj+K0YutPBxYGbpQrkbV2FfBNKfvmAe+XI894oLMxpmnRisAfWrfAaxLfIt7erLW5ge8zxBgTfOLjRNwBxu7a0WCgCu6+IokPFfncGI8bHemkoH3TgP8DPrDW7gisDrXtzsUNr13Sdtm4kaEkvsVSeyvNqbj7IhfsZjuJfdFobyGJ0HGf+OZ73O9YWYDRwHbgCqAn8DhQCBxbbLupwA/F1vULbPtkYN/LA+91d7Ht6uIOSE/E3QC6Nuh566DtquDOWCzADSc5EPfH9xNQ1ffvSkvctLf2uIOBd3CTAg/DHYy+WWy7D4HrA+3sKNwkiFuAr4FKvn9XWkJuUyF9bgD74ubuuLHY/q/hLj05N9Be3gq0qw7Ftgu17V4QWH9bYLtbAs9H+v5daUms9oab8HxS4DOuF3A8MA7XM3St79+Vlrhqbx1xx2RDAu3nDXYep2UFbRfS/8Na4mfxHiBWFiAVN8vwStylRt8CJ5aw3QxKmFA18AH8TWDfVbgJulKLbdMz8AdW0nJTsW0b4WZB3gRsxt3g3tj370lL/LS3wHY9cGfrtwNrgAeCP9QD2zwALA60s1zc5Xb3ANV9/560lLld7fZzA2hcymdOJnAf8HugvcwDepbwPUJqu4Ftz8eN+rUDN0nwhb5/R1oSr73hhk1+H1gd2GYL7t7OU3z/jrTEXXt7jtKP04p/r5D+H9YSH4sJ/KOKiIiIiIgkHd1DJCIiIiIiSUsFkYiIiIiIJC0VRCIiIiIikrRUEImIiIiISNJSQSQiIiIiIklLBZGIiIiIiCQtFUQiIiIiIpK0VBCJiIiIiEjS+n8brJfFwWC5ZQAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 1008x504 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
     }
    ],
    "source": [
     "###\n",
-    "# Synthetic data 1\n",
+    "# Synthetic data 1 - without unobservables Z.\n",
     "\n",
     "# R ~ U(0,1)\n",
-    "# X ~ N(0,1)\n",
+    "# X ~ U(0,1)\n",
     "# T ~ Bin(1, sigmoid(a * r + b * x))\n",
-    "# Y ~ Bin(1, sigmoid(c * t + d * x))\n",
+    "# Y ~ Bin(1, sigmoid(c * t + d * x)), but if T = 0 then Y = 1\n",
     "\n",
-    "# Weights:\n",
+    "# Weights for edges:\n",
     "# R -> T: a\n",
     "# X -> T: b\n",
     "# T -> Y: c\n",
@@ -115,6 +142,7 @@
     "\n",
     "    y = npr.binomial(n=1, p=sigmoid(c * t + d * x), size=N)\n",
     "\n",
+    "    # If decision is negative, we set Y to positive.\n",
     "    y[t == 0] = 1\n",
     "\n",
     "    return r, x, t, y\n",
@@ -129,44 +157,48 @@
     "        1)\n",
     "\n",
     "    # P(y=0|t=0, x) * p(t=0|r, x) * p(x)\n",
-    "    # Now equal to 0, hence commented out. See data generation.\n",
+    "    # Now equal to 0, hence commented out. Y=0 and T=0 mutually exclusive.\n",
     "    t_0 = si.quad(\n",
     "        lambda x: (1 - sigmoid(c * 0 + d * x)) * (1 - sigmoid(a * r + b * x)) *\n",
     "        1, 0, 1)\n",
     "\n",
-    "    return t_1[0]  #+ t_0[0]\n",
+    "    return t_1[0]  # + t_0[0]\n",
+    "\n",
+    "# Calculate the difference between the analytic value and an estimate obtained \n",
+    "# from synthetic data.\n",
     "\n",
+    "diffs = np.zeros(0)\n",
     "\n",
-    "r, x, t, y = generateData()\n",
+    "for i in range(1000):\n",
     "\n",
-    "# Fit predictive models.\n",
-    "lr_t = LogisticRegression(solver='lbfgs')\n",
-    "lr_y = LogisticRegression(solver='lbfgs')\n",
+    "    r, x, t, y = generateData()\n",
     "\n",
-    "lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
-    "lr_y = lr_y.fit(np.array([t[t == 1], x[t == 1]]).T, y[t == 1])\n",
+    "    # Fit predictive models.\n",
+    "    lr_t = LogisticRegression(solver='lbfgs')\n",
+    "    lr_y = LogisticRegression(solver='lbfgs')\n",
     "\n",
+    "    lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
+    "    lr_y = lr_y.fit(np.array([t, x]).T, y)\n",
     "\n",
-    "def causal_effect_of_R_on_Y(r):\n",
+    "    def causal_effect_of_R_on_Y(r):\n",
     "\n",
-    "    # Equal to integration of P(Y=0|T=0, X=x) * P(T=0|R=r, X=x) * f(x)\n",
-    "    # from 0 to 1\n",
-    "    print(\"Bias:\",\n",
-    "        si.quad(\n",
-    "            lambda x: lr_y.predict_proba(np.array([[0, x]]))[0, 0] * lr_t.\n",
-    "            predict_proba(np.array([[r, x]]))[0, 0], 0, 1))\n",
+    "        # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
+    "        return si.quad(\n",
+    "            lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
+    "            predict_proba(np.array([[r, x]]))[0, 1], 0, 1)\n",
     "\n",
-    "    # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
-    "    return (si.quad(\n",
-    "        lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
-    "        predict_proba(np.array([[r, x]]))[0, 1], 0, 1))\n",
+    "    r0 = causal_effect_of_R_on_Y(0)\n",
+    "    r1 = causal_effect_of_R_on_Y(1)\n",
     "\n",
+    "    analytical = analytic_R_on_Y(1, a, b, c, d) - \\\n",
+    "        analytic_R_on_Y(0, a, b, c, d)\n",
     "\n",
-    "r0 = causal_effect_of_R_on_Y(0)\n",
-    "r1 = causal_effect_of_R_on_Y(1)\n",
+    "    diffs = np.append(diffs, analytical - r1[0] + r0[0])\n",
     "\n",
-    "analytical = analytic_R_on_Y(1, a, b, c, d) - analytic_R_on_Y(0, a, b, c, d)\n",
+    "    if i % 100 == 0:\n",
+    "        print(i / 10, \"%\", end=\" \")\n",
     "\n",
+    "print()\n",
     "print(\"Analytical:\", analytical)\n",
     "print(\"Estimated: \", r1[0] - r0[0])\n",
     "print(\"Difference:\", analytical - r1[0] + r0[0])\n",
@@ -174,29 +206,69 @@
     "print(\"Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\\n\")\n",
     "print(\"Analytical:\", analytic_R_on_Y(1, a, b, c, d),\n",
     "      analytic_R_on_Y(0, a, b, c, d))\n",
-    "print(\"Estimated: \", r1[0], r0[0])"
+    "print(\"Estimated: \", r1[0], r0[0])\n",
+    "print()\n",
+    "print(\"Average difference:\", np.mean(diffs))\n",
+    "plt.hist(diffs, bins=50)\n",
+    "plt.title(\"Histogram of differences over 2000 simulations (analytic - estimate)\")\n",
+    "sns.kdeplot(diffs)\n",
+    "plt.axvline(x=0, c='r')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It can be observed that when there are no unobservables, the causal effect $$P(Y=0|do(R=r))=\\int_x P(Y=0|T=1, X=x)P(T=1|R=r, X=x)f_x(x)~dx$$ can be estimated without bias.\n",
+    "\n",
+    "Next we generate data and estimate the causal effect from a causal model with unobservables Z. \n",
+    "\n",
+    "Z influences both T and Y."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 91,
+   "execution_count": 319,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Bias: (0.048298659980495415, 5.362228436888762e-16)\n",
-      "Bias: (0.02685101277580984, 2.981061261820832e-16)\n",
-      "Analytical: 0.030762705619782782\n",
-      "Estimated:  0.021455460896022127\n",
-      "Difference: 0.009307244723760655\n",
+      "0.0 % 10.0 % 20.0 % 30.0 % 40.0 % 50.0 % 60.0 % 70.0 % 80.0 % 90.0 % \n",
+      "Analytical: 0.02047548901524536\n",
+      "Estimated:  0.01934600091154902\n",
       "\n",
       "Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\n",
       "\n",
-      "Analytical: 0.16344054312353337 0.1326778375037506\n",
-      "Estimated:  0.15963912402042202 0.1381836631243999\n"
+      "Analytical: 0.10799301700601992 0.08751752799077456\n",
+      "Estimated:  0.11028400555152502 0.090938004639976\n",
+      "\n",
+      "Average difference: 0.0013544947748748567\n",
+      "Std of differences: 0.0016880435339749395\n"
      ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.lines.Line2D at 0x1a20fb7668>"
+      ]
+     },
+     "execution_count": 319,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0QAAAG1CAYAAADZSDS7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XecFEX+//HXZ3dhQZQkCCICoiIKCKiAKII5/8wZjHdiulPvzvP0kuH8enqnd2bMooKnZzoDJkRARKIIIipGgoqCkuOyM/X7o3plHGZ2Z3dnt3Zn3s/HYx4L3TU9n5nuru5PV3W1OecQERERERHJRwWhAxAREREREQlFCZGIiIiIiOQtJUQiIiIiIpK3lBCJiIiIiEjeUkIkIiIiIiJ5SwmRiIiIiIjkLSVEecTM5pnZvNBxyObMrNjMbjSzL8xso5k5M+tVxWWdE73/nIRpnaJpw1OUP8rMppjZyqjMbQnzzjazD8xsTTTv8qrEJJKrzGx4tG90Ch0L1Hw85dUl9YGZPWBmX5tZcehY0qnpY7WZjTOzvH3mSl3bZ7PJzF43s2lmZqFjqW+UENUzCQej/5VT5riozLVZ/NycrUDqiN8DVwPzgZuB64DvavpDzawz8CywPXB/9LmvRfMGAMOBYuDOaN7kmo5J6g8z287MfmNmb5rZQjMrMbNvzOwJM+tezvt2NbNnzexHM1trZu+Z2VkVfM5wM/vOzNab2Zzoc1Mew8ysuZndEcW0IbrQcL2ZNcrG9851uVrfm1k34FzgRufchtDx1BQz2z/b5wD1SX1L2qNYx2VpcdcCewGnZml5eaModABSqw4KHYCkdQSwGjjMObexBpb/DbArsCJp+kH4hOe3zrknU8QEcLZzTomQpPJr4A/AZ8CrwFKgO3AacIKZHe6cG5f4BjPrAUwEGgJPAd8DxwKPmtkOzrnrkspvB0wF2uKT9y+Bg4F/AbvjT3ATy28FvA30wCf3s4C+wF+A/lFMsSx9f/AXMm7C72P5IF1dUh9cA6wEHgodSGBnAVuEDiKgnN1nnXOTzOxdfGKUfEyXcighyiPOuS9CxyBpbQv8WEPJENFyP0nzuZC6Naq8eSLgE5X9nHPvJE40s5OB/wL3ALslvWcYsCU++R8dlb8GeAf4i5k96Zybm1D+ZqAd8Evn3ENR+ULgBeAcMxvhnBuTUP4P+GToBufcXxJiGgZcCJxDFk+InXOLgEXZWl5dV05dUqeZWRvgOGB4LrcOZcI5tyB0DCHlwT47ErjbzAY6594OHUy94ZzTqx69gE6AA/5XTpnjojLXJk2fB8xLmtYCuBF/gFsLLAM+xJ/IbJnwPpfiNTxpWb8ApgFr8Ffh3gaOTRPjjvirvSuism8APfFdtBzQKaHsOdG0c4Dj8d221gDjovnNgKuACfiT9xJgAXAf0DbFZ5d9Rmf8ydMXwDpgJnB4VKYpcDe+0lwHvAXsUsl1NRB/hXpZtIzZwBVAUUKZa9P8tuMyWH4T4Fbg22j5M4ATE3+vFNvN8KT/p3qdk25e0ucfALwC/AisBz6K1kNRUrly119UpgA4H5iCbylbDbwLnFDO+tsBuBT4FNgQrcfL0vxWjfDdEmdEy16Jbzm4AWiQVLY38DS+5aJsuX8n2h+Syp6CP5H/IVoH84H/AQMqsZ1UuN/gr2w74KQ0y/hjNP/0mlpHVair5kbLbJUwrWs0bXSK8sdG8/6eMK1ptA4+TVG+Z1T+PwnTDL8/rAC2SCrfCl83TMow/grrxqTtMV2ddSwwPVrGfOB3CbH+Ntp+1wNzgKNTxDGPpHq7vHlp4sm4jqSC+p6kuiTpvRXWeVG5/aNlXAv0Ad7E75dL8SdzrVMs+2BgdBT/euBr4HXgmAzX52+izzw0xbx2wPX4BH9JtM19DtwCbJWi/LhoWQ3wXYnnR++ZQ9I+GJXvAvwTf4xZlrC+/0xS/ZNqvQIPR5/XO813e6xsPumPKS45/hTLybiezOYLfz7wSLROS6K/91RlO6Cc4xfRPkHN77N74s8f5kS/4ZroN70EsBT7QarX/gnliqP1MiuKaQV+nxmU5vdsEy3joZpaZ7n4UgtRHotuunsd39/0DeBFfDeWHfGVwo34SvG26P89gduB5dEiZiYs607gV/jK475oOScD/zOz3znn/pVQdnt8l5ltos/8COiFPxH8aZkpnIrv4vUCvkIvjabvij8ovQU8gz8Q9wSGAoea2R7OuWUplvdvfMX1Er61dDDwopnti6+Mi/BNzp3wSebLZtbVZdDdxsxOAZ7AV15P4Q+CR+IPigPM7Hjna65x0VvKBisoG9BgXgXLLwBexleoM/AHxHbRZ46uKD78Orwuev8g4NGEz5wZzTuOzdd52ef/CrgDnwi8EH2/AfjEoS9wQorPTLn+ou3wP/jk4qMoFoCjgGfN7HLn3O0plndL9Jkv47ffk4HbzGyDc+7ehFi3AMYAe0fLfwB/sOgKXBktZ3lU9nj8Oi/BJzbfAXvgTyYPiK64lURlLwHuwidMT+L3lXbR77k/PlEqVyX2mxH4E53B+G082RnR57+QsOysraMqKmvtTFzGoOhvqm30TSCeUAagP/43eTO5sHNulpl9n1S+C75l81Xn3Nqk8j+Y2Qygr5k1cs6tTxd4JerGipyA/z2fx28PJwC3mNlafL11An77LQCGAM+Z2W7Ouc8zWHZlVKaOrLC+T6USdV6iPvh98E3gXmBf/Lbc2cz2KStvZkfj18F3+G1zKb4LZT/gmGheRQ7A7/dTUswbiE+YxuAvxLho2b8DBprZvi516/2T+PphFP54cTrwhJktd869mlDuBHzXzrei71qM327/ht/Gjqsg9vuj95+H76b6EzNrir8Q9r5z7n0za4Y/Zp0NjGfTMaZclakns8nM+uOT6Eb49fhV9JkXAoeZWR/n3NKobCbbwUz8dnsZPoFIvN86k/izsc+eDxyNP6cZBWwFHIo/XuzMpuP9PPx+eQ3+GDA8YRnzou/cCF8H7Ye/cHY//mLoscAYMzvFOfdc4hdwzn0fDcpxYAbfV8qEzsj0qtyLTVfnPsGfIKV6PUkGLUT4/vcO+FeKz2kKNEz4/3CSrqgkzBsUzZvFz6+ctsX30d0IdE6YPjIqf0nScv5M0pWcaPo50bRSfPec5M9vBrRMMX1I9L4/J00v+y4fA1snTD8xmr4Mf4JemDDvzmjeZi0WaX675fgTpl0TphfhT7IccFZ56yaDzzgvWs7/gIKE6Qcl/IbnpNhuhict51qSrkZVtM6BbtE6nQw0S5hu+ArfkdCSkcH6uyCaf3fSb94Ef/KyAWiXIq7PgDYJ03eK4pqbtPx/ReWHkXB1LprXhujqNb4FYSX+HpV2SeV+Hy3jioRpM/BXJ5NbIizV9piF/WZy9Fs0T1pOWUvJYzW1jir7wl9ocMC0pOm3UM5+BCwEFif8/1dR+d+mKT8hmt8k+v/R0f/vSFP+8Wj+bhXEX626MeH33EDCVX1gO/yV5eX4k85U9c8dSZ83j+y0EFWljuyU4j2dSKpLqGSdx8+vjJ+YML0Af1LugP4J05+LfstULQZbJ09L81stAT5PM2+bsm0oaXrZMWlI0vRx0fR3+fm+W7ZPv55UfrvEbSZhXyxLOgYkzUu1XmfjE4DipOll9eclCdPKft9r03zfcWze4p9RPZnNF/4iw/zoe+2aNO/kKJ67KrsdpNpGM9hHziF7+2wHEo7LCfvCa0AM6Jg0z5GmNR5/AcsBVyVNbx1tJ0uAxine90z0vnaplqvX5i+NMld/7YK/qpDqVdnRRdYlT3DOrXTRlfAMnB39/Ytz7qerps657/BXB8taX4iGOj0Rf+JzX9JybsVXjOk875ybkCLWFS66gpRkJL5p+eA0y7vROfdj4vLxLQPNgSvdz1uCnor+9iwnvjLH4U9A7nfOfZwQZym+ix5s+s2qakj094/OuXjCZ4whsxai6rgAv05/5Zz76cZq52vhsq5bp6V4X8r1h+9GsAz4TeJv7pxbg7+C2pDUrRn/55z7PqH85/grel2iG+sxsyLgl/hWkiujGEl4z/fRegF/o/FWUblvkz7rVvyBJ/l7lZDUiuK88rbjMhnvN5GR+N/ipKTlDE6YXybb6yhjZrYlm044/pA0u2n0d2Wat6/E7zuVKZ9YLtPyzdLMT1bdunGEc+79hPd+g99Gm+G7BqaqfzKpYyqlGnVkpqpa5413zj2bUD6Ob+0G33KSaCObWh1/kvQbphQdd1rhu8Fuxjm3OKpvkt0T/U33+/wxad8djz9J/VnszrlvkreZaF8cVsHyEz2A78Z5fNL08/An7E9ksIyUKllPZtPR+OThxsTtJvrMp4H32LyeqvJ2kKFq77POuQWJx+VoWim+dacA31pZoagnyIXAHOfcTUnLW4K/wNSK1ANmLY7+bpfJZ4kGVajPXnDOpWxmN7Pj8DtqRT7C94m/2vwzb0bhr7h+mFwhVqDseTnjUswbl1SmK767wNTkCtY5t87MZpK+mXd6ugDM7CB8M3RfYGugMGH2tinf5K/MJ35+3MyW4K/4L0wqWzawQLt0MSRI+3s452aa2YqEMlXVE1jqnPsoxbx3gEOqufzy9MOf7B4TdWFItg6/npNttv6ibhrd8fcz/NE2f3RC6+hvquW9n2Ja2ahBzYFV0fu2wl99W5WifKJ+0d8BaYaM3pgUx1P4kYo+NLOn8N1TJqU5sUqlMvsN+Jbff+EToAfhp65dp+O3z8RuZVlbR5VhZg3wgyl0B65xzr2VXCT6m2n9UtPl08lW3TgrxbTvUs1LqH8yqWMqrYp1ZKaqWudVtA+XeQqfCHxoZk9Gn/OOcy7TLlxbR39TdZ0GfhoI5IIozhb8/LEk6X6fdPF3SFp2Af5ewXPwrbdN2bStlrf8RI/jBxj5BdHoYVE91Rd4wqXuFp6pytSTKZl/7l2npMnDnXPzynlbWZ3bPc0Q4Y2Brc2slXPuB6q/HWSi2vtslIBfir84vQt+EJlEme5vu+D3g/lpfp+do79d8d34EpVtD60y/Ky8p4QojznnSs3sQPzNpCfg+3sDLDSzG5xz92e4qKbAeudcqquy3yWUAV/pgr/ansriNNPTzov6rj+JPwF+DX+FruzK7uX4BCyVVPGWljMd/E20FSn7rimvRuJ/kx0zWE55muG7jKWS7nOzpSX+YP6Xcso0STEt1fprES2rI751szLLSzXsb9l6KjvZK2sJSG7xSaVl9PeyDMoC/AN/0LkI37Xmz8D66ED92wxOUCqz3+CcW2Jmo/H96reLrlwOAtoDtyW1aGZzHWUkusr8H/xw7bc6565PUaxsnaVroWnKz9drJuVh0z6baflyh4zOYt1YXl2Sbl4mdUylVKOOzFRV67xM9mGcc0+ZWSn+hvbf4O/tKTWzl4DLXcWjppV918apZprZ7/H782L88PHf4FtdwNdLKX+fxNbXpPiTe9/cCVyM7x72HJsGtmiOr28q/P2dc8vM7FngDDPr6Jybj0+OILpAUg2VqSfTOYef388HPmGZV857yurcinpMNAF+yMJ2kIls7LPP4u+B/QTfcrckKtcJ/10z3d/Kfp+elN9ynKouL9vW16aYJykoIcpzUbPrRdEN4t3xLQuXA/eZ2RLnXCYtTSuBHc2saYqTuzYJZcAfkGHTlf9k25QXbprp1+APeHu4hKHFo6vnV5azvJpS9l3bpJnfhvRdejK1gvS/VbrPzZaV+H7QTVzlhq9Ntf7KfoeJzrkB1Y5sc2VXDjO56l4Wy84ug5vao5aC+4H7oyF9B+G7r5yDP5Adm8HnZbrflBmJTzhOx3eXSNVdrux92VpHFYqGwR6B7w57p3PuijRFy5L4nVIsowl+PU3JpHxkZ2BRQqtcJuXj+PvEypWlujEb4qRPkpqSWV1S03Vkjdd5Ude6Z82sBf4G89Px3ak6m1nvClruluNbeFsmz4gS+T/jk4GeUUtE2bw2lH+hpkLRMi7Cty70d86tS5jXj8wvwIDvNjcYONfMbsR3nf6SDAdOKEdl6smUnHP7V+FtZdvEIc65zQZOSfM51dkOapyZ9cEnQ68BRyV2nTOzU6lcd/my32ekc25IuSU31yL6m+7isyTRPUQC+KZf59wHzrlb8RUM+FFbypRdfS5kc2WjDw1MMW9QUplP8Dct9o0ORD8xs8ZUrf/8jsDHbvPnLPUmzRXBGpb29zCz3fFXBcsdsSkDs4CWZpb8jBfwI4nVpKn47aBPdRcUdc/4BN9lIrlbQTbMxSfh/TNY/tTo796V/ZCoj/1/8S0JnwFHJG/fKVRmvynzP/wQrmdE3TJOwg8ikdzVLWvrqCJRMvQ4vnvIvc65S8spPj76m6pL58H4Y1LiczMm46+kb3aPhZn1xJ9oJ5b/FD9U/r5RfZJYvhV+RLBprpwR5pJlUDfWtOVAm+h3/omZdWTTSU9FKltHllffp1IbdR7gW0qccy86507HD8DQk6Quaine4/BDIO8YdV9L1AqfWE5KTIYi+2Yh5B3wrbVvJiZDVVl+dI/SXPyIc8fhY384RRJQ2fVXmXoym6pT55a3HVT2+2dTWUvoqOT7iEi/vuOkjvVj/Hrpk7z/Z6AL/iJItkeszFlKiPKYme1gZqnuIyi7ypdYeZfdkJvqBr2ym2Cvi+4JKVv+Nvirj6VEN3xGV6ufB7bH99dO9Fs29fWujAXAztHnlX12U/zQmyG8gL+yM9TMfrpSHVVoN0f/fSzVGythRPT3/xIP8NF9AjV5/xD4G41j+Ae/bdYX2szamNmulVjenfguG/dEQ4wmL69b4rqtjOg+tQfwJw7/tKSblMxsm4TE5RH8KFk3m9nOJDGz5mbWO+H/h6Y4SG2B7y9egj/IlSfj/Sbh+6zBJ0W98d1FmrNpW0iU7XWUUrTtPYJPFB7CdwtKyzn3CX5kroPM7KftNEpe/hrF/EhC+RX40ZJ2NrNfJJQvxD8bBRK6C0Unho/gT3D/mPTx1+NbWirsXlTJurGmvYePuywZK7tX65ZKLKOydWR59X0qNVrnmdmB0QWAxGlFbGrxySTBnYDvsp28XhcTtZ4lJtHRfnNjVWNOUNaNq39i/WNmXYCrq7C8B/En/v/G7y/DU5Sp1PqrZD2ZTf/DD7B0lZn1TZ5pZo2jVrSy/2e6HZR1Vw4xoEDZ+v5Z8mNme+OHuU9lKSlijdbLvfjk5oZUSZGZ9Us8fkTTivDHiMnJg3lIeuoyl996As+b2WT81bPF+KtZx+H7nd6bUHYs/gF795rZM9H82c65Uc65ceafAn8R/mbH5/GjYZ2C79Z1ZdKVyavwo6LcaWYH46+C9MJXIBPwzeAVnUwmugt/YJ8R9bEuxncr+pbq9YmuEufcCjO7EH+i+p75e0qW41sPuuNvfqxuQjQcOBO/rqaZv7ekHf4q/Sh8k32NcM7NNrNf43/3T83sFXw/8Rb4LkkD8PeufJx2IT83DNgH/332N7O38H3st8UPf9wL/zyaqt7f8hf8tnUhsJ+ZvYHfvroAh+FPcpc75xab2WD8jbtzou/1Gb5/dmc2Pa/pwmi5/wVWm9k7+HsDtsD/7tsCN6S4OvgzVdhvyozEd5u5Nvr/ZqNL1cA6Suev+PW2HL+vXWObD4yRfGP1hfjnkL0U7RuL8d0LuwDXOefmJr3/Svwwwveb2WH4LkIH44f2fjRFV5ub8C04fzazPfGtqf3wIzuNISHhKkdl6saadje+G+bDURK5HF9/rsK3hmWisnVk2vo+1cJroc77F9DezMbht+NC/DbQHT8qWCb3Tb6Af4bPwfhBM8pij5vZvfh7Ut43s1H4E+yy58jsUo24cc59G+3bx+Pr6rH4uvoY/JDkJ1ZykY8C/xctY1R0L2Gyufht4zQzW0M0UIVLGqksSUb1ZCVjLZdzboP5wSxeBSZHn/kR/ty0E77OnQQcHr0lo+3AObfazKYBg8zsQfyz4hwwLM19X9k0BT84zWlm1hb/7KDObHpOUqr1PRY42fzgPB/gE90nonui/ooftfAq4Hgzm4BPoNpH08uevZZ4r9A++JbfF5DMuTow9rdemb/YNL7+/8opcxwpnkHA5s8hao8/eZiCP+Cvx59sPErSMwGi8lfjK5aNbP4cCsM/jOw9/I65Gp/cpHveyE74lqKV0esN/EnIS9GyWySUPSeadk6aZRl+6OaP8Ff6FuIP/lslf+eo/HDSP2Njs/JJv/vwVDGkiWt//AFvOZueav0HMngyeYbL3xJ/gFgUfe/38ZXtZr9XuvipwnOIEub3B56OPr8En8RMxlfgHTJdfwnlBuMPDMvw3SoXRL/fRSQ8I6SC9ZdyHv7gcDX+WR7r8PdgzSRqNUgqu1u0nIXR9/oB/8yhm4CuCeUuirbX+dH6XYw/gTq1EuuwUvtN9J4i/M3rDn/vVXnLz+o6Kuf3Lu+VatvaDX9z+dJofcwAzi7nc9rj66XF0bbxMb5FuTBN+Rb4lsevo/Jf4odwb5Th98q4bky1zZX3e1aw/c4jdf1zKP4ka0O07u+iEvUblawjo/ekrO8ppy4kwzqPcp6Tk2oe/kLPf6N41gI/4rtbXUCGz8eJfoPP8COcJs9riN8nPo/i/gJfNzYkxTNiSPEcn/LmRb/zbWyqK+bgRyHbIdVvmW6dJMx/MXrf8eWU2Qdfl6yOyrqK4qcS9WQ2X/gWr7KHXG/AHwNm4/fhPlXZDvAtga9H36GsLupUG/ssPnkcjr/YsBZfvw9OtW1H5dvhW8J/xCehP6s38XX+Jfi6e2W0br7Et7CdleK73xdtZxU+D0+vTS+LfjyR4KLuN1/ih72uUhcpERGRVMwPkHEXfvCED0LHUxXRcfILfIt0e+fcZs/kkfxl/hl884HnnHO/DB1PfaJ7iCSIVPc14LvGdETNvCIikn3341uB/ho6kGo4Ht9K95CSIUnhV0AjNnWplgzpHiIJ5WMzm4rvwlGAf7hcP3x3kGsDxiUiIjnIObfR/ANEDzKzYle5IemDMrOL8F3LhuK7TYUaNEjqtlX47sdfhw6kvlGXOQnCzG7G34C+Pb7p/zt8f9/rXXYeriYiIpITzGwefiSyj/EPfs7ouT0ikhklRCIiIiIikrfqZZe5Vq1auU6dOoUOQ0QkjLnRyNS7VGtEYBERkZz23nvv/eCca11RuXqZEHXq1Inp05MfzC4ikif239//HTcuZBQiIiJ1mpnNz6ScRpkTEREREZG8pYRIRERERETylhIiERERERHJW0qIREREREQkbykhEhERERGRvKWESERERERE8pYSIhERERERyVtKiEREREREJG8pIRIRERERkbylhEhERERERPKWEiIREREREclbSohERERERCRvKSESEREREZG8pYRIRERERETylhIiERERERHJW0qIREREREQkbxWFDkBERKS2dbpqVI0uf95NR9Xo8kVEJHvUQiQiIiIiInlLCZGIiIiIiOQtJUQiIiIiIpK3lBCJiIiIiEjeUkIkIiIiIiJ5S6PMiYhI/bVhNfwwFxZ/Akui1/oVULxVwqspNGoGHfeFDv2hQNcCRURkEyVEIiJSfzgHi2bCB0/D3FGwbN6meYUNoVUX2GJrWLcMli+A9SthwyrYuMaXadoeepxIV9uWT9z2gIX4FiIiUocoIRIRkbrvxy9g9jMw+2n48TOf/Ox4EPQeAq27QutdoUUnKExzWNuwGua+6t8/6W5eKy5lbrw9/40N4tHYYZTqcCgikrd0BBARkbrr6+kw+hqY/w5g0GkA7PNr2O0YaNwi8+UUbwm7n+xfa37kzzf+jWMLJ/KXBiM5qnAKl278NV+71jX2NUREpO5SQiQiInXPsvkw5jr48Flo0hoOvg56nAzNtqv+sptszYjYIYyIHcJRBZP5e4MHeKXh1fxh4/m8Gu9X/eWLiEi9ooRIRETqjnXLYcKtMOVesEIY+HvY9zI/OEINGBXfm1klnbmrwZ0Ma3g7I0sP4vrSM9lAwxr5PBERqXuUEImISN3w0Yvw0mV+QISep8OBf85Oi1AFvnbbcFLJNVxR9DQXFr3EXgVzuXjjZXzhav6zRUQkPI09KiIiYcXjMPZG+O+Z0HIHuGA8HD+sVpKhMqUUcVPp6ZxV8ge2tpWMbHgjbfmx1j5fRETCUUIkIiLhbFjlE6HxN0OvIXDuq7Btz2DhvB3vyZCSP7Il63i44S00YV2wWEREpHYoIRIRkTCWfgkPHuKHwz78Zjj2LigqDh0Vn7gOXLLxMrrYQu5qcAeFxEKHJCIiNUgJkYiI1L4vx8H9B8CqRXDmc7D3hWB15yGp4+M9+WvpuRxQOItrix4FXOiQRESkhmhQBRERqV1fjIWRJ8HWO8PpT0DLzqEjSumJ2EF0sO+5sOhl5rk2PBQ7KnRIIiJSA5QQiYhI7fluNjx1JrTq4u8Xatw8dETlurn0NDrYYv5U9ARfu214Pd4ndEgiIpJl6jInIiK1Y8XXMPJk/0yhwc/U+WQIwFHAbzZezCy3I7c1uJud7evQIYmISJYpIRIRkZq3bjmMOAlK1sCQZ2p1SO3q2kBDzi/5HetpyN8bPIgRDx2SiIhkkRIiERGpWaUb4Kkh8OPncOoIaNMtdESV9gPNuGHjEPYq+JTBhWNChyMiIlmkhEhERGpOPA7/uwjmTYDjhkHnQaEjqrJn4/sxIdadPxQ9SRuWhg5HRESyRAmRiIjUnPE3w4fPwsHXwu4nh46mmow/lf6CBpRyfYPhoYMREZEsUUIkIiI145v34O1/wu6nwb6Xh44mKxa4NtxWeiKHFU7nsIJpocMREZEsUEIkIiLZt3E9PH8RbNUWjvxHnXroanU9GDuSj+Idub7BI2zF2tDhiIhINSkhEhGR7Bt7A/wwF465Exo1Cx1NVpVSxFUbf0krVnBl0ZOhwxERkWpSQiQiItm1YAq8exfseS7sdFDoaGrEB25HhscO58yiN9nT5oYOR0REqqEodAAiIpJDStb6UeWabw+H/q3Ki+l01agsBlUzbi09mcMKp3FDg4c5quTvxHWNUUSkXlLtLSIi2TPmOlj6BRx7NxRvFTqaGrWWRty08XR2LVjI0QWTQocjIiJVpIRIRESy46sJMOVe6HsB7DAwdDS1YlS8Hx/FO/LbomcoojQqoAHnAAAgAElEQVR0OCIiUgVKiEREpPo2rIYXLoaWneHga0JHU2scBdxSejKdCr7npMK3Q4cjIiJVoIRIRESq7907YPkC31WuYZPQ0dSqt+K9mRHfiUuLnqOYktDhiIhIJSkhEhGR6lm5CN69E7odDx33CR1NAMY/S0+lnS1lcOGY0MGIiEglKSESEZHqGft/ENsIB+VPV7lkk+LdmBDrzsVFL9CEdaHDERGRSlBCJCIiVff9HJg5EvoOhZY7hI4mqFtKT6GVreTcwtdChyIiIpWghEhERKpu9F/98NoDrwgdSXCz3E68EduToUWjYO3S0OGIiEiGlBCJiEjVfPEWfP4mDPw9bNEydDR1wq2lJ7Ml6/wgEyIiUi8oIRIRkcqLx+CNv0LzDr67nAAw13XgxXh/mHIfrPo+dDgiIpIBJUQiIlJ5HzwF38/2AykUFYeOpk75d+lJULoBJt0ZOhQREclAUegARESknilZC2P+xsx4Z44b2QgYFTqiOmW+a+uHIJ/+COz3O2jcInRIIiJSDrUQiYhI5Uy+B1Z9y40bBwMWOpq6acDlULIapj0UOhIREamAEiIREcnc+hUw8Q7Y5Uimul1DR1N3te0BOx0MU+6FjXoukYhIXaaESEREMjflftiwAgb9IXQkdd++l8OaJTDzidCRiIhIOZQQiYhIZjashsl3w86HQbteoaOp+zoNgO32hHfv9KPyiYhInaSESEREMjP9IVi3DAZdGTqS+sHMtxIt+wo+eiF0NCIikoYSIhERqVjJWt/S0fkAaL9X6Gjqj65HwdY7wcTbwLnQ0YiISApKiEREpGIzHvX3w6h1qHIKCmGfS2HRLPhybOhoREQkBSVEIiJSvo3rYeLt0HEAdNwndDT1T8/TYMu28M5toSMREZEUlBCJiEj5Zo6AVYtg4BWhI6mfioqh/8Xw1Xj49v3Q0YiISBIlRCIikl5piW/ZaN8HOu8fOpr6a89zobiZWolEROqgChMiMzvJzJ41s/lmts7M5prZ381sq6RyLczsQTP7wczWmNmbZtYjxfIamdk/zWxRtLxJZjYwm19KRESy5IOnYMVCGHilHzVNqqZRU+hznh9tbulXoaMREZEEmbQQXQHEgD8ChwPDgIuA0WZWAGBmBrwYzf81cCLQABhrZu2TlvcQcD7wV+BoYBHwupnpoRYiInVJrBQm3Arb9oKdDwkdTf3Xd6gfZGHq/aEjERGRBEUZlPl/zrklCf8fb2ZLgUeB/YG3gGOAAcCBzrmxAGY2CfgKuBK4NJrWEzgDOM8590g0bTwwB7g+Wo6IiNQFc57zz9A5daRah7KhaTvodjzMeBz2v9q3GomISHAVthAlJUNlpkV/t4v+HgN8W5YMRe9bAbwEHJvwvmOAjcBTCeVKgSeBw8ysuFLRi4hIzXDO3+/SuivscmToaHLH3hdDySp4f0ToSEREJFLVQRUGRX8/jv52Az5MUW4O0MHMtkwo95Vzbm2Kcg2BnaoYj4iIZNPnb8LiObDvZVCg8XeyZrs9oEN/mHIvxGOhoxEREaqQEJnZdvjubW8656ZHk1sCy1IUXxr9bZFhuZblfO5QM5tuZtOXLEnVaCUiIlnzzm3QdDvoflLoSHLP3hfB8vkw95XQkYiICJVMiKKWnheAUuDcxFmAS/WWFP/PpNxmnHP3O+f2cs7t1bp16wwjFhGRSvt6Osx/x3fvKmoYOprc0/VoaN4BJg8LHYmIiFCJhMjMGuFHkusMHOac+zph9lJSt+6UtQwty7Dc0hTzRESkNk28DRo1gz3PDh1JbioohL4XwPyJ8O3M0NGIiOS9jBIiM2sAPAv0BY50zs1OKjIHf39Qst2ABc651QnldjCzLVKUKwE+zzRwERGpAT98Bh+/DH3Oh+KtKi4vVbPHmdBwS5h8T+hIRETyXiYPZi0ARgIHAcc65yanKPYisJ2ZDUp4X1Pg/0XzEss1AE5OKFcEnAq84ZzbUJUvISIiWfLuHVDYEPpdGDqS3NaoGfQ+Ez58DlYuCh2NiEhey6SF6G58AnMLsMbM9k54lT109UVgEjDCzE4zs8OiaQb8o2xBzrmZ+CG3bzOzX5rZQfght3cArsnatxIRkcpb9R3MehJ6D4Ytda9mjes3FOKlMO3B0JGIiOS1TBKiI6K/f8InPYmvXwI45+LA0cBo4B7geSAGHOCcW5i0vHOBR4AbgFHA9sDhzrkZ1fomIiJSPZOH+RP0fX4dOpL80LIzdD0Kpj8MG9eFjkZEJG8VVVTAOdcpkwU555YC50Wv8sqtA34bvUREpC5Yv8KfmO92rD9Rl9qx90XwycvwwVOw5zmhoxERyUt62p6IiMB7w2HDSv8gVqk9HfeFtj1gyn3gUj2VQkREapoSIhGRfFe6ASbdAzsMgna9Q0eTX8z8ABaLP4J574SORkQkLykhEhHJdx88Bau/gwGXh44kP3U/ERq3hCn3ho5ERCQvKSESEcln8ThMvAPa7g6dDwgdTX5q0Ng/BHfuK7B8QehoRETyjhIiEZF8NvcV+PEz3zpkFjqa/LXXL/zfaQ+FjUNEJA8pIRIRyVfOwcTboHlH2PXY0NHkt+bbQ9ejYcajGoJbRKSWKSESEclX89+Fr6f55w4VVvgUBqlp/S6Adctg9tOhIxERyStKiERE8tXE22GLVtB7SOhIBPwQ3Nt0gyn3awhuEZFapIRIRCQfff8RfPa6H/K5QePQ0QhEQ3BfAN/PhgWTQkcjIpI3lBCJiOSjibdDgybQ5xehI5FEPU6GRs01BLeISC1Sp3ERkRzV6apRKae34wfGFz/NY7FD+dv1aomoUxpuAXucBZPuhhVfQ7P2oSMSEcl5aiESEckzvyh6FYAHS48MHImk1OeXgIPpD4eOREQkLyghEhHJI81YzWmFb/FifB8WsXXocCSVFh1hlyPhveGwcX3oaEREcp4SIhGRPHJW4Rs0sQ3cV3p06FCkPH2Hwtof4cNnQ0ciIpLzlBCJiOSJJqzjvKLXGB3bg0/d9qHDkfLsMBBa7wpT79MQ3CIiNUwJkYhInjizcDQtbDV3lh4fOhSpiBn0PR8WzYKFU0NHIyKS05QQiYjkgcas55dFrzAu1pMP3I6hw5FM9DwNiptpCG4RkRqmhEhEJA+cUTiGVraSO9Q6VH80bAJ7nAkfvwgrF4WORkQkZykhEhHJccWUcEHRKCbGujHDdQkdjlRGn19CPKYhuEVEapAezCoikuNOLRzLNracS2O/Ch1K3kj3UNyqeKBBb3qNv5d93+hGCQ0AmHfTUVlbvohIvlMLkYhIDmvIRi4qeokp8a5Mju8aOhypgkdjh9HaVnJUweTQoYiI5CQlRCIiOeykwrfZ1pZGI8tZ6HCkCt6Jd+fzeDvOLno9dCgiIjlJCZGISI4qopSLi17g/fhOvBPvHjocqTLj0dih9Cr4kl72eehgRERyjhIiEZEcdXzhO7S3H6KR5dQ6VJ89F9uPla4x5xS9FjoUEZGco4RIRCQXxUq5pPAFZsc7MTbeK3Q0Uk1raMwzsUEcWTCF1iwPHY6ISE5RQiQikos+eJJOBd/r3qEc8ljsEBpajDMKx4QORUQkpyghEhHJNaUbYNxNzIp35o34XqGjkSyZ57ZlbKwng4vGQGlJ6HBERHKGEiIRkVzz3qOwYiG3lJ6CWodyy/DY4Wxjy+GjF0KHIiKSM5QQiYjkkpI18PY/oeMAJsR7hI5GsuzteA++jLeFqfeFDkVEJGcoIRIRySVT7oM1i+Ggv6DWodzjKOCx2KHw9TT45r3Q4YiI5AQlRCIiuWLdcph4O+x8KHTYO3Q0UkOeiQ2EhlvClPtDhyIikhOUEImI5IpJd8H65XDgn0NHIjVoNVtArzNgznOweknocERE6j0lRCIiuWD1Eph0D+x2HGzbM3Q0UtP6DoVYCbw3PHQkIiL1nhIiEZFc8M6/oXQdHPCn0JFIbWi1M+x4IEx/CGIbQ0cjIlKvKSESEanvVnwN0x6EnqdD6y6ho5Ha0u9CWLUIPn4xdCQiIvWaEiIRkfpu/D/AxWHQH0JHIrVpp0OgxQ4aXEFEpJqUEImI1GdLPoX3R8Be50KLjqGjkdpUUAB9z4eFk+HbmaGjERGpt5QQiYjUZ2OugwaNYeCVoSOREHoNhgZbwFS1EomIVJUSIhGR+mrBZPjkZdj3MtiydehoJITGzf29Y7OfgTU/hI5GRKReUkIkIlJfjf4rbNkG+l8SOhIJqe9QiG2AGY+GjkREpF5SQiQiUh+t/REWToH9r4KGTUJHIyFt0xV2GATTHoJYaehoRETqHSVEIiL1jXOwbB5svRP0PjN0NFIX9LsAVn7ju1CKiEilKCESEalvVn8PG9fBQddAYYPQ0Uhd0OVwaN5BgyuIiFSBEiIRkfqkZA0sXwDFW8Gu/y90NFJXFBRCn/Nh/kT47sPQ0YiI1CtKiERE6pPJwyBW4h/IaRY6GqlLeg+BosYw9b7QkYiI1CtKiERE6os1P8LE22GLltCoaehopK7ZoiXsfgp88F9YuzR0NCIi9YYSIhGR+uLtf0LJamjRKXQkUlf1uwBK18OMx0JHIiJSbyghEhGpD5Z+BdMehF6DocEWoaORuqpNN+i0n99WNAS3iEhGlBCJiNQHb90ABUVwwB9DRyJ1Xd+hsGIhfPpq6EhEROoFJUQiInXdtzPhw2dg74ugabvQ0Uhdt8uR0Gx7mKLBFUREMqGESESkrnvzGmjcEgZcHjoSqQ8Ki6DPL2DeBPj+o9DRiIjUeUqIRETqss/HwJfjYODvoVGz0NFIfbHH2VDUSA9qFRHJgBIiEZG6Kh73rUPNO/gr/iKZ2qIl9DgJPngK1i0LHY2ISJ2mhEhEpK6a/TR8NxsO/AsUFYeORuqbvhfAxrXw/ojQkYiI1GlKiERE6qLSDX5kuba7Q/eTQkcj9dG2u0OHfWDqAxCPhY5GRKTOUkIkIlIXTXsQViyAQ66DAlXVUkX9hsLy+fDp66EjERGps3SUFRGpa9avhLdvgc4HwI4Hho5G6rOuR8NW7WCqhuAWEUlHCZGISF0z5V5YtxQO+mvoSKS+K2zgB+T4chwsmRs6GhGROkkJkYhIXbJuObx7l3+45nZ7hI5GcsGe50BhsYbgFhFJoyh0ACIi+arTVaM2m/aboqe5rGgFR3ywHx/P2nw+wJNf/gjAaSneL7KZJq2g+4kw8z++1VHPsxIR+Rm1EImI1BHNWcV5ha8xKtaXj13H0OFILuk3FDaugfdHho5ERKTOUUIkIlJHXFD0Mk1Yz22lGmZbsqxdb9i+n+82F4+HjkZEpE5RQiQiUgdszQrOLnyDl+L9+cy1Dx2O5KK+Q2HZV/D56NCRiIjUKRklRGbW3szuNLNJZrbWzJyZdUpRzqV59UoqV2BmV5vZPDNbb2azzOzE7HwlEZH658KilyimhNtLTwgdiuSq3Y6FLdvCFA3BLSKSKNMWop2AU4BlwIQKyg4H+ie9Pk0q8zfgWuAu4AhgMvC0mR2ZYTwiIjljG5ZxZuFono/vx5euXehwJFcVNoC9zoMvxsAPn4WORkSkzsg0IXrbOdfGOXck8HQFZb9xzk1Oeq0tm2lm2wBXADc5525xzo11zl0AjAVuqtK3EBGpxy4ueoEiYtxeenzoUCTX7XUuFDbUENwiIgkySoicc9m8A/MwoCEwImn6CKCHme2Qxc8SEanTtuVHTi98i6djg1jo2oQOR3Ldltv4IbjfH+mfeSUiIjUyqMJFZrYhutfoLTPbL2l+N2AD8HnS9DnR391qICYRkTrpoqIXMRx3lR4XOhTJF/0ujIbgTr4uKSKSn7KdEI0ALgYOBoYCWwNvmdn+CWVaAsudcy7pvUsT5m/GzIaa2XQzm75kyZLsRi0iEkArVnBq4TiejQ3kG1qHDkfyRbte0GEfmHofxGOhoxERCS6rCZFz7kzn3FPOuQnOuRHAAOBb4IaEYgYkJ0Nl08tb9v3Oub2cc3u1bq0TBxGp/84repUGlHJf7OjQoUi+2fsiWL4A5r4SOhIRkeBq9DlEzrlVwCigT8LkpUALM0tOgFokzBcRyW3rVzCkcDSvxPsyz20bOhrJN12PgmYdYPKw0JGIiARXGw9mTW4RmgMUAzsmlSu7d+ijWohJRCSsaQ/S1NYxrPTY0JFIPioohH5DYf5EWDQrdDQiIkEV1eTCzawpcBQwJWHya0AJMBi4LmH6EOBD59xXNRmTiEhwG9fB5GGMj+3OHNcpdDRSD3W6alS1l9GUNkwqLubVe/7MFRsv/Gn6vJuOqvayRUTqk4wTIjM7KfrnntHfI8xsCbDEOTfezK4AdsE/T+hboCP+eUNt8ckPAM65xWb2b+BqM1sFzABOBQ4EdKlURHLf+yNgzRLuKb2w4rIiNWQlTXgmNpDTCsdyE6fzA81ChyQiEkRlWoiSH8h6T/R3PLA/MBc4Pno1A1YCE4FfOOemJr33T8Bq4DJ8wjQXOMU591JlghcRqXdiG2HiHdC+D1M+7xo6Gslzw2OHc3bRaAYXvsntsRNDhyMiEkTGCZFzrqJR4F4CMkponHMx/MhzN1RUVkQkp3z4HKxYAEfcDJ+nGnBTpPZ85bZlTKw3Q4pGMyx2DCU0CB2SiEitq41BFUREBCAeh3f+Da13hS6Hh45GBIBHYofT2lZydMGk0KGIiAShhEhEpLZ8+hos+RgG/AYKVP1K3fBOvDufxrfjvKLXSP2YQBGR3KYjsohIbXAO3vkXNO8A3XWvhtQlxsOxI+heMI9+9knoYEREap0SIhGR2rBgMnw9Dfa5FApr9IkHIpX2fGwAP7im/LKo+sN5i4jUN0qIRERqw+R7oFFz6HVG6EhENrOBhoyIHcwhhTPgh89ChyMiUquUEImI1LRl8+CTl2HPc6Bhk9DRiKT0eOkhbHANYNLdoUMREalVSohERGra1AcAg75DQ0ciktaPNOO52ACY9R9Y80PocEREao0SIhGRmrRhFcx4DLodB822Cx2NSLkejB0Jpeth2kOhQxERqTVKiEREatL7I2HDStj7ktCRiFToC7cd7HwoTHsANq4PHY6ISK1QQiQiUlPiMZhyL7TvC+33DB2NSGb6/wrWLIHZ/w0diYhIrVBCJCJSUz59DZZ9BXtfFDoSkcztMBDa9vCDKzg9qFVEcp8SIhGRmjJ5GDRtD7seEzoSkcyZQf9fw5JP4PM3Q0cjIlLjlBCJiNSERR/AvAnQb6gexCr1T/cTYKt28O6doSMREalxSohERGrC5GHQoAnscVboSEQqr7AB9LsAvhrvk3sRkRymhEhEJNtWfQ8fPgO9zoDGLUJHI1I1e57tk/pJd4WORESkRikhEhHJtukPQ6xEgylI/da4hU+KPnwWli8MHY2ISI1RQiQikk2lJT4h2vkw2HrH0NGIVE//6PlZk+4OG4eISA1SQiQikk0fvwhrFkPf80NHIlJ9zdpDj5NhxqOwdmnoaEREaoQSIhGRbJr6ALTYAXY8KHQkItmx72Wwca3ftkVEcpASIhGRbPluNiycDH1+AQWqXiVHbLMrdDkCptwLJWtCRyMiknU6YouIZMvUB6CoEfQaHDoSkewacDmsWwrvjwgdiYhI1ikhEhHJhnXLYfbT0OMk2KJl6GhEsqvD3rD93vDuXRDbGDoaEZGsUkIkIpINM5/w91n00WAKkqMGXA4rFsCHz4WOREQkq5QQiYhUVzwO0x6E9n2gXa/Q0YjUjJ0Pg9a7wsTbwbnQ0YiIZI0SIhGR6vpyLCz9Qq1DktsKCvyIc4vnwGejQ0cjIpI1SohERKpr2oOwRSvodlzoSERqVo+ToGl7eOffoSMREckaJUQiItWxfAF8+hrscRYUFYeORqRmFTaA/pfAgndhwZTQ0YiIZIUSIhGR6pj+sP+713lh4xCpLXueDY1bwoRbQ0ciIpIVRaEDEBGpyzpdNSrtvGJKeLf4Qd6L92boTR8AH9ReYCKhNGwC/S+Gt26ARbNg256hIxIRqRa1EImIVNGRBVPY2lbxWOzQ0KGI1K4+50NxU7USiUhOUEIkIlJFZxWN5ov4tkyMdwsdikjtatwc+g6Fj16EJXNDRyMiUi1KiEREqqC7fUnvgs8ZETsYp6pU8tHeF0ODxjDhX6EjERGpFh3FRUSq4KzC0ax1xTwbGxg6FJEwmmztBxOZ/TQs/TJ0NCIiVaaESESkkpqzimMK3+V/sX1ZSZPQ4YiE0/9XUFAE79wWOhIRkSpTQiQiUkknF46nkW3UYAoiTbeF3kNg5hOw4pvQ0YiIVIkSIhGRSjDiDCl8k6nxXfjEdQgdjkh4+14GOHj3jtCRiIhUiRIiEZFKGFQwi44Fi3m89JDQoYjUDS06wu6nwnuPwurFoaMREak0JUQiIpVwVuFolrhmvBbvGzoUkbpjwG+hdD1Mujt0JCIilaaESEQkQ9vb9+xfMIsnYgeykaLQ4YjUHa12gu4nwLQHYe3S0NGIiFSKEiIRkQwNKXyTOMYTpQeFDkWk7tnvCihZDZPvCR2JiEilKCESEclAMSWcUjieN+J78T0tQ4cjUve02Q12Oxam3AfrloWORkQkY0qIREQycEzhu7Sw1Twe02AKImkN/D1sWAmT7w0diYhIxpQQiYhUyHFm4Wg+jW/HpPhuoYMRqbva9oCuR8PkYbB+RehoREQyooRIRKQCe9hn7F7wVfQgVgsdjkjdNuhK2LDCd50TEakHlBCJiFTg7KI3WOka81xsv9ChiNR92/aELkf4IbjXrwwdjYhIhTRurIhIOVqzjCMLpvB47BDW0ih0OCI1rtNVo6q9jB62Ly8Vv8o//vY77okd97N58246qtrLFxHJJrUQiYiUY3DRGBpYjMc0mIJIxma7zrwV68X5Ra/QhHWhwxERKZcSIhGRdEpLOKPwLcbGejLPbRs6GpF65Y7SE2hhqzmzcHToUEREyqWESEQknY9eYBtbzqOxw0JHIlLvzHQ7MT62O+cXjaIx60OHIyKSlhIiEZF0pt7Hl/G2jI/vHjoSkXrp9tIT2NpWMaTwzdChiIikpYRIRCSVb2bA19N4PHYITlWlSJXMcF2YEOvO0KKXacSG0OGIiKSko7yISCpT74cGTXgmNih0JCL12u2lJ9DaVjK4cEzoUEREUlJCJCKSbPUS+PBZ6HU6q9gidDQi9dp015V3Y7txQdHLFFMSOhwRkc0oIRIRSTZjOMRKoO/Q0JGI5IQ7YiewjS3ntMKxoUMREdmMEiIRkUSxjTDtYei8P7TeJXQ0Ijlhcnw3psS7clHRi7BRI86JSN2ihEhEJNEnL8Oqb6HvBaEjEckpt5eeQFtbBu8/HjoUEZGfUUIkIpJo8jBo3hG66NlDItn0brwb0+Jd4J1/Q6lGnBORukMJkYhImYXTYOEU2PtiKCgMHY1IjjHuKD0BVn4DM0eGDkZE5CdKiEREyky+G4qbQe8hoSMRyUkT4j2gfR+Y8C8o1YhzIlI3KCESEQFYNh8+egH2OgeKtwwdjUiOMhj0B1ixEGb9J3QwIiKAEiIREW/KfWAFGkxBpKbtdDC02wMm3OJHdRQRCUwJkYjI+hUw4zHodjw02y50NCK5zaJWouUL4IOnQkcjIqKESESEGY9DySrof0noSETyQ5fDYNue8PYtECsNHY2I5DklRCKS32KlMOVe6DgA2vUOHY1IfihrJVr2Fcx+OnQ0IpLnMkqIzKy9md1pZpPMbK2ZOTPrlKJcIzP7p5ktMrN1UfmBKcoVmNnVZjbPzNab2SwzO7H6X0dEpJI+fsHf4K3WIZHatcuR0KYHvP1PiMdCRyMieSzTFqKdgFOAZcCEcso9BJwP/BU4GlgEvG5mvZLK/Q24FrgLOAKYDDxtZkdmHLmISHU5B+/eBS07Q5fDQ0cjkl/MYNCVsPQL+PDZ0NGISB7LNCF62znXxjl3JJCybdvMegJnAL9xzj3gnBuDT6IWANcnlNsGuAK4yTl3i3NurHPuAmAscFM1vouISOUsnALfzogexKoexCK1ruvRsM1uaiUSkaAyOgNwzsUzKHYMsBH4acgY51wp8CRwmJkVR5MPAxoCI5LePwLoYWY7ZBKTiEi1TboLGreAXmeEjkQkPxUU+FaiHz6FOc+HjkZE8lQ2L4l2A75yzq1Nmj4HnwDtlFBuA/B5inIAu2UxJhGR1H78Aj5+GfY8Fxo2CR2NSP7a9Vho3TVqJcrk+quISHZlMyFqib/HKNnShPllf5c751wF5X7GzIaa2XQzm75kyZJqBysieW7i7VDYEPpdGDoSkfxWUAADfw9LPvGDnIiI1LJsJkQGJCc5ZdOrUu5nnHP3O+f2cs7t1bp16yqGKCICrFwEs/4DvQfDVm1CRyMi3Y6HrXeG8WolEpHal82EaCmpW3daJMwv+9vCzJIToORyIiI1Y/LdEC+FfS4NHYmIABQU+laixXPgk5dDRyMieSabCdEcYAcz2yJp+m5ACZvuGZoDFAM7pigH8FEWYxIR+bl1y2D6I9DtBGipMVxE6ozuJ0LLHWH8P/yQ+CIitSSbCdGLQAPg5LIJZlYEnAq84ZzbEE1+DZ8gDU56/xDgQ+fcV1mMSUTk56Y+CCWrYcBvQkciIokKi2DgFfD9bJj7SuhoRCSPFGVa0MxOiv65Z/T3CDNbAixxzo13zs00s6eA28ysAfAVcBGwAwnJj3NusZn9G7jazFYBM/BJ04HAsdX+RiIi6ZSshSnDYOdDoW330NGISLIep8D4m/1rlyP9w1tFRGpYxgkRmz+Q9Z7o73hg/+jf5wL/B9wANAdmAYc752YkvfdPwGrgMqAtMBc4xTn3UiXiERGpnPcfh7U/woDfho5ERFIpLIL9roAXfwWfvQFdDgsdkYjkgYy7zDnnLM1r/4Qy65xzv3XOtXXONXLO9XPOjUuxrJhz7gbnXEfnXLFzbnfn3DPZ+UoiIinENsK7d0KH/tCxf+hoRCSdnqdB8w4w7ibdSyQitSKb9xCJiG0rfpkAACAASURBVNRds5+BFQt175BIXVfYAPb7HXw7Az4fEzoaEckDSohEJPfF4zDxNtimm79/SETqtp5nQLPtYbxaiUSk5ikhEpHc9+mrsOQT3zqkm7RF6r6ihn5//XoafDk2dDQikuOUEIlIbnMOJtwKzTtCt+NDRyMimeo9BJpuB+NuViuRiNQoJUQiktu+GAPfvOfvSSiszMCaIhJUUbFvJVo4Gb56O3Q0IpLDlBCJSO5yzl9dbrY99Dw9dDQiUlm9z4SttvXPJRIRqSG6XCoiuevLcfD1VDjqX/6eBBEJrtNVoypV/pzCQ7h21WOcevUtTHG7Vlh+3k1HVTU0EclTaiESkdzknL+qvFU7fy+CiNRL/4kdyGLXnMuKng0diojkKCVEIpKb5r0DCyb5exCKikNHIyJVtIGG3Fd6NPsUfkQf+yR0OCKSg5QQiUhuGn8zbNkW9jgrdCQiUk0jYwexxDXl10XPhw5FRHKQEiIRyT3z34V5E2Dfy6BBo9DRiEg1raeY+0uPZmDhbPawT0OHIyI5RgmRiOSe8f+AJq1hz3NCRyIiWTIydjA/uq24VK1EIpJlSohEJLcsnOqfbL/PpdBwi9DRiEiWrKURD5Yexf6Fs+hpn4cO5/+3d99xUtT3H8dfn929Ozh6FVQUREBpoqIGURCxRcFu9Bd7wxaNJRrTLNEY1GgsscQuir2iookNNBQFBEFQEATpCNLr3e5+f3/MnBznwR1wd9/Z3ffz8ZjH3M7O7L53vnt3+9nvzHdEJIuoIBKR7DL8DihsAt3P851ERKrYoNThLHN11UskIlVKBZGIZI+542D6+9DjN1BQ13caEalia6jNY8mj6RsfTxf7znccEckSujCriGS00hd5fDLvdvaK1eXgd3ZhzTtbd/FHEckMT6eOYEDiba5IvM6Fxdf4jiMiWUA9RCKSFfa2b+kT/5JHk/1YQ23fcUSkmqymkMeTR3N4fBydbJbvOCKSBVQQiUhWuDLxKj+6ejydOsJ3FBGpZk+ljmSlK+SKxGu+o4hIFlBBJCIZbx+bRu/4RB5J9mMtuu6QSLZbSR2eTB3FkfGx7GGzfccRkQyngkhEMt6ViVdZ4uozKHW47ygiUkOeSB7FKleb3yTe8B1FRDKcCiIRyWj72lR6xSfx72Q/1ql3SCRnrKAug1KHc3TsM9raPN9xRCSDqSASkYx2VeIVFrv6PJs6zHcUEalhjyWPZh356iUSke2igkhEMtf3IzkoPpmHk/3VOySSg5YRfBlybGwkrW2B7zgikqFUEIlI5vr4Nha7BgxW75BIzno02Y9iElwWf9N3FBHJUCqIRCQzzfofzPqUh5LHsp4C32lExJMlNOC5VF9OiP+PVrbIdxwRyUAqiEQkMw0bCHVbMDjV13cSEfHs4WR/UsS5JD7EdxQRyUAqiEQk84S9Qxx0JRvI951GRDz7gUa8mDqEk+OfwPI5vuOISIZRQSQimWfYQKi7A+x7ju8kIhIRDyf7Bz+MuMdvEBHJOCqIRCSzzBoR9A71vBLyavtOIyIRMZ+mvJLqBV8MgpXzfccRkQyigkhEMsvwsHeo+7m+k4hIxDyYOg7SKRhxn+8oIpJBVBCJSOb4fiTM/AR6/la9QyLyM3Ndc9jr/2Dck7BKI86JSOWoIBKRzDFsINRpDvuqd0hENuPgqyFVBKPu951ERDKECiIRyQzfj4KZw4PeofxC32lEJKqatIUup8CYx2HNEt9pRCQDqCASkcwwfCDUaQbdz/OdRESi7uDfQfE6GPWA7yQikgFUEIlI9M0eDd8NU++QiFROs/bQ6QT4/BFYu9R3GhGJOBVEIhJ9w9Q7JCJbqdfvoGg1fPaw7yQiEnEqiEQk2mZ/Bt99DAdeAfl1fKcRkUyxQyfYsz+MfhjWr/CdRkQiTAWRiETb8IFQ2BT2O993EhHJNL2uhQ0r4LNHfCcRkQhTQSQi0TXnc5jxEfRU75CIbIOWe0H7X8LoB2DDKt9pRCSiVBCJSHQNGwiFTWC/C3wnEZFM1ftaWLcMxjzmO4mIRJQKIhGJpjljYMaHOndIRLbPTvvC7ofByPthw2rfaUQkglQQiUg0DVfvkIhUkd7Xw9of4fN/+04iIhGkgkhEomfuWJj+ARx4ORTU9Z1GRDJdq/2g/VEw4l5Yt9x3GhGJGBVEIhI9wwZC7caw34W+k4hItujzp2D47VEP+E4iIhGjgkhEomXuOJj+vnqHRKRqtewKHY+H0Q/Cmh99pxGRCFFBJCLRMnwg1G4E+6t3SESqWJ8/QvFaGPFP30lEJEISvgOISHZrff07lV63q81gSMF/uaP4Vzx44yfVmEpEclKzDtD1VPj8UejxG6jXwnciEYkA9RCJSGT8NvEay1xdBqWO8B1FRLJV7+sgnYRP/uE7iYhEhAoiEYmETjaTvvHxPJY8mtUU+o4jItmq8W6w9xkw7ilYPtt3GhGJABVEIhIJv0m8wUpXqN4hEal+va4Di8Hw230nEZEIUEEkIt61s7n8Mj6GJ1NHskq9QyJS3RrsBN3PgwnPw5LpvtOIiGcqiETEu8sSb7DGFfBk8ijfUUQkVxx8NSQK4OO/+U4iIp6pIBIRr3a1hfSPjeLZ1GEsp57vOCKSK+o2hx6XweTXYN4432lExCMVRCLi1SXxISRJ8FjyGN9RRCTXHHgFFDaF//4FnPOdRkQ8UUEkIt7syBJOin/K86k+LKah7zgikmtq1YdDrofvR8C093ynERFPVBCJiDcXJd7CAY8k+/mOIiK5at9zoMnu8P4NkEr6TiMiHqggEhEvmrGM0+LDeDXVi/k09R1HRHJVPA8OuxmWTIPxg3ynEREPVBCJiBcXJIaSIMnDqf6+o4hIrtvjGNilB3z8d9iwyncaEalhCd8BRCT3NGIlZ8Q/YEj6QL53LXzHEZEs0vr6d7Zpu73tKF4vuJF7brmCe5Inb3a9WQM1AIxItlEPkYjUuHMT71HHNvBA8jjfUUREABjv2vF26gAGxN+hGct8xxGRGqSCSERqVH3WcE78PwxN7c90t7PvOCIiP7kjeRoJklyVeMV3FBGpQSqIRKRGnRl/n/q2jgeSx/uOIiKyidluB55NHc6p8WG0tzm+44hIDVFBJCI1pjbrOT8xlI9S3ZjsWvuOIyLyM/clT2Aldbg58TSgi7WK5AIVRCJSY34d/5DGtpp/qXdIRCJqOfW4M3kqPeJT6B8b5TuOiNSAKi2IzOwQM3PlTMvLrNfIzB4zsyVmtsbMPjCzLlWZRUSipYAiLkq8w4hUJ75w7X3HERHZrBdSfZiUbs2f8gZTyHrfcUSkmlVXD9EVQI9S02Eld5iZAUOAo4DLgZOAPOBjM9MZ1iJZ6pT4cJrbcv6VUu+QiERbmhg3Fp9DC1vG5YnXfccRkWpWXdch+to5N3oz9x0LHAQc6pz7GMDMRgEzgesIiikRySIJklyceItx6XaMSnf0HUdEpEJfuPa8lOzN+fGhvJzqzXduR9+RRKSa+DiH6FhgfkkxBOCcWwG8BeiiJCJZ6IT4/9jZlnB/8njAfMcREamUO5KnsZ4CbtIACyJZrboKosFmljKzH83sOTPbpdR9nYCvytlmMrCLmdWtpkwi4kGMNJfG3+SrdGuGpbv5jiMiUmlLaMDdyZPpFZ/EkbExvuOISDWp6oJoBXAXcAFwKHALwflDo8ysebhOYyj3EtBLw3mj8h7YzAaY2VgzG7t48eKqTS0i1aZfbDRtYovUOyQiGemZ1OF8nW7FX/KepRYbfMcRkWpQpQWRc268c+53zrm3nHPDnXP3EAyesAMbzw0yyu933uInJefcI8657s657s2aNavK2CJSXdJpLku8wbT0Tvw33d13GhGRrZYizo3F57CzLeHSxJu+44hINaj2c4icc18A04D9wkVLCXqJyirpGSqv90hEMtE3b9MhNpcHksfhdNkzEclQn7s9eT3Vk4vjb8GiKb7jiEgVq6lPKKV7hSYTnEdUVkdgtnNudQ1lEpHqlE7D8NuZkW7J2+kevtOIiGyXW4rPZBWFMOQ3kE75jiMiVajaCyIz6w60Bz4LFw0BdjKz3qXWqQ/0D+8TkWww9R1Y9BX/Sh5PirjvNCIi22Up9bm5+GyYNw5GP+g7johUoSotiMxssJndamYnmtmhZnYN8B4wD7g/XG0IMAp41sxOM7Mjw2UG3FGVeUTEk3Qaht0OjdsyJH2g7zQiIlViSLoHdDgaProVfpzhO46IVJGq7iH6iuA6Q08C/wGuBF4DDnDOLQFwzqWBfsD7wIPA60AK6OOcm1PFeUTEh6lDYdEk6HWteodEJIsYHHM3xAtgyBXBlz8ikvGqepS5vzvnujrnGjjn8pxzrZxzA5xzC8qst9Q5d55zrrFzrtA519c592VVZhERT5yD4QOh8W7Q5RTfaUREqlb9lnDkrfD9/2Dck77TiEgV0LBPIlK1pg6FhZOg13UQT/hOIyJS9fY+E3Y7BN6/EVbM9Z1GRLaTCiIRqTrOwTD1DolIljOD/veCS8FbVwZ/+0QkY6kgEpGqM/VdWDgRel2r3iERyW6NWkPfG2H6+zBhsO80IrIdVBCJSNVwDob9HRq1gS6/8p1GRKT67T8AWh8MQ6/TqHMiGUwFkYhUjWnvqXdIRHJLLAYn/BsS+fDKeZAs8p1IRLaBCiIR2X7p9Mbeoa6n+k4jIlJzGuwEx/4LFkyAj27xnUZEtoEKIhHZflNehwVfwiHXq3dIRHLPnv2g+3kw8j6Y8ZHvNCKylVQQicj2SRbBh7fADp01spyI5K4j/gbN9oDXL4Y1S3ynEZGtoIJIRLbPF0/DspnBaEuxuO80IiJ+5BfCSY/DuuXwxqUailskg6ggEpFtt2E1DL8ddu0J7Q73nUZExK8WneGIW+Db/8Dnj/hOIyKVpIJIRLbd6AdhzWI47ObgQoUiIrlu/wHQ/ij4759h7ljfaUSkElQQici2WbMERtwHe/SDVvv5TiMiEg1mcPxDUK8lvHgGrFrkO5GIVEAFkYhsm0/+AcVroO8NvpOIiERLYWM47TlYvwJeOlPXJxKJOBVEIrL1ls2CMY9Bt9OhWQffaUREoqdFZzjuAZjzGbx7ne80IrIFumCIiGy9j28LRpQ75A++k4iIRFfnE4NrtI24B1ruBd3P9Z1IRMqhHiIR2ToLJ8HEl+CAi4IrtIuIyOb1vQHa9oWh18Lsz3ynEZFyqCASkcpzLvinXrshHHSV7zQiItEXi8PJj0ODnYPziVbO951IRMpQQSQilffl8zB7VDDMdu1GvtOIiGSG2o3g/54Prt323KmwYZXvRCJSigoiEamcdcvgv3+BnfeDvc/0nUZEJLM03xNOeRIWTYaXzoJUse9EIhJSQSQilfPhLbBuKRxzN8T0p0NEZKu1PxL63wMzPoIhlweHIYuIdxplTkQqNm8cjH0CDrgYWnb1nUZEJHPtcxasXADDboP6O+pabiIRoIJIRLYsnYK3r4a6O0CfP/pOIyKS+XpfByvnwad3Qb2WsP+FvhOJ5DQVRCKyZWOfgAUT4KTHoVZ932lERDKfWXD48epFwcid9VrCnv18pxLJWToRQEQ2b/UPwblDbXpD55N8pxERyR7xBJz8BOy0L7x6Psz81HcikZylHiIR2bz//gWK18IxdwXfaIqI5LjW179TpY/XiAt4Mf8WdnrqJOqc9ybs2qNKH19EKqYeIhEp39dvwcQXoOdvoWk732lERLLSMupzetGfWOgaw+BTYO5Y35FEco4KIhH5ueVz4M3LYMe9offvfacREclqi2nIr4v+BHWawjMnwvzxviOJ5BQVRCKyqVQSXr0A0ulgIIVEvu9EIiJZbxGN4ey3oHYDGHQ8LJzkO5JIzlBBJCKbGn47zBkN/f4JTdr6TiMikjsatgqKovw6MOg4WDTFdyKRnKCCSEQ2mvkJfHIndDsDup7iO42ISO5p1DooimJ58HQ/WDDRdyKRrKdR5kRyXMmISY1ZydCCP7DGtaDf6ENZN7pqR1ISEZFKatIWzh0KTx8bFEWnvwKt9vedSiRrqYdIRADHnXn/phGruLz4ctZRy3cgEZHc1qQtnPcuFDYJzin6brjvRCJZSwWRiHBx/C36xsfzt+TpTHGtfccRERGAhrvAue8G88GnwLT/+E4kkpVUEInkuFPiw7g+7wXeTv2CQakjfMcREZHS6rUIDp/boSO88Gv46jXfiUSyjgoikVw2+Q0GJh7lk1QXri6+BDDfiUREpKzCxnDWENh5P3j1fBj7pO9EIllFgyqI5KrpH8CrFzDeteOi4qsoIs93IhGRnFYyyM3m1GIAD+Wtpc/bV3LvG8P5Z/JktuaLrFkDj9nOhCLZST1EIrlo9mfw4pnQbA/OK7pWgyiIiGSA9RRwYfE1vJg8hN8mXueOxCMkSPqOJZLxVBCJ5JqFk4KTc+u1hDNfYyV1fCcSEZFKSpLg98kLuTd5Ir9KDOfRvLsoZL3vWCIZTQWRSC6ZOw6eOQEK6sJZb0Ld5r4TiYjIVjP+mTyZ64svoFdsIi/k30JTVvgOJZKxVBCJ5IpxT8GTR0GiNpz5BjRs5TuRiIhshxdSh3Jh8TW0s3m8ln8D7Wyu70giGUkFkUi2K14PQy6Ht34LrQ+Ci4ZDs/a+U4mISBX4KL0P/1f0Z2pbEa/l38ihsS98RxLJOCqIRLLZ8jlBr9AXg+Dga+D0V4LhW0VEJGtMcLtz7IZbmeV24LG8uxgQfwtwvmOJZAwVRCLZ6tsP4JHesGQ6nDoY+t4AsbjvVCIiUg0W0IRTim5kaPoA/pj3PHflPUQBRb5jiWQEXYdIJNv88A28fwN8+x9o2gFOGwxN2/lOJSIi1Ww9Bfym+HKmpnfmmrxXaGMLuajoKhbTyHc0kUhTD5FItli1EIZcAQ/1gNmj4bCb4aJPVAyJiOQU4/7UiVxcdCV72ByGFvyRHrHJvkOJRJp6iEQy3dql8PkjMOI+SG2A/S+CXtdCnSa+k4mIiCfvpffnu6KWPJh3L4PzbuPe5ImQPkqHTouUQwWRSCZKJWHGRzBhMEwdCqki6Hgc9L0RmrT1nU5ERCJgmmvFsUW3ckveE1yV9yo8sxhOfAzq7eA7mkikqCASyRTOwQ9fw8QX4csXYPVCKGwC3c+Hvc+AFp19JxQRkYhZSy2uKb6E0emO3DnnGXj4IDjpUdjtEN/RRCJDBZFIlK3+Ab4bBjM+Duar5oPFod0R0O3X0P4oSOT7TikiIpFmvJw6hDsvPBdePhsGHQ8HXBSMPppfx3c4Ee9UEIlERao46AFaMAHmTwgGRvghPBG2diNo0zv4Rq/D0TrcQUREtt4OHWHAMPjgJvjsYZj6Lhx7n3qLJOepIBLxIVkEi78OCp+SAmjR5GBQBICC+rDj3sE5QW37QIuuOhFWRES2X34dOPpO6Hg8DLkcBh0H+5wNR9wCtRr4TifihQoikeqWLIIfpmwsfBaUFD/hBfMKGkDLrnDAAGjZLSiEGrWBmEbFFxGRatK6J1wyAj6+DUb9C759H46+A/boB2a+04nUKBVEIlUplQyKn3njNhZAP0zZtPjZcS844GLYsVtQAFWi+Gl9/Ts1EF5ERHJKXu2gZ6jT8fDm5fDiGdDqF8GyVvv7TidSY1QQiWwr52DFXJg3lkeef5m9YjPoYjMptOCwtxWukEnpNnzljmRSug2TXBtmr28OKwy+LnmQb8JJRETEk532DS7kPX4QDBsIjx8Oe/aHvjdB0919pxOpdiqIRCpr/UqY/wXMHRv0AM0bB6sXAXB2PI/JbldeSPVhQrotX7q2fO92AHTYgYiIZIB4ArqfB11+BaMegJH3wTdDYd+z4aCroOEuvhOKVBsVRJITtuWQszqsY7/YNxwYm0LP2FfsabOJmQNgRrolE1w7vkwfxYT07nztdqVYv04iIpLpCurCIb+H7ufC8Nth3FPBtGd/+MWl0OoAnWMkWUef4ERCBRSxT+xbDoxN5sDYZPayGSQszQaXYFy6PfemTmS8250v021ZQV3fcUVERKpP3eZwzF3Q80oY82hQFE15Mxj45xeXBqPU6Tp4kiVUEEnOSpCkq31Hj9gUDoxNpntsGgVWTNLFmOh24+FUf0amOzEu3Z4N6I++iIjkoIat4PC/Qu/fw5fPw+iH4bUL4b0/QKcToPNJQa+RRkaVDKaCSHKGkWZPm02P2GR6xiazf+wb6tp6AKakd+WZ1GGMTHfi8/QerKbQc1oREZEIya8D+10A+54HMz6C8c8E05hHof7O0Dksjlp20yF1knFUEEn2cg5+nA7fDePBvJfoEZtCI1sNBOcAvZ46iJHpToxO78ky6nsOKyIikgFiMWh3WDBtWAVT34OvXg16jkbeD/V2hLaHwu6Hwm59oLCx78QiFVJBJNll+WyY+cnGadUCAPaKNeGD1D6MTHdiVLojC2niOaiIiEiGK6gHXU8JprVL4Zt3YPr78M1bMOFZwIJzjtoeCrv0gFb7Qa0GvlOL/IwKIslsq3/YtABaNjNYXtgU2vT6aep559doCGwREZFqUtiY1i81Bk4lxinsZTM4ODaJXnMn0m3e3SQsTdoZU93OjEu3Z0y6A+NdO2a75lTm//OsgcdU+0uQ3KWCSDLLinkwe1QwzRoBi8MrnBbUh9YHwQEXQZve0HzPMscw6+KnIiIiNSFNjPGuHeNT7bgvdSKFrKdbbDrdbRrdY1M5Lj6SMxIfArDc1QkvYt6GiendmOTaMNc1Q19iSk1SQSTRVbwefpgM88fD7M9g9mhYMTu4L78utNof9jo16AVqsVdwUTkREREp17Zck68qrKUWI9OdGUlnSEGMNB1sDnvFZtDVvqNzbCbnx4aSn0gBsNLVZpprxTfpVnzjdmFquhWsOxBqN/KSX7Kft0+QZtYK+CdwOMHXAB8AVzrnZvvKJJ44xz5/eIG2Np8OsTl0sZl0ic2knc0lz4I/jj+4hnye7sDY9CGMSXfgm/W7kJoch8kAC8NJREREoi5NjK/drnyd2pUXOBSAfIrpYHPoEptJB5vNHrE59I+P4gwLepK4/a/B4fBN20GT3YOpaTto1Boa7Kxzk2S7eCmIzKwQ+AjYAJwNOOBW4GMz6+qcW+Mjl1Sj5IZggIMV82Dl/KCn58cZsGQaLPmWL2ot/2nVH109vkq34aN0N74Ku9Hnuqao+1xERCQ7FZHHJLcbk1K7lVrqaMlSOsRm89Qx9YKRY5dMh2nvwZrFmz5AQYOgMGrYCurvCHWaQ52m4dQsKKZqN4KCupBXqKHBZRO+eoguBHYDOjjnpgOY2UTgW+Ai4G5PuaQ86VRQ0CTXb5wXrYGi1bBhdTAvWg3rV8K6pcFIM6Xnqxb+/A8XQN0Wwbc7nU/kr6OKmeF25Nv0TsynCSp+REREcp2xgCYsSDeBnmUGVVi3PCiQls+GFXNhxZxgvnwOzPk8+Pyx2YeNBSPk5dcLCqR4PuPnr6WIBMUuTjEJiklQRIIkcYrdxtvF4bIiEiRdnCQJigjmwXbB+kniG2+T4JkLe0I8H+J5kKgVTHm1N52rSPPGV0F0LDC6pBgCcM7NNLMRwHFkQ0HkHLj0xjmuzM/l3Vfm/lQxpIuDeaoYUkWQTpb5uWjj/emy6xWVKWTKFDWl56kt3J9OVvplp52xnDosc/VYTl2Wubosdp1Z4JqwgMbB3AXztetrwZJq2v8iIiKSvWo3hJ27B1N5UsmgKFqzGNYsCebrVwTXTipaHcw3rIYNKyFVzMp588gjSYEVU5d15JMiQZI8kuTFUuSTJJ9iEqTC9Sr/2QiApyuxzuYKpZ/Na0GiduXnsURw/ahYAiwe3o4HU+nbFt+4fJN1Y1v3WjOQr4KoE/BmOcsnA6fUcJZt98HNMOqBnxc5ON/JNmWx4BciUVBqqrXpvFb9jbfj5a9z07sz2EAeRS6PNdQKJhfOqc0qV5uVFOLI/l8cERERibB4Auo2D6ZKOHurB5xwxEmTICiWguIpKJbyLFlqebDs1QH7bfwSO7k+mIrXVXK+HtYvD464Sa4Lbpee1xgr1YtV8nN4u+RnM7hkJDRpW4O5tp85V/Mf3s2sCLjbOXd9meW3Atc7535WqJnZAGBAeLMDMLWa4jVF/RZRovaIDrVFtKg9okXtES1qj+hQW0RLrrXHrs65ZhWt5HOc4vIqsc0ePOmcewR4pPrihAHMxjrnNtP/KjVN7REdaotoUXtEi9ojWtQe0aG2iBa1R/l8Hdu0DGhczvJG4X0iIiIiIiLVzldBNJngPKKyOgJTajiLiIiIiIjkKF8F0RDgF2b202DzZtYa6Bne51O1H5YnW0XtER1qi2hRe0SL2iNa1B7RobaIFrVHOXwNqlAH+BJYB/yZ4HyiW4B6QFfn3OoaDyUiIiIiIjnHSw+Rc24NcCgwDXgGGAzMBA5VMSQiIiIiIjXFSw+RiIiIiIhIFGTlFTTNLGZmfzCzWWa23sy+NLOTtmL7481sfLjt92b2ZzOLl7PeQWY20szWmdlCM7vbzGpv5jHPMrMxZrbWzJaZ2f/MrMv2vM5MEcX2CNfPM7NJZubM7IJtfX2ZJkrtYWYXmtlQM5tnZmvM7Cszu9bM8qvitUaFmbUys1fMbIWZrTSz18xsl0puW8vM7jSzBeG+HGVmvcpZr9LtGu73b8xsg5lNNbOLt/c1ZpKotIeZtTSzv5vZ2DDLYjP7sLzHy2ZRaY9ytjnQzNLh/wiflympUVFrDzNrZGb3mNns8G/WXDN7ajtfZkaIUluYWaGZ3Wxm08LHm2NmgywYAyDzOeeybgL+BmwAfgf0Af4NpIGjK7HtkUCK4KSzPsDVwHrg9jLrdSU4B+oNoC9wAcGQ4S+W85i3AWuBm8LHPAa4GTjA977KxfYotc0fgAUE57Bd4Hs/5WJ7AHOAJ4CTw8f7S/h4L/veT1W4vwuBb4GvgOOB44BJwAygTiW2HwwsBy4M9+Vr4b7tti3tGj5OOly/r66oBwAACMxJREFUD3BrePsS3/sq19oD6Ad8F77vDw//N7wTrtfP977KtfYos35emKPkf0TC977KxfYguBzL5DDPWUAv4DTgft/7Kgfb4jmCz7LXhuudDcwK89T1vb+2e3/7DlANb6DmYcPeXGb5h8DESmw/HhheZtkNQBHQotSy18M3al6pZWeFfzj3KbWsR/jGOt73vlF7bPIYuwFrgNPJoYIoau0BNCvnOW4I19vN9/6qon3+W4IicvdSy9oASeDqCrbdK9wX55ZalgCmAkO2tl3DbX8Ani6z3hMEVy7P25rXlolTxNqjIWU+aJd6vE9876tca48y9/2R4IPo38itgihS7QE8DHwP1Pe9b3K5LYDa4fPeVma9o8LnOdL3/treKRsPmTsSyAeeLbP8WaCLmbXZ3IZm1groVs62zxB8W/TLcL08gjfBS8654lLrvUTwwfC4UssuAWY6597Y+peSFaLWHiUeAl4ARlT6lWSHSLWHc25xOU81JpzvVNGLyRDHAqOdc9NLFjjnZhK898p7b5bdthh4sdS2SYL37pFmVhAurmy79gCalbPeM0AT4KBKvqZMFpn2cM4tD7f/SXh7Atnz/q9IZNqjhJm1Bf4EXBo+fi6JTHtYMCLxWcBjzrmV2/yKMldk2oKgmIoDZdtheTjP+Hoi419AOToRVLvTyyyfHM47VrAtBN8K/SR8A64ttW1boFY5660n6Dos/RwHAV+a2XUWnCeRtOA8iVMq+XoyXdTaAzM7HegO/L5SryC7RK49ytGboFd1WgXrZYpOlNkXoclUvC86EXyhsracbfOB3UutV5l2LbcNy1kvm0WpPX7GgvPnegBfV5AlW0SxPR4CXnHOfVLB82ejKLXHvgQ9E4vC82jWmdlqM3tjS1/eZZHItIVzbhXBF2dXmFkfM6trZp2AOwkuo/NhxS8n2rKxIGoMLHdhX14pS0vdv6VtITjXoaxlpe7f0npLyzzHjsBhBD1F1xJ8i/418JKZVVThZ4NItYeZNQLuBn7vnFuyhefOVpFqj7LMrCvBYQJPOOcWbSFLJmnM5vdFo+3YtuT+knll2nVzbVOZ9s8WUWqP8twE7AzcXkGWbBGp9jCzMwi+MLu2gufOVlFqjx3D+T8IDh07FhgA7A0MM7N6FeTJdFFqC4BzCQ6H/whYRVCs5QGHO+eKKsgTeZEviMzssHCEl4qmYSWbEBzP+LOHqszThfOKtq/sehDs43rASc6555xz7wOnElTff6xEpkjJgva4k6CX4vFKPH/kZUF7bLzDrCXwJkH7XF2JPJlke/Z5ZbbdmvU2lyeXRKU9Nr3T7NfA9cAtzrlPK5EnW0SiPcysMXAX8Efn3A+VeP5sFYn2YONn1JnAac65951zzwG/AnYBzqhEpkwXlbaAYACeMwgGX+gNnElwqPW74eGNGS0ThpEcCexZifVKugWXAo3MzMpUvI1K3b85W/r2rmGp+7e0XsmIKCV+BPKdc1+ULHDOpc3sQyATh7rN2PYwswOAcwhGW2lgZgD1w/Vqm1lDYEU535REWca2R2lm1gR4n+CP8JFh93y2KN17Vlojyv8Gr7SlBP/4y9u25P6SeWXatXTbLCi1XuMy92ezKLXHT8ysP/AU8Lhz7sYKcmSTKLXHrcAigiM4GobLaoXzBma23gUXls9mUWqPH8P5B6XXc859ZmYrCXqKsllk2iI8PO56ggGofvpC2cw+Izi8/QLg3goyRVrkC6Lw+MdvtmKTyUABwXkMpY+JLDneckoF20JwTOWokoXhGOuFpbadQXDMZadS22JmtQhGL3u5zGN2K+e5NleVR1qGt8eeBCcFDivnue4Lp0ZsPEkw8jK8PUqW1wf+Q/BN08HOuXmVfTEZYjJl9kWoI1ve3yXbnmBmhWWOBe9IMEDF9FLrVaZdS7fhgi2sl82i1B4AmFlfgt+L14GLKvEaskmU2qMj0IWNH8RLW0LQg318BZkyXZTao+Tv1eY+K6UryJPpotQWJdfNHFNqHZxz35rZcir3xWy0VecQdj4mNg4heGOZ5R8Akyqx/QTg4zLL/szPhxV+g6AqTpRadgbBL+6+pZZdFi7rXmpZjOBN9qnv/ZVL7QG0AA4pM50WrnNneDurh1aNUnuEywqBTwk+gHTxvX+qaZ9fSTBc6W6llrUmGAHomgq27Rbus7NLLUsQnIf41ta2K8Hx3ouBJ8us9xgbe7O977NcaY9wWQ9gNcH1h7J+2PMot0f4eGX/RzwVPkdfoLPv/ZVL7REuGwNMBKzUsh7h85zve3/lSluEvwuOMpcoAdqHy//ie39t9/72HaCa3kQDCS7ueHXYiA8RfJPQv8x6HwLTyyw7Olz33+G2V4WPdWc5b7Z1BBe66gucT9C1+HKZ9WoRFD/fExxv+UuCbwFTQF/f+yrX2qOcbK3L+yXP5ilK7QG8Gz7e5cAvykw/u0ZRJk5AHYJv3iYRDJV6LMGoPN9R6mJ2wK4E//xuKLP9CwSHR1wQ7stXwn2+T5n1KtuuF4fLbw3X+2t4+zLf+yrX2gPYI/y9mBWus8nvgO99lWvtsZl8N5Fb1yGKVHuEj5EEXiX4/HQWwQW9vwZq+95fudIWBEfXTABWsOmFWb8lOKpmF9/7a7v3t+8A1fQmihN8a/09QeU7ETi5nPWGAbPKWX5i+KbbAMwmuFBkvJz1ehEcOrSe4Ljje4DCctZrSTCm+9Jw3VHAEb73U662R5ltWpN7BVFk2iPc95ubzvG9r6pwn+9C8A99JcHoPG8ArTfzXrypzPLaBCMjLgz35WfAIdvaruG6FxH04G0g+Id2qe99lIvtQXBO42Z/B3zvp1xrj81ku4kcKoii2B4EhdCY8PF+BAYBO/jeT7nWFgSHtd9F8D9jHUFh+iLQwfd+qorJwhcpIiIiIiKScyI/7LaIiIiIiEh1UUEkIiIiIiI5SwWRiIiIiIjkLBVEIiIiIiKSs1QQiYiIiIhIzlJBJCIiIiIiOUsFkYiIiIiI5CwVRCIiIiIikrP+Hzsq/3RX6fKRAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 1008x504 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
     }
    ],
    "source": [
@@ -204,11 +276,12 @@
     "# Synthetic data 2 - with unobservable Z\n",
     "\n",
     "# R ~ U(0,1)\n",
-    "# X ~ N(0,1)\n",
-    "# T ~ Bin(1, sigmoid(a * r + b * x))\n",
-    "# Y ~ Bin(1, sigmoid(c * t + d * x))\n",
+    "# X ~ U(0,1)\n",
+    "# Z ~ U(0,1)\n",
+    "# T ~ Bin(1, sigmoid(a * r + b * x + e * z))\n",
+    "# Y ~ Bin(1, sigmoid(c * t + d * x + f * z))\n",
     "\n",
-    "# Weights:\n",
+    "# Weights for edges:\n",
     "# R -> T: a\n",
     "# X -> T: b\n",
     "# T -> Y: c\n",
@@ -219,22 +292,23 @@
     "N = 12000\n",
     "a = 1\n",
     "b = 1\n",
-    "c = .5\n",
+    "c = 1\n",
     "d = 1\n",
     "e = 1\n",
     "f = 1\n",
     "\n",
+    "\n",
     "def generateDataWithZ(N=N, a=a, b=b, c=c, d=d, e=e, f=f):\n",
     "\n",
     "    r = npr.uniform(size=N)\n",
     "    x = npr.uniform(size=N)\n",
-    "    \n",
     "    z = npr.uniform(size=N)\n",
     "\n",
     "    t = npr.binomial(n=1, p=sigmoid(a * r + b * x + e * z), size=N)\n",
     "\n",
     "    y = npr.binomial(n=1, p=sigmoid(c * t + d * x + f * z), size=N)\n",
     "\n",
+    "    # If decision is negative, we set Y to positive.\n",
     "    y[t == 0] = 1\n",
     "\n",
     "    return r, x, z, t, y\n",
@@ -243,66 +317,88 @@
     "def analytic_R_on_Y(r, a, b, c, d, e, f):\n",
     "    '''Return probability of Y=0 with given parameters.'''\n",
     "\n",
-    "    # P(y=0|t=1, x) * p(t=1|r, x) * p(x)\n",
-    "    def integrable(t, x, z):\n",
+    "    # P(y=0|T=t, x, z) * p(T=t|r, x, z) * p(x) * p(z)\n",
+    "    def integrand(t, x, z):\n",
     "        p_y0 = 1 - sigmoid(c * t + d * x + f * z)\n",
-    "        \n",
+    "\n",
     "        if t == 1:\n",
     "            p_t = sigmoid(a * r + b * x + e * z)\n",
     "        elif t == 0:\n",
     "            p_t = 1 - sigmoid(a * r + b * x + e * z)\n",
-    "        \n",
+    "\n",
     "        p_x = 1\n",
     "        \n",
-    "        return p_y0 * p_t * p_x\n",
-    "    \n",
-    "    t_1 = si.dblquad(lambda x, z: integrable(1, x, z), 0, 1, lambda x: 0, lambda x: 1)\n",
+    "        p_z = 1\n",
+    "\n",
+    "        return p_y0 * p_t * p_x * p_z\n",
+    "\n",
+    "    t_1 = si.dblquad(lambda x, z: integrand(1, x, z), 0,\n",
+    "                     1, lambda x: 0, lambda x: 1)\n",
     "\n",
     "    # P(y=0|t=0, x) * p(t=0|r, x) * p(x)\n",
     "    # Now equal to 0, hence commented out. See data generation.\n",
-    "    t_0 = si.dblquad(lambda x, z: integrable(0, x, z), 0, 1, lambda x: 0, lambda x: 1)\n",
+    "    t_0 = si.dblquad(lambda x, z: integrand(0, x, z), 0,\n",
+    "                     1, lambda x: 0, lambda x: 1)\n",
     "\n",
-    "    return t_1[0]  #+ t_0[0]\n",
+    "    return t_1[0]  # + t_0[0]\n",
     "\n",
     "\n",
-    "r, x, z, t, y = generateDataWithZ()\n",
+    "diffs = np.zeros(0)\n",
     "\n",
-    "# Fit predictive models.\n",
-    "lr_t = LogisticRegression(solver='lbfgs')\n",
-    "lr_y = LogisticRegression(solver='lbfgs')\n",
+    "for i in range(2000):\n",
     "\n",
-    "lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
-    "lr_y = lr_y.fit(np.array([t[t == 1], x[t == 1]]).T, y[t == 1])\n",
+    "    r, x, z, t, y = generateDataWithZ()\n",
     "\n",
+    "    # Fit predictive models.\n",
+    "    lr_t = LogisticRegression(solver='lbfgs')\n",
+    "    lr_y = LogisticRegression(solver='lbfgs')\n",
     "\n",
-    "def causal_effect_of_R_on_Y(r):\n",
+    "    lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
+    "    lr_y = lr_y.fit(np.array([t, x]).T, y)\n",
     "\n",
-    "    # Equal to integration of P(Y=0|T=0, X=x) * P(T=0|R=r, X=x) * f(x)\n",
-    "    # from 0 to 1\n",
-    "    print(\"Bias:\",\n",
-    "        si.quad(\n",
-    "            lambda x: lr_y.predict_proba(np.array([[0, x]]))[0, 0] * lr_t.\n",
-    "            predict_proba(np.array([[r, x]]))[0, 0], 0, 1))\n",
+    "    def causal_effect_of_R_on_Y(r):\n",
     "\n",
-    "    # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
-    "    return (si.quad(\n",
-    "        lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
-    "        predict_proba(np.array([[r, x]]))[0, 1], 0, 1))\n",
+    "        # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
+    "        return (si.quad(\n",
+    "            lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
+    "            predict_proba(np.array([[r, x]]))[0, 1], 0, 1))\n",
     "\n",
+    "    r0 = causal_effect_of_R_on_Y(0)\n",
+    "    r1 = causal_effect_of_R_on_Y(1)\n",
     "\n",
-    "r0 = causal_effect_of_R_on_Y(0)\n",
-    "r1 = causal_effect_of_R_on_Y(1)\n",
+    "    analytical = analytic_R_on_Y(1, a, b, c, d, e, f) - analytic_R_on_Y(\n",
+    "        0, a, b, c, d, e, f)\n",
     "\n",
-    "analytical = analytic_R_on_Y(1, a, b, c, d, e, f) - analytic_R_on_Y(0, a, b, c, d, e, f)\n",
+    "    diffs = np.append(diffs, analytical - r1[0] + r0[0])\n",
     "\n",
+    "    if i % 200 == 0:\n",
+    "        print(i / 20, \"%\", end=\" \")\n",
+    "\n",
+    "print()\n",
     "print(\"Analytical:\", analytical)\n",
     "print(\"Estimated: \", r1[0] - r0[0])\n",
-    "print(\"Difference:\", analytical - r1[0] + r0[0])\n",
     "print()\n",
     "print(\"Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\\n\")\n",
     "print(\"Analytical:\", analytic_R_on_Y(1, a, b, c, d, e, f),\n",
     "      analytic_R_on_Y(0, a, b, c, d, e, f))\n",
-    "print(\"Estimated: \", r1[0], r0[0])"
+    "print(\"Estimated: \", r1[0], r0[0])\n",
+    "print()\n",
+    "print(\"Average difference:\", np.mean(diffs))\n",
+    "print(\"Std of differences:\", np.std(diffs))\n",
+    "plt.hist(diffs, bins=25)\n",
+    "plt.title(\"Histogram of differences over 2000 simulations (analytic - estimate)\")\n",
+    "sns.kdeplot(diffs)\n",
+    "plt.axvline(x=0, c='r')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now it is visible that the quantity $P(Y=0|do(R=r))$ is slightly biased when estimated with equation \\ref{int}.\n",
+    "\n",
+    "\n",
+    "According to [Pearl's Theorem 3](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2836213/): \"A *sufficient* condition for identifying the causal effect P(y|do(x)) is that every path between X and any of its children traces at least one arrow emanating from a measured variable.\""
    ]
   }
  ],
-- 
GitLab