diff --git a/CompareIncomeDist.ipynb b/CompareIncomeDist.ipynb
new file mode 100644
index 0000000..cb7da69
--- /dev/null
+++ b/CompareIncomeDist.ipynb
@@ -0,0 +1,232 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1afb4980",
+   "metadata": {},
+   "source": [
+    "# Comparing Income Distribution Estimates\n",
+    "\n",
+    "This notebook estimates the income distribution using various distributional assumptions, and compares the resulting estimates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "16ed9b3c",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/var/folders/2m/db9b_l2950xckts1f65_rc2m0000gn/T/ipykernel_99466/1679099148.py:3: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display\n",
+      "  from IPython.core.display import display, HTML\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "        <script type=\"text/javascript\">\n",
+       "        window.PlotlyConfig = {MathJaxConfig: 'local'};\n",
+       "        if (window.MathJax && window.MathJax.Hub && window.MathJax.Hub.Config) {window.MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}\n",
+       "        if (typeof require !== 'undefined') {\n",
+       "        require.undef(\"plotly\");\n",
+       "        requirejs.config({\n",
+       "            paths: {\n",
+       "                'plotly': ['https://cdn.plot.ly/plotly-2.27.0.min']\n",
+       "            }\n",
+       "        });\n",
+       "        require(['plotly'], function(Plotly) {\n",
+       "            window._Plotly = Plotly;\n",
+       "        });\n",
+       "        }\n",
+       "        </script>\n",
+       "        "
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# imports\n",
+    "\n",
+    "from IPython.core.display import display, HTML\n",
+    "from plotly.offline import init_notebook_mode, plot\n",
+    "import plotly.express as px\n",
+    "init_notebook_mode(connected=True)\n",
+    "from iot.generate_data import gen_microdata\n",
+    "from iot.inverse_optimal_tax import IOT\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "9b7c81a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get data\n",
+    "data = gen_microdata()\n",
+    "data = data[data[\"expanded_income\"] > 0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "ac8fa99d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "iot_ln = IOT(data=data, income_measure=\"expanded_income\", dist_type=\"log_normal\") \n",
+    "iot_Pln = IOT(data=data, income_measure=\"expanded_income\", dist_type=\"Pln\")\n",
+    "\n",
+    "df_ln = iot_ln.df()\n",
+    "df_Pln = iot_Pln.df()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "e7dc7302",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.legend.Legend at 0x2aa6bc710>"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHACAYAAACCrj2lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqIUlEQVR4nO3deXhU5d3/8fc5s2VfICRhCYR931FERKGiiEirPrVUfURpa9VKq6Kt0qq4VWir1taqWH8q1qeKS92qqFUEF6QgS5B9kxC2hD37Msv5/TGZSQIBQ0gyM5nP67rOlcyZMzN3ZqTz6ff+nvsYlmVZiIiIiMgpMUM9ABEREZFIpBAlIiIi0ggKUSIiIiKNoBAlIiIi0ggKUSIiIiKNoBAlIiIi0ggKUSIiIiKNoBAlIiIi0ggKUSIiIiKNoBAlIiIi0ghRHaI+//xzJk+eTIcOHTAMg7fffrtZX+++++7DMIw6W58+fZr1NUVERKR5RHWIKi0tZfDgwTz55JMt9pr9+/dn3759we3LL79ssdcWERGRpmMP9QBCaeLEiUycOPGE91dWVvK73/2OV155haNHjzJgwAD+8Ic/MHbs2Ea/pt1uJzMzs9GPFxERkfAQ1ZWo7zJ9+nSWLl3K/Pnz+eabb7jiiiu46KKL2Lp1a6Ofc+vWrXTo0IFu3bpx9dVXk5eX14QjFhERkZZiWJZlhXoQ4cAwDN566y0uvfRSAPLy8ujWrRt5eXl06NAheNz48eM588wzefjhh0/5NT744ANKSkro3bs3+/bt4/7772fPnj2sW7eOxMTEpvpTREREpAVE9XTeyaxduxav10uvXr3q7K+srKRt27YAbNq0ib59+570ee68807mzJkDUGfqcNCgQYwcOZIuXbrw2muv8dOf/rSJ/wIRERFpTgpRJ1BSUoLNZmPlypXYbLY69yUkJADQrVs3Nm7ceNLnCQSu+qSkpNCrVy+2bdt2+gMWERGRFqUQdQJDhw7F6/Wyf/9+xowZU+8xTqfztJYoKCkpYfv27VxzzTWNfg4REREJjagOUSUlJXWqQDt27CAnJ4c2bdrQq1cvrr76aqZOncqjjz7K0KFDOXDgAAsXLmTQoEFMmjTplF/vjjvuYPLkyXTp0oW9e/cya9YsbDYbV155ZVP+WSIiItICorqxfPHixYwbN+64/ddeey3z5s3D7Xbz0EMP8Y9//IM9e/aQlpbGWWedxf3338/AgQNP+fV+/OMf8/nnn3Po0CHatWvHOeecw+9//3u6d+/eFH+OiIiItKCoDlEiIiIijaV1okREREQaQSFKREREpBGirrHc5/Oxd+9eEhMTMQwj1MMRERGRBrAsi+LiYjp06IBphkcNKOpC1N69e8nKygr1MERERKQRdu3aRadOnUI9DCAKQ1Tg8iq7du0iKSkpxKMRERGRhigqKiIrKyusLpMWdSEqMIWXlJSkECUiIhJhwqkVJzwmFUVEREQijEKUiIiISCMoRImIiIg0QtT1RImISGTyer243e5QD0OakdPpDJvlCxpCIUpERMKaZVnk5+dz9OjRUA9FmplpmnTt2hWn0xnqoTSIQpSIiIS1QIBKT08nLi4urM7OkqYTWAx73759dO7cOSI+Z4UoEREJW16vNxig2rZtG+rhSDNr164de/fuxePx4HA4Qj2c7xQ5E48iIhJ1Aj1QcXFxIR6JtITANJ7X6w3xSBpGIUpERMJeJEztyOmLtM9ZIUpERESkERSiREREmtjYsWO59dZbQz2MsLN48WIMw2g1Z1oqRImIiIg0gkKUiIiIBFVVVYV6CBFDIeo7FBevorR0faiHISIiEezIkSNMnTqV1NRU4uLimDhxIlu3bq1zzLPPPktWVhZxcXFcdtllPPbYY6SkpJzwOXNzczEMgzfffJNx48YRFxfH4MGDWbp0aZ3j/vWvf9G/f39cLhfZ2dk8+uijde7Pzs7mwQcfZOrUqSQlJfHzn/+cefPmkZKSwnvvvUfv3r2Ji4vjhz/8IWVlZbz44otkZ2eTmprKr371qzpn0r300kuMGDGCxMREMjMzueqqq9i/f//pv4FhSiHqJHJz72flyuF8/fUAdu36c6iHIyIi+Fcw93pLW3yzLKvRY77uuutYsWIF7777LkuXLsWyLC6++OLgEg5Llizhxhtv5JZbbiEnJ4cLLriA3//+9w167t/97nfccccd5OTk0KtXL6688ko8Hg8AK1eu5Ec/+hE//vGPWbt2Lffddx/33HMP8+bNq/McjzzyCIMHD2b16tXcc889AJSVlfHXv/6V+fPn8+GHH7J48WIuu+wyFixYwIIFC3jppZd45plneOONN4LP43a7efDBB1mzZg1vv/02ubm5XHfddY1+38KdFts8gfLyXHbufCh4+9tvf0Na2g+Ije0WwlGJiIjPV8YXXyS0+OuOGVOCzRZ/yo/bunUr7777LkuWLOHss88G4J///CdZWVm8/fbbXHHFFTzxxBNMnDiRO+64A4BevXrx1Vdf8d57733n899xxx1MmjQJgPvvv5/+/fuzbds2+vTpw2OPPcb5558fDEa9evViw4YN/OlPf6oTbr73ve9x++23B29/8cUXuN1unn76abp37w7AD3/4Q1566SUKCgpISEigX79+jBs3jkWLFjFlyhQAfvKTnwSfo1u3bvz1r3/ljDPOoKSkhISElv/MmpsqUSewf//LWJaHlJRxpKZOwLI87Nr1SKiHJSIiEWbjxo3Y7XZGjhwZ3Ne2bVt69+7Nxo0bAdi8eTNnnnlmnccde/tEBg0aFPy9ffv2AMEptI0bNzJ69Og6x48ePZqtW7fWmYYbMWLEcc8bFxcXDFAAGRkZZGdn1wlDGRkZdabrVq5cyeTJk+ncuTOJiYmcd955AOTl5TXob4k0qkSdwMGDbwOQkXE1MTHZHDnyEQUFL9O9+6PYbLGhHZyISBQzzTjGjCkJyeuGo9qXRwksVunz+U7pOeLjj6+wHXvZFcMw6t0XeK3S0lImTJjAhAkT+Oc//0m7du3Iy8tjwoQJrbZZXSGqHl5vOSUlqwFITb0Al6sTLlcXKit3cvjwB7Rrd3mIRygiEr0Mw2jUtFqo9O3bF4/Hw7Jly4LTeYcOHWLz5s3069cPgN69e/P111/Xedyxtxv72kuWLKmzb8mSJfTq1QubzXbaz1/bpk2bOHToEHPmzCErKwuAFStWNOlrhBtN59WjpGQVluXB6czE5crCMEzatbsMgEOHvnt+WkREJKBnz5784Ac/4Prrr+fLL79kzZo1/O///i8dO3bkBz/4AQC//OUvWbBgAY899hhbt27lmWee4YMPPjjty6DcfvvtLFy4kAcffJAtW7bw4osv8re//S3Ye9WUOnfujNPp5IknnuDbb7/l3Xff5cEHH2zy1wknClH1KCryp//ExDOD/wG3bXsJAIcOvY9lnVqZVEREotsLL7zA8OHDueSSSxg1ahSWZbFgwYLg9Njo0aOZO3cujz32GIMHD+bDDz/ktttuIyYm5rRed9iwYbz22mvMnz+fAQMGcO+99/LAAw80yxlz7dq1Y968ebz++uv069ePOXPm8MgjrbuX2LBO55zNCFRUVERycjKFhYUkJSXVe8yWLTexd+9cOnf+Ld26+U8x9fmqWLKkLV5vCcOHryIxcWhLDltEJCpVVFSwY8cOunbtetqBItJcf/31bNq0iS+++CLUQ2kxJ/u8G/L93dJUiapHWdkWAOLiegX3maaT5ORzACgs/Dwk4xIRkdbrkUceYc2aNWzbto0nnniCF198kWuvvTbUw5KTUIiqR3m5P0TFxvausz852X+q5tGjClEiItK0li9fzgUXXMDAgQOZO3cuf/3rX/nZz34W6mHJSejsvGN4vaVUVu4G6laiAFJSzgX8lSjLsk674U9ERCTgtddeC/UQ5BSpEnWMiopcAOz2VByONnXuS0wcgWnG4nYfpKxsYwhGJyIiIuFClahjBKpQLlen4+4rXl6O7fEH8XVYSmH3ZcTH92vp4YmIiEiYUIg6RmXlXgBcro519hevLibnvBysquHAcHbt20GHf4VggCIiIhIWNJ13jMrKPQA4nR3q7N9++3asKgtnTw+YXsrf7MqhDw6FYogiIiISBhSijlFV5Q9RtStRZVvKOLroKJjQ771suPxNAHLv20GULbMlIiIi1RSijhGoRNUOUQfePABAmwltSO7ZHfvUj8BZSfHyEoq/Lg7JOEVERCS0FKKOUTOdVxOijnx8BIA2F7fBMAySuvSCc/1rReXPy2/5QYqISESbN28eKSkpoR6GnCaFqGNUVe0DwOXy90R5K7wUflkIQJsL/EseJCaeCRd9CMD+V/bjc+taeiIiUtd1112HYRgYhoHT6aRHjx488MADeDyeUA9NmojOzqvFsizc7oMAOBztAChdW4pVZWFvaye2VywAiYnDYMiDkFKC52gChV8Ukvq91JCNW0REwtNFF13ECy+8QGVlJQsWLODmm2/G4XDQvn37UA9NmoAqUbV4vcVYlhsAh6MtACU5JQAkDk0MrlCekDAYbD44awkAh/6ts/REROR4LpeLzMxMunTpwk033cT48eN59913jzvuvvvuY8iQIbz00ktkZ2eTnJzMj3/8Y4qL1XcbzhSianG7/WHINGOx2eKAmhCVMDQheJzL1Rm7PQXO/hKAg+8c1Fl6IiItxLIsvKXeFt+a4n/nY2Njqaqqqve+7du38/bbb/Pee+/x3nvv8dlnnzFnzpzTfk1pPprOqyUQogJVKKgVoobUhCjDMIiPH0zhiGUYTouKHRWUbysnrmdcyw5YRCQK+cp8fJHwRYu/7piSMdjibY16rGVZLFy4kI8++ohf/vKX9R7j8/mYN28eiYmJAFxzzTUsXLiQ3//+940eszQvVaJqCfRD2e01IapsUxkA8f3j6xybkDAYYitwDPYHr6OfHm2ZQYqISMR47733SEhIICYmhokTJzJlyhTuu+++eo/Nzs4OBiiA9u3bs3///hYaqTSGKlG1eDyBSlQaAO4jbjyH/WdRxHSPqXNsQsJgAMzh6+Hr8ziy6Agdbqi7yrmIiDQ9M85kTMmYkLzuqRo3bhxPP/00TqeTDh06YLef+GvX4XDUuW0YBj6fzv4OZwpRtRw7nVe+vdx/O8OBPaHuW5WQMMT/mAEfA+dx9NOjWJYVbD4XEZHmYRhGo6fVWlp8fDw9evQI9TCkmWg6r5aa5Q38IapiewUAsd1jjzs2Lq4fYMPbczlmrIH7gJvS9aUtNlYREREJrZCGqM8//5zJkyfToUMHDMPg7bffPunxb775JhdccAHt2rUjKSmJUaNG8dFHHzXZeGoqUf7pvEAlqr4QZbPFEBfXB5xu4kb6p/wKPytssrGIiIhIeAtpiCotLWXw4ME8+eSTDTr+888/54ILLmDBggWsXLmScePGMXnyZFavXt0k4znRdF5sj+NDFNT0RdmH+C8VU/iVQpSIiPjNmzfvhMWB6667jqNHjwZv33fffeTk5NQ55tZbbyU3N7fZxienL6Q9URMnTmTixIkNPv7xxx+vc/vhhx/mnXfe4d///jdDhw497fF4PP5r5Nnt/tXHK/MqAYjpElPv8QkJg9m//2WsfmuAzhQtLTrtMYiIiEhkiOieKJ/PR3FxMW3atDnhMZWVlRQVFdXZTsTr9d9ntyf7H7vXH6KcHZ31Hh8fPxCAqu6fgQEVOyqozK9s1N8iIiIikSWiQ9QjjzxCSUkJP/rRj054zOzZs0lOTg5uWVlZJzzW4/GHKJstCYDKPf5A5Oroqvf4+PgBAFTY1hA/wL/QpqpRIiIi0SFiQ9TLL7/M/fffz2uvvUZ6evoJj5s5cyaFhYXBbdeuXSc81uPx9zTZ7Un+Zf4LvQC4OtQfolyuTthsSViWh9gR/rU8FKJERESiQ0SuEzV//nx+9rOf8frrrzN+/PiTHutyuXC56g9BxwpM59lsycGpPDPOxJZU/3ok/su/9KeoaCmOwQVAsprLRUSaga5PGh0i7XOOuErUK6+8wrRp03jllVeYNGlSkz2vZfnwev1Xy7bbk6ja679ApKuD66QLaAam9BiwHoDiFcX4qrTCrIhIUwis4l1WVhbikUhLCFyc2WaLjMVUQ1qJKikpYdu2bcHbO3bsICcnhzZt2tC5c2dmzpzJnj17+Mc//gH4p/CuvfZa/vKXvzBy5Ejy8/MB/1Wxk5OTT2ssgQAF/p6oyj3+itKJmsoDAiGqKv1r7Knn4jnioXRtKYnDE0/6OBER+W42m42UlJTgNeTi4uJ0ZYhWyufzceDAAeLi4k56eZxwEtJRrlixgnHjxgVvz5gxA4Brr72WefPmsW/fPvLy8oL3//3vf8fj8XDzzTdz8803B/cHjj8dgaZyw3Bis8VQtfcAcOJ+qIBAiCorW0fiiESOfHyE4hXFClEiIk0kMzMTQBfjjQKmadK5c+eICcohDVFjx4496fznscFo8eLFzTaWmuUNqs/MCyxv0L5hlajy8u2kDYvhyMf+KT1uaLahiohEFcMwaN++Penp6bjd7lAPR5qR0+nENCOn0ygy6mUtoGZ5A/+0oPuA/x+qM+PkIcrpTMfhSMPtPohjoH+xzuKVxSd9jIiInDqbzRYxvTISHSIn7jWz2ssbALgP+kOUI83xnY8NVKOM3v7+rtK1pXgrvM0xTBEREQkTClHVapY3qA5Rh049RFWmrMGR5sDyWJR+U9pMIxUREZFwoBBVLTCdF7jkS2MqUWVl60kc4W8oL16hKT0REZHWTCGqmtd7+tN5paXrFKJERESihEJUtdrXzfNV+fAW+XuaGhKi4uL6A1BZuYu4of63VCFKRESkdVOIqub1lgBgsyUE+6EwwZ7y3ScwOhwpOJ0d/Q/puxuA0vWleMvUXC4iItJaKURV8/n8lxSw2eJqpvLaOjDMhi34FVy5PGk9zkwn+KAkp6R5BisiIiIhpxBVzev1n01nmvGn1A8VULu5PGF4AqD1okRERFozhahqgRBls51eiCotXUfiUH9zuSpRIiIirZdCVLV6p/MaFaLWkjDUX4kqWa0QJSIi0lopRFWrM513qKYnqqHi4/sC4HYfwNXPf9290vWl+Kp8TTxSERERCQcKUdVqT+d5jnqAhp2ZF2CzxRMT0w0AT/pWbMk2rCqLso1lTT9YERERCTmFqGq1p/MaE6LgmObyIdXN5avVXC4iItIaKURVqz2d5y30r+9kSz61q4XXbi4PhCg1l4uIiLROClHVTnc6DyA+3r9yeZ0z9NRcLiIi0iopRFWrCVFxeAqrQ1Ry46bzSkvXET84HvBXoizLasKRioiISDhQiAIsy4tl+c+oM83GV6Li4noDNrzeQuw9CjGcBt4iLxU7Kpp6yCIiIhJiClGA11tzBp3NFt/oSpRpuoiL6wVAedV64gfUVKNERESkdVGIomYqDwxMM6bRlSio2xcVbC5XX5SIiEiroxBFzfIGphmHr9KHVeXvYTrVShQcc4beUJ2hJyIi0lopRFH/mXkYYEs8tSUOoHaI0lpRIiIirZlCFHXPzAuuEZVkwzCNU36u2iEqflAcGFC1p4qqA1VNN2AREREJuVOfr2qFaqbzap2Z14ipPICYmO4YhgufrwyPYzexPWIp31pOSU4JwxYubdBz5M6Z1KjXFhERkZajShTHTOcVNr6pHMA07cTF9QHUXC4iItKaKUTRNAtt1qbmchERkdZPIQrw+coB/9l5p7O8QUB9zeWqRImIiLQuClGAz+dfUdw0Y4OVKFvSqZ+ZF1BfJapscxlO9ZaLiIi0GgpR1A5RMXhL/Gfn2RNPpxLlX3CzrGwTjnQTZ6YTLMg6oLdbRESktdC3OvWHKFtC4ytRMTFdMM14LKuK8vJtwSm9zvv1douIiLQW+lan6UOUYZi1Lv+yPjil16VAb7eIiEhroW91ajeWN02Igvr7olSJEhERaT30rU7TV6Kg/hDV6YCJ6TutpxUREZEwoRDFCUJUI66bV1vNdN46YrvFYku04fQYtD906peSERERkfCjEEXzVqLKy7fhsypJGBzoizq95xUREZHwoBBF84Qop7M9dnsq4KW8fHNNX5Say0VERFoFfaPTPCHKMIx6+6K6qLlcRESkVdA3Os0ToqBuX1RwragCE6zTfmoREREJMYUomjNE1VxDL75/PB7TIr7SIK1IzeUiIiKRTiGKmhBleGOwKv1loqYNUeswnSZ70vzrG6gvSkREJPLp25yaEEVFTHBfU4SouDj/dF5FxQ48nhJ2ZvhDlFYuFxERiXz6NqcmRFllLgAMh4HpPP23xulMw+HIAKCsbAN5ClEiIiKthr7NqVWJKveHqKaoQgXU7ovKS/eHqCydoSciIhLxQvpt/vnnnzN58mQ6dOiAYRi8/fbb3/mYxYsXM2zYMFwuFz169GDevHmnPY6aSpQDaK4QtS4YotoWmySUNdlLiIiISAiENESVlpYyePBgnnzyyQYdv2PHDiZNmsS4cePIycnh1ltv5Wc/+xkfffTRaY2jpUJUhQvyUzWlJyIi0hrYQ/niEydOZOLEiQ0+fu7cuXTt2pVHH30UgL59+/Lll1/y5z//mQkTJjR6HC0VogDy0n1kHjHpvN/G+q71X404+673G/TcuXMmNc0gRURE5JRFVDlk6dKljB8/vs6+CRMmsHTp0hM+prKykqKiojpbbZblxbLc/t/L/OGpaUNUPwCqqvYS7yjWGXoiIiKtRER9k+fn55ORkVFnX0ZGBkVFRZSXl9f7mNmzZ5OcnBzcsrKy6tzv81XW/F7a9CHKbk/C5eoCQKeE3GBflNaKEhERiWyt/pt85syZFBYWBrddu3bVuT94Zh7gK/O/Hbb4pgtRAAkJgwHIStoRXOYg87CBs6pJX0ZERERaUESFqMzMTAoKCursKygoICkpidjY2Hof43K5SEpKqrPVVhOibFjVxSwzvmnflkCI6py4g8IEi6PxPkwMsg5E1NsvIiIitUTUt/ioUaNYuHBhnX0ff/wxo0aNavRz1rluXln1dfPimqkSlbgDQH1RIiIirUBIv8VLSkrIyckhJycH8C9hkJOTQ15eHuCfips6dWrw+BtvvJFvv/2W3/zmN2zatImnnnqK1157jdtuu63RY6gdonzl/nBjxjbt2xIf7w9RnRJ2YhreYF+UQpSIiEjkCum3+IoVKxg6dChDhw4FYMaMGQwdOpR7770XgH379gUDFUDXrl15//33+fjjjxk8eDCPPvoo/+///b8mWd7ANGPwlfnDTVNXomJju2GzJeCwucmM2xPsi9LK5SIiIpErpOtEjR07FsuyTnh/fauRjx07ltWrVzfZGHw+fyNU7em8pq5EGYZJfPxAioqWkpW0gx0Z/rP1sg6Y2LzgbdrMJiIiIi0g6kshgSUOTNNZM50X1/RvS+3m8gMpFuVOC4fXoP1ho8lfS0RERJqfQpTPv86AabpqGstjm740FOiLykrcgWVQa70olaFEREQiUdSHKMvyhyjDaJlKVOAMvTydoSciIhLRov4bvKYS5Qw2ljd1TxRAfPxAfJZBasxhEp1H2akz9ERERCJa1H+D11SiXHjLm2edKAC7PYH9ZZmAvy8qUInqvN+EE/fWi4iISJiK+hDVUpUogF3FXQH/lN6eNB8e0yKu0iCtUM3lIiIikSbqQ1R9PVHNUYmCmhDVOXEHXhvsbqcpPRERkUgV9d/etStRwXWimqGxHCCvqBvgvxAxoJXLRUREIljUf3vXe3ZeM03n5VVXotrH78JuuMnN9L9edn7UfwwiIiIRJ+q/vQOVqJaYzjtc0Y5Sdzx200uHhF3sCIYom5rLRUREIkzUh6hgJaoqNrivuSpRYNRpLt+d7m8uTyo3aFuk5nIREZFIEvUhKlCJorIlQhTkFVWHqKRvcdtrmsu7akpPREQkokT9N3egEkVlDACGw8C0N2OIKvY3l3dJ/BZAfVEiIiIRKuq/uQMXIDaq/CGquc7MC8gt7AFAl6TtGPhq+qL26Rp6IiIikSTqQ1SwJ6q6EtUcFx+ubV9pFm6vgzhHGelx+4KVqK75WrlcREQkkkR9iAr0RFkVLqD5K1Feyx5c6qBL0nZ2t/PhtlnEVxq0O6rmchERkUgR9SGqpieqOkQ1Y1N5wM6i7gBkJ23Da6tZdFPN5SIiIpEj6r+1a87OcwLNt0ZUbblFNX1RQN0pPREREYkIUf+tHaxEVfhDVEtWovwhyqq76KaIiIhEhKgPUcGeqEoH0DKVqN3FXXD77CQ4S2gXW0Bue/81+7LzTQw1l4uIiESEqA9RNZUof4hq7sZyAK/lYHdxNgBdkraxt61Fpd0itsog87Cay0VERCJB1IeomkqUHWiZ6Tyo1VyevA2fCXkZWnRTREQkkkT9N3agEmWV+0NUc68TFXBsc3mgL6qrFt0UERGJCFEfooKVKLc/vLRYJaowsMyBv7lcl38RERGJLFH/jV2zTlT1dJ6rZd6S3SXZeHw2Ep1FtI05wI5Mf3N5l/0mpq9FhiAiIiKnIepDVE1PlP+tMGNa5i1x+5zsLekM+JvL89tYlDstXG6DDofUXC4iIhLuoj5EBXuiqvxvheFquQAT6IvKTt6OZcLO6ubybnvVFyUiIhLuoj5EhaoSBZBbGGgu3wbA9g7+Kb1u+6L+YxEREQl7Uf9t7fNVAjWVqBYNUdWVqK7JWwGL7R38lajuqkSJiIiEvagPUcHpvEr/NF5LNZYD7CruisdnJ8lZRFpsQTBEdTxgEFPZYsMQERGRRoj6EFUznVcdolqwEuX2Ockr7gpA9+TNFCZYHEzyYWJoqQMREZEwF9Xf1JZl1apE+fe1ZCUK4NujvQDomrIFQFN6IiIiESLKQ5QX8F/x1xeCShTAt4X+ENUteStQ01zefW9UfzQiIiJhL6q/qYMLbQJWpT9MtXyI6g1AdtI2bIaHb9vXqkRZLToUEREROQVRHaIC/VAAvhBN5xWUdqDMHY/TVkXHhJ3szPDhMS2SywzSirTopoiISLiK6hBVpxJVEZpKlIUZnNLrnrIZtwPyqhfd7L4nqj8eERGRsBbV39KBSpRhOPFV+oNLS65YHhAIUV2T1VwuIiISKaI6RAUqUYbhwFfhDy4tXYmCmjP0uleHqG/ba+VyERGRcBfV39KW5QbANB3BSlRL90RBTSWqfcIuYmxlbKuuRHUpMLF7Wnw4IiIi0gBRHaJ8Pnf1b6GtRBVVpXKwPB3TsOiavJUDKRZFsRYOr0Hn/VH9EYmIiIStqP6GDlSiDF8M+DNUSEIU1EzpdUveAgZ8G7gYsdaLEhERCUtR/Q1tWf65MsMdE9wXiuk8qLXoZspmoKa5vIeay0VERMJSyEPUk08+SXZ2NjExMYwcOZLly5ef9PjHH3+c3r17ExsbS1ZWFrfddhsVFRWNeu1gJcodF9wX8hCVvAWw2NbRX4nquTvkH5GIiIjUI6Tf0K+++iozZsxg1qxZrFq1isGDBzNhwgT2799f7/Evv/wyd911F7NmzWLjxo0899xzvPrqq/z2t79t1OsHK1EefyXKsBsYttAscJlb2AOPz0ZqzGHSYvezvb0Pr2HRttikjRbdFBERCTshDVGPPfYY119/PdOmTaNfv37MnTuXuLg4nn/++XqP/+qrrxg9ejRXXXUV2dnZXHjhhVx55ZXfWb06kZpKVCwQun4ogCpfDHlF3QDombKBKifszPRP6akaJSIiEn5C9u1cVVXFypUrGT9+fM1gTJPx48ezdOnSeh9z9tlns3LlymBo+vbbb1mwYAEXX3xxo8YQPDuvyl+JCmWIAth6tB8APVI3ArCloz9E9dqtvigREZFwYw/VCx88eBCv10tGRkad/RkZGWzatKnex1x11VUcPHiQc845B8uy8Hg83HjjjSedzqusrKSysjJ4u6ioKPj7cdN5IVitvLatR/oyIfsdeqb4Q9S2Tl5Y4ThhJSr7rvcb9Ly5cyY12RhFRETEL6LmiRYvXszDDz/MU089xapVq3jzzTd5//33efDBB0/4mNmzZ5OcnBzcsrKygvcFpvPCrRLVKXEnMbYytlZXojodMIlrXO+8iIiINJOQpYa0tDRsNhsFBQV19hcUFJCZmVnvY+655x6uueYafvaznzFw4EAuu+wyHn74YWbPno3P56v3MTNnzqSwsDC47dq1K3hfTU+UCwjdmXkBhZVt2F+WgWn46J6yicIEi/xUHyYG3fdoSk9ERCSchCw1OJ1Ohg8fzsKFC4P7fD4fCxcuZNSoUfU+pqysDNOsO2SbzR8uLMuq9zEul4ukpKQ6W8Cx60SFuhIFsO2IvxrVs7ovamtgqYM9oR+biIiI1AjpN/OMGTN49tlnefHFF9m4cSM33XQTpaWlTJs2DYCpU6cyc+bM4PGTJ0/m6aefZv78+ezYsYOPP/6Ye+65h8mTJwfD1KkITue5nUB4hKitR/sC/jP0ALZ2UnO5iIhIOApZYznAlClTOHDgAPfeey/5+fkMGTKEDz/8MNhsnpeXV6fydPfdd2MYBnfffTd79uyhXbt2TJ48md///veNev2as/PCYzoP/M3lAN1TNmMaXrZ28je7d9tnYvOCV1lKREQkLDQqRC1atIhx48Y1yQCmT5/O9OnT671v8eLFdW7b7XZmzZrFrFmzmuS1A9N5VJ1eJaqhZ8k1xJ6SzpS544hzlJGVuIOdvh4Ux1oklhtk55ts71h/75eIiIi0rEalhosuuoju3bvz0EMP1WnUjjTHTeeFQSXKwsb2o30A/EsdGLC1ky4BIyIiEm4a9a28Z88epk+fzhtvvEG3bt2YMGECr732GlVVVU09vmZVs8RB+PREQU1fVM2im4Hmcs3liYiIhItGpYa0tDRuu+02cnJyWLZsGb169eIXv/gFHTp04Fe/+hVr1qxp6nE2i+B0ntsBhFGICpyhV73oZqC5vOduG9R/EqKIiIi0sNNODcOGDWPmzJlMnz6dkpISnn/+eYYPH86YMWNYv359U4yx2dRUovwhKtQrlgdsL+yNx2ejbewB0mILyM30UemwSCo36HAwPMYoIiIS7RodotxuN2+88QYXX3wxXbp04aOPPuJvf/sbBQUFbNu2jS5dunDFFVc05VibXM3ZeeE1nVfljSG3qAcAfdqsxWsjuHp53zxN6YmIiISDRqWGX/7yl7Rv354bbriBXr16sXr1apYuXcrPfvYz4uPjyc7O5pFHHjnhNfDCRc3Zef6TFMOhsTxg0+GBAPROXee/3dnfF9VHIUpERCQsNGqJgw0bNvDEE09w+eWX43K56j0mLS2NRYsWndbgmltgOs8KwxC1+fBALun2Bn3arAVgY3WI6r3LhmGBpVk9ERGRkGpUapg1axZXXHHFcQHK4/Hw+eefA/41nc4777zTH2EzqlniIPxC1NYjffH6TNrFFdA2Zn+dvqiO6osSEREJuUalhnHjxnH48OHj9hcWFjbZIpwtITid5/FPkRnO8AknFd44cot6AtC7ui8qsNSBpvRERERCr1EhyrIsDOP4wHHo0CHi4+NPe1AtJViJCoQoR/iEKIBNhwcABKf0NnXxN5crRImIiITeKfVEXX755QAYhsF1111XZzrP6/XyzTffcPbZZzftCJtR4Ow8qzpEmc7wmc4Df3P5pG7/ok+b6ubyrJpKlPqiREREQuuUQlRycjLgr0QlJiYSGxsbvM/pdHLWWWdx/fXXN+0Im1HNYpvhWYnaeqQfPsskPS6fNjEHyM1sR4XDIqHCoOMBg93pWnlTREQkVE4pRL3wwgsAZGdnc8cdd0TU1F19aqbz/BWoYytRTXlh4cao8MaRW9idbilb6Z26jqUV49jaycvAHXb65NnYne4J6fhERESiWaPPzov0AAW1ljgI00oUwOYj/vWiapY60KKbIiIi4aDBlahhw4axcOFCUlNTGTp0aL2N5QGrVq1qksE1pwGzPuJXI3czIhN25JeTSQI3v7aaFd94Qz20OjYdHsjErm/WNJdrvSgREZGw0OAQ9YMf/CDYSH7ppZc213halN30T4eZ1Y3l3jAs7myp7ovKiN9Hm5j97MxIp9zp74vK2m+Sl+EL9RBFRESiUoND1KxZs+r9PZLZDH9Vx/D6ZzU94XVyHgDlnni+LexJj5TNDGibw+cVF7I5y8uQ7Xb659oUokREREKkUbFh165d7N69O3h7+fLl3Hrrrfz9739vsoG1BLM6RJnV6clrC8+z3dYfHApAv7QcANZl+8fdf0cYls5ERESiRKNC1FVXXRW8Ll5+fj7jx49n+fLl/O53v+OBBx5o0gE2J5t5TCUqTDPJ+kNDAOjXNgcDH+u7VvdF7TZxuEM4MBERkSjWqBC1bt06zjzzTABee+01Bg4cyFdffcU///lP5s2b15Tja1Y2w98TZasOUd4wnM4D2H60NxWeGJKcRWQl7mBfG4tDiT4cXoNeu8M0+YmIiLRyjYoNbrc72GT+ySef8P3vfx+APn36sG/fvqYbXTM7ricqTPOI13Kw6bB/qYP+bXPAgPXVU3oDcsM0+YmIiLRyjfoG7t+/P3PnzuWLL77g448/5qKLLgJg7969tG3btkkH2JwC03m2QE+UGZ49UVAzpde/ui9qvfqiREREQuqUViwP+MMf/sBll13Gn/70J6699loGDx4MwLvvvhuc5osEgek80+dfbCkclzgICDSX90pdj8OsYkO2E4DOB2wklULRSdY+bejK67lzJp32OEVERKJFo0LU2LFjOXjwIEVFRaSmpgb3//znPycuLq7JBtfcbIZ/eYDA2XnhuMRBwN7SLI5UtCE15jA9Ujaw0TeE3Awv2QU2+ufaWNo/vBYJFRERae0aHRtsNludAAX+a+qlp6ef9qBaiq16sU1bBFSiwAhO6Q04dkovN6wHLiIi0io1KkQVFBRwzTXX0KFDB+x2Ozabrc4WKQLrRNm81SEqjCtRUDOl17/tav/t2iEqfNu5REREWqVGTeddd9115OXlcc8999C+ffuTXkcvnNkML3hNjOoL0HnCdLHNgA3VlajOSd+S6DzK1k4pVNktUktMOh402NMuvMcvIiLSmjQqRH355Zd88cUXDBkypImH07Jspgc8NW9BuFeiCqtSyS3sTnbydgamreKrvd9jc5aXgTvsDPzWxp52nlAPUUREJGo0KjZkZWVhWZFf9bAZ3johKlzXiartmwMjABjc7msA1nTzT+kN+rZReVhEREQaqVEh6vHHH+euu+4iNze3iYfTso4NUeHdWO635qA/RA1IW4VpePmmuz9E9dptElsZypGJiIhEl0aVL6ZMmUJZWRndu3cnLi4Oh8NR5/7Dhw83yeCam386z78kg9ewsCKgtevbo70orkoi0VlE9+RNbLX6s6+Nj/aHTfrvsLGij5Y6EBERaQmNClGPP/54Ew8jNGyGD9z+ABgJVSgACxtrDw7j7A6LGZz+NVuP9mdNNy/tD5sM+lYhSkREpKU0KkRde+21TT2OFmfgwzR8wfQUzgttHuubAyM4u8NiBqWt4I0t17Gmu4eLVjgY/K0NwyIiKmoiIiKRrtHRYfv27dx9991ceeWV7N+/H4APPviA9evXN9ngmlPg4sORVokCWHtwGD7LoHNSLqmug2zJ8lHutEguNemSH0FpUEREJII16hv3s88+Y+DAgSxbtow333yTkpISANasWcOsWbOadIDNxai++HCgsTzc14iqrdSdxPajvQEY1G4FXlvNwpuDvo2gNCgiIhLBGhWi7rrrLh566CE+/vhjnE5ncP/3vvc9/vvf/zbZ4JqTvfriw4EQFe5rRB3rmwNnADC43Qr/7eqlDgZvV4gSERFpCY2KDmvXruWyyy47bn96ejoHDx487UG1hMAlX2oqUSEcTCOsqV4vql/bHBxmVTBEdd1nklgaypGJiIhEh0aFqJSUFPbt23fc/tWrV9OxY8fTHlRLsBn+6TtfoCcqwipRecXdOFzRlhh7BX3brOFookVuhhcTQwtvioiItIBGRYcf//jH3HnnneTn52MYBj6fjyVLlnDHHXcwderUph5jszCqK1GW2z8d6Y2gnig/g1UFZwEwLMM/hZpTvfDm0G0RVlYTERGJQI0KUQ8//DB9+vQhKyuLkpIS+vXrx5gxYzj77LO5++67m3qMzcJm+ADwVYeoSFriIGDV/lEADE1fhoGXVb38IWrgtzac7lCOTEREpPVr1LyP0+nk2Wef5d5772Xt2rWUlJQwdOhQevbs2dTjazYm1SHKE6hEhXI0jbP58ABK3fEku47SPWUz26x+HEzykVZk0j/XxuqeWnhTRESkuTQ4RM2YMeOk99c+K++xxx5r/IhaiHnsOlERWInyWnbW7D+TszsuYnjGf9l2tB+renq5cKXJsC0KUSIiIs2pwSFq9erVdW6vWrUKj8dD797+9Yq2bNmCzWZj+PDhTTvCZmIYdStRkXZ2XsCq/WdxdsdFDEtfyqubp7Gyl4cLVzoYst2O6avCF4HhUEREJBI0+Ct20aJFwW3y5Mmcd9557N69m1WrVrFq1Sp27drFuHHjmDRp0ikN4MknnyQ7O5uYmBhGjhzJ8uXLT3r80aNHufnmm2nfvj0ul4tevXqxYMGCU3pNqFmx3ApWoiKtsdxv7cFhuL0OMuL30TFhJ1s7+SiOtUgsN+i1WwlKRESkuTTqW/bRRx9l9uzZpKamBvelpqby0EMP8eijjzb4eV599VVmzJjBrFmzWLVqFYMHD2bChAnBy8gcq6qqigsuuIDc3FzeeOMNNm/ezLPPPtuoZRXM6kpUIERFaiWq0hvL+kNDABiW/l98Jqzu4V9IdNgWLXUgIiLSXBoVooqKijhw4MBx+w8cOEBxcXGDn+exxx7j+uuvZ9q0afTr14+5c+cSFxfH888/X+/xzz//PIcPH+btt99m9OjRZGdnc9555zF48OBT/huCIcoTuT1RASsL/GfpDc9YCsCq6l6oYVttEJkFNhERkbDXqOhw2WWXMW3aNN588012797N7t27+de//sVPf/pTLr/88gY9R1VVFStXrmT8+PE1gzFNxo8fz9KlS+t9zLvvvsuoUaO4+eabycjIYMCAATz88MN4vSduoK6srKSoqKjOBsdP50VqJQogZ/9IfJZJdvJ20mILWJ/tpcJhkVZk0qUggtOhiIhIGGvUN+zcuXOZOHEiV111FV26dKFLly5cddVVXHTRRTz11FMNeo6DBw/i9XrJyMiosz8jI4P8/Px6H/Ptt9/yxhtv4PV6WbBgAffccw+PPvooDz300AlfZ/bs2SQnJwe3rKwsoKaxHE9k90QBFLuT2XR4AABnZH6B2wFru/pD4vAtEZwORUREwlijQlRcXBxPPfUUhw4dYvXq1axevZrDhw/z1FNPER8f39RjDPL5fKSnp/P3v/+d4cOHM2XKFH73u98xd+7cEz5m5syZFBYWBrddu3YBtSpRnsivRAEszx8DwMjMLwBY2dvfF3XGZrum9ERERJrBaXUex8fHM2jQoEY9Ni0tDZvNRkFBQZ39BQUFZGZm1vuY9u3b43A4sNlqEk/fvn3Jz8+nqqoKp9N53GNcLhcul+u4/eaxZ+dFeIhamX821/R9muzk7aTH7SWneweq7BbtD5t03m+Sl+EL9RBFRERalZA1zDidToYPH87ChQuD+3w+HwsXLmTUqFH1Pmb06NFs27YNn68mEGzZsoX27dvXG6BOpmY6z58jI/GyL7UVu5PZeNjfYH9m5pdUuOCbbv6geMamCE+IIiIiYSik0WHGjBk8++yzvPjii2zcuJGbbrqJ0tJSpk2bBsDUqVOZOXNm8PibbrqJw4cPc8stt7Blyxbef/99Hn74YW6++eZTfu3AtfOCK5ZH3AWIj7d8n39K78zMzwFY1sc/pTdyk6b0REREmlpIFxKaMmUKBw4c4N577yU/P58hQ4bw4YcfBpvN8/LyMM2anJeVlcVHH33EbbfdxqBBg+jYsSO33HILd9555ym/dnCJA6//LYjkJQ4CVu4fxVTfk3ROyqV9/C7WdM+i0mGRftQkO98kt72m9ERERJpKyFdjnD59OtOnT6/3vsWLFx+3b9SoUXWu09dYgZ4ooxUscRBQ6k5k/cGhDE5fwZmZX/BO6VXkdPcycpOdkZtsClEiIiJNqBXUXxqnZrHN6kpUKwhRAMvyA1N6XwKwvHpK7wxN6YmIiDSpqA1RgSUOAo3lkbxOVG2r95+F22enY2IenRJ28E03L+VO/8Kb3fdG7cctIiLS5KL2W9U89uy8VlKJKvfEs2b/mQCM7rgItwNW9/AHxpEbQz57KyIi0mooRHlaT2N5wJK93wNgVPtFmIaXZX0DZ+nZMNUWJSIi0iSitjRhUt1YXl2Cai2VKIBvDgynuCqJlJgj9GubwzrPcIpiLZJLTfrn2ljbrf5rDWbf9X6Dnj93zqSmHK6IiEhEakX1l1NjGNU9UK2wEuW1HCzbdy4Aozt8itcGy/r5q1Fnr4va3CwiItKkWlF0ODWBxvKaSlTraCwPWLJ3HADDMv5LjK2MJf39IWr4VhsxlaEcmYiISOsQtSHKNAMhqvVVogB2FPZiX0knXLZKRmQuITfTx942PpwegxGbVY0SERE5Xa0sOjScSXWHtbf19UT5GcEG89EdPgUDlgzwV6NGr1eIEhEROV1RG6JqpvNaZyUK4Ku9YwHo23YtabEFLO3vwYdF3zwbbQuN0A5OREQkwrXC6NAwgSUOWmtPFMDhinQ2HBoEwDkdP+FwksWmLv6/e9QGVaNERERORxSHqLqN5a3lsi/H+mz3BADGdPwYAy9fVTeYj16ny8CIiIicjqgNUYElDozqeTxPK30nVhWMorgqibaxBxnUbhUrenmocFi0P2zSa3cr/aNFRERaQNR+ix67xEFrrUS5fc7gcgfndfqQChfBFczPW6MpPRERkcaK2hAVXOIgcHZeK7kAcX0+23URAIPbfU2K6xCLB/tD1Bmb7cRVhHJkIiIikSt6Q1T1EgdmK69EAewrzWLz4X7YTB9jOn7CjvY+8tp5cXoMRmm5AxERkUaJ3hAVrET534LWHKKgpsH83E7/wTB8fFZdjRq7xqEGcxERkUaI3hAVqES18sbygK/zz6HUHU+7uAL6p61maX8PVXaLrAMm3fa18j9eRESkGUTtt6fN8IEFpjcwnde6yzFun4sle/wrmJ/f+X3KYuDr3mowFxERaayoDVGm6QNPTXhojSuWH2th3iWAv8G8XWx+cEpv5Ea7LkosIiJyiqIgOtTPoG6Ian3XzjteQVlHvjkwHNOwOL/ze2zp5GNPWx8xboNz1qkaJSIiciqiNkTZDG/UVaIAPt45GYBzO32My17OwmFuAM5f5cBo3TOaIiIiTSpKosPxTKJvOg9g3cFh5Jd2IM5RytkdFrFkgIcyp38F8/47oqAcJyIi0kSiJDoczzRrKlEe0wIjxANqIRYmn1T3Ro3v8m8qnRZfDvT3Rl2wUlN6IiIiDRW9IapWJSoa+qFq+3L3eMo9sXRM2EW/tmuCU3oDv7WRfiRK0qSIiMhpit4QZXjB7QCiZyovoMIbx5I95wMwocvbFLSxWNPNg4nB+ascIR6diIhIZIiy+FDDZvqCy5RHWyUK4D87v4/PMhmcvoJOCbl8Msw/pTfmGzuuqhAPTkREJAJEbROMia9WJSr6TkvbX9aBFflnc2b7L7mo65s8VzyD/FQfmUdMzllrZ+Fwzwkfm33X+w16jdw5k5pquCIiImEnaitRRu3G8iisRAF8sONyAM5q/xmpsQf4zwh/b9SErx2YvlCOTEREJPxFbYiq3Vje2i8+fCI7inqx4dAg7KaXC7u8w5cDPRTFWqQXmozYHKVvioiISANFbYiy1VniIMSDCaEPdvwPAOdlfYQ9tiR4pt7FyxwQfbOcIiIiDRa18cGoU4mK3rSw9uAw8oqyibWX873OC1g4zE2l3SK7wEbfvKj9z0NEROQ7Re23ZO1KVLQtcVCXEaxGXdjlHaoSK/hikL+p/OJlWu5ARETkRKI2PphGzdl50dpYHrAs/1wKStuT5Crke1kL+OgMNz7DYuAOO1n7o/Y/ERERkZOK2m9IEyvYUR6tjeUBPsvGv7+dAsDFXf9FYdtyvu7t9d/+r6pRIiIi9YnaEGUzvVG9TtSxvto7jv1lmcFq1Ptn+RvMR2600f6QLgUjIiJyrKgNUUYUXzuvPj7Lxrvb/dWoiV3fJL9DGat6+C8FM/krZ4hHJyIiEn6iNkSZaiw/ztLqalSy6yjfy1rAO6P91aizNtrIOKxqlIiISG1RGx9shlYsP5bXsvPv7T8C/NWofR3LyOnuwbQMvv+VeqNERERqi9oQZRq11olST1TQV3u/F6xGTejyNm9XV6NGbbCTfkTVKBERkYDoDVGoElUfr2XnX1uvAeDibv/iYOcjrOnmr0ZNVjVKREQkKIpDlKWeqBNYvm8MuYXdibWXM7n7a8Fq1Oj1dp2pJyIiUi0s4sOTTz5JdnY2MTExjBw5kuXLlzfocfPnz8cwDC699NJTfs3ajeWqRNVlYfL6lusAOL/z+xR338uqnv5q1P98rjP1REREAOyhHsCrr77KjBkzmDt3LiNHjuTxxx9nwoQJbN68mfT09BM+Ljc3lzvuuIMxY8Y06nVNfLXWiWrUU7Rq6w8NZd3BIQxIy+HyHv/kjXPvYMg2GyO22Om21823HXzf+RzZd73foNfKnTPpdIcrIiLS4kIeHx577DGuv/56pk2bRr9+/Zg7dy5xcXE8//zzJ3yM1+vl6quv5v7776dbt26Nel2b6Q0uVe6J4gsQn0ygGnVWh8XYu27nywH+a+pd8ZkT9JaJiEiUC2mIqqqqYuXKlYwfPz64zzRNxo8fz9KlS0/4uAceeID09HR++tOfNvq1DUOVqO+ys6gH/913LqZhcWWfZ3n7nCrcNou+eTYG7tAcqIiIRLeQxoeDBw/i9XrJyMiosz8jI4P8/Px6H/Pll1/y3HPP8eyzzzboNSorKykqKqqzAZi11omK9mvnncwbW66lyuukb9u1dOu5hE+GBapRDgxVo0REJIpFVA2muLiYa665hmeffZa0tLQGPWb27NkkJycHt6ysLABshi770hAHyzP4YMflAPy4z3N8dHYJZU6LzvttjFof8pY6ERGRkAlpiEpLS8Nms1FQUFBnf0FBAZmZmccdv337dnJzc5k8eTJ2ux273c4//vEP3n33Xex2O9u3bz/uMTNnzqSwsDC47dq1CwDTsLTYZgO9v+OHHCpPIy12P+f2/xfvjfIveXDFZw5cVSEenIiISIiENEQ5nU6GDx/OwoULg/t8Ph8LFy5k1KhRxx3fp08f1q5dS05OTnD7/ve/z7hx48jJyQlWmWpzuVwkJSXV2YI0ndcgVd4YXts8DYBJ3d7g63P2UpDiI7XE5JL/agFOERGJTiGfzpsxYwbPPvssL774Ihs3buSmm26itLSUadP8X9pTp05l5syZAMTExDBgwIA6W0pKComJiQwYMACn8xTXMApM54X8XQh/y/LPZfPhfrhslfyw//PM/56/BHXRcgftjmoBThERiT4hjw9TpkzhkUce4d5772XIkCHk5OTw4YcfBpvN8/Ly2LdvX/O8uCpRp8DgnxtvwGeZjGz/BVUjv2Z9Fy8Or8GURVqAU0REok9YdAZPnz6d6dOn13vf4sWLT/rYefPmNf6Fg43l6olqiLzi7nyycxIXZv+bawc8xVMXDKTvc8mM2GKn7043G7t89wKcIiIirUXIK1EhpXWiTtmbW6/hcEVb0uPyGTHqn3w61L/kwf9+7MLmDfHgREREWlB0x4fgiuUhHkcEqfDG8X8bbgTgouy3WDZhO0VxFh0PmVy0XE3mIiISPaI6RFmqRDXKqv2jWFlwFnbTy4/O+Avzx1UA8IOvHLQ7oiZzERGJDtEdH9QT1Wj/t+FGyj2x9EzdRNzEd9jQxYvTYzD1Y11XT0REokNYNJaHiuV2YKBKVGMcqUzjjS3Xck2/uVzRZx5//f5Iej6VzcAddkZu9LCsX8MbpLLver9Bx+XOmdTY4YqIiDS56I4PXl325XR8mncxGw4NwmWr5NJxf+K9UZUAXPWpk7iKEA9ORESkmUV3iFJP1GmxMHlu7a3V03ob8fzoDfa28ZFcanL1J1o7SkREWrfojg+ewNl5auJprEMV6byy8XoALu33D975n534DIvR6x0M3aoSn4iItF5RHqJUiWoKn++5gDX7R+CwuZk46fd8dKZ/Wu+6j5zEl4d4cCIiIs0kuuODR+tENQ2DF9b/kpKqRLKTt2P+5Dn2tPVP6/3vx65QD05ERKRZRHWIMnTtvCZztLItz627BYAJvd5g0dU5eA2LURvtjNisN1hERFqfqA5RgbPzNJ3XNFbvP4tPdl4CwPcvfoCFo4sBuO5DFynFWoRTRERal+iND14Dw+f/89VY3nRe3fwT8oq6kuQqpNMv7ic3w0tChcEN77kwdH1iERFpRaI4RNVMMakS1XTcPidPr/kNlR4X/dqvYvPNb1LhsOibZ2PSMl1bT0REWo/ojQ+emi909UQ1rX2lWby08SYALjjnaRZdvgWAy75w0H1P9P4nJyIirUv0fqN5aq5444ned6HZfLlnPIvyLsI0LMZMu4NVA4uxWQY3vevSauYiItIqRG98qA5RPsPCit53oVn9c+MNbD/amwRXCe1+eycHUrykFZn87H0XhtrQREQkwkVvfNAaUc3OYzn42+qZFFamkJW5kb23P4fbZjFsm51Llqo/SkREIlsUhygtb9ASjlSm8VTOnXh9JoPPeYXVV30F+Puj+u9QghURkcgVvRGieo0oVaKa3+YjA3llk//6eiOm3cO6s3Zj4u+PSivU+lEiIhKZ7N99SCvlDlw3T805LeGTvMlkxu9hfJf36HnvTey5fj4d98Qz/S0XD19dQVUDZvey73q/Qa+VO2fSaY5WRETku0VxJSqw0GaIxxFFXt50PWsODMcVW0LCw7dREuclu8DG9e+p0VxERCJP9IaoYCUqxOOIIj7LxtM5d7KrOJvkzlspvXcObpvFGVvsXP65Gs1FRCSyRG+EUE9USFR44/jzyns5WpFKxshPOPiLfwAw+b9OzlkbvbPLIiISeaI3RAXPztM8Uks7XJHOoyvvp8wdT/vL55F/2UIArvvQSe+86P1PUkREIkv0fmNpnaiQ2lXcjcdX3UOV10nm9N9z8Mz12H0Gv3ozho4HdMaeiIiEvygOUVonKtS2HBnAkzl34cUg7YEZHO2+l/hKgztei9HSByIiEvaiN0JU90Tp4sOhtebAmTy37hZwVZHy2I0UdzhEaonJHa/GkFga6tGJiIicWPSGqOqz83Tx4dD7au/5/GP9TZBUTOLjN1LWtojMIya3vx5DTGWoRyciIlK/6I0Q1SUor02N5eHg012TeGnDjdDuIHF/vpnKhHKyC2zc+q8YnFWhHp2IiMjxojdEuTWdF24W5l3CSxtugKzduB65FXdMFX122bj1zRic7lCPTkREpK7oDVHVjeWazgsvC/Mm838bboDeW3A8chuemCr67bTxq3/F4FCQEhGRMBK9qxuqsTxsfZI3GZ9l8r/95mL/4+14f/0oA3Y6+dWbLv76P5WBIuIJ6Rp7IiLSEqK3DhOsRKknKhx9umsSz3xzB57+G7H98Q68rkoG5tq5RT1SIiISJqI4RAUay0M8DjmhZfvO44nVd1PVfzO2OXfii6lkQK6NX78aQ1xFqEcnIiLRLnpDVGCJA4WosLbmwBk8suIByvpvxXz0NnwJJfTca2Pmy7Ekl2hBThERCZ3oDVG+6kpU9L4DEWPLkQE8vOwPHOp6APMvv8KXepisAya//adWNhcRkdCJ3ggRrESpJyoS7C7pyoP/fZTctibmE7/EytxHxlGTu1+KoUt+9P5nLCIioRO93z66dl7EOVrZloeX/4FV9iyMJ34J3baTUmry25djGLJV87IiItKyojdCVDdDqScqslR5Y3hi9W/5qOQc+OuvYMTXuNwGv3rTxfgV0btih4iItLwoDlGqREUqCxuvbLqeZ7ffgPuhe2HSe5gY/O9CF1d/4sT0hXqEIiISDaL3/7oHF9tUT1SkWrL3fHaXdOGXv/g9aR33wN9v4IKVDjoeNHn6+xUUx5388Q1dlBO0MKeIiBwveuswuuxLq7CzqAf3L32cDRdugPvvhZhy+u20cd+LajgXEZHmFRbfMk8++STZ2dnExMQwcuRIli9ffsJjn332WcaMGUNqaiqpqamMHz/+pMefkBbbbDWK3ck8suJB3uuUge/Jm6HTLtoW2bj7ny5Gr43eYquIiDSvkIeoV199lRkzZjBr1ixWrVrF4MGDmTBhAvv376/3+MWLF3PllVeyaNEili5dSlZWFhdeeCF79uw5tRdWT1Sr4rNsvLHlOh499BMKH78LRn2Fw2Ny/QIX133oxKmLF4uISBMzLMsKaVPQyJEjOeOMM/jb3/4GgM/nIysri1/+8pfcdddd3/l4r9dLamoqf/vb35g6dep3Hl9UVERycjLvDX2A+NVj+PukSr4a4Dntv0PCR5LzCNf3f5yBHwyAf0wFy2Rfmpu/fd/NnnaN+89dPVEiIqEV+P4uLCwkKSkp1MMBQlyJqqqqYuXKlYwfPz64zzRNxo8fz9KlSxv0HGVlZbjdbtq0aVPv/ZWVlRQVFdXZgOA8nlcXIG51iqpSeWz1LOafZcP7x7ugzSHaH3Rw/0suxubYQB+5iIg0gZCGqIMHD+L1esnIyKizPyMjg/z8/AY9x5133kmHDh3qBLHaZs+eTXJycnDLysry3+GubixXT1SrZGHyYe7lzKq4lrzHHoIzl2F327juoxhufdcksSzUIxQRkUgX0R1Bc+bMYf78+bz11lvExMTUe8zMmTMpLCwMbrt27fLfEeiJUohq1XaXZHP/pvt46yff4L1hLtg8DNkUyx9esDNsS0T/5y8iIiEW0lOX0tLSsNlsFBQU1NlfUFBAZmbmSR/7yCOPMGfOHD755BMGDRp0wuNcLhcul+v4O7xa4iBaeC077+y4ipzO27jpkYfI/Mu1xOV25VdvwapBRTw31kZp7Mmfo6FrSql3SkQkeoQ0QjidToYPH87ChQuD+3w+HwsXLmTUqFEnfNwf//hHHnzwQT788ENGjBjRuBf3aLHNaLOzqAe/K7iVN+/4HO+V88H0MuybJB6dZzBim/47EBGRUxPyOsyMGTN49tlnefHFF9m4cSM33XQTpaWlTJs2DYCpU6cyc+bM4PF/+MMfuOeee3j++efJzs4mPz+f/Px8SkpKTu2FPeqJikZey8G7u37IXb3PYOuDT0JWHjFFcUz/VwJ3v1dOm6JQj1BERCJFyEPUlClTeOSRR7j33nsZMmQIOTk5fPjhh8Fm87y8PPbt2xc8/umnn6aqqoof/vCHtG/fPrg98sgjp/bCbk3nRbMD5Zn8vugnzL15M+U/ehNsHnqsT+NPz9v50ZpiXX9PRES+U8jXiWppwXWi2swj/nAX7r2unLwMfWNGsxhbGVe6PmHM6wMx1w0AoLhzAc+Os/FN5nc0Sx1DPVEiIs1D60SFE4/D/0M9UVGvwhvHC2Xf567LY9nx8zchqZDEvAxmvJjGgwsLaF9SGeohiohIGIreSlTcq8SXpfObn5exPzWq3gL5DkMcG/npsjISPzkHLBPLUcW6721lbr+OlDodTfIaqliJiJwaVaLCiRbblBPIcfflV8OG8MbtX1I1aD2G28nAj/rzxDyLn+3cgsPnDfUQRUQkDERviNIFiOUkLGy8ZxvOTRM7sfD6r/C234d5JJVz5g/l6f8r47r89TjQVY1FRKJZ9EYIy1+CUk+UnIwXOy+1GczN17hY+cMV+JKPYt+XydgXz+TpVwuZdmQlTrMi1MMUEZEQiN4QVU2VKGmICsPFE937csv1XnK+n4MvoQR7bhfO+/u5PP2vw9xQtJR42ymuVSYiIhEtpJd9CQfqiZJTUWyL4fG+PUnq5uGnGzcy8LOu2Lb0ZNSWnozsspONF67knx26s7eiU6iHKiIizSzqQ5QuQCyNUeSy8+chnUnu7ebajbsZsqQD5s4u9H+2Cw9n7mP3pI94vVsGa4oGA0aohysiIs0gqkOU17Sw9P0mp6Ew1uSvw9oT29/HlE07Gf1VWxz57en03OXcllRI8YRP+WiIxadlZ1HmSQj1cEVEpAlF7zpRvIfDEcvPZ5SHekjSijjdMHFjGRctdxJ7KNm/0+bBe+6XrDtvB/92DWDb0b40tDql9aRERPzCcZ2o6K5E6cw8aWJVDnhnUBz/HgAjtxdy6ZpKMranY1s0lsGLxjK4z0aOTHiZj7om88XB8yh1J4Z6yCIi0kgKUSLNwGfC0p52lva007mgjMu/KWPQNymYm/qSuqkvP44v4YrzP2XTWTv5NKYPa/afgcdqmtXQRUSkZUR1iNIaUdIS8jIsHr8glsSzKxm/roLvrYXEQwnY3v0+/d+F/r02U3nR/7Gsv5svDo9m69F+qBldRCT8RXWI0hpR0pKK4+GtkQZvnwl9d5YzcX0V/TcmYG7pjWtLb851VnLu2V9xdMzTfNUpgeUHR2FZFoahQCUiEo6iurG8pK2TO39WFeohSRRLLIPR62xcsM5H2wNxNXckFcLYxTguXkv694aQnv5DkpLOxDCU/EUkOoVjY3lUh6jCdAe/nabrn0kYsKBLgck5G+DsTQbxxTE192Xug/M+o3DkNyxLTmPVwVFsPdIPr1VTSNZZfCLS2oVjiIrq6TytVi5hw4CdmT52ZsLLY6HfznJGbzQYscWOM789vPpjkl/9MRemHeDCMV9QcfY7rGnvYs2hEXxzYHioRy8iEpWiuhJV0MnBA1erEiXhy+mGQdttnLHVZMh2E1dlrf/fk3IERi/Bd9Z/SR5rJy3rQtq0uZiEhEGa9hORVkeVqDDj0feMhLkqB6zo42VFHy8OD/TL9TBis40R2w1ij6bC+5dgvn8JxY4qigevYcfIh7Cfs5U2g/qTmjqe1NTxxMR0DvWfISLSKkV1JSov28EfpqgSJZHH5oXeu0yGbbUz9Ftoe/SYNaY67YKRy+CMr4k5s5w2HceQmjqelJSxOBypoRm0iMhpCMdKVFSHqO3dHTz2Q4UoiXAWrJ82lsPvH+bQgoMc/bwQPLXut3mg3wYYtgqGrSZuhJ3UdmeTnDyG5OQxuFztQzVyEZEGU4gKA7VD1OaeDv56uUKUtC4xlTAg18ag7Tb655m0LTzmDIqYchj0jT9UDV6Da4BFatrZJCefS3LyGGJju2ttKhEJO+EYoqK7J8qmLwppfSpcsKK3lxW9vWBBu6MG/Xba/FueQWJZLCwf6d+Ayphy8vtuJH/AFzDwKYp67CK7+3CSkkaSlDSSxMQzNQUoIlIPhSiR1syAA6kWn6V6+GyIB8OCjgf8oapvno2ee0wSymNh9TD/BiSZXg53+5bDA9dC/z9Cn03EdI8nOXkkiYn+YJWQMAjTdIX4jxMRCa3oDlFmVM1kimAZsDvdYne6h/+c4Q9V7Q8Z9Nxto9duk167TdoV2mBbT//21uUAVCQWUdF7MwV9VkOfV6DvNuKzMklIGEJCwlASE4eSkDAEuz05xH+hiEjLieoQ5dVimxLlLAP2plnsTfPw2RD/vpRig/+MHUHhF4UULS+iZHUJVnESrDjDv1UrTTtAad+NFPRaAd1fhe7bcWUlVAeqodUBayAuV2f1WIlIqxTVIcqtECVynKOJFmeu/BrigLFgGwOdDph022fStXrreMjEPNgOvmgHX5wbfGxlUiGV3bdzsPt26P4l9NiG2fUQ8Sm9iI/vT3x8f+Li+hMfPwCXq6PClYhEtKgOUVUOTeeJfBevLXBJGh+Lhvr3uar81/p7btgAilcXU7qmlNKNpVCUXKe/CsBn81DcZSfF2bnQZQ10eRe67MTsXExCSu/qYNWPuLhexMb2IiYmG9N01D8YEZEwEtUhSpUokcapdMKWLB9jDnwDnYBO4JgAHQ6aZO036Vy9Ze03ia+0w7fd/VstPpuHok67KeqcB9lrofP7kJ0LnfcRm9wpGKpq/3Q6O6h6JSJhI7pDVFT/9SJNy22vqVgFWdCm2OA/E8+kbEMZpRtLKdtQRtnGUrwldtiZ7d9qTQkClLfbT3nHPdBxD3RcDJ3+Dzruweh0lLjUzsTGdiUmphsxMV2rf+9KTEw2NltcS/7JIhLlojpGuO2azhNpVgYcTrJIm5wGk2t2W5ZF5e7KOsGqdIP/p+eIBw6k+7ecoXWezqK6ob3T7uqAtQQy/wWZ+ZCZjyM9JhiqYmO7VYerrsTEdMbl6qRlGUSkSUV5iAr1CESiQ/Zd75/8gDbAObDji4txH3JTvrWc8m3ldX6WbS3DW+iFg+382zEBC8DtqsCdUUBRdagi45NgwCIzH3u6g5iYLGJisnC5arbAbaezA6ap/2EQkYaJ6v+1qIrqv14k/BiGgTPNiTPNSfKoumtOWZblD1jb6gasip0VVORWULW3CipjIK+Lf6uHx1FFSbsDlLQ7AGkHod02SPsvBG8fwpHhIjahIy5XJ5zO9sHN5ar53eFoi2GYLfGWiEgYi+oY4bZpOk8knHxnxapa7pxJx+3zVfmo3FVJRa4/VAXCVWCr3FMJbifs7ejfTsBtenG3OVwdqg5A2j5o9w2kHoE2h/0/2xbhbBeDKy6j3pDldGbgcKTjdLbDZotv9PshIuEtukNUVP/1IpGrwWHrHzVhy+f2Ubm7smbbU0nVnqrg75W7K6naV4nlsdVMGW7qe8LnrjJ8VCUX+kNVMGDlQJtP6wQuo20FjjQXrri2OBztcDja4XSm1/lZ+3c1x4tEjqiOEQpRItHDdJjEdo0ltmvsCY+xfBZV+/3Bqk7A2luJu8BNVX4VVfmVVO13g8+Eo6n+bceJX9cCqoCq+BJILqy1HYLkb4/ZV4iR4saRZsfRxoUzpg12exscjto/2x5z2//TNJ1N/p6JyMlFdYzQ2XkirdupTg8apoEr04Ur0wUjTny85fX3Z1UVVPmDVeBnflVN2CrwBy73QU/1aYUJ/u0kU4lQK3QZPkoTi/3hKrEYEkog4Sgk7qp1u8T/e2IxZrIXe4oDRxsn9qQ4nK422O2p2O3J2O0p2GzJ1b/7b9vtydX7UrDbEzEMLZwncqqiO0TpfzNEhIaHre9kBzpC7ks104iW18Jz1IP7oPukW9XBKtzVm/eoBZbpXwG+qGEXdfZRHb4ATG9NyAoErYSDkJALcWUQX+r/Wet3M9HClmTDnuzAnujEnhSLw5V0TOBKwmZLrN4SsNv9P2vv02rzEk2iO0SpEiUizcywGTjaOnC0dUDvhj3G5/bhOVwreB1x4zni8W9H/T9r73MfqcRzxI3nqA+rEvDZTimAgT+E+QB37Z2x9QWug3X3xVT4t9hy/884D7Z4A1uCiS3eji3egS3RgT3BhS0mDoejJnDVDl82Wzw2WxymGYfNFl/9M/B7rM6GlLAU1SGqQi0EIhKGTIeJM8OJM+PU/0fKW+49YdjyFHnwFnn9P4u91b+78RRW4Sly4y3y4i2ysNzVl9Ypj/NvpzqG6u04Nk9N4AqErthCiCmoG8QCP12Vwc2IsTBiwBZrYsaZmLEmtjg7Zqwde5wDW7wTW6wLW0xMMJDVDmOmGY/NFotpxlRvscf8DPwe1V+Lcoqi+r+WCqcqUSLS9JpsevAU5c6ZhC3Whi3WhqtD41dn91X6jgtcdW7X+uktrd6KPXhKq/AWe/CWevCWevGVWvjKwKqsDmVee01v2Cmyqjffdx1oesFZVRPCnFXgKgPXkZpQduz9tTeHF8MFpsvCcBmYThMjBkyXDdNlYsbYMGPsmC4btlgHthj/ZsY4scU4sTmODWguTNOJYfh/mqbrmN+d1cfU/G4YDl0jMkJEeYgK9QhERJrO6ayzVZvpMnG2c0K7phgV+Dw+fKU+vCXVgaukVvgq8dZ/X4kXX7kXb5nbv5V78JZ58JX7arYK/CGtwgCrOnT4bFAR698ayeIElbSGsHnA4a4VytzgLPX/DOyv/dNefbzNW3Pb7gGHD8MJht3y/3SA4TAwnAamEwy7ieEyMZ0GhsOG6TT9Qc9hw3TaMF12TIfd/9PpD32mM/C7A9N0VIe1ur/X3Hae5D5HdeCrfV90hr6wCFFPPvkkf/rTn8jPz2fw4ME88cQTnHnmmSc8/vXXX+eee+4hNzeXnj178oc//IGLL774lF7TZ/PgsYfFny8i0qKaulL2naHMbmImm9iTm+d/cy3Lwqqy8JZ78ZX5A5a33OsPWmW1fq/evGXeur9X+PBWuPFVuPFWePBVevBVePFVev0/Kyx8VT6sSgtfBVhV1VuFiVVVK8CBv9rmtZ9WiAv+XdVbs7AdG94qwV7qv23z+oNc4Pdjt2Pvs3vA5gObheGw/D/tYNgBe83v/s0AO5h2wx8KbQYEfrcbmA7DHxADt+0mhsPEsJuUVlU117vRaCFPEa+++iozZsxg7ty5jBw5kscff5wJEyawefNm0tPTjzv+q6++4sorr2T27NlccsklvPzyy1x66aWsWrWKAQMGNPh1PTFuwuDPFxGJeKGavmyo7wp5p8OyLCyPha/Ch6+yOmhV+oK3A7+fcL/bwnJb+Cq9eCs9+Krc+Kqqg5zbi1VVHebcXqwqnz/MuatDXVX1Y6uofh6w3IY/2HnAcleHvPpORQ+EvaZ+P4752ZRKKW2GZz09hmVZIW0MGjlyJGeccQZ/+9vfAPD5fGRlZfHLX/6Su+6667jjp0yZQmlpKe+9915w31lnncWQIUOYO3fud75eUVERycnJvNbu//j1T1Ka7O8QEZHI1pxhK5Qsy8LyWrVClz+A+dzVPwOBzFO9uWt+97l9x+93W/6w5/Hic3uw3NWBz+MPe74qL5bH43+s2+ff7wn8bmF5/D99VdW3vdUh0IM//B2z4TGwPFBaVc7E/VdQWFhIUlJSqN9WIMSlmKqqKlauXMnMmTOD+0zTZPz48SxdurTexyxdupQZM2bU2TdhwgTefvvtU3rt8qRyIOUURywiIq1VuFfUIoYJOKu3ejQ0rB77efgqLXj8tEbW5EIaog4ePIjX6yUjI6PO/oyMDDZt2lTvY/Lz8+s9Pj8/v97jKysrqaysDN4uLCz0P0/KEXyVDV9DRURERE5f59teb9TjfJVlgL+yFi5afVPQ7Nmzuf/++4/bf+fK38DKEAxIREREGu3QoUMkJ4dHESSkISotLQ2bzUZBQUGd/QUFBWRmZtb7mMzMzFM6fubMmXWm/44ePUqXLl3Iy8sLmw8hWhUVFZGVlcWuXbvCZn47munzCB/6LMKHPovwUVhYSOfOnWnTpk2ohxIU0hDldDoZPnw4Cxcu5NJLLwX8jeULFy5k+vTp9T5m1KhRLFy4kFtvvTW47+OPP2bUqFH1Hu9yuXC5jl90Ljk5Wf8gwkRSUpI+izCizyN86LMIH/oswodphs8lgEI+nTdjxgyuvfZaRowYwZlnnsnjjz9OaWkp06ZNA2Dq1Kl07NiR2bNnA3DLLbdw3nnn8eijjzJp0iTmz5/PihUr+Pvf/x7KP0NERESiTMhD1JQpUzhw4AD33nsv+fn5DBkyhA8//DDYPJ6Xl1cndZ599tm8/PLL3H333fz2t7+lZ8+evP3226e0RpSIiIjI6Qp5iAKYPn36CafvFi9efNy+K664giuuuKJRr+VyuZg1a1a9U3zSsvRZhBd9HuFDn0X40GcRPsLxswj5YpsiIiIikSh8urNEREREIohClIiIiEgjKESJiIiINELUhagnn3yS7OxsYmJiGDlyJMuXLw/1kMLa559/zuTJk+nQoQOGYRx3jULLsrj33ntp3749sbGxjB8/nq1bt9Y55vDhw1x99dUkJSWRkpLCT3/6U0pKSuoc88033zBmzBhiYmLIysrij3/843Fjef311+nTpw8xMTEMHDiQBQsWnPJYItXs2bM544wzSExMJD09nUsvvZTNmzfXOaaiooKbb76Ztm3bkpCQwP/8z/8ctzBtXl4ekyZNIi4ujvT0dH7961/j8XjqHLN48WKGDRuGy+WiR48ezJs377jxfNe/o4aMJZI9/fTTDBo0KLh20KhRo/jggw+C9+uzCJ05c+ZgGEadtQT1ebSM++67D8Mw6mx9+vQJ3t8qPwcrisyfP99yOp3W888/b61fv966/vrrrZSUFKugoCDUQwtbCxYssH73u99Zb775pgVYb731Vp3758yZYyUnJ1tvv/22tWbNGuv73/++1bVrV6u8vDx4zEUXXWQNHjzY+u9//2t98cUXVo8ePawrr7wyeH9hYaGVkZFhXX311da6deusV155xYqNjbWeeeaZ4DFLliyxbDab9cc//tHasGGDdffdd1sOh8Nau3btKY0lUk2YMMF64YUXrHXr1lk5OTnWxRdfbHXu3NkqKSkJHnPjjTdaWVlZ1sKFC60VK1ZYZ511lnX22WcH7/d4PNaAAQOs8ePHW6tXr7YWLFhgpaWlWTNnzgwe8+2331pxcXHWjBkzrA0bNlhPPPGEZbPZrA8//DB4TEP+HX3XWCLdu+++a73//vvWli1brM2bN1u//e1vLYfDYa1bt86yLH0WobJ8+XIrOzvbGjRokHXLLbcE9+vzaBmzZs2y+vfvb+3bty+4HThwIHh/a/wcoipEnXnmmdbNN98cvO31eq0OHTpYs2fPDuGoIsexIcrn81mZmZnWn/70p+C+o0ePWi6Xy3rllVcsy7KsDRs2WID19ddfB4/54IMPLMMwrD179liWZVlPPfWUlZqaalVWVgaPufPOO63evXsHb//oRz+yJk2aVGc8I0eOtG644YYGj6U12b9/vwVYn332mWVZ/r/V4XBYr7/+evCYjRs3WoC1dOlSy7L8gdg0TSs/Pz94zNNPP20lJSUF3/vf/OY3Vv/+/eu81pQpU6wJEyYEb3/Xv6OGjKU1Sk1Ntf7f//t/+ixCpLi42OrZs6f18ccfW+edd14wROnzaDmzZs2yBg8eXO99rfVziJrpvKqqKlauXMn48eOD+0zTZPz48SxdujSEI4tcO3bsID8/v857mpyczMiRI4Pv6dKlS0lJSWHEiBHBY8aPH49pmixbtix4zLnnnovT6QweM2HCBDZv3syRI0eCx9R+ncAxgddpyFhak8LCQoDgNaRWrlyJ2+2u8/f36dOHzp071/ksBg4cGFzIFvzvYVFREevXrw8ec7L3uSH/jhoyltbE6/Uyf/58SktLGTVqlD6LELn55puZNGnSce+ZPo+WtXXrVjp06EC3bt24+uqrycvLA1rv5xA1IergwYN4vd46Hw5ARkYG+fn5IRpVZAu8byd7T/Pz80lPT69zv91up02bNnWOqe85ar/GiY6pff93jaW18Pl83HrrrYwePTq4Un9+fj5Op5OUlJQ6xx77HjX2fS4qKqK8vLxB/44aMpbWYO3atSQkJOByubjxxht566236Nevnz6LEJg/fz6rVq0KXh6sNn0eLWfkyJHMmzePDz/8kKeffpodO3YwZswYiouLW+3nEBYrlotIw918882sW7eOL7/8MtRDiWq9e/cmJyeHwsJC3njjDa699lo+++yzUA8r6uzatYtbbrmFjz/+mJiYmFAPJ6pNnDgx+PugQYMYOXIkXbp04bXXXiM2NjaEI2s+UVOJSktLw2azHdd9X1BQQGZmZohGFdkC79vJ3tPMzEz2799f536Px8Phw4frHFPfc9R+jRMdU/v+7xpLazB9+nTee+89Fi1aRKdOnYL7MzMzqaqq4ujRo3WOP/Y9auz7nJSURGxsbIP+HTVkLK2B0+mkR48eDB8+nNmzZzN48GD+8pe/6LNoYStXrmT//v0MGzYMu92O3W7ns88+469//St2u52MjAx9HiGSkpJCr1692LZtW6v9dxE1IcrpdDJ8+HAWLlwY3Ofz+Vi4cCGjRo0K4cgiV9euXcnMzKzznhYVFbFs2bLgezpq1CiOHj3KypUrg8d8+umn+Hw+Ro4cGTzm888/x+12B4/5+OOP6d27N6mpqcFjar9O4JjA6zRkLJHMsiymT5/OW2+9xaeffkrXrl3r3D98+HAcDkedv3/z5s3k5eXV+SzWrl1bJ9R+/PHHJCUl0a9fv+AxJ3ufG/LvqCFjaY18Ph+VlZX6LFrY+eefz9q1a8nJyQluI0aM4Oqrrw7+rs8jNEpKSti+fTvt27dvvf8uTqkNPcLNnz/fcrlc1rx586wNGzZYP//5z62UlJQ6ZwJIXcXFxdbq1aut1atXW4D12GOPWatXr7Z27txpWZZ/WYGUlBTrnXfesb755hvrBz/4Qb1LHAwdOtRatmyZ9eWXX1o9e/ass8TB0aNHrYyMDOuaa66x1q1bZ82fP9+Ki4s7bokDu91uPfLII9bGjRutWbNm1bvEwXeNJVLddNNNVnJysrV48eI6pw+XlZUFj7nxxhutzp07W59++qm1YsUKa9SoUdaoUaOC9wdOH77wwgutnJwc68MPP7TatWtX7+nDv/71r62NGzdaTz75ZL2nD3/Xv6PvGkuku+uuu6zPPvvM2rFjh/XNN99Yd911l2UYhvWf//zHsix9FqFW++w8y9Ln0VJuv/12a/HixdaOHTusJUuWWOPHj7fS0tKs/fv3W5bVOj+HqApRlmVZTzzxhNW5c2fL6XRaZ555pvXf//431EMKa4sWLbKA47Zrr73Wsiz/0gL33HOPlZGRYblcLuv888+3Nm/eXOc5Dh06ZF155ZVWQkKClZSUZE2bNs0qLi6uc8yaNWusc845x3K5XFbHjh2tOXPmHDeW1157zerVq5fldDqt/v37W++//36d+xsylkhV32cAWC+88ELwmPLycusXv/iFlZqaasXFxVmXXXaZtW/fvjrPk5uba02cONGKjY210tLSrNtvv91yu911jlm0aJE1ZMgQy+l0Wt26davzGgHf9e+oIWOJZD/5yU+sLl26WE6n02rXrp11/vnnBwOUZemzCLVjQ5Q+j5YxZcoUq3379pbT6bQ6duxoTZkyxdq2bVvw/tb4ORiWZVmnVrsSERERkajpiRIRERFpSgpRIiIiIo2gECUiIiLSCApRIiIiIo2gECUiIiLSCApRIiIiIo2gECUiIiLSCApRIiIiIo2gECUiLWrs2LHceuutoR6GiMhp04rlItKiDh8+jMPhIDExMdRDERE5LQpRIiIiIo2g6TwRaVG1p/Oys7N5+OGH+clPfkJiYiKdO3fm73//e53jd+/ezZVXXkmbNm2Ij49nxIgRLFu2LHj/008/Tffu3XE6nfTu3ZuXXnqpzuMNw+CZZ57hkksuIS4ujr59+7J06VK2bdvG2LFjiY+P5+yzz2b79u11HvfOO+8wbNgwYmJi6NatG/fffz8ej6d53hQRiUgKUSISUo8++igjRoxg9erV/OIXv+Cmm25i8+bNAJSUlHDeeeexZ88e3n33XdasWcNvfvMbfD4fAG+99Ra33HILt99+O+vWreOGG25g2rRpLFq0qM5rPPjgg0ydOpWcnBz69OnDVVddxQ033MDMmTNZsWIFlmUxffr04PFffPEFU6dO5ZZbbmHDhg0888wzzJs3j9///vct98aISPizRERa0HnnnWfdcsstlmVZVpcuXaz//d//Dd7n8/ms9PR06+mnn7Ysy7KeeeYZKzEx0Tp06FC9z3X22Wdb119/fZ19V1xxhXXxxRcHbwPW3XffHby9dOlSC7Cee+654L5XXnnFiomJCd4+//zzrYcffrjO87700ktW+/btT/GvFZHWTJUoEQmpQYMGBX83DIPMzEz2798PQE5ODkOHDqVNmzb1Pnbjxo2MHj26zr7Ro0ezcePGE75GRkYGAAMHDqyzr6KigqKiIgDWrFnDAw88QEJCQnC7/vrr2bdvH2VlZafx14pIa2IP9QBEJLo5HI46tw3DCE7XxcbGNvlrGIZxwn2B1y0pKeH+++/n8ssvP+65YmJimmRMIhL5VIkSkbA1aNAgcnJyOHz4cL339+3blyVLltTZt2TJEvr163darzts2DA2b95Mjx49jttMU/+zKSJ+qkSJSNi68sorefjhh7n00kuZPXs27du3Z/Xq1XTo0IFRo0bx61//mh/96EcMHTqU8ePH8+9//5s333yTTz755LRe99577+WSSy6hc+fO/PCHP8Q0TdasWcO6det46KGHmuivE5FIp/9LJSJhy+l08p///If09HQuvvhiBg4cyJw5c7DZbABceuml/OUvf+GRRx6hf//+PPPMM7zwwguMHTv2tF53woQJvPfee/znP//hjDPO4KyzzuLPf/4zXbp0aYK/SkRaCy22KSIiItIIqkSJiIiINIJClIiIiEgjKESJiIiINIJClIiIiEgjKESJiIiINIJClIiIiEgjKESJiIiINIJClIiIiEgjKESJiIiINIJClIiIiEgjKESJiIiINIJClIiIiEgj/H/it0DjbN+TpQAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "# truncate data.expanded_income to 99.5th percentile \n",
+    "top99 = data[\"expanded_income\"].quantile(0.99)\n",
+    "data = data[data[\"expanded_income\"] < top99]\n",
+    "\n",
+    "\n",
+    "# plot weighted histogram of income from data.expanded_income, with weights s006\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.hist(data[\"expanded_income\"], bins=50, density=True, weights=data[\"s006\"])\n",
+    "ax.set_xlabel(\"income\")\n",
+    "ax.set_ylabel(\"density\")\n",
+    "ax.set_xlim(left=0, right=500_000)\n",
+    "ax.plot(df_ln[\"z\"], df_ln[\"f\"], 'y-', label='log normal')\n",
+    "ax.plot(df_Pln[\"z\"], df_Pln[\"f\"], 'm-', label='Pln')\n",
+    "\n",
+    "ax.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "7a0b71a4",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.legend.Legend at 0x2aa6ce090>"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAGzCAYAAADjbSfcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8SUlEQVR4nO3deXxU1f3/8fdkmySQTBISQgIJq4CssigCUkBRRKVov60LomDrjopS28qXr4pSClbkZ2sV1xK0tqitIkVEFEUFUSESRQjIKmEngcxkgZBk7u+PazaSQBJyZ3Izr+fjcR6ZuXPn3k+8tnl7zrnnOgzDMAQAANDEBfm7AAAAgLogtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFsgtAAAAFuwPLTs27dPEyZMUKtWrRQREaHevXtr/fr1Vp8WAAA0MyFWHvzYsWMaOnSoRo4cqffff18JCQnatm2bYmNjrTwtAABohhxWPjDxoYce0po1a/T555836Pter1f79+9XVFSUHA5HI1cHAACsYBiG8vLylJycrKCgxhvUsTS09OjRQ6NHj9bevXv16aefqm3btrr77rt122231bh/UVGRioqKyt/v27dPPXr0sKo8AABgoaysLLVr167RjmdpaAkPD5ckTZ06Vb/61a+0bt06TZkyRc8//7wmTpxYbf8ZM2boscceq7Y9KytL0dHRVpUJAAAakcfjUUpKinJzc+VyuRrtuJaGlrCwMA0cOFBffPFF+bb77rtP69at09q1a6vtf2pPS9kv7Xa7CS0AANiEx+ORy+Vq9L/flt49lJSUVG1459xzz9WePXtq3N/pdCo6OrpKAwAAkCwOLUOHDtXWrVurbPvhhx/Uvn17K08LAACaIUtveX7ggQc0ZMgQ/elPf9K1116rr7/+Wi+++KJefPFFK08LAPCz0tJSFRcX+7sMWCgsLKxR7wyqC0vntEjS0qVLNW3aNG3btk0dO3bU1KlTa7176FRWjYkBAKxhGIYOHjyo3Nxcf5cCiwUFBaljx44KCwur9plVf78tDy1ng9ACAPZy4MAB5ebmqnXr1oqMjGSNrWaqbB210NBQpaamVrvOVv39tnR4CAAQOEpLS8sDS6tWrfxdDiyWkJCg/fv3q6SkRKGhoT45Jw9MBAA0irI5LJGRkX6uBL5QNixUWlrqs3MSWgAAjYohocDgj+tMaAEAALZAaAEABLQRI0bo/vvv93cZTc6qVavkcDia1J1ghBYAAGALhBYAAGzq5MmT/i7BpwgtAABUcuzYMd18882KjY1VZGSkxowZo23btlXZ56WXXlJKSooiIyN1zTXXaN68eYqJian1mLt375bD4dDbb7+tkSNHKjIyUn379q328OD//Oc/6tmzp5xOpzp06KCnnnqqyucdOnTQzJkzdfPNNys6Olq333670tLSFBMTo6VLl6pbt26KjIzUL3/5SxUWFmrhwoXq0KGDYmNjdd9991W50+e1117TwIEDFRUVpTZt2mj8+PE6fPjw2f8DtBChBQBgDcOQCgr8085i3dRJkyZp/fr1WrJkidauXSvDMHTFFVeU39K9Zs0a3XnnnZoyZYoyMjJ06aWXatasWXU69vTp0/Xggw8qIyNDXbt21Q033KCSkhJJUnp6uq699lpdf/312rhxo2bMmKGHH35YaWlpVY4xd+5c9e3bVxs2bNDDDz8sSSosLNRf//pXLVq0SMuXL9eqVat0zTXXaNmyZVq2bJlee+01vfDCC/r3v/9dfpzi4mLNnDlT3377rRYvXqzdu3dr0qRJDf7n5hNGE+Z2uw1Jhtvt9ncpAIAzOH78uLF582bj+PHj5ob8fMMw44PvW35+nesePny4MWXKFMMwDOOHH34wJBlr1qwp/zw7O9uIiIgw3nzzTcMwDOO6664zrrzyyirHuPHGGw2Xy1XrOXbt2mVIMl5++eXybZs2bTIkGZmZmYZhGMb48eONSy+9tMr3fve73xk9evQof9++fXvj6quvrrLPggULDEnG9u3by7fdcccdRmRkpJGXl1e+bfTo0cYdd9xRa43r1q0zJJV/55NPPjEkGceOHatx/2rXuxKr/n7T0wIAwE8yMzMVEhKiQYMGlW9r1aqVunXrpszMTEnS1q1bdcEFF1T53qnva9OnT5/y10lJSZJUPiSTmZmpoUOHVtl/6NCh2rZtW5VhnYEDB1Y7bmRkpDp37lz+PjExUR06dFDLli2rbKs8/JOenq6xY8cqNTVVUVFRGj58uCRpz549dfpd/IFl/AEA1oiMlPLz/XfuJqjycvdli7N5vd56HaNFixanPW7ZsWvaVnaugoICjR49WqNHj9brr7+uhIQE7dmzR6NHj27Sk3sJLQAAazgcUg1/YJuyc889VyUlJfrqq680ZMgQSVJOTo62bt2qHj16SJK6deumdevWVfneqe8beu41a9ZU2bZmzRp17dpVwcHBZ338yrZs2aKcnBzNmTNHKSkpkqT169c36jmswPAQAAA/OeecczRu3DjddtttWr16tb799ltNmDBBbdu21bhx4yRJ9957r5YtW6Z58+Zp27ZteuGFF/T++++f9bL2v/3tb7Vy5UrNnDlTP/zwgxYuXKi//e1vevDBBxvjV6siNTVVYWFheuaZZ7Rz504tWbJEM2fObPTzNDZCCwAAlSxYsEADBgzQVVddpcGDB8swDC1btqx8uGXo0KF6/vnnNW/ePPXt21fLly/XAw88oPDw8LM6b//+/fXmm29q0aJF6tWrlx555BE9/vjjltzRk5CQoLS0NL311lvq0aOH5syZo7lz5zb6eRqbwzDO4r4wi3k8HrlcLrndbkVHR/u7HADAaZw4cUK7du1Sx44dz/oPuN3cdttt2rJliz7//HN/l+Izp7veVv39Zk4LAAD1NHfuXF166aVq0aKF3n//fS1cuFDPPfecv8tq9ggtAADU09dff60///nPysvLU6dOnfTXv/5Vt956q7/LavYILQAA1NObb77p7xICEhNxAQCALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAAA4jbS0NMXExPi7DIjQAgCAJk2aJIfDIYfDobCwMHXp0kWPP/64SkpK/F0aKmFxOQAAJF1++eVasGCBioqKtGzZMk2ePFmhoaFKSkryd2n4CT0tAABIcjqdatOmjdq3b6+77rpLo0aN0pIlS6rtN2PGDJ133nl67bXX1KFDB7lcLl1//fXKy8vzQ9WBhZ4WAIAlDEMqLPTPuSMjJYfj7I4RERGhnJycGj/bsWOHFi9erKVLl+rYsWO69tprNWfOHM2aNevsTorTIrQAACxRWCi1bOmfc+fnSy1aNOy7hmFo5cqV+uCDD3TvvffWuI/X61VaWpqioqIkSTfddJNWrlxJaLEYw0MAAEhaunSpWrZsqfDwcI0ZM0bXXXedZsyYUeO+HTp0KA8skpSUlKTDhw/7qNLARU8LAMASkZFmj4e/zl1fI0eO1Pz58xUWFqbk5GSFhNT+JzI0NLTKe4fDIa/XW/+Tol4sDS0zZszQY489VmVbt27dtGXLFitPCwBoAhyOhg/R+EOLFi3UpUsXf5eB07C8p6Vnz5766KOPKk54muQKAABQG8sTREhIiNq0aWP1aQAAQDPnMAzDsOrgM2bM0JNPPimXy6Xw8HANHjxYs2fPVmpqao37FxUVqaioqPy9x+NRSkqK3G63oqOjrSoTANAITpw4oV27dqljx44KDw/3dzmw2Omut8fjkcvlavS/35bePTRo0CClpaVp+fLlmj9/vnbt2qVhw4bVugDP7Nmz5XK5yltKSoqV5QEAABuxtKflVLm5uWrfvr3mzZun3/zmN9U+p6cFAOyLnpbA4o+eFp/Oio2JiVHXrl21ffv2Gj93Op1yOp2+LAkAANiETxeXy8/P144dO3j4FAAAqDdLQ8uDDz6oTz/9VLt379YXX3yha665RsHBwbrhhhusPC0AAGiGLB0e2rt3r2644Qbl5OQoISFBF110kb788kslJCRYeVoAgB/5cKok/Mgf19nS0LJo0SIrDw8AaELKlrYvLCxURESEn6uB1U6ePClJCg4O9tk5WZ4WANAogoODFRMTU/7gwMjISDkcDj9XBSt4vV4dOXJEkZGRPl3pntACAGg0ZSug88Tj5i8oKEipqak+DaaEFgBAo3E4HEpKSlLr1q1VXFzs73JgobCwMAUF+fQmZEILAKDxBQcH+3SuAwKDbyMSAABAAxFaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALRBaAACALfgstMyZM0cOh0P333+/r04JAACaEZ+ElnXr1umFF15Qnz59fHE6AADQDFkeWvLz83XjjTfqpZdeUmxsrNWnAwAAzZTloWXy5Mm68sorNWrUqDPuW1RUJI/HU6UBAABIUoiVB1+0aJG++eYbrVu3rk77z549W4899piVJQEAAJuyrKclKytLU6ZM0euvv67w8PA6fWfatGlyu93lLSsry6ryAACAzTgMwzCsOPDixYt1zTXXKDg4uHxbaWmpHA6HgoKCVFRUVOWzmng8HrlcLrndbkVHR1tRJgAAaGRW/f22bHjokksu0caNG6tsu+WWW9S9e3f94Q9/OGNgAQAAqMyy0BIVFaVevXpV2daiRQu1atWq2nYAAIAzYUVcAABgC5bePXSqVatW+fJ0AACgGaGnBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2IItQsuBA/6uAAAA+JstQsvWrf6uAAAA+JstQktJib8rAAAA/maL0FJa6u8KAACAv9kitBQX+7sCAADgb7YILQwPAQAAQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFW4QWVsQFAAC2CC30tAAAAFuEFh6YCAAAbBFaGB4CAACWhpb58+erT58+io6OVnR0tAYPHqz333+/3sdheAgAAFgaWtq1a6c5c+YoPT1d69ev18UXX6xx48Zp06ZN9ToOoQUAAIRYefCxY8dWeT9r1izNnz9fX375pXr27Fnn4xBaAACApaGlstLSUr311lsqKCjQ4MGDa9ynqKhIRUVF5e89Ho8kQgsAAPDBRNyNGzeqZcuWcjqduvPOO/XOO++oR48eNe47e/ZsuVyu8paSkiKJ0AIAAHwQWrp166aMjAx99dVXuuuuuzRx4kRt3ry5xn2nTZsmt9td3rKysiRx9xAAAPDB8FBYWJi6dOkiSRowYIDWrVunv/zlL3rhhReq7et0OuV0OqttZ50WAADg83VavF5vlXkrdUFPCwAAsLSnZdq0aRozZoxSU1OVl5enf/7zn1q1apU++OCDeh3n5EmLCgQAALZhaWg5fPiwbr75Zh04cEAul0t9+vTRBx98oEsvvbRexyG0AAAAS0PLK6+80ijHqedoEgAAaIZs8eyhEyf8XQEAAPA3W4QWhocAAIAtQgvDQwAAgNACAABsgdACAABsgdACAABswRahhYm4AADAFqHlxAnD3yUAAAA/s0VoKSlxyOv1dxUAAMCfbBFaJIaIAAAIdLYJLayKCwBAYLNNaOEOIgAAAhuhBQAA2AKhBQAA2IJtQgtzWgAACGy2CS3Hj/u7AgAA4E+2CS0FBf6uAAAA+BOhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2AKhBQAA2IJtQkt+vr8rAAAA/mSb0EJPCwAAgY3QAgAAbIHQAgAAbIHQAgAAbIHQAgAAbIHQAgAAbME2ocXj8XcFAADAn2wTWnJz/V0BAADwJ9uEFrdb8nr9XQUAAPAXS0PL7Nmzdf755ysqKkqtW7fW1Vdfra1btzboWIbBqrgAAAQyS0PLp59+qsmTJ+vLL7/Uhx9+qOLiYl122WUqqOes2lAVSWKICACAQBZi5cGXL19e5X1aWppat26t9PR0/exnP6vzcaLlUY4S5HY3doUAAMAuLA0tp3L/lDri4uJq/LyoqEhFRUXl7z0/3TLkkls5oqcFAIBA5rOJuF6vV/fff7+GDh2qXr161bjP7Nmz5XK5yltKSookKUa5kkRPCwAAAcxnoWXy5Mn6/vvvtWjRolr3mTZtmtxud3nLysqSZA4PSfS0AAAQyHwyPHTPPfdo6dKl+uyzz9SuXbta93M6nXI6ndW2u2R2sdDTAgBA4LI0tBiGoXvvvVfvvPOOVq1apY4dOzboOPS0AAAAS0PL5MmT9c9//lPvvvuuoqKidPDgQUmSy+VSREREnY9T1tNCaAEAIHBZOqdl/vz5crvdGjFihJKSksrbG2+8Ua/jlE3EPXrUgiIBAIAtWD481BjilS1JOnKkUQ4HAABsyBbPHiK0AAAAQgsAALAFQgsAALAFW4SWVsqRJHk8UqVV/gEAQACxRWiJUa6CHaWSpOxsPxcDAAD8whahJUhSqxBzrRZCCwAAgckWoUWSEoLMRVqY1wIAQGCyT2gxDksitAAAEKhsE1pal+yTJB0+7OdCAACAX9gmtCR5zdCyf7+fCwEAAH5hj9ASEaG2MkPLvn1+rgUAAPiFPUJLQgKhBQCAAGeP0BIfT2gBACDA2SO0tGpVJbQ00sOjAQCAjdgjtFTqaSkslNxuP9cDAAB8zjahJUInFBdeIIkhIgAAApE9QktioiSprdN8cCKhBQCAwGOP0JKcLElqG2Qu0rJ3rz+LAQAA/mCP0NK2rSSpQ8kOSdLu3X6sBQAA+IU9QstPPS2dCr+XJO3c6c9iAACAP9gjtCQlSQ6HOpX+IInQAgBAILJHaAkNldq0USeZaYXQAgBA4LFHaJGkdu3KQ8uhQ1JBgZ/rAQAAPmWr0OKSR3GRJyRJu3b5uR4AAOBT9gkt7dtLkjpFHZHEEBEAAIHGPqGlSxdJUqfgHyVJO3b4sxgAAOBr9gkt55xj/ji5SZK0das/iwEAAL5mn9DyU09Lj9wvJEmbN/uzGAAA4Gv2CS2pqVJIiM4t+U4SoQUAgEBjn9ASEiJ17Khu2iqHw1BOjnTkiL+LAgAAvmKf0CJJ55yjSB1Xx1YeSfS2AAAQSOwVWrp1kySd29J8zDOhBQCAwGGv0NKnjySpR+lGSYQWAAACiaWh5bPPPtPYsWOVnJwsh8OhxYsXn90Be/eWJPU8+rkk6bvvzrJAAABgG5aGloKCAvXt21fPPvts4xywRw8pKEj9Cz6TJG3YIHm9jXNoAADQtIVYefAxY8ZozJgxjXfAiAjpnHN07tZMhYeVKi8vWDt2lK87BwAAmrEmNaelqKhIHo+nSqumd2+FqFR92xySJKWn+7hIAADgF00qtMyePVsul6u8paSkVN/pvPMkSf3DzOX8v/nGhwUCAAC/aVKhZdq0aXK73eUtKyur+k4XXCBJ6p/7sSRCCwAAgcLSOS315XQ65XQ6T7/TBRdIDof6Z38gabbS083JuEFNKn4BAIDGZr8/9S6X1L27eul7RYSVKjdX2rLF30UBAACrWRpa8vPzlZGRoYyMDEnSrl27lJGRoT179pzdgS+8UGEq1qCkHyVJa9acZaEAAKDJszS0rF+/Xv369VO/fv0kSVOnTlW/fv30yCOPnN2BBw2SJA3VF5Kk1avP7nAAAKDps3ROy4gRI2QYRuMf+KKLzB/735Q0gdACAEAAsN+cFslcGTcxUYOLP5XDYWjnTunAAX8XBQAArGTP0OJwSBdfLJc86tP6oCTp88/9XBMAALCUPUOLJF1yiSRpRJCZVj76yJ/FAAAAq9k+tFx26DVJ0ooVkhXTZwAAQNNg39DSoYPUpYuGez9WaLBXP/4obdvm76IAAIBV7BtaJGncOLVQoS6KN1eXW7HCz/UAAADL2Du0/PznkqTLPP+WRGgBAKA5s3doGTJEatVKlx1fLElauVI6fty/JQEAAGvYO7SEhEhXXaV+2qCUlkdVWCh9+KG/iwIAAFawd2iRpF/9Sg5J15T+R5L0zjv+LQcAAFjD/qHlssuk+Hhdc/x1SdKSJVJJiZ9rAgAAjc7+oSU0VLr+el2k1WrlzNPRo9Jnn/m7KAAA0NjsH1okacIEhahU40reliS99Zaf6wEAAI2ueYSWCy6QzjlHN5Saq+O+8YZUVOTnmgAAQKNqHqHF4ZDuvFMj9Ynahh7SsWPSe+/5uygAANCYmkdokaRbblFwhFMTihdIkl57zc/1AACARtV8QktsrDR+vG6SmVbee0/KzvZzTQAAoNE0n9AiSXffrZ7arP76RsXF0quv+rsgAADQWJpXaOnfX7rkEt2p+ZKk556TvF4/1wQAABpF8wotkjR9usbrn3IpVzt28BBFAACai+YXWkaMUIsL+2iS0iRJzz7r33IAAEDjaH6hxeGQpk/X3XpOkvTee4Z27PBzTQAA4Kw1v9AiSVdeqa5DEjRGy2QYDj35pL8LAgAAZ6t5hhaHQ/rznzVNsyVJC/7u1f79fq4JAACcleYZWiRp6FANG9dKF+lznSwO0v/7f/4uCAAAnI3mG1okac4cTQv6syRp/t9KlJPj53oAAECDNe/Q0r27xvyul/rpGxWcCNHsx4v9XREAAGig5h1aJDkeeVizWv9VkvS3Z6U9e/xcEAAAaJBmH1oUGanLX/6lhmuVikpDNeOeI/6uCAAANEDzDy2SHGOv0pxRKyVJC/8bp+/XHfdzRQAAoL4CIrRI0oVvPKBfhL8nr4J17zVZMgx/VwQAAOojYEKL4uL01AtRilChVu3rqkV/2ODvigAAQD0ETmiR1OHmn2n64I8lSb+d20ae73b7tyAAAFBnPgktzz77rDp06KDw8HANGjRIX3/9tS9OW6MHl1+qLs4sHTCS9H+jvpROnPBbLQAAoO4sDy1vvPGGpk6dqkcffVTffPON+vbtq9GjR+vw4cNWn7pGzminnn3FKUl65sj1+vQXT4sJLgAANH2Wh5Z58+bptttu0y233KIePXro+eefV2RkpP7+979bfepaXXZja906Zq8k6dfv/0r5M1njHwCAps7S0HLy5Emlp6dr1KhRFScMCtKoUaO0du3aavsXFRXJ4/FUaVZ5alE7pcZ6tFOd9YdHndKbb1p2LgAAcPYsDS3Z2dkqLS1VYmJile2JiYk6ePBgtf1nz54tl8tV3lJSUiyrLTpaeuXNaEnSc5qsJTe+Ia1ebdn5AADA2WlSdw9NmzZNbre7vGVlZVl6vlGjpPvv80qSJpW8pB/H3Cmlp1t6TgAA0DCWhpb4+HgFBwfr0KFDVbYfOnRIbdq0qba/0+lUdHR0lWa1J54M0vkDSnVMcbou/2WdvPRKaeNGy88LAADqx9LQEhYWpgEDBmjlypXl27xer1auXKnBgwdbeeo6CwuT3ngrWDExhr7ShfrdsWlmF8yWLf4uDQAAVGL58NDUqVP10ksvaeHChcrMzNRdd92lgoIC3XLLLVafus46dpTS0hySpL9qil48PE762c+kDayaCwBAU2F5aLnuuus0d+5cPfLIIzrvvPOUkZGh5cuXV5uc62/jxkmPP26+nqzn9MmRntLIkdKaNf4tDAAASJIchtF0V1bzeDxyuVxyu90+md9iGNL48dKiRVJciFtflgzUORH7pLffli6/3PLzAwDQHFj197tJ3T3kbw6H9Pe/S+efLx0tcenSiNXadzxWuuoqaf58f5cHAEBAI7ScIiJC+u9/pS5dpB+PJ2p0zFc6Whot3X239MADUmmpv0sEACAgEVpqkJgoffihlJwsbcptpytSNipfLaSnn5auvlpyu/1dIgAAAYfQUosOHaQVK6S4OOmrrLa6vNtueZwJ0tKl0sCB0rff+rtEAAACCqHlNHr2lJYvl2JipDVb43XZOTuVm9Jb2r5duvBCKS3N3yUCABAwCC1ncP750sqVP/W4fN9So+K+0dFLfiWdOCHdcov0619L+fn+LhMAgGaP0FIH/ftLH38sxcdL6d+G6KL9b+jH3/7VvN1owQLpvPOkGp5aDQAAGg+hpY769pVWrZLatZMyMx268PV7lfHi11JqqrRjh3TRRdIjj0jFxf4uFQCAZonQUg89e5odKr17SwcPSsMeGKgVT2+WJkyQvF5p5kxp0CDpm2/8XSoAAM0OoaWe2rWTPv/cXOE/P18a88sWmtfvNRn/WiTFxprPK7rgAun3v5cKC/1dLgAAzQahpQFcLvOuokmTzA6W3/5WmvDf61S4frN07bXmAnRPPml2yXz4ob/LBQCgWSC0NFBYmLnk/zPPSMHB0j//KQ39nzbaOfsNc0nddu2knTulyy6TfvlLafduf5cMAICtEVrOgsMh3XOPeUt0QoKUkSH16yctyr9K2rxZuu8+M9H85z9S9+7Sww9LBQX+LhsAAFsitDSC4cOl9euloUMlj0e64Qbp11OiVPCnv5hJ5uKLpaIi6Y9/lLp1k157jWcYAQBQT4SWRpKaat4S/fDDFcu3DBggfVXQS/roI+ntt81nA+zbJ918s9kls3SpZBj+Lh0AAFsgtDSikBDp8cfNhejatpW2bpWGDJF+/weHjl9+jZSZKc2ebT4XYONGaexYc32Xzz7zd+kAADR5hBYLjBghffddxfItTz5pLpr7xTfh0kMPmRN0H3pIioiQvvjCHF8aPVpavdrfpQMA0GQRWiwSF2dOXVmyREpOln74wexUuesu6agRa/a4bN8u3Xmn2UWzYoU0bJiZeD76iGEjAABOQWix2Nix0qZN5rMVDUN6/nmpa1fplVckb5tkaf58cxzp9tul0FDp00+lSy+VBg9mzgsAAJUQWnwgJsZc02XVKvNRADk50q23mvNd1q+X1KmT9MIL5jOM7rtPCg+XvvrKTDw9e0ovvigdP+7n3wIAAP8itPjQ8OHmKv/z5klRUWYuOf98afx4adcuSSkp0l/+Yi5E9/vfmztlZkp33GF+9n//Jx044O9fAwAAvyC0+FhoqPTAA9KWLeZEXUn617/M5VseeMDshVFiovTEE9LevWbC6dDB/GDWLKl9e/OW6a++YugIABBQCC1+kpxsTtT95htp1CipuFh6+mmpc2dzDTq3W1J0tJlktm2T/v1vc/W64mLzixdeKPXvb06Sycvz968DAIDlHIbRdP9z3ePxyOVyye12Kzo62t/lWGrFCnNE6NtvzfcxMdL990tTppivy339tfS3v0lvvmmusitJLVuaY0x33GEGGQAA/Miqv9+ElibE65XeeEOaOdOcyiKZT5SeMsVscXGVdj56VHr1VbOnZevWiu39+0sTJ5rPEkhI8Gn9AABIhJaACC1lSkvN0aCZM83bpSUpMlL69a/N3pfOnSvtbBjmirrPP28+mLG42NweEiJdeaU5/+XKKyWn09e/BgAgQBFaAii0lPF6zUcW/fGPFcNGQUHSNddIv/2tuZRLFdnZ5qzehQul9PSK7XFx0vXXSzfdJA0aZD4cCQAAixBaAjC0lDEM83lGc+dKy5dXbB88WJo8Wfqf/zGXdqli0yZz+Ogf/5D276/Y3r69dO210nXXmUNJBBgAQCMjtARwaKls0ybzLuh//EM6edLc1qqVOXR0xx2nDB1J5ljTypVm78u770oFBRWfde5cEWD69CHAAAAaBaGF0FLFwYPSSy+Zi+Xu3Vux/bLLzPBy1VVSWNgpXzp+XFq2zJztu3Rp1VV2u3WTfvEL6eqrpYEDzXEoAAAagNBCaKlRSYmZQ+bPlz74oGK9ufh48waiiRNrGQXKz5fee88MMMuWVdw+LUlJSdK4cWYbOZJJvACAeiG0EFrOaOdOs+dl4UKzJ6ZMz57mTUQTJpiL2lXj8Zg9L+++awaY/PyKz6KipDFjzABzxRWnLBoDAEB1hBZCS52VlEgffmiGl8WLKzpRgoKkESPMaSy/+EUty7gUFUmffGIGmHffrfqso+Bgc/bv5ZebQea88xhGAgBUY7vQMmvWLL333nvKyMhQWFiYcnNz630MQsvZy82V3nrLDDBr1lRsDw42R36uvda8hTo+voYve73mY6jffddMP5s3V/08MVEaPdoMMJdeas4IBgAEPNuFlkcffVQxMTHau3evXnnlFUJLE7Bzp7lo3ZtvVl3GJThYuvhiM7yMHSu1a1fLAXbvNu+5Xr7cvCOp8jBSUJB0wQVmiLnkEnM9mGozgQEAgcB2oaVMWlqa7r//fkJLE7Njh9kD89Zb5kMbK+vfX/r5z80A069fLXdCnzwprV5tBpj335e+/77q5y1aSMOGmWnokksYSgKAABIQoaWoqEhFle5i8Xg8SklJIbRYbPt28wkA//2v9MUXFXcgSWavy9ix5i3Uw4ebWaRGWVlmgPnoI3MlvOzsqp/HxZkTai65xGxdu7IuDAA0UwERWmbMmKHHHnus2nZCi+8cPmzeQLRkiXkLdWFhxWdhYdJFF5lrwVx2mdS3by2dJ16v2fOycqXZPv206lCSJLVpY/bElLXevc1xKgCA7TWJ0PLQQw/piSeeOO0+mZmZ6t69e/l7elrs68SJihuJli+Xfvyx6uetW5vzby+7TBo1qpbbqSXzIY7r15sB5uOPzRnBZcv5lnG5pCFDpJ/9zAwxAweyPgwA2FSTCC1HjhxRTk7Oaffp1KmTwipNwGROS/NgGNK2bdKKFWb7+OOqTwSQpHPOMUeAhg83W60Teo8fl77+Wvr8c7N98UX1nhin05zMO2yYGWYGDeLuJACwiSYRWhqC0NI8nTwprV1rBpgPPjAn8576b1LnzhUhZsQIKSWlloOVlJiPsS4LMZ9/Lh05Un2/c86RLrywovXuLYWGNvJvBgA4W7YLLXv27NHRo0e1ZMkSPfnkk/r8888lSV26dFHLli3rdAxCi33k5po3E61aZbYNG8ypLZWlpppr0w0ebHaenHdeLZnDMKQffjDDy+rV0pdfSlu3Vt8vIsIcRqocZGodowIA+IrtQsukSZO0cOHCats/+eQTjRgxok7HILTYl9tt5o1PPzVDTHp69RATHm5mjiFDKsJMYmItBzx61BxS+vJLs331lZmUTtWunXT++dKAARWtxqV/AQBWsV1oaQyEluYjL8/MHGvXVrRjx6rv17GjGWQGDDB/9u8vxcbWcECv1+yNKQsxX34pbdxYPRlJZhdPWYApO3iNSwADABoDoYXQ0qyUZY7KIWbTpurzYiSpU6fqQabG5zbm55tdOmVt/XrzJDUpCzJlB+zb17wNm7VjAOCsEVoILc2e223mjLK8kZ5uPnqgJp07m3Ni+vQxW9++Uvv2Nawb4/GYE2zKDpieXnuQSUgwD1TWzjtP6t6dyb4AUE+EFkJLQDp61LwzqXKQ2bWr5n2joswbispCTJ8+5vuoqFN2dLvNIFN20IwMM8jUNLQUFib16FE1zPTty+3XAHAahBZCC36Sk2Nmju++M++U/u478wHUp65XV6ZjR6lnTzN7nHtuRavyr1RhoTk+9e23ZogpO7DHU/NB27Y1U1HPnhWtR4/TPOcAAAIHoYXQgtMoLjY7S8qyRlnbt6/277Rta4aXymGmR49KNxsZhvlk62+/rdpqG7OSKhJSWevVyxxiiohozF8XAJo0QguhBQ2QnW3eVJSZafbGZGaa7cCB2r/TqpXUrZu5ll3l1qXLT0NNHo950O+/N3tnyn4ePlzzAYOCzNnElYNMjx7mQSMjLfm9AcCfCC2EFjSi3NyKAFM5zOzeXfMdTGXatKkeZsoCTWRhthleKgeZTZvM8azatG9vJqTu3c1W9jopiTuZANgWoYXQAh8oLDQX3922rXqr6ckClSUlmaNDVVoHQx2jstUu93uFbPm+IsxkZpqzjGsTFVVzmOnSxVyVDwCaMEILoQV+lpsrbd9ec6A5Xf6QpOBg89lLVQJNfJ46GjvVsXCTEg9kKGhrppmYduyo+U4myRxq6tChYvyqS5eKn+3bc3s2gCaB0EJoQROWk2Peil1T27279jubyoSGmhODU1OllORSpbQ8phTHXqWc2KbUY98q5cDXit32tRwed+0HCQkxA03lIFP2ukMHAg0AnyG0EFpgU16vOfG3LMCcGmqysmrvWKksMtJQSlKpUmLylBp+SClGllIKtyrl2HdKOvCNkk/uUpyOqsaZMMHBZk/Mqb0zZT00DDkBaESEFkILmqmSEjPUZGVJe/aYP8ta2fszzacpExbqVXJ0gZKcR5Ws/Uou2qVkzxYlFf9ovv+pxSi3arhJTjbvcCobu6r8um3bGpYaBoDaEVoILQhgJ05Ie/fWHGz27jVDT3Z23Y8XHnRSSSGHlVySpWRvlpJ0QIk6VN5a63D56/Aww+yNqRxmKoea2FjudAJQBaGF0AKcVlGRdPCgtH+/2Q4cqHhd+f2ZJg2fKlruKiHm1FCT2KJArVPDldi5paK6JMrRPtUMOak//WzVilADBBhCC6EFaBQnTlQEmMrB5vBh6dChinb48JknEJ8qXMfVWocVr+yKFuJWfEyJ4hMcim8bpvjUForvEqP47vFq1StJYR3bmpOIATQbhBZCC+BThmE+W/LUIFP5/aEDpTp8oFSHjgSp4ETDgke03IoPzlV8RIHio08qPs6rVokhim/nVHyHKMWfE6tW7SIUGyvFxZmjUZGRdN4ATZlVf7/5zxsANXI4pJgYs3XrVttewT81qaCgItjk5JhzbLIPlih7V56y9xQq+2CJud0dquzCCOUUR8urYHnkkqfUpZ35kvIl7Zf0/elrCw0qVVzLIsW6vIpr5VBs6zDFtg5RXJyjSrip/LrsZ1hYI/0DAuBzhBYAjaJFC3N+bqdOlbeGSIr9qVXl9Uq5R73K3pKt7Mwjyt52TNm785S9r0jZh7zKPhqknLwwHSmJ0VHF6ZhidVRxKlGoir3BOuSJ1CGPpKz611k50LhcZouJqXh9ukYvD+A/DA8BaNry883Hde/dKyNrrwp2HNTRnbk6lpWvY/uP6+ihYh3LDykPNZV/Vn6dqxgZOvtbt0NCpOjouoccl8vcPyqq4mdUFNN40Lwxp4XQAqA2x4+XB5sqreye8L175T18RG65qoUbt1xVWq5i5A6NN1tQnNyKlrukhdxF4fIajbdeTXh4RYA529aiBUvpoGlhTgsA1CYiomKF31oEFRUpdv9+xe7dq05lgab8FqrMituoTpyQimW2SgxJBWpRNeREpcgd0165USlyR7aRO6y13CFxcjti5C6Nkrs4Qu7jYfLkBSkvT8rLM29Nl8zTnDhR94UDz6Rly4oQ07KlGWQqt5q21eUznv6ApoSeFgAoYxjmkzFPXeSmpoVv6nM/eFyc1KaNlJio4tZtlRfXXnnRbZUXlay8yETlhScoLzROeUEu5RVUBByPR+Wva2t1eQTE2QgLq3/QiYgw5/5U/lnTtrKfwcHW/g7wPYaHCC0AmgrDMFfpO1O4OXBAKi4+8/HKBAVJ8fHlAUdt2tT+Oi5OhiNIx49XDzL5+ebdXDW1031WeZ/SUuv+8Z0qLOzMwaYu4efUfcLDa26EJOsRWggtAOzG65WOHTOXKj540LwnvPLPyq8PHzbDUF2FhEitW1cPM61bSwkJ5s+y1wkJ9RrnMQyzI6khwaew0GzHj5ut7PWpP48fb8A/z0YSElJ7oLG6OZ2BMQmb0EJoAdCclZaai9ucGmZqep2TU//jx8ZWDzS1vY+Ls7w7wjDMOT21hZq6BJ8z/SybN3TihPlg0qYiONgML2Fh5s+6vK7Pvg09Rmho493Oz0RcAGjOgoPN3pLERKlv39Pve/KkOYO3plBz5IjZa1P5Z1mPz7Fj0tatZ66lbJjqdAGn8muXq95/7RyOiqEcXygpMSdBVw4ytbW67lefVnmUsLS0okeqqakt7NTWQkNr3m5Vdwg9LQDQnJUFlsOHK1pZoDn19eHD9X+ipmSOd8TH195ataq+rUWLgFqlr7S0ahgqKjKzZ1FR/V835vesm7vkkcTwEADASsXF5vBTTYGmprDj8TTsPE5n3QNOWfNVt0wAKS2tW/ApLjbfl7VT35/a8vI8euYZQgsAoCk5ceKnB01lV3ro1GnakSMVi9XUV2Rk7SGnVStzLk5cXNXXDRi6wtljTgsAoOkJD5fatTNbXRiGOZnj1DBzpsBTXGx+b88es9VVcHDFg6ZqCzY1fRYdTdhpgggtAADfcTgqVqFr375u3zEMcwGa04Wao0crWk6O+bOwsOKurOzs+tVZFnbqE3Ti4gg7FiO0AACaNofDDAPR0ac+Rvz0TpyoOczU9L7y67MNOzExFY8SL2t12eZy8RCpMyC0AACap/BwKTnZbPVROezUNejk5JiLw5SWmq8bspZOWTirb9gpex8AD4qyLLTs3r1bM2fO1Mcff6yDBw8qOTlZEyZM0PTp0xUWFmbVaQEAODsNDTvHj5sB5tgx8xlWZWvjlLVTt1V+X1hoDoO53Wbbvbv+dbdoceZgExNjNper6uvoaFv08lgWWrZs2SKv16sXXnhBXbp00ffff6/bbrtNBQUFmjt3rlWnBQDAPyIipLZtzVZfJ0+eOdjUtq3stvOyZyns3Vv/8zsc5iPCaws1p76uaZsPOiR8esvzk08+qfnz52vnzp112p9bngEAOIPSUrN3pi69Om63+To3t+L1iRONU0dERHmA8bRsKdf69fa+5dntdisuLq7Wz4uKilRU6f59T0MXLQIAIFAEB1fcxdQQRUU1h5mattX0ednf6rKHRR08eNa/Um18Flq2b9+uZ5555rRDQ7Nnz9Zjjz3mq5IAAIDTWfEcqYYoLTVvSa8cavbvl8aPb8wqJTVgeOihhx7SE088cdp9MjMz1b179/L3+/bt0/DhwzVixAi9/PLLtX6vpp6WlJQUhocAALARq6Z31Du0HDlyRDlnuJWrU6dO5XcI7d+/XyNGjNCFF16otLQ0BdVjdjJzWgAAsJ8ms4x/QkKCEhIS6rTvvn37NHLkSA0YMEALFiyoV2ABAACozLI5Lfv27dOIESPUvn17zZ07V0eOHCn/rE2bNladFgAANFOWhZYPP/xQ27dv1/bt29XulAdpNeEHSwMAgCbKsvGaSZMmyTCMGhsAAEB9MckEAADYAqEFAADYAqEFAADYAqEFAADYAqEFAADYAqEFAADYAqEFAADYAqEFAADYgmUr4jaGsoXoPB6PnysBAAB1VfZ3u7EXlG3SoSUvL0+SlJKS4udKAABAfeXk5MjlcjXa8RxGE15X3+v1av/+/YqKipLD4fB3OQHN4/EoJSVFWVlZjfqYcdQf16Jp4Xo0HVyLpsPtdis1NVXHjh1TTExMox23Sfe0BAUFVXvYIvwrOjqa/zNoIrgWTQvXo+ngWjQdQUGNO3WWibgAAMAWCC0AAMAWCC2oE6fTqUcffVROp9PfpQQ8rkXTwvVoOrgWTYdV16JJT8QFAAAoQ08LAACwBUILAACwBUILAACwBUILAACwBUILyj377LPq0KGDwsPDNWjQIH399de17vvSSy9p2LBhio2NVWxsrEaNGnXa/VE/9bkWlS1atEgOh0NXX321tQUGkPpei9zcXE2ePFlJSUlyOp3q2rWrli1b5qNqm7/6Xo+nn35a3bp1U0REhFJSUvTAAw/oxIkTPqq2+frss880duxYJScny+FwaPHixWf8zqpVq9S/f385nU516dJFaWlp9T+xARiGsWjRIiMsLMz4+9//bmzatMm47bbbjJiYGOPQoUM17j9+/Hjj2WefNTZs2GBkZmYakyZNMlwul7F3714fV9781PdalNm1a5fRtm1bY9iwYca4ceN8U2wzV99rUVRUZAwcONC44oorjNWrVxu7du0yVq1aZWRkZPi48uapvtfj9ddfN5xOp/H6668bu3btMj744AMjKSnJeOCBB3xcefOzbNkyY/r06cbbb79tSDLeeeed0+6/c+dOIzIy0pg6daqxefNm45lnnjGCg4ON5cuX1+u8hBYYhmEYF1xwgTF58uTy96WlpUZycrIxe/bsOn2/pKTEiIqKMhYuXGhViQGjIdeipKTEGDJkiPHyyy8bEydOJLQ0kvpei/nz5xudOnUyTp486asSA0p9r8fkyZONiy++uMq2qVOnGkOHDrW0zkBTl9Dy+9//3ujZs2eVbdddd50xevToep2L4SHo5MmTSk9P16hRo8q3BQUFadSoUVq7dm2djlFYWKji4mLFxcVZVWZAaOi1ePzxx9W6dWv95je/8UWZAaEh12LJkiUaPHiwJk+erMTERPXq1Ut/+tOfVFpa6quym62GXI8hQ4YoPT29fAhp586dWrZsma644gqf1IwKa9eurXLtJGn06NF1/htTpkk/MBG+kZ2drdLSUiUmJlbZnpiYqC1bttTpGH/4wx+UnJxc7V9K1E9DrsXq1av1yiuvKCMjwwcVBo6GXIudO3fq448/1o033qhly5Zp+/btuvvuu1VcXKxHH33UF2U3Ww25HuPHj1d2drYuuugiGYahkpIS3Xnnnfrf//1fX5SMSg4ePFjjtfN4PDp+/LgiIiLqdBx6WnDW5syZo0WLFumdd95ReHi4v8sJKHl5ebrpppv00ksvKT4+3t/lBDyv16vWrVvrxRdf1IABA3Tddddp+vTpev755/1dWkBatWqV/vSnP+m5557TN998o7ffflvvvfeeZs6c6e/S0ED0tEDx8fEKDg7WoUOHqmw/dOiQ2rRpc9rvzp07V3PmzNFHH32kPn36WFlmQKjvtdixY4d2796tsWPHlm/zer2SpJCQEG3dulWdO3e2tuhmqiH/u0hKSlJoaKiCg4PLt5177rk6ePCgTp48qbCwMEtrbs4acj0efvhh3XTTTbr11lslSb1791ZBQYFuv/12TZ8+XUFB/He7r7Rp06bGaxcdHV3nXhaJnhZICgsL04ABA7Ry5crybV6vVytXrtTgwYNr/d6f//xnzZw5U8uXL9fAgQN9UWqzV99r0b17d23cuFEZGRnl7ec//7lGjhypjIwMpaSk+LL8ZqUh/7sYOnSotm/fXh4cJemHH35QUlISgeUsNeR6FBYWVgsmZYHS4LF7PjV48OAq106SPvzww9P+jalR/eYIo7latGiR4XQ6jbS0NGPz5s3G7bffbsTExBgHDx40DMMwbrrpJuOhhx4q33/OnDlGWFiY8e9//9s4cOBAecvLy/PXr9Bs1PdanIq7hxpPfa/Fnj17jKioKOOee+4xtm7daixdutRo3bq18cc//tFfv0KzUt/r8eijjxpRUVHGv/71L2Pnzp3GihUrjM6dOxvXXnutv36FZiMvL8/YsGGDsWHDBkOSMW/ePGPDhg3Gjz/+aBiGYTz00EPGTTfdVL5/2S3Pv/vd74zMzEzj2Wef5ZZnnJ1nnnnGSE1NNcLCwowLLrjA+PLLL8s/Gz58uDFx4sTy9+3btzckVWuPPvqo7wtvhupzLU5FaGlc9b0WX3zxhTFo0CDD6XQanTp1MmbNmmWUlJT4uOrmqz7Xo7i42JgxY4bRuXNnIzw83EhJSTHuvvtu49ixY74vvJn55JNPavwbUPbPf+LEicbw4cOrfee8884zwsLCjE6dOhkLFiyo93kdhkEfGQAAaPqY0wIAAGyB0AIAAGyB0AIAAGyB0AIAAGyB0AIAAGyB0AIAAGyB0AIAAGyB0AIAQAD67LPPNHbsWCUnJ8vhcGjx4sX1PoZhGJo7d666du0qp9Optm3batasWY1f7E94YCIAAAGooKBAffv21a9//Wv94he/aNAxpkyZohUrVmju3Lnq3bu3jh49qqNHjzZypRVYERcAgADncDj0zjvv6Oqrry7fVlRUpOnTp+tf//qXcnNz1atXLz3xxBMaMWKEJCkzM1N9+vTR999/r27duvmkToaHAABANffcc4/Wrl2rRYsW6bvvvtOvfvUrXX755dq2bZsk6b///a86deqkpUuXqmPHjurQoYNuvfVWS3taCC0AAKCKPXv2aMGCBXrrrbc0bNgwde7cWQ8++KAuuugiLViwQJK0c+dO/fjjj3rrrbf06quvKi0tTenp6frlL39pWV3MaQEAAFVs3LhRpaWl6tq1a5XtRUVFatWqlSTJ6/WqqKhIr776avl+r7zyigYMGKCtW7daMmREaAEAAFXk5+crODhY6enpCg4OrvJZy5YtJUlJSUkKCQmpEmzOPfdcSWZPDaEFAABYrl+/fiotLdXhw4c1bNiwGvcZOnSoSkpKtGPHDnXu3FmS9MMPP0iS2rdvb0ld3D0EAEAAys/P1/bt2yWZIWXevHkaOXKk4uLilJqaqgkTJmjNmjV66qmn1K9fPx05ckQrV65Unz59dOWVV8rr9er8889Xy5Yt9fTTT8vr9Wry5MmKjo7WihUrLKmZ0AIAQABatWqVRo4cWW37xIkTlZaWpuLiYv3xj3/Uq6++qn379ik+Pl4XXnihHnvsMfXu3VuStH//ft17771asWKFWrRooTFjxuipp55SXFycJTUTWgAAgC1wyzMAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALAFQgsAALCF/w/NFt0Mk6EvuwAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(df_ln[\"z\"], df_ln[\"theta_z\"], 'r-')\n",
+    "ax.plot(df_Pln[\"z\"], df_Pln[\"theta_z\"], 'b-')\n",
+    "ax.set_xlim(left=1000, right=1_000_000)\n",
+    "ax.legend([\"log normal\", \"Pln\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "e90ef79d",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.legend.Legend at 0x2aa7617d0>"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHACAYAAABeV0mSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEGElEQVR4nO3dd3hUZd7/8c/MJJn0hJAEEgiEgKFIMciiiAUURUBW1hVUVGCxrI+oKLqPsq4FC+Iu8nN1WbuAPqysBRARsYAsiqwKEkRK6IYeWjppM+f3x0kGQgnpkznzfl3XXJlz5pTvnEDmM/d9n3NshmEYAgAAsAi7twsAAACoT4QbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKX4dbpYvX66hQ4cqMTFRNptN8+fPb/B97tmzR7fccouaN2+ukJAQdevWTatWrWrw/QIA4C/8OtwUFBSoR48emj59eqPs7+jRo+rbt68CAwP12WefacOGDXrhhRfUrFmzRtk/AAD+wMaNM002m03z5s3TsGHDPPOKi4v16KOP6r333lN2dra6du2q559/Xv369avVPh555BGtWLFC33zzTf0UDQAATuHXLTdnc88992jlypWaM2eOfv75Zw0fPlxXX321tmzZUqvtLViwQL169dLw4cMVHx+vtLQ0vfHGG/VcNQAA/o2Wm3Int9xkZmYqJSVFmZmZSkxM9Cw3YMAA9e7dW5MnT67xPoKDgyVJEyZM0PDhw/Xjjz9q/PjxevXVVzV69Oh6eR8AAPi7AG8X0FStW7dOLpdLqampleYXFxerefPmkqRNmzapc+fOVW7n4Ycf1pQpUyRJbrdbvXr18gSjtLQ0/fLLL4QbAADqEeHmDPLz8+VwOLR69Wo5HI5Kr4WHh0uSUlJStHHjxiq3UxGEJCkhIUFdunSp9Hrnzp310Ucf1VPVAACAcHMGaWlpcrlcysrK0iWXXHLaZYKCgtSpU6dqb7Nv377KyMioNG/z5s1q27ZtnWoFAADH+XW4yc/P19atWz3TO3bsUHp6umJiYpSamqqbb75Zo0aN0gsvvKC0tDQdPHhQS5YsUffu3TVkyJAa7++BBx7QRRddpMmTJ2vEiBH64Ycf9Prrr+v111+vz7cFAIBf8+sBxcuWLVP//v1PmT969GjNnDlTpaWleuaZZ/TOO+9oz549io2N1YUXXqhJkyapW7dutdrnwoULNXHiRG3ZskXt2rXThAkTdMcdd9T1rQAAgHJ+HW4AAID1cJ0bAABgKYQbAABgKX43oNjtdmvv3r2KiIiQzWbzdjkAAKAaDMNQXl6eEhMTZbdX3Tbjd+Fm7969SkpK8nYZAACgFnbt2qXWrVtXuYzfhZuIiAhJ5sGJjIz0cjUAAKA6cnNzlZSU5Pkcr4rfhZuKrqjIyEjCDQAAPqY6Q0oYUAwAACyFcAMAACyFcAMAACzF78bcAAB8k8vlUmlpqbfLQAMKCgo662ne1UG4AQA0aYZhaP/+/crOzvZ2KWhgdrtd7dq1U1BQUJ22Q7gBADRpFcEmPj5eoaGhXIDVoiousrtv3z61adOmTr9nwg0AoMlyuVyeYNO8eXNvl4MGFhcXp71796qsrEyBgYG13g4DigEATVbFGJvQ0FAvV4LGUNEd5XK56rQdwg0AoMmjK8o/1NfvmXADAAAshXADAEAD6Nevn+6//35vl9HkLFu2TDabrUHPfiPcAAAASyHcAACASkpKSrxdQp0QbgAAaARHjx7VqFGj1KxZM4WGhmrQoEHasmVLpWXeeOMNJSUlKTQ0VL/73e80bdo0RUdHn3GbO3fulM1m09y5c9W/f3+FhoaqR48eWrlyZaXlPvroI5177rlyOp1KTk7WCy+8UOn15ORkPf300xo1apQiIyN15513aubMmYqOjtbChQvVsWNHhYaG6vrrr1dhYaFmzZql5ORkNWvWTPfdd1+ls5veffdd9erVSxEREWrZsqVGjhyprKysuh/AGiDcAAB8imEYKiwpa/SHYRh1qnvMmDFatWqVFixYoJUrV8owDA0ePNhzuvuKFSt01113afz48UpPT9eVV16pZ599tlrbfvTRR/XQQw8pPT1dqampuummm1RWViZJWr16tUaMGKEbb7xR69at05NPPqnHHntMM2fOrLSNqVOnqkePHlqzZo0ee+wxSVJhYaFeeuklzZkzR4sXL9ayZcv0u9/9TosWLdKiRYv07rvv6rXXXtOHH37o2U5paamefvpprV27VvPnz9fOnTs1ZsyYOh27muIifgAAn3Ks1KUuj3/e6Pvd8NRAhQbV7mNzy5YtWrBggVasWKGLLrpIkjR79mwlJSVp/vz5Gj58uF5++WUNGjRIDz30kCQpNTVV3333nRYuXHjW7T/00EMaMmSIJGnSpEk699xztXXrVnXq1EnTpk3TFVdc4Qksqamp2rBhg/72t79VCh2XX365HnzwQc/0N998o9LSUr3yyitq3769JOn666/Xu+++qwMHDig8PFxdunRR//799fXXX+uGG26QJI0dO9azjZSUFL300kv6zW9+o/z8fIWHh9fq+NUULTcAADSwjRs3KiAgQBdccIFnXvPmzdWxY0dt3LhRkpSRkaHevXtXWu/k6TPp3r2753lCQoIkebqCNm7cqL59+1Zavm/fvtqyZUul7qRevXqdst3Q0FBPsJGkFi1aKDk5uVJIadGiRaVup9WrV2vo0KFq06aNIiIidNlll0mSMjMzq/Ve6gMtNwAAnxIS6NCGpwZ6Zb9N1Ym3Kqi4EJ7b7a7RNsLCwqrcbsW2TzevYl8FBQUaOHCgBg4cqNmzZysuLk6ZmZkaOHBgow5SJtwAAHyKzWardfeQt3Tu3FllZWX6/vvvPd1Shw8fVkZGhrp06SJJ6tixo3788cdK6508Xdt9r1ixotK8FStWKDU1VQ5H/Qa2TZs26fDhw5oyZYqSkpIkSatWrarXfVQH3VIAADSwc845R9dee63uuOMOffvtt1q7dq1uueUWtWrVStdee60k6d5779WiRYs0bdo0bdmyRa+99po+++yzOt+S4MEHH9SSJUv09NNPa/PmzZo1a5b+8Y9/eMb21Kc2bdooKChIL7/8srZv364FCxbo6aefrvf9nA3hBgCARjBjxgydf/75uuaaa9SnTx8ZhqFFixZ5unn69u2rV199VdOmTVOPHj20ePFiPfDAAwoODq7Tfnv27Kn3339fc+bMUdeuXfX444/rqaeeapAzmOLi4jRz5kx98MEH6tKli6ZMmaKpU6fW+37OxmbU9dw2H5Obm6uoqCjl5OQoMjLS2+UAAKpQVFSkHTt2qF27dnX+kPdFd9xxhzZt2qRvvvnG26U0iqp+3zX5/PZqy83y5cs1dOhQJSYmymazaf78+WddZ/bs2erRo4dCQ0OVkJCgsWPH6vDhww1fLAAADWzq1Klau3attm7dqpdfflmzZs3S6NGjvV2Wz/FquCkoKFCPHj00ffr0ai2/YsUKjRo1SrfddpvWr1+vDz74QD/88IPuuOOOBq4UAICG98MPP+jKK69Ut27d9Oqrr+qll17S7bff7u2yfI5Xh5sPGjRIgwYNqvbyK1euVHJysu677z5JUrt27fTHP/5Rzz//fEOVCABAo3n//fe9XYIl+NSA4j59+mjXrl1atGiRDMPQgQMH9OGHH2rw4MHeLg0AADQRPhVu+vbtq9mzZ+uGG25QUFCQWrZsqaioqCq7tYqLi5Wbm1vpAQAArMunws2GDRs0fvx4Pf7441q9erUWL16snTt36q677jrjOs8995yioqI8j4qLCgEAAGtqMqeC22w2zZs3T8OGDTvjMrfeequKior0wQcfeOZ9++23uuSSS7R3717P/TROVFxcrOLiYs90bm6ukpKSOBUcAHyAv58K7m/q61Rwn7p+dWFhoQICKpdccenoM2U0p9Mpp9PZ4LUBAICmwavdUvn5+UpPT1d6erokaceOHUpPT/fcOXTixIkaNWqUZ/mhQ4dq7ty5euWVV7R9+3atWLFC9913n3r37q3ExERvvAUAANDEeDXcrFq1SmlpaUpLS5MkTZgwQWlpaXr88cclSfv27at0i/QxY8Zo2rRp+sc//qGuXbtq+PDh6tixo+bOneuV+gEAqI2ZM2cqOjra22VYlle7pfr163fG7iTJ/OWf7N5779W9997bgFUBAFB3Y8aM0axZsyRJgYGBatOmjUaNGqU///nPXq7M+nxqzA0AAL7k6quv1owZM1RcXKxFixZp3LhxCgwMPO0JMKg/PnUqOAAAvsTpdKply5Zq27at/ud//kcDBgzQggULTlnuySef1Hnnnad3331XycnJioqK0o033qi8vDwvVO37aLkBAPgWw5BKCxt/v4Ghks1Wp02EhISc8WbP27Zt0/z587Vw4UIdPXpUI0aM0JQpU/Tss8/WaZ/+iHADAPAtpYXSZC+cIfvnvVJQWK1WNQxDS5Ys0eeff37GcaNut1szZ85URESEJPPabkuWLCHc1ALdUgAANJCFCxcqPDxcwcHBGjRokG644QY9+eSTp102OTnZE2wkKSEhQVlZWY1UqbXQcgMA8C2BoWYrijf2W0P9+/fXK6+8oqCgICUmJp5yIdpKmw8MrDRts9nkdrtrvE8QbgAAvsZmq3X3UGMLCwtThw4dvF2G36FbCgAAWArhBgAAWEqTuSt4Y6nJXUUBAN7FXcH9S33dFZyWGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwBAk+dn5774rfr6PRNuAABNVsVVewsLvXCjTDS6kpISSZLD4ajTdrhCMQCgyXI4HIqOjvbcYyk0NFS2Ot6ZG02T2+3WwYMHFRoaWuVtKqqDcAMAaNJatmwpSdxE0g/Y7Xa1adOmzgGWcAMAaNJsNpsSEhIUHx+v0tJSb5eDBhQUFCS7ve4jZgg3AACf4HA46jwWA/6BAcUAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSvBpuli9frqFDhyoxMVE2m03z588/6zrFxcV69NFH1bZtWzmdTiUnJ+vtt99u+GIBAIBPCPDmzgsKCtSjRw+NHTtW1113XbXWGTFihA4cOKC33npLHTp00L59++R2uxu4UgAA4Cu8Gm4GDRqkQYMGVXv5xYsX6z//+Y+2b9+umJgYSVJycnIDVQcAAHyRT425WbBggXr16qW//vWvatWqlVJTU/XQQw/p2LFj3i4NAAA0EV5tuamp7du369tvv1VwcLDmzZunQ4cO6e6779bhw4c1Y8aM065TXFys4uJiz3Rubm5jlQsAALzAp1pu3G63bDabZs+erd69e2vw4MGaNm2aZs2adcbWm+eee05RUVGeR1JSUiNXDQAAGpNPhZuEhAS1atVKUVFRnnmdO3eWYRjavXv3adeZOHGicnJyPI9du3Y1VrkAAMALfCrc9O3bV3v37lV+fr5n3ubNm2W329W6devTruN0OhUZGVnpAQAArMur4SY/P1/p6elKT0+XJO3YsUPp6enKzMyUZLa6jBo1yrP8yJEj1bx5c/3hD3/Qhg0btHz5cv3pT3/S2LFjFRIS4o23AAAAmhivhptVq1YpLS1NaWlpkqQJEyYoLS1Njz/+uCRp3759nqAjSeHh4fryyy+VnZ2tXr166eabb9bQoUP10ksveaV+AADQ9NgMwzC8XURjys3NVVRUlHJycuiiAgDAR9Tk89unxtwAAACcDeEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYiv+GmyM7vV0BAABoAP4bbrZ+4e0KAABAA/DfcON2e7sCAADQAPw33Bgub1cAAAAagP+GG3eZtysAAAANwI/DDd1SAABYkf+GG7qlAACwJP8NN3RLAQBgSX4cbuiWAgDAivw43NByAwCAFflvuDFouQEAwIr8ONwwoBgAACvy33BDtxQAAJbkv+GGlhsAACzJf8ONm3ADAIAV+XG4YUAxAABW5L/hhm4pAAAsyX/DDQOKAQCwJD8ON7TcAABgRYQbAABgKf4bbrhCMQAAluS/4YaWGwAALMl/w43BgGIAAKzIj8MNLTcAAFiR/4YbuqUAALAkPw43DCgGAMCK/Dfc0C0FAIAl+W+44QrFAABYkh+HG1puAACwIsINAACwFP8NN4y5AQDAkvw33HC2FAAAluS/4cYo9XYFAACgAfhvuHFxthQAAFbkx+GmxNsVAACABuC/4Ybr3AAAYEn+G25cjLkBAMCK/Djc0C0FAIAV+XG4oeUGAAAr8t9wIzdXKQYAwIL8ONyI1hsAACzIz8MN424AALCaWoWbsWPHatasWafMz83N1dixY+tcVKOh5QYAAMupVbiZOXOm7r77bt13331yn3CPpmPHjp029DRZtNwAAGA5te6W+vTTT7Vo0SINHDhQR48erc+aGo+blhsAAKym1uGmS5cu+v7771VaWqrevXtr48aNNd7G8uXLNXToUCUmJspms2n+/PnVXnfFihUKCAjQeeedV+P9etAtBQCA5dQq3NhsNklS8+bN9dVXX+myyy5Tnz59tGDBghptp6CgQD169ND06dNrtF52drZGjRqlK664okbrnYJuKQAALCegNisZhnF8AwEBevPNN9WlSxfdfffdNdrOoEGDNGjQoBrv/6677tLIkSPlcDhq1NpzCsINAACWU6tw8/XXXysmJqbSvAkTJqh79+5asWJFvRR2JjNmzND27dv1f//3f3rmmWfqtjEXN88EAMBqahVuLrvsstPOHzBggAYMGOCZjoyMVHp6ulJSUmpX3Um2bNmiRx55RN98840CAqpXenFxsYqLiz3Tubm5x1+k5QYAAMtp0Iv4ndh9VVcul0sjR47UpEmTlJqaWu31nnvuOUVFRXkeSUlJJ2yUcAMAgNX4zBWK8/LytGrVKt1zzz0KCAhQQECAnnrqKa1du1YBAQFaunTpadebOHGicnJyPI9du3Ydf5GzpQAAsJxadUt5Q2RkpNatW1dp3j//+U8tXbpUH374odq1a3fa9ZxOp5xO5+k3ynVuAACwHK+Gm/z8fG3dutUzvWPHDqWnpysmJkZt2rTRxIkTtWfPHr3zzjuy2+3q2rVrpfXj4+MVHBx8yvxqo1sKAADLadBwU3E9nDNZtWqV+vfv75meMGGCJGn06NGaOXOm9u3bp8zMzIYrsIxwAwCA1diM+hz1e5KIiAitXbu23s6Wqg+5ubmKiopSziMRihzxTyntFm+XBAAAzsLz+Z2To8jIyCqXrVXLTUULy8lsNpuCg4N1zjnn6Le//a0+++wztWrVqja7aBylx7xdAQAAqGe1Cjdr1qzRTz/9JJfLpY4dO0qSNm/eLIfDoU6dOumf//ynJkyYoG+++ebMg3mbgrIib1cAAADqWa1OBb/22ms1YMAA7d27V6tXr9bq1au1e/duXXnllbrpppu0Z88eXXrppWds4WkyCDcAAFhOrcbctGrVSl9++aW6dOlSaf769et11VVXac+ePfrpp5901VVX6dChQ/VWbH2oNOZmwJ+kKx7zdkkAAOAsajLmplYtNzk5OcrKyjpl/sGDBz23N4iOjlZJSRM/G4mWGwAALKfW3VJjx47VvHnztHv3bu3evVvz5s3TbbfdpmHDhkmSfvjhhxrdJsErGFAMAIDl1GpA8WuvvaYHHnhAN954o8rKzDtrBwQEaPTo0fp//+//SZI6deqkN998s/4qbQi03AAAYDl1us5Nfn6+tm/fLklKSUlReHh4vRXWUCqNuen5e2n4DG+XBAAAzqLBr3NTITw8XN27d6/LJryLlhsAACzHZ+4K3iAYcwMAgOX4d7ih5QYAAMvx73BDyw0AAJbj3+GmrNjbFQAAgHrm5+GGlhsAAKzGr8ONUcqYGwAArMavww0tNwAAWI9/hxtabgAAsBy/Djc2V7HkKvN2GQAAoB75dbiRJJUWeLsCAABQj/w23JQaDvNJcb53CwEAAPXKb8NNgZzmkxLCDQAAVuLH4SbEfEK4AQDAUvw23BwzgswndEsBAGApfhtuaLkBAMCa/DbcFFa03JRwthQAAFbit+HG03JTnOfdQgAAQL3y43DD2VIAAFiR34abY0aw+YRuKQAALMVvw83xbilabgAAsBK/DTfHBxQz5gYAACvx23Bz/FRwuqUAALASPw435QOK6ZYCAMBS/Dbc5Buh5pOiHO8WAgAA6pXfhptcEW4AALAivw03OZ6Wm2yv1gEAAOqX34abvIqWm2PZXq0DAADUL78NN7kVLTdlx6SyYu8WAwAA6o3fhpt8hcotmzlB6w0AAJbht+HGkF35YtwNAABW47fhRpJyjDDzCS03AABYhl+Hm2zOmAIAwHL8Ntw47DZabgAAsCC/DTcRTodyVB5uaLkBAMAy/DbcRIYE0nIDAIAF+W+4CQ5UtiLMiWNHvFsMAACoN34bbiJCAnXYiDQnCg55txgAAFBv/DbcRAUHnBBuDnq3GAAAUG/8NtxEhgTqsGi5AQDAavw33AQH0nIDAIAF+W+4CQnQISPKnCg8JLnd3i0IAADUC/8NN8GBOlpxtpThlo4d9W5BAACgXvhtuIkKDVKpApRnCzdnFDLuBgAAK/DbcBMbHihJOqLyrinG3QAAYAl+G26ahzslSVnu8q4pwg0AAJbgv+EmzAw3B93lZ0zlE24AALACvw03Yc4AhQQ6dLDijKn8/d4tCAAA1Au/DTeSFBsRpANGjDmRu8+7xQAAgHrh3+Em3Kl9nnCzx7vFAACAeuH34Wa/ysNNHi03AABYgd+HmwNGM3Mid693iwEAAPXCr8NNXIRT+yu6pUrypaJc7xYEAADqzK/DTWJUsAoVrEJbmDmDrikAAHyeX4ebhOgQSdJBG4OKAQCwCr8ON62igyVJu90V4YZxNwAA+Dq/DjcJUWbLza9l5eEme5cXqwEAAPXBr8NNmDNA0aGB2mW0MGcc3enVegAAQN35dbiRpMSoEGUa8eYE4QYAAJ/n1XCzfPlyDR06VImJibLZbJo/f36Vy8+dO1dXXnml4uLiFBkZqT59+ujzzz+vUw2J0YQbAACsxKvhpqCgQD169ND06dOrtfzy5ct15ZVXatGiRVq9erX69++voUOHas2aNbWuoVV08PFwk79fKims9bYAAID3BXhz54MGDdKgQYOqvfyLL75YaXry5Mn6+OOP9cknnygtLa1WNbRuFqochanQHqZQd4GU/asU37lW2wIAAN7n02Nu3G638vLyFBMTc8ZliouLlZubW+lxouTYMEk27bW1NGfQNQUAgE/z6XAzdepU5efna8SIEWdc5rnnnlNUVJTnkZSUVOn1drHm1Ym3lcWaMwg3AAD4NJ8NN//61780adIkvf/++4qPjz/jchMnTlROTo7nsWtX5WvZtIkJld0m7XDFmTMINwAA+DSvjrmprTlz5uj222/XBx98oAEDBlS5rNPplNPpPOPrQQF2JcWEKjO7/Fo3h7fVZ6kAAKCR+VzLzXvvvac//OEPeu+99zRkyJB62WZy8zBtcbcyJw5m1Ms2AQCAd3g13OTn5ys9PV3p6emSpB07dig9PV2ZmZmSzC6lUaNGeZb/17/+pVGjRumFF17QBRdcoP3792v//v3KycmpUx3tYsO0xSgPNzmZUnF+nbYHAAC8x6vhZtWqVUpLS/Ocxj1hwgSlpaXp8ccflyTt27fPE3Qk6fXXX1dZWZnGjRunhIQEz2P8+PF1qqN9fLiyFaEcezNzxiFabwAA8FVeHXPTr18/GYZxxtdnzpxZaXrZsmUNUkeXhAhJ0hajlXrpqNk11er8BtkXAABoWD435qYhdGwZKZtNWl+aYM7I2ujdggAAQK0RbiSFOwPUNiZUW4zW5gwGFQMA4LMIN+U6J0Rqi7si3NByAwCAryLclOuSEKnNFWdMZXPGFAAAvopwU65zQqSOKlJHbOVnTDHuBgAAn0S4Kdc5MVKStNbV1pyxL917xQAAgFoj3JRLjApW87AgrXMnmzP2rfVqPQAAoHYIN+VsNpt6tm2m9e525gxabgAA8EmEmxOc37aZfqloucnaKJUWebUeAABQc4SbE/Rs00x7FKsjipTcZXRNAQDggwg3J+jeOkoBdrt+cHU0Z2R+592CAABAjRFuThAc6FCPpGj96E41Z2T+17sFAQCAGiPcnOTiDrH60d3JnMj8r+R2e7cgAABQI4Sbk1xyTqw2GG11TE6pKFs6uMnbJQEAgBog3JykR1K0gp3B+snVwZyRudK7BQEAgBoh3Jwk0GHXhSkx+tGoGFTMuBsAAHwJ4eY0zHE3FeGGlhsAAHwJ4eY0rujcQmvc56jMsEs5u6TsXd4uCQAAVBPh5jSSYkLVoXUL/WIkmzN2LPdqPQAAoPoIN2cwuFuClru7mxPblni3GAAAUG2EmzMY1LWl/uPqIUlyb10quV1erggAAFQH4eYM2jYPU0nLnso1QmUvOirtXePtkgAAQDUQbqpwdffWWu7uZk5kfObdYgAAQLUQbqowpFuCvnSdL0kq27DQy9UAAIDqINxUITk2TNmt+qnMsCvg8CbpyA5vlwQAAM6CcHMW1/Xtqu/dnSVJpes/9nI1AADgbAg3ZzG4W4K+c14sScpd9b6XqwEAAGdDuDmLQIddCX1GyGXY1DxnvdyHtnu7JAAAUAXCTTUMu/g8/WDrKknatmyWl6sBAABVIdxUQ7gzQIdTrpMkRWz8t+R2e7kiAABwJoSbavrNkNHKN0LU0rVPm39c7O1yAADAGRBuqqlF8+ZaFzNAknT027e8XA0AADgTwk0NtOp/hySpR+5/lLlnr5erAQAAp0O4qYE23S7V7sBkBdtK9eMnr3u7HAAAcBqEm5qw2WTrOUqSdN7eOfp51xEvFwQAAE5GuKmhVv3v0DF7mNrb92nRhzNlGIa3SwIAACcg3NRUcKRcPcdIki4/OkcfpzP2BgCApoRwUwvhl94rly1Ave0Z+vDjudqXc8zbJQEAgHKEm9qITJCt+w2SpFtcH+vB99fK7aZ7CgCApoBwU0v2vvdJkq6yr9KB7T/r7RU7vFwRAACQCDe1F99J6jhYdpuhhwPm6K+LM7Rpf663qwIAwO8Rburiiidk2By6yrFavYyfdf+cdBWVurxdFQAAfo1wUxfxnWT7zW2SpKeC/k9b9mfrhS8yvFwUAAD+jXBTV/0mSsHR6qBM3eRYqje+2aGlmw54uyoAAPwW4aauQmOk/n+WJE0MnqtI5eu+99K1cR/jbwAA8AbCTX3oNVaK7agwV46mNP9M+cVl+sOMH7n+DQAAXkC4qQ+OQOnqyZKkQcc+Ub+Yo9qfW6RRb/2gowUlXi4OAAD/QripLx0GSOcMlM1dpteavavEiEBtycrXrW9/r5zCUm9XBwCA3yDc1KdBU6SgcDn3/Fef9PhOzcOC9MueXAIOAACNiHBTn2JSpGtelCQ1X/Wi5g12KSYsSD/vztENr69UVl6Rd+sDAMAPEG7qW/fhUtotkgy1+Xq83r+lg2LDndq0P0/Xv7JSOw4VeLtCAAAsjXDTEAb9VYrtKOXvV4cVD+mjuy5Qm5hQZR4p1O/+uUL/3X7Y2xUCAGBZhJuGEBQmDZ8pBQRLW79S201v6aP/uUg9kqKVXViqW978XrO+2ynD4E7iAADUN8JNQ2nRRbp6ivl86dOKy16rf995oYb2SFSZ29ATC9brwffXqqC4zLt1AgBgMYSbhnT+GOnc30nuMunftyi4YI9euvE8/WVIZznsNs1ds0dD//Gt1u/N8XalAABYBuGmIdls0tCXpPhzpfwD0uzhshVl6/ZLUjT79gvUItKp7QcLNGz6Ck3/eqvKXG5vVwwAgM+zGX428CM3N1dRUVHKyclRZGRk4+w0Z7f05gApb5+UdIF0y1zJGa4jBSV65KOf9cUG80ab3VtHacp13dUlsZHqAgCgJgxDKsqRCg9LhUfKfx6WCg9JXz5uLtPnHmngs/W+65p8fhNuGsuB9dKMQeY/irZ9pZHvS85wGYahj37ao0mfrFdeUZkcdpvG9k3W+AGpCncGNF59AAD/U1Io7f9Zkk3a9V8pLO6E8HL4NCHmsDnU4myerP/hFoSbKngt3EjS7tXSu8Ok4lypTR/p5g8kZ4QkKSu3SE8sWK/PftnvWXz8FefovivOkcNua9w6AQC+ye2Wjh2VCg6WP7KkgkMnTB+SNi2s+36CwqXQGCm0ufkIiZHWvW++dtMcqeOguu/jJISbKng13EjlAed3UnGOlHShNPLfUki05+WvM7I0acF67TxcKEnq1DJC/3t1R/XvGC+bjZADAH6npPCEcHKwclApOCjlnxBgCg9Lhqv2+3JGSR0ul0JjjweXE0NMWKwZZAKD6+/9VRPhpgpeDzeStOcnswWnKEdq0dXsoopq5Xm5uMyl6V9v00tLtnjm9WrbTOMHnKOLO8QScgDAl7ldZgg5OaRUPPJPCjCltbiyfUgzs4spLM4MJGHxx5+HNJOWPiOdP1rqMdKc5wOfK4SbKjSJcCNJ+9dJ//d78yyqiASzBSehR6VFsgtL9MqybZr53U4Vl5lnUvVIitYfL03RwHNb0l0FAE2BYUgl+dULKhWtK6rhR6/DKYXHnxBYKkLLCc8rXg9tLjkCG+StehPhpgpNJtxI0tFfpX+NkA5ukgLDpN+9InW59pTFDuQW6dX/bNN7P2SqqNQMOW1iQjW2b7J+f35rRQRb7x8xAHiVq/TU1pX8rNMEmPLnZcdquANbeTfPySElTgo/TYAJCveJ1pWGRLipQpMKN5J0LFv6YIy0/Wtz+sJx0oAnpYCgUxY9lF+sWd/t1Lv//VXZhaWSpNAgh4altdLI3m3UtVVUo5UNAD7FMMyTOSqNUzldUCmff+xozfcRGHbm1pSTQ0xIjOTgjNiaINxUocmFG0lylUlfPSGt/Ic5nXCedP3bUvP2p128sKRMH63erVkrf9XWrHzP/C4Jkfr9+a117XmJig13NkLhAOBFrjLp8Fbzy+HiRyq/1mPkqWcLuUpqtn2b3RxYe6bWlJO7iILC6u+94RSEmyo0yXBTYdMi6eO7zW8MQeHS4KlSjxvP2BRpGIZWbj+sf32fqS/WH1BJ+RWOHXabLjknVkO6Jeiqc1sqKoRuKwA+wu0y/wYe3Sn9ukL67mUzmNSXoIgqgsqJA2/jzIG3di7k31T4TLhZvny5/va3v2n16tXat2+f5s2bp2HDhlW5zrJlyzRhwgStX79eSUlJ+stf/qIxY8ZUe59NOtxIUs4eae4d5n9qSeo4RBryghSZUOVq2YUl+jh9r+b+tFtrdx+/eFKA3aaebZupX8c49UuNV+eECM62AtD4ivOkn96RAoKlgxmS4ZZ+fKPh9jfgycpnCFX8DAxpuH2iQflMuPnss8+0YsUKnX/++bruuuvOGm527Nihrl276q677tLtt9+uJUuW6P7779enn36qgQMHVmufTT7cSOY3l2+nScuel9zm2Bp1/q00fKZkd5x19W0H87Vw7T59um6vNh/Ir/Ray8hgXZYap34d49T3nFhFMhgZQG25SqUtX0p7VkvfTDW/jGV8eupygWG1O535TDoOls4baV5KI7pNtf4uwvf5TLg5kc1mO2u4efjhh/Xpp5/ql19+8cy78cYblZ2drcWLF1drPz4RbiocWC+9ctHx6ZbdpKufl5L7VnsTmYcLtWxzlpZlHNR32w55zraSKrfqXHpOnDonRHJ6OQDzTM4NH0vblh4/2aFNHylzZcPut/cfpda/kc79HYNtcQrLhptLL71UPXv21IsvvuiZN2PGDN1///3KyTn9fSyKi4tVXFzsmc7NzVVSUpJvhBvJHDD3yXhp/bzj33w6DjabXOM61mhTRaUufb/jiJZlZOk/GQe1/VDlb1LhzgCltYnWb5Jj1LNNM/VsG63QIP7AAJbgdkm7vpd+fl+KaCkte65x9nvDbKlFF7OLyBneOPv0I4ZhqNRlyOU2VOZ2y+WuPF3mMlTmNrT5QJ4++2W/hnZP0Nyf9ujzDft1pk9/u01yn/BaRHCA2sSE6qbebeQ2zG273IYMQ3KVT+/NPqbZ32dKksb1b68/DexU7+/VsuEmNTVVf/jDHzRx4kTPvEWLFmnIkCEqLCxUSMipfalPPvmkJk2adMp8nwk3FQoOSV8/K62eZV5a22aXut8oXfrQGc+qOpuKVp2vN2Xpx51HlV986s3Q2seFqUfraHVtFaXUFhHqEB+uFpFOxu0ATUXpMWn3j2b3UEyKtPD+htlPbEfpUIbUYYA5biYyUUodKKVc7pODbg3DUInLrZIy85FbVKb1e3OUX2T+Hfx03T7FhAXp4/S9ah4WpGOlLhWWuNS/Y5wKil36YecR3XfFOcouLNE7K3+tdR0XtW+uMrehUpdbazKzPfPbx4WVBxTjhMDi9kybocVdKYQ0JTunDKn3bRJuTuDzLTcnO7hZWjLp+I3PbHap6/XSxQ+Y345qyeU2lLE/T6t+PaIfdx7Vd1sP6XDB6U+bjHAGqH18uM6JD1eH+HCd0yJcHeIi1LpZiOx0awF1Zxjm1ct3fCOtnnH8BIOG1Po35rVfzh0m9RorRbett4vGVbQulLjcKi1ze0JFcXmwODFklLhcZ37t5HVP81qp64R1T9l25edWFmC3yWG3KcBuk81mO+2X15q6vFO8nAF22e02OWzm9u02mxx2aWtWvn4qD2d/v/E8XXteq6o3Vgs1CTc+1efQsmVLHThwoNK8AwcOKDIy8rTBRpKcTqecTgtd8yUuVbpxtnkDzuV/lTYvNu/Euu59KaWf1Os2826sNbz0tsNuU5fESHVJjNSoPsmSzCsjr92VrQ37crV+b662ZeXr1yOFyisuU/qubKXvyq60DWeAXe3jygPPCcGnbfMwBTp875sdUO9Ki6S8vVLm9+YVbf/zN3O6vkS3lbJPbUXIPu8uHWg7VEuOxikoMFAd4sMVERyg379ywhiareU/90tasl7Ses9LYy5KVpn7DAGjqgBR5lZxeeBoGl+jm5YJV6aqbfNQBdjtctilLQfyNT99j8b0badz4sMVYLcpwGGvFFQcdpsCHfZK06dbxt9b132q5ebhhx/WokWLtG7dOs+8kSNH6siRI9YcUFwde9dI30wzW3KM8m8i4S2lnrdKPUdL0Un1urviMpd2HirU1qx8bcnK09asfG3Nytf2gwVVfhNqGRms6NBApcSFqUN8hFxut+LCnerWOkopseGKDg30+/+M8FGGYV6mf/86adXb5o1xc3c36C7XuZO1y4hXiQI0tewG7TZiJfnO/x+H3aYgh11BAeUPh13OgMrTp3vuLJ8OPPH1067vOGV95+m2XTHtsNPq7AN8plsqPz9fW7eaXxfS0tI0bdo09e/fXzExMWrTpo0mTpyoPXv26J133pF0/FTwcePGaezYsVq6dKnuu+8+650KXhtHf5V+mmVeR6Liglc2u3TOVdL5Y6T2V5z2lg71pczl1q6jxzxhZ0tWnrZl5WtLVr4KS1w12taQ7uY1fVpGBmvVziNqGRWsOy9tr+jQQEUEByg2zMkfIlSLYRgqLnOroLhMNptNGfvzlF1YorbNw7Ry+2FN/3qrWkQGa+O+XDULDdTRwlK1iQlVTFiQ1u46ohTbPl1t/1G3BSzSBndbtYoMULuCtQ1W73Z3S/3o7qQ5rv7aZCTpmIIbbF9BDnu1u2b+eFmKggMcZw4J1XjuPCFwcFYmasNnws2yZcvUv3//U+aPHj1aM2fO1JgxY7Rz504tW7as0joPPPCANmzYoNatW+uxxx6z1kX86qqsxLzOxKq3pR3Lj88PjpbaXWremLPLtY12x1jDMHQwr1hbD+Zrb3aRth/M19HCUmXsz/X0z9aHId0T9NWGA567p3dJiFRidLBaNwtVz7bNFGC3qWVUsJqHBSncGaDw4AAFOey0FtWzijBR4nKruNQMFdsP5SsiOFCLf9mvw/nFGtClhWJCg1Rc3n1RekI3RqnLrZ8yszVvzR7PNn+X1qrS9JmkxIapuMyt4jKXDuWf6TL7hiJVoDhbjuJsOYpVjuJs2epk26UbApbVz0E4jQx3a60zUvSvssu1xuggQ5W7aTu1jFBchFPfbDnkmXdxh1h9u/XQyZuqJCokUDnHStUlIVIJUcH67XmJSogK0dJNWbqic7zax4UrwGFTaKBDAXQNw8f5TLjxBsuHmxMd2mqGnF8+kvL3H58fHG2eTp56lZTSXwqJ9laFksx7ZR3MK9aazGyt35uj5uFObd6fp6OFJfo6w2yFqvhW3ZgigwP0m+QYBQWYzeCrfz2qPdnH1L9jnFJbRph92zab7OX93J6ftuP93l9sOKCdhwv0P5d1UKDDpqy8Yv3t8wwN7tZSbrc0sGsLrdudq9bNQpRbVKr4iGD9d/thLVi7V6P7tFXPts20N7tIyzKyFBEcqEtTY7X4l/36cecRPTH0XP1lvnnNpwtTYvT9jiMyDGlItwR1bRUlt2FoycYDZwyR0aGBGtQ1QWUut+dsjYU/72vEI1x/IlSo6x3/UbwtW5lGvALkUpwtW3HKUazNDDCxtlzFKVtOW+0HVu6POFct89braFh7HY3uqpKwBO1uNUgzNkpXdWurzKPHFOYM0EXtmysk0KHm4UEKDnQoKiTQM7ATQO0QbqrgV+Gmgtsl7fxWWj9XWj2z8ms2h3mWRGwHKW2U1Or8Jn3xrIpv5Vm5Rdq4L0+frturizvEqajUpb8v2aLEqGDtzSnydpk4i9QW4Z7QeOKYiECHXf/dflhHC0sVqQK1tR3QbbEbNCzvX6dsY7cRq9a2qls2zsgZZd6tueIRECKtLd9Hi27SeTeZA/Oj2jTp/w+APyHcVMEvw82J3C7p1++kjM+kLV9Ih7ecukybi8xr53QfYYYdi93pttTlVnZhqY4UlCi/uExFpS7tzT6mH3YcUeeESNls0tpd2erdrrkcdpWfAWJoTeZRLfx5n265sI2CAxyei1dVepw077NfzBazS86JlcNu07KMyjcAjA4NVHZ5i1RSTIhiw52VrnVxsuZhQZ5T9Lu2itQve3JPu9zve7aWYRiaW0V3TlyEU9f1bKWokEAF2s1xEHuzj+nNb3fo4g6x6t46Sks3ZWnT/jzPOi/dlKaI4AB9ueGAUuPD1T0pWiVlbq3cdlipLSLUKSFCzUKDFBrkqHQaqnngi6StX0nbl0m5e80r35YWnu3XVTudf1seXFqU39G5xfEgExYvBTbcWBYADYNwUwW/DzcnO/qrtG2JlLFY2vL5qa87gqTYVPN6F+Etpcv/IsV3lgIsdHo9aqf0mHQs2zxjb85Np18m6QJzgHv+Qakk7/TL1NaId8x7GxVlS11/b3a30u0DWBbhpgqEmyq4XdLedLNFZ+c35vPT3ezOHijFnmPetC6uoxl+4jtLzZIbbaAy6llxnpS7T8rOlD7/szmv9JgUHmfeFLExdBkmdfmtlLNHat3LvJeaM6Jx9g2gySPcVIFwUwOGIeXsMsfrLH1Gyt1jfjsuyj798vYAswugxblSTHsz7DRrK0W2Mu/cG9KMb9YNyTDMK8xmZ5rHuSRfOnbUnPfZ/zZ+Pb+5XWp3mflvIizODErOSP4NAKgVwk0VCDd1VBF4DmyQDvwiZW2U9qVL2bskV3HV6wYEmzfsi0gwH5GJZreCzSYlX2x+AAZHmUEoMNT6H4Jut9kyVlJoBpHSQqmkwHxeUv68tKB8XvkyGYukI9sbt864ztLBjdK515k1xneRul5ndlOGxVr/9wSgSSDcVIFw00DcbjP0ZP9qfvge3iYd3Vk+L9O8gmtNhbeUgiPNQFUx8Dn+XPN6PXaHOe4nIPiEn8EnTZf/PLpDat5BCgwxzw4zXGbXWlGOFBRqtm44gsyAFpNi1lqUbXbLtOhqLu8qldyl5k9XqeQqMcPG1q+koHCpZffj4aS0PIiUnC64nPAoO1afv4HKnFFmeAyNMcfFZJVfSj+hh7RvrdmCcskEszuqVS9zfkSCT94AEYB/INxUgXDjJWXFUt5+KW+feaZMxc9tX5sfvBWhwy/ZzIAUFGqemRYYZv48ZTrM7Br85SMzTGVtlC6+X4rrZHYFBjilqCTGPQGwJMveOBM+LMBpjr9p1rbq5QzjeLdM7h6zxSM/S9q4wGxhaf0bqTjf/AAvK5bKiqr+WVpottwERZhhIf+A2fVVlGPeniI4ytxuVZolmy09jkBzXJEj0LwSdEi0OfBaMm9zERp7PIRUPAJDzxBcTngEBNesa+f6t6u/LAD4IcINmhabzTxDxhkhRbQ4Pr/rdd6rCQDgU+hgBwAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlhLg7QIam2EYkqTc3FwvVwIAAKqr4nO74nO8Kn4XbvLy8iRJSUlJXq4EAADUVF5enqKioqpcxmZUJwJZiNvt1t69exURESGbzebtcvxGbm6ukpKStGvXLkVGRnq7HL/Csfcejr13cNy9pyGPvWEYysvLU2Jiouz2qkfV+F3Ljd1uV+vWrb1dht+KjIzkj42XcOy9h2PvHRx372moY3+2FpsKDCgGAACWQrgBAACWQrhBo3A6nXriiSfkdDq9XYrf4dh7D8feOzju3tNUjr3fDSgGAADWRssNAACwFMINAACwFMINAACwFMIN6s306dOVnJys4OBgXXDBBfrhhx/OuOwbb7yhSy65RM2aNVOzZs00YMCAKpdH1Wpy7E80Z84c2Ww2DRs2rGELtLCaHvvs7GyNGzdOCQkJcjqdSk1N1aJFixqpWuuo6XF/8cUX1bFjR4WEhCgpKUkPPPCAioqKGqlaa1i+fLmGDh2qxMRE2Ww2zZ8//6zrLFu2TD179pTT6VSHDh00c+bMBq9TkmQA9WDOnDlGUFCQ8fbbbxvr16837rjjDiM6Oto4cODAaZcfOXKkMX36dGPNmjXGxo0bjTFjxhhRUVHG7t27G7ly31fTY19hx44dRqtWrYxLLrnEuPbaaxunWIup6bEvLi42evXqZQwePNj49ttvjR07dhjLli0z0tPTG7ly31bT4z579mzD6XQas2fPNnbs2GF8/vnnRkJCgvHAAw80cuW+bdGiRcajjz5qzJ0715BkzJs3r8rlt2/fboSGhhoTJkwwNmzYYLz88suGw+EwFi9e3OC1Em5QL3r37m2MGzfOM+1yuYzExETjueeeq9b6ZWVlRkREhDFr1qyGKtGyanPsy8rKjIsuush48803jdGjRxNuaqmmx/6VV14xUlJSjJKSksYq0ZJqetzHjRtnXH755ZXmTZgwwejbt2+D1mll1Qk3//u//2uce+65lebdcMMNxsCBAxuwMhPdUqizkpISrV69WgMGDPDMs9vtGjBggFauXFmtbRQWFqq0tFQxMTENVaYl1fbYP/XUU4qPj9dtt93WGGVaUm2O/YIFC9SnTx+NGzdOLVq0UNeuXTV58mS5XK7GKtvn1ea4X3TRRVq9erWn62r79u1atGiRBg8e3Cg1+6uVK1dW+j1J0sCBA6v9uVAXfndvKdS/Q4cOyeVyqUWLFpXmt2jRQps2barWNh5++GElJiae8h8BVavNsf/222/11ltvKT09vREqtK7aHPvt27dr6dKluvnmm7Vo0SJt3bpVd999t0pLS/XEE080Rtk+rzbHfeTIkTp06JAuvvhiGYahsrIy3XXXXfrzn//cGCX7rf3795/295Sbm6tjx44pJCSkwfZNyw28bsqUKZozZ47mzZun4OBgb5djaXl5ebr11lv1xhtvKDY21tvl+B232634+Hi9/vrrOv/883XDDTfo0Ucf1auvvurt0ixt2bJlmjx5sv75z3/qp59+0ty5c/Xpp5/q6aef9nZpaCC03KDOYmNj5XA4dODAgUrzDxw4oJYtW1a57tSpUzVlyhR99dVX6t69e0OWaUk1Pfbbtm3Tzp07NXToUM88t9stSQoICFBGRobat2/fsEVbRG3+3SckJCgwMFAOh8Mzr3Pnztq/f79KSkoUFBTUoDVbQW2O+2OPPaZbb71Vt99+uySpW7duKigo0J133qlHH31Udjvf8xtCy5YtT/t7ioyMbNBWG4mWG9SDoKAgnX/++VqyZIlnntvt1pIlS9SnT58zrvfXv/5VTz/9tBYvXqxevXo1RqmWU9Nj36lTJ61bt07p6emex29/+1v1799f6enpSkpKaszyfVpt/t337dtXW7du9QRKSdq8ebMSEhIINtVUm+NeWFh4SoCpCJgGdyBqMH369Kn0e5KkL7/8ssrPhXrT4EOW4RfmzJljOJ1OY+bMmcaGDRuMO++804iOjjb2799vGIZh3HrrrcYjjzziWX7KlClGUFCQ8eGHHxr79u3zPPLy8rz1FnxWTY/9yThbqvZqeuwzMzONiIgI45577jEyMjKMhQsXGvHx8cYzzzzjrbfgk2p63J944gkjIiLCeO+994zt27cbX3zxhdG+fXtjxIgR3noLPikvL89Ys2aNsWbNGkOSMW3aNGPNmjXGr7/+ahiGYTzyyCPGrbfe6lm+4lTwP/3pT8bGjRuN6dOncyo4fM/LL79stGnTxggKCjJ69+5t/Pe///W8dtlllxmjR4/2TLdt29aQdMrjiSeeaPzCLaAmx/5khJu6qemx/+6774wLLrjAcDqdRkpKivHss88aZWVljVy176vJcS8tLTWefPJJo3379kZwcLCRlJRk3H333cbRo0cbv3Af9vXXX5/273bFsR49erRx2WWXnbLOeeedZwQFBRkpKSnGjBkzGqVW7goOAAAshTE3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3ABpVv379dP/993u7DAANYPny5Ro6dKgSExNls9k0f/78Gm/DMAxNnTpVqampcjqdatWqlZ599tkabYO7ggNoVHPnzlVgYKC3ywDQAAoKCtSjRw+NHTtW1113Xa22MX78eH3xxReaOnWqunXrpiNHjujIkSM12ga3XwAAAPXOZrNp3rx5GjZsmGdecXGxHn30Ub333nvKzs5W165d9fzzz6tfv36SpI0bN6p79+765Zdf1LFjx1rvm24pAI3qxG6p5ORkTZ48WWPHjlVERITatGmj119/vdLyu3fv1k033aSYmBiFhYWpV69e+v777z2vv/LKK2rfvr2CgoLUsWNHvfvuu5XWt9lseu2113TNNdcoNDRUnTt31sqVK7V161b169dPYWFhuuiii7Rt27ZK63388cfq2bOngoODlZKSokmTJqmsrKxhDgrgJ+655x6tXLlSc+bM0c8//6zhw4fr6quv1pYtWyRJn3zyiVJSUrRw4UK1a9dOycnJuv3222vccsNdwQE0qssuu8wYP368YRjm3eFjYmKM6dOnG1u2bDGee+45w263G5s2bTIMwzDy8vKMlJQU45JLLjG++eYbY8uWLca///1v47vvvjMMwzDmzp1rBAYGGtOnTzcyMjKMF154wXA4HMbSpUs9+5NktGrVyvj3v/9tZGRkGMOGDTOSk5ONyy+/3Fi8eLGxYcMG48ILLzSuvvpqzzrLly83IiMjjZkzZxrbtm0zvvjiCyM5Odl48sknG+9AAT5OkjFv3jzP9K+//mo4HA5jz549lZa74oorjIkTJxqGYRh//OMfDafTaVxwwQXG8uXLPXcV79+/f832XefqAaAGTg43t9xyi+c1t9ttxMfHG6+88ophGIbx2muvGREREcbhw4dPu62LLrrIuOOOOyrNGz58uDF48GDPtCTjL3/5i2d65cqVhiTjrbfe8sx77733jODgYM/0FVdcYUyePLnSdt99910jISGhhu8W8F8nh5uFCxcakoywsLBKj4CAAGPEiBGGYRjGHXfcYUgyMjIyPOutXr3akOT50lMdDCgG4FXdu3f3PLfZbGrZsqWysrIkSenp6UpLS1NMTMxp1924caPuvPPOSvP69u2rv//972fcR4sWLSRJ3bp1qzSvqKhIubm5ioyM1Nq1a7VixYpKZ2i4XC4VFRWpsLBQoaGhtXy3gP/Kz8+Xw+HQ6tWr5XA4Kr0WHh4uSUpISFBAQIBSU1M9r3Xu3FmSlJmZWe1xOIQbAF518plTNptNbrdbkhQSElLv+7DZbGecV7Hf/Px8TZo06bRnewQHB9dLTYC/SUtLk8vlUlZWli655JLTLtO3b1+VlZVp27Ztat++vSRp8+bNkqS2bdtWe1+EGwBNVvfu3fXmm2/qyJEjp2296dy5s1asWKHRo0d75q1YsUJdunSp03579uypjIwMdejQoU7bAfxNfn6+tm7d6pnesWOH0tPTFRMTo9TUVN18880aNWqUXnjhBaWlpengwYNasmSJunfvriFDhmjAgAHq2bOnxo4dqxdffFFut1vjxo3TlVdeWak152wINwCarJtuukmTJ0/WsGHD9NxzzykhIUFr1qxRYmKi+vTpoz/96U8aMWKE0tLSNGDAAH3yySeaO3euvvrqqzrt9/HHH9c111yjNm3a6Prrr5fdbtfatWv1yy+/6JlnnqmndwdYz6pVq9S/f3/P9IQJEyRJo0eP1syZMzVjxgw988wzevDBB7Vnzx7Fxsbqwgsv1DXXXCNJstvt+uSTT3Tvvffq0ksvVVhYmAYNGqQXXnihRnUQbgA0WUFBQfriiy/04IMPavDgwSorK1OXLl00ffp0SdKwYcP097//XVOnTtX48ePVrl07zZgxw3PNjNoaOHCgFi5cqKeeekrPP/+8AgMD1alTJ91+++318K4A6+rXr5+MKi6fFxgYqEmTJmnSpElnXCYxMVEfffRRnergIn4AAMBSuIgfAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwlP8Pg+vQq49dlR8AAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# plot g_z across z\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(df_ln[\"z\"], df_ln[\"g_z\"])\n",
+    "ax.plot(df_Pln[\"z\"], df_Pln[\"g_z\"])\n",
+    "\n",
+    "ax.set_xlim(left=1000)\n",
+    "ax.set_xlabel(\"income\")\n",
+    "ax.set_ylabel(\"g_z\")\n",
+    "ax.legend([\"log normal\", \"Pln\"])"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Model.md b/Model.md
new file mode 100644
index 0000000..c2c8936
--- /dev/null
+++ b/Model.md
@@ -0,0 +1,181 @@
+# Model
+
+Based Saez (2001), Diamond (1998), Saez, Slemrod, Giertz (2012), lecture notes by Rishabh Kirpalani.
+
+We derive the Diamond-Saez-Mirrlees optimal tax formula, and invert it following Lockwood and Weinzierl (2016) and Jacobs, Jongen and Zoutman (2017)
+
+## Environment
+
+Suppose households have an ability type, $\theta \in \Theta$ with distribution $f(\theta)$, which is private information. The production technology is such that type $\theta$ has production (or taxable income) $z(\theta) = \theta * l(\theta)$. The planner or policy maker can observe income / production $z$ for each type, but not labor $l$. The government has exogenous expenditures $E$, which must be funded by labor of households or (in the decentralized problem) taxes on labor.
+
+ An allocation is a set of consumption and output $\{ c(\theta), z(\theta) \}_{\forall \theta \in \Theta}$. The household has preferences given by 
+
+ $$ U(\theta) = u(c(\theta)) - v(l(\theta)) = u(c(\theta)) - v\left (\frac{z(\theta)}{\theta}\right ),$$
+
+where $u', v' > 0$, and $u'' \leq 0 \leq v''$. The marginal rate of substitution for the household between $c$ and $z$ is 
+
+ $$MRS_{c, z} = -\theta\frac{u'(c)}{v'(\frac{z}{\theta})}.$$
+
+
+Now, note that 
+
+$$\frac{\partial}{\partial \theta} MRS_{c, z} = -u'(c) \left[\frac{1}{v'(\frac{z}{\theta})}+ \frac{zv''(\frac{z}{\theta})}{\theta v'(\frac{z}{\theta})}\right] < 0 $$
+
+by assumptions on $u$ and $v$. Thus, the single crossing property holds. 
+## Incentive Compatibility
+
+Since labor effort and productivity type are unobservable, the planner or policy maker cannot achieve the first best outcome in which policy can be a function of ability directly. By the revelation principle, a decentralized equilibrium in which agents truthfully reveal their type is equivalent to planner's problem with incentive compatibility constraints. Thus, policy makers must design the tax system to ensure that agents truthfully reveal their productivity through their income choice. An allocation is incentive compatible (globally) if
+
+$$ u(c(\theta)) - v\left (\frac{z(\theta)}{\theta}\right ) > c(\hat \theta) - v\left (\frac{z(\hat \theta)}{\theta}\right ), \forall \theta, \hat{\theta}.$$
+
+However, it can be shown that this global incentive compatibility constraint can be replaced by the following conditions for local incentive compatibility:
+
+1. $U'(\theta) = \frac{z(\theta)}{\theta^2}v'\left (\frac{z(\theta)}{\theta}\right )$
+2. $z(\theta)$ increasing
+
+Mechanism design literature drops condition 2, to be verified ex-post.
+
+## Constrained Social Planner's Problem
+
+The constrained social planner wants to maximize welfare subject to incentive compatibility and resource constraints. The constrained social planner's problem is:
+
+$$\begin{align*}
+\max_{c(\theta),z(\theta)} & \int_{\Theta} G(U(\theta)) dF(\theta) \\
+\text{subject to} & \\
+& \int_{\Theta} (z(\theta) - c(\theta)) dF(\theta) \geq E \quad (\text{RC}) \\
+& U'(\theta) = \frac{z(\theta)}{\theta^2} v'\left(\frac{z(\theta)}{\theta}\right) \quad \forall \theta \in \Theta \quad (\text{LIC}) \\
+& U(\theta) = u(c(\theta)) - v\left(\frac{z(\theta)}{\theta}\right) \quad \forall \theta \in \Theta \\
+& z(\theta) \text{  increasing in } \theta
+\end{align*}$$
+
+Where G is a weighting function for the social planner which determines their redistributive preferences. We will drop the monotonicity condition to be verified later. The Lagrangian for the social planner's problem is:
+
+$$
+\begin{align*}
+L =& \int_{\Theta} G(U(\theta)) + \lambda[z(\theta) - c(\theta) - E] dF(\theta) \\
+& + \int_{\Theta} \gamma(\theta) \left [u(c(\theta)) - v\left(\frac{y(\theta)}{\theta}\right) - U(\theta)\right ] + \mu(\theta)\left [U'(\theta)- \frac{z(\theta)}{\theta^2} v'\left(\frac{z(\theta)}{\theta}\right)\right]d\theta\\
+\end{align*}
+$$
+
+
+Where $\lambda, \gamma(\theta), \mu(\theta)$ are Lagrange multipliers. There is an equivalent social planner's problem in which the planner has an exogenous value of raising public funds; another interpretation of $\lambda$ is the marginal value of public funds. 
+
+
+Using integration by parts on $\int_{\Theta} \mu(\theta)U'(\theta)d\theta$, we have:
+
+$$
+\begin{align*}
+L =& \int_{\Theta} G(U(\theta)) + \lambda[z(\theta) - c(\theta)- E] dF(\theta) \\
+& + \int_{\Theta} \gamma(\theta) \left [u(c(\theta)) - v\left(\frac{y(\theta)}{\theta}\right) - U(\theta)\right ] - \mu(\theta)\left [\frac{z(\theta)}{\theta^2} v'\left(\frac{z(\theta)}{\theta}\right)\right]d\theta\\ \\
+&+ \mu(\bar{\theta}) U(\bar{\theta}) - \mu(\underline{\theta}) U(\underline{\theta})
+\end{align*}
+$$
+
+It must be the case that $\mu(\underline\theta) = \mu(\bar\theta) = 0$. Otherwise, the planner would like to set $U(\barθ)  = ∞$ $(U (\underlineθ) = −∞)$ which would violate incentive constraints. 
+
+Taking first order conditions:
+
+$$
+\begin{align*}
+U(\theta): \quad & G'(U(\theta)) f(\theta) - \gamma(\theta) - \mu'(\theta) = 0 \\
+c(\theta): \quad & \gamma(\theta) u'(c(\theta)) - \lambda f(\theta) = 0 \\  
+z(\theta): \quad & -\gamma(\theta) \frac{1}{\theta} v'\left(\frac{z(\theta)}{\theta}\right) + \lambda f(\theta) - \mu(\theta) \left[\frac{1}{\theta^2} v'\left(\frac{z(\theta)}{\theta}\right) + \frac{z(\theta)}{\theta^3} v''\left(\frac{z(\theta)}{\theta}\right)\right] = 0
+\end{align*}
+$$
+
+Boundary conditions: $\mu(\bar{\theta}) = \mu(\underline{\theta}) = 0$.
+
+Using the first order conditions for $U(\theta)$ and $c(\theta)$, we have:
+
+$$
+\mu(\theta) = \int_\theta^{\bar\theta} \left[\frac{\lambda f(y)}{u'(c(y))} - G'(U(y)) f(y)\right] dy
+$$
+
+Substituting this into the first order condition for $z(\theta)$ yields:
+
+$$
+\frac{1}{\frac{1}{\theta} v'\left(\frac{z(\theta)}{\theta}\right)} - \frac{1}{u'(c(\theta))} = \frac{1-F(\theta)}{\theta f(\theta)} \left[1 + \frac{z(\theta)}{\theta} \frac{v''\left(\frac{z(\theta)}{\theta}\right)}{v'\left(\frac{z(\theta)}{\theta}\right)}\right] \int_\theta^{\bar{\theta}} \left[\frac{1}{u'(c(y))} - \frac{G'(U(y))}{\lambda}\right] \frac{dF(y)}{1-F(\theta)}
+$$
+
+This, togethether with the resource constraint characterizes optimal allocations. 
+
+## The Diamond-Saez-Mirrlees Optimal Tax Formula
+
+The policy maker wants to implement the constrained planner's problem with a tax policy $T(z)$. The household's problem is:
+
+$$\begin{align*}
+\max_{c, z} \quad &u(c) - v(\frac{z}{\theta}) \\
+\text{subject to} & \\
+& c = z - T(z)
+\end{align*}$$
+
+Assuming households do not optimize with respect to a highly nonlinear tax code, the first order condtion of the household is $\frac{1}{\theta} v'(\frac{z}{\theta}) = u'(c)(1-T'(z))$. 
+Letting $T(z)$ be the tax function that implements the efficient allocation and substituting the household FOC, the formula becomes:
+
+$$
+\frac{T'(\theta)}{1-T'(\theta)} = u'(c(\theta)) \frac{1-F(\theta)}{\theta f(\theta)} \left[1 + \frac{z(\theta)}{\theta} \frac{v''\left(\frac{z(\theta)}{\theta}\right)}{v'\left(\frac{z(\theta)}{\theta}\right)}\right] \int_\theta^{\bar{\theta}} \left[\frac{1}{u'(c(y))} - \frac{G'(U(y))}{\lambda}\right] \frac{dF(y)}{1-F(\theta)}
+$$
+
+As per Diamond (1998), we take preferences to be GHH, which has the interpretation of no income effects, so that the formula does not depend on both consumption and income simultaneously. This also has the benefit of being able to interpret utility in dollars.
+
+$$U(c, l) = c - \psi (\varepsilon l^{\frac{1}{\varepsilon}})$$
+
+where $\varepsilon$ is the elasticity of taxable income. Thus, $u'(c) = 1$, and $x \frac{v''(x)}{v'(x)} = \frac{1}{\epsilon}-1$. Also note that there is a 1 to 1 mapping between $\theta$ and $z$ under incentive compatibility. This gives us the intuition that we can replace $\theta$ with $z$. While this is an oversimplification, Saez (2001) shows that the formula with $z$ does indeed hold. However, one should be careful that we now are using $f$ as the virtual density, which makes the assumption that taxes are linearized around $T(z)$. This is fine as long as individuals are not optimizing with respect to a highly nonlinear tax code. Therefore, the formula becomes:
+
+
+$$
+\frac{T'(z)}{1-T'(z)} =  \frac{1-F(z)}{\varepsilon z f(z)}  \int_z^{\bar{z}} \left[1 - \frac{G'(U(y))}{\lambda}\right] \frac{dF(y)}{1-F(z)}
+$$
+
+
+Let $g(z) = G'(U(z))/\lambda$ be the marginal social welfare weight. The interpretation is that $g(z)$ is the social welfare or value to the policy maker from giving an additional dollar of income or consumption to an agent earning $z$. Therefore, we get the formula used in Lockwood and Weinzierl (2015) and Jacobs, Jongen and Zoutman (2017):
+
+$$
+\frac{T'(z)}{1-T'(z)} =  \frac{1-F(z)}{\varepsilon z f(z)}  \int_z^{\bar{z}} \left[1 - g(y)\right] \frac{dF(y)}{1-F(z)}
+$$
+
+
+We assume there is an unbounded distribution of income, so $\bar{z} = \infty$. The key components are:
+
+1. A hazard ratio $\frac{1-F(z)}{zf(z)}$. For a thin-tailed distribution such as the lognormal distribution, this converges to 0, giving us the "no distortion at the top" result of a 0 marginal tax rate at the top of the income distribution. For a Pareto distribution, this term is $\alpha$, which represents the thinness of the tail.
+2. The elasticity of taxable income $\varepsilon$. This has estimates ranging from .12 to .4, with .25 as a middle of the road estimate. See Saez, Slemrod, Giertz (2012).
+3. The planner's redistributionary motives, captured by $g(z)$. Note in general that these are endogenous; they depend on the equilibrium allocation. 
+
+
+## Inverting the optimal tax formula
+
+As per Lockwood and Weinzierl (2016), we can invert this formula:
+
+$$\bar{g}(z) = 1-F(z) - \frac{\varepsilon z f(z) T'(z)}{1-T'(z)}$$
+
+Where $\bar{g}(z) \equiv \int_z^\infty g(y)dF(y)$. By the fundamental theorem of calculus, $\frac{d}{dz}\bar{g}(z) = - g(z)f(z)$. Thus,
+
+$$
+\begin{align*}
+g(z) &= -\frac{1}{f(z)}\frac{d}{dz}\left[ 1-F(z) - \varepsilon z f(z)\frac{ T'(z)}{1-T'(z)}\right]\\
+&=1 + \frac{1}{f(z)}\frac{d}{dz}\left[\varepsilon z f(z)\frac{ T'(z)}{1-T'(z)}\right] \\
+& = 1 + \theta_z \varepsilon \frac{T'(z)}{1-T'(z)} + \varepsilon \frac{zT''(z)}{(1-T'(z))^2}
+\end{align*}
+$$
+
+Where $\theta_z \equiv 1 + \frac{zf'(z)}{f(z)}$ is the elasticity of the local tax base with respect to income. In the case of constant ETI, we get the formula used by Lockwood & Weinzierl (2016) on the first line, and Jacobs, Jongen and Zoutman (2017) on the 3rd line.
+
+In the case of variable ETI, we get an extra term for the JJZ formula:
+
+$$g(z) = 1 + \theta_z \varepsilon(z) \frac{T'(z)}{1-T'(z)} + \varepsilon(z) \frac{zT''(z)}{(1-T'(z))^2}+ \varepsilon'(z)\frac{zT'(z)}{1-T'(z)}$$
+
+However, note that this is subject to the constraint that 
+
+$$\int_0^\infty g(z) dF(z) = 1$$
+
+To see why this is true, consider a reform in which the government collects an additional dollar from everyone. Since we have GHH preferences and utility is in terms of dollars, the welfare loss from this reform is 
+
+$$ \int_0^\infty G(U(z)) - G(U(z)-1) dF(z) \approx \int_0^\infty G'(U(z))dF(z).$$
+
+The gain to the social planner of collecting this dollar is $\lambda$, the value (or shadow price) of relaxing the government budget constraint. At optimum, the marginal benefit equals the marginal cost, so
+
+$$\int_0^\infty G'(U(z))dF(z) = \lambda$$
+
+$$\implies \int_0^\infty g(z) dF(z) = 1$$
+
+by definition of $g$. 
diff --git a/iot/inverse_optimal_tax.py b/iot/inverse_optimal_tax.py
index 6e30f4d..7288897 100644
--- a/iot/inverse_optimal_tax.py
+++ b/iot/inverse_optimal_tax.py
@@ -40,9 +40,6 @@ def __init__(
         income_measure="e00200",
         weight_var="s006",
         eti=0.25,
-        bandwidth=1000,
-        lower_bound=0,
-        upper_bound=500000,
         dist_type="log_normal",
         kde_bw=None,
         mtr_smoother="kreg",
@@ -187,7 +184,7 @@ def compute_income_dist(
                 * f_prime (array_like): slope of the density function for
                     income bin z
         """
-        z_line = np.linspace(1, 1000000, 100000)
+        z_line = np.linspace(100, 1000000, 100000)
         # drop zero income observations
         data = data[data[income_measure] > 0]
         if dist_type == "log_normal":
@@ -234,10 +231,81 @@ def compute_income_dist(
                 # bw_method=kde_bw,
                 weights=data[weight_var],
             )
-            F = f_function.cdf(z_line)
             f = f_function.pdf(z_line)
-            f = f / np.sum(f)
+            F = np.cumsum(f)
             f_prime = np.gradient(f, edge_order=2)
+        elif dist_type == "Pln":
+
+            def pln_pdf(y, mu, sigma, alpha):
+                x1 = alpha * sigma - (np.log(y) - mu) / sigma
+                phi = st.norm.pdf((np.log(y) - mu) / sigma)
+                R = (1 - st.norm.cdf(x1)) / (st.norm.pdf(x1) + 1e-15)
+                # 1e-15 to avoid division by zero
+                pdf = alpha / y * phi * R
+                return pdf
+
+            def neg_weighted_log_likelihood(params, data, weights):
+                mu, sigma, alpha = params
+                likelihood = np.sum(
+                    weights * np.log(pln_pdf(data, mu, sigma, alpha) + 1e-15)
+                )
+                # 1e-15 to avoid log(0)
+                return -likelihood
+
+            def fit_pln(data, weights, initial_guess):
+                bounds = [(None, None), (0.01, None), (0.01, None)]
+                result = scipy.optimize.minimize(
+                    neg_weighted_log_likelihood,
+                    initial_guess,
+                    args=(data, weights),
+                    method="L-BFGS-B",
+                    bounds=bounds,
+                )
+                return result.x
+
+            mu_initial = (
+                np.log(data[income_measure]) * data[weight_var]
+            ).sum() / data[weight_var].sum()
+            sigmasq = (
+                (
+                    ((np.log(data[income_measure]) - mu_initial) ** 2)
+                    * data[weight_var]
+                ).values
+                / data[weight_var].sum()
+            ).sum()
+            sigma_initial = np.sqrt(sigmasq)
+            # Initial guess for m, sigma, alpha
+            initial_guess = np.array([mu_initial, sigma_initial, 1.5])
+            mu, sigma, alpha = fit_pln(
+                data[income_measure], data[weight_var], initial_guess
+            )
+
+            def pln_cdf(y, mu, sigma, alpha):
+                x1 = alpha * sigma - (np.log(y) - mu) / sigma
+                R = (1 - st.norm.cdf(x1)) / (st.norm.pdf(x1) + 1e-12)
+                CDF = (
+                    st.norm.cdf((np.log(y) - mu) / sigma)
+                    - st.norm.pdf((np.log(y) - mu) / sigma) * R
+                )
+                return CDF
+
+            def pln_dpdf(y, mu, sigma, alpha):
+                x = (np.log(y) - mu) / sigma
+                R = (1 - st.norm.cdf(alpha * sigma - x)) / (
+                    st.norm.pdf(alpha * sigma - x) + 1e-15
+                )
+                left = (1 + x / sigma) * pln_pdf(y, mu, sigma, alpha)
+                right = (
+                    alpha
+                    * st.norm.pdf(x)
+                    * ((alpha * sigma - x) * R - 1)
+                    / (sigma * y)
+                )
+                return -(left + right) / y
+
+            f = pln_pdf(z_line, mu, sigma, alpha)
+            F = pln_cdf(z_line, mu, sigma, alpha)
+            f_prime = pln_dpdf(z_line, mu, sigma, alpha)
         else:
             print("Please enter a valid value for dist_type")
             assert False
@@ -267,6 +335,8 @@ def sw_weights(self):
             + ((self.theta_z * self.eti * self.mtr) / (1 - self.mtr))
             + ((self.eti * self.z * self.mtr_prime) / (1 - self.mtr) ** 2)
         )
+        integral = np.trapz(g_z, self.z)
+        g_z = g_z / integral
         # use Lockwood and Weinzierl formula, which should be equivalent but using numerical differentiation
         bracket_term = (
             1
@@ -277,6 +347,8 @@ def sw_weights(self):
         d_dz_bracket = np.diff(bracket_term) / np.diff(self.z)
         d_dz_bracket = np.append(d_dz_bracket, d_dz_bracket[-1])
         g_z_numerical = -(1 / self.f) * d_dz_bracket
+        integral = np.trapz(g_z_numerical, self.z)
+        g_z_numerical = g_z_numerical / integral
         return g_z, g_z_numerical