12_GLM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**pySpark Basics: GLM and Data Preparation**\n",
    "\n",
    "*by Jeff Levy (jlevy@urban.org)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last Updated: 29 Aug 2016, Spark v2.0\n",
    "\n",
    "_**Abstract**: We will explore preparing your data for analysis in pySpark, in particular configuring the independent variables.  Then we will test it on the GLM function and show how to view the results._\n",
    "\n",
    "**_Main operations used:_** *Pipeline, StringIndexer, OneHotEncoder, VectorAssembler, GeneralizedLinearRegression, fit, build_indep_vars, summarizer*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To look at pySpark's GLM capabilities, we will use the simplest example of a gausian error distribution with an identity link function, which is equivalent to OLS.  In our model we'll try to answer a simple question: how does a change in size (carat) affect a dimaond's price when we control for clarity?  We will test the equation:\n",
    "\n",
    "    price = a + B*carat + clarity_dummy_vars + e\n",
    "    \n",
    "Don't worry about the accuracy of this model specification, or the assumptions behind it; our goal is just to compare pySpark's GLM function to something most users will be more familiar with.  \n",
    "\n",
    "# Stata Output\n",
    "\n",
    "Here are the results from Stata, which we can use for comparison because the diamonds dataset we're using for testing isn't actually that large.  Note that `clarity` has to be converted to numeric because Stata can't work with strings - this is also true of pySpark.  In Stata this was done using:\n",
    "\n",
    "    encode clarity, gen(clarity_)\n",
    "    \n",
    "Using the `i.` prefix on a column name in the regression command makes Stata split a categorical column out into dummy variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAIMCAYAAAB47eFeAAAAAXNSR0IArs4c6QAAAARnQU1BAACx\njwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAFpxSURBVHhe7d2Jceu4EkDRiUsBKR5F42ReMB40\nFhForNwkgr6nijXfgkACja0tzx/99wsAAIBbILEDAAC4CRI7AACAmyCxAwAAuAkSOwAAgJsgsQMA\nALiJ+yZ2P8/f//77z1zP3x//0tn+vR7+meZ6nvXUn9+nuf/j9c//nLJtOO3ZV+LisKartx8fO+c/\nN98BANdz70/svnTQ/TzPTBzavvnsT3JJ2raxve/4tJNKAMD9HZzYuYPl/alIdPC+Py0xh9ryycnj\nNz2D/v2+HnF9/R5V/niZV2LR86UsS+zq7RvTa59TO7zt63HdpP1R20zd5b3L/dufOOm+Re/79/p9\nvF/z93t/oinXQcmvfs4rPCP0Ye/8CNx98gSG8XFtOGg8AQDTOS6x84dTfJ65Q2Y5+Oxh+Hj8Pv0L\ntvx9ePpDOT5M/eHm7ukPxugB7nANh5g67N+HpS8faF9br32LWuKQcvfTycmSMLh22zaqe7XuXy/L\nkyH73iR52cHHd7l/SGTG4+/aU5sfERt3PW6Mj1WIMwDg7zgssXOHdDj00iscMtlBZQ/e+ODXh7Uc\ndv61+L1v7jCU+xeTgKjOSPuaeu2LVA9vn2jEVzlx0P1M1ZODdpl7five2x0R/+b8iLh7qdcZHy9P\nEAEAf8fBid3KA2/NQVY85Ncmdu32Na04aIuHt60fJ5GNT4RqB7+3OXHwz/zv+fL/bD9njSPi35wf\nkeK9GB+PxA4A/rJCYucPF31I9/iDMTmM7MGsPpHJysMBHQ61+DCLD1d3YMXl9n6hfuVg1uWt9rX1\n2rcoHt46SfHPPjpxsEnPu6zQPv/c1p+gXeK08tMiH9/lWX68VsS/PT8i9nXdvr8zPk3ZOgAA/CXH\nJXZWOMyjy9/HJQv+NTl1wkEfvef97PB6UiZ0uTr043vKwdj9l/fNtaqfrfYV7m2v5YC2B354PXwq\nY/+378X7UHeXTiryZ5QOf/WerH+FZEIZSV6KkvibsdHJUiP+Y/MjcPfJ+/A3xqfFxbGQDAMA/oTj\n/s8TmEL/4JfEo5SQrGcTpSxxOcZdE5h9/aolvACAv4LE7k8ofFpVS7jkU6ktn9aJ5BM7uc5MvFyf\ntjb1WlaMT0v2CSkA4K8hsQMAALgJEjvg4pJ/969y8edXAIAgsQMAALgJEjsAAICbILEDAAC4CRI7\nAACAmyCxAwAAuAkSOwAAgJsgsQMAALgJEjsAAICbILEDAAC4CRI7AACAmyCxAwAAuAkSOwAAgJsg\nsQMAALgJEjsAAICbILEDAAC4CRI7AACAmyCxAwAAuAkSOwAAgJsgsQMAALgJEjsAAICbILEDAAC4\nCRI7AACAmyCxAwAAuAkSOwAAgJu4VWL37/X4/e+//9z1/PGvHu3n92nu/3j98z/P4zPxCVycjn/M\nifH/eZrYPM0TAACY0y0/sft5fiJxmdcn4uOSyNmSpHmTdgAAxLGJ3b/X7yN8IvTf4/f1kk9A/P+2\nZ6U7ON+fGkUH//vTJJNwLJ8shXrrZIlL1K7nz7/f18P97zTxiNpm6tp72J+XNvQ/8YrvHa64D/X+\nf1I9sVPtf7zMK5Hu+AYqQZoo/nMmpAAAOMcldv7wXj7tCIeoPyR9eXweu0N0OXjtYf54/D79C7Zc\nJxcDyomLTzYeT5UopEnJklC4dtv3qHuV7++Tiri99k97vs8D/a8rJSzRtTJGrfjEr7tYpONXHd+Y\n7bfu1yTxL7wPAIBZHJbY2UNSJxjRv7PkDtFwaKdXOETtgZ0dzIXEoaOVuKR/ZnPJQPxaksxUFO9v\nE4JSMuNeG+n/pxTbX4y1i5m8tTe+MdfX8r2uH/9SOwEAmMOHE7uVB/YJiV36ciWxyOqmxhOLxUj/\n6z7wid2HErvrx5/EDgAwr0Ji55OIlcmCO1jjA9EdkO/D1JcnB7JNDJbDPjuwT0jskvvZ+6fJwObE\nIsQteT1KXAb6/ynN+ESv2/ep8auOb6wQ1+L7rxh//75OEwAAuKTjEjsRDs9wgNuDMz74w+EeXf45\n7hMV/5qcqvG9htpSuLe9QuLgyp8/8ftU8uEP+nDln9roZ+jkxccufk/S9nr/z9eLj9DtV0lbd3wD\n96w0fu61q8d/3yerAAB813H/54kC+8nKxxKXnpBY+B+xW2t88wRphvi7NuYJJQAAczg2sUs+0ZHr\nIp98XLVds1kVxyiRmyX+1U8gAQCYw6mf2AEAAOBzSOwAAABugsQOAADgJkjsAAAAboLEDgAA4CZI\n7AAAAG6CxA4AAOAmSOwAAABugsQOAADgJkjsAAAAboLEDgAA4CZI7AAAAG6CxA4AAOAmSOwAAABu\ngsQOAADgJkjsAAAAboLEDgAA4CZI7AAAAG6CxA4AAOAmSOwAAABugsQOAADgJkjsAAAAboLEDgAA\n4CZI7AAAAG6CxA4AAOAmSOwAAABu4laJ3b/X4/e///5z1/PHv3q0n9+nuf/j9c///LfYGO+K7Xnx\nu8L4748PAADb3fITu5/nmQf73zZDbL/ZRuYeAOCbDkzs3CcZ4dMSe8DZnx+/y4cb0Xvs9TSvRP69\nfh/vMlPv9Vz+94oPeLLDNbrv8+ff7+sRnhE/v9/+/idC8b3DtaL/Xe37L2321+NlagT9/g3Xj69i\nHMr2f6LWi6+Tjb+3Nz7t9jfio+e13O8nzG251s4DAADKDv/EbjkQ3WH1/tOUP9zi89AdlP6g8+XL\nn7jCQbn+0Csf7O5+j8dTHdRpYlBtf6R8f590xMmCP7ztW3v97+rcP+Per/9kONI/p1G//MBh2+4x\n3v+x+2+PT+v+9TI//6Ln2fcmySUAAPuclNjlyZhLYqQsv+QctOX6kLMH98GJXXKQ54d7rf2x4v1t\n4qaTNHmme63X/67O/d2P8adA7ionLpX+jdYfanDdpnuM9N+r3n9vfLxW+5t9i+dzsT8AAOxzTmJX\nONhcYlM/MD+V2KUvVxK72sHsFd/TOah7/e/qJQK2PO7fyk/c9tZfYdM9ViRC9fE5pn+t97Tru2f+\n93z5f7afAwDAWoXEzh8+G/9EVD3Y/MGalPlPUOxLvnw5aF0idnRil9zPPj9NFrYf7OHQjl+Pkode\n/7s699dJsL/3cOIyWN8mqO/60fNXaMfY9zObf53+R4r33xufSOs93fj4547/CR4AgHHHJnbvQ8td\n+YEfkqvoip8Tkh97mUNYH8ZNhXvbKxyg4RO7+H3qcF3dfn04+9jF70ni2Ol/V/v+NuEIr4dPhez/\n9hHs9K9b31J9WNX+XvxEa/61+l+Irb2WZ+yNz1j7e/FxfcjvDQDAfof/KfZI9iBelTi0uAO38kEL\n8BG7/yQPAEDDtRK75BM7uQ46AM+6LzCk8GniYb+wAACwuPQndgAAABhHYgd8WfLv/lUu/p08AMAI\nEjsAAICbILEDAAC4CRI7AACAmyCxAwAAuAkSOwAAgJsgsQMAALgJEjsAAICbILEDAAC4CRI7AACA\nmyCxAwAAuAkSOwAAgJsgsQMAALgJEjsAAICbILEDAAC4CRI7AACAmyCxAwAAuAkSOwAAgJsgsQMA\nALgJEjsAAICbILEDAAC4CRK7Sf33H0MHAABSZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKDdLjt4\nPf6zSY+9ns/f5+PlS+7lryZ2vfH9K+Pf9/P7lBgU+v/zjGJkruePL7gtH4v/Hr+vf/4l8fN0MYhi\nxPy5pzDnH34C/Hs9lnE+YIzd/dX8uqgz1v9M/d8s7BfxlQQv7DPheppXUnreFcP17/X7KNbv3z+4\nVXYgQQsLV9jJdtONWQb2r+mN718a/zFmI9D9l03jL8YkbJbRRvzzNPPlsWyOzJ970+N99Fp4PSZI\nbE5c/1P0fw+b2NWTKbtfJPuL2z/eIVH19fsDef1hYpk+65/7pbN1/8itE7uS5DfyKHBSV14LcbNB\nMz/H9wt15bXwfj3QoZ69TNCTyf7OxJfyraT+X9Mb35HxvzU1v54/aWKXzM1w7ZiDU5HYmMX9Cge7\n/tm+9Mfnz81JYvd6mcM1bPJRktPc/8O6Mnu5+8TEzJn3pzfL/LF7/U+8BlWio9anXnsj50tP7XwT\nR6z/1v27/TeS+ub9o0nmEo/8+tia1Z/Ymba3nmzjHb3H9iFMMGHngxrfkPxlSaRL7JK+ynsqbbhd\ndpBOnHzhlAITSODjuEt5NmlswKN7m5+TzSC6gZuMy+TWv9HYctXGUdKGv6g1vqJXfl/yMX08v/zH\n9joG0WH2p0i/ZW36NS1rVX6MEzvxd+fP/dnEzqyP8E+9Fpr7v0/KbLk9A9y8ieu4uROtQXVwD+3/\njfOlp3e+WTvWf+/+vf5Lf+P6IfGZRjTu7kfX/igiRvzn0nRvcfVDfPwncMmeHSVv6llJWeDjnz7f\nuXV2kCwcmWSFFfJe5EZzYQel16zCn71ifmN4Hxq1wR8kdf+64sYY6ZXfSmleljbxHRv71N7rP014\ndWIX+1Pz5w/IEjq1Fpr7f/xeOVD9G+M6cvAm9Y13udRP9v1w5Yd/+XzpeM/vVHy+WarPwwbu3+y/\nFSc9hb432LWY1F2uTfE6gsyDah984mZinYZ/6cfzJXMizT/e98vuTWIXkYnkgzMwMdNJaJQWWXXh\ndRK7bvk6MjEQjW9Rr/xGSvOytIlv3dhnF6//n5/3nGgldn9q/vwB2V4vB2u0Fpr7f7xu5ED1b4zr\n6E/kxFI+uP9Xz5eOgfPN2rr+B+7f7n+BxHHmvUhi0tofbHkek8B+4veuL3vNkqy+LxMfV/3PJnYS\nmMIkjmZVPvHSxSaTMFnIJrDZImssPBmotMwNVmhCXr7d30vseuPbH/97k/6nm0zxTx0Sk5k3060q\nc2FJ7P76/Lm/NMlxe3O8Ppr7f7xuqomdWm/2HsuaHNr/tyZ2Ru98s+J+rNS7f6//Wf0ojpdn+5K2\nP+yv9qVCucyNd7liy6LcICOxSfZz/wlgVCF5vnKzxM50NLnURLLBT9+TBDYpN3XlX7Q1/9sttNL9\n84GxwV9RXv8ot03q/i298R0Y/7uzm0HUf/nPdcg//Wabzz1zVXeWG9HrWuZE8pqsQebPndlD0I/r\nO3GS9ZIlItHYh/3frKP367Je/DqT+4QDOtR5muQx/lnPn/r+P3a+NHXOt93rv3F/F992/486+74l\nHWtzqaQqnmPuSvuXltf7nrwveYaeI/V78Pe8ScnAAgAAxMgOJkViBwAANLKDSZHYAQAAjexgUiR2\nAABAIzuYFIkdAADQyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0AANDIDiZFYgcAADSy\ng0mR2AHAuOJXaslV+Vqt0pfa75V8XdSq7wsbkH3lV/6VU8nXYsVfp+b12jdav/h9t9K+Qp3DveOg\n+9/6Sq7yV7pt+kqv5KsdV8wh/ZWQcpXmSLV/C7KDScmgAwDGvR7qMJRDsnR4yiF7cOJlk6LonpIE\nHfqIWl88m3RFiZW0J0nA1Hfnyvvj8l79d//MfYqJnXF4nwtsu01SrhMf234V/yVxc0nbUlz50v1q\nfSckvpv6aBO7erIW1PoXIzuYFIkdAKwTErssqVHkUN90ODfIM+N79tqwWiex0/TzdftsshMlclq1\n/Y3Ezraxcc8aedb7Uyx1Jc8KydFAklRKzN78p2KtcGb17TN3zBtf/32V2jbYP7KDScnAAwDGxX+O\nrSdV8ulNOynYRv0pb0OC0+STkff9i32I2qCfL8lC9JpLpvQ9GvUDc592bI//E7fjPmWzz64mPvEY\n1Me43HdRr2/rmJi8JOGz5Sv7qdqcJ54j/XPe2YFrCBcXFxfXpy98xtAndjZBqh+am5nDOP3E7nlS\nguPJ4d/4+Mj+SU81IPlkTJIU/afrSKm+ZZ5bja2x5dPQpF3qCs9KkrFO4mPebduRJk5B/mfYXF4/\ntLHYni1UH9b0jx1lUjKBAADjWonK20mJ3c9TfYLTSYD2+2n+KfW08hMSu774k7ToihKvjB3nwqdq\nA3+GtVR9m3jFz6vdf1QyD9f1j+xgUjKoAIBxQ4mdPUSPT+zk4E+ShVICZD+JMQf2hsxHEqb4fjbR\niO6jEyr7p77Gc6Q8Lh6u30zsJLY7kp1RNo7RGBaSrCwR85JPxoKR+vo9tg26r5VPCgv3t/GtJG5Z\n/xSyg0mR2AHAOHughk86sgM3pZOYY+hPXfKD2R7mnba1JH00SUFKPV+Xh6QyXFkA2vVd29WlnyEJ\nTNauYyXtMM8KoczbV06M7PsKgz9UX8Uwv039T8AuoVzqlt4jav2LkR1MSgYVAHACOaCPz+y65NA/\n98+z3yVJyRfC+ueQHUyKxA4AznPGf6C4yf45rv7ntel94NM6OGQHkyKxAwAAGtnBpEjsAACARnYw\nKRI7AACgkR1MisQOAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyKxA4AAGh/JDsIX4Wi\n/oOT4es/Nv1HE/09V9SNvzJk739d/KuJnf0PaUo/9H9Ms/+VOe+Y1/7z471yUX1+i26bv5KvZOm3\nP/nal8rXuVzXwPgcIflanXjNHTA/xKbxF73n7y339vQfp0m/UsxdR33Lg/uapw//B42VUv/sxVc9\nDNq7Ps/eX8bb93c+9gmHQTTJf56P38fQl0LXmEBvSAolOZg5sZNN7PHQX5TsvwMvia/5+Z38+PIw\nKbPNple+KD+/xy2K5ba6vb32G/bAXp5pyxvtvJaB/h0gJL55WHrPP3v88/HS/W+Xj8Wv3v/+83E2\n2QPSOXPkt0vsu5fbn/aeCy99nsm5V5qMUPbuj4P722nlqb+V2JmgvCe+/tmzwQuHiz443p8UuOv5\noxI7VS5BL5k6sQvJjUpyTK+GJ559PXqf1iyvPn8lP1bLY1z7k3GRZ0Ttt4d23C57jx1t+Kh1G8Mm\ndkzKSc1IfIOPjL/R639aPhC/Zv9zh8cfHXliJ2s6jJcdXzN+MkdDgq7nWHhPucwkdj/xGbA+0XP3\n354ghvPsiDPmSpbxyK9j+jm+P5X16p9dnvpziZ0EQ4Ijm6r8GCd2tcD5H8ymEC849xvWUu4XdlTd\nTsaoPDhi0cmE/rxocklssoPVx+R96XLHHmiN069e3nv+OLdRxPUHFo59Zhhj9/5v//llnbHx2SrM\n95eMn71/HJvxjenc8Y9jUKrfKm/Hr93/oPd8nEdin8Y8+5TNziszNmHfNj+HqVibv0G2H2z9xS98\nQBDde5Rrg7uStqJjfH8q69U/uzz19xK7sLGagAj9CZ4mf661sTRBTIIqpI6/z3sxZle+sGdN7JJk\nSCaV6lv6CYSbiKUDemti13v+uFLbxhaOa4Opa67nS8a8dHhf0+j4bBViE2KYjNdgfMX54y98/wvP\nd/LyXvza/dd6z8fxdGKu5qMw8yp7TfTOB0PGU79F5sDWJabn0wg+sdtqfH8q69U/uzz1BxM74+fn\nvdkeltjJpvH+321HLDqZ0J+Vb4r2akw8uxgLE29bYtd7/go+CU8fsW7hCNvOXcnFJ42Pz1bZ/Wyc\nw8E3Ht/Txz9I2lfQab/ub7v/Bb1yHEzmUGe9lvZ5IWPVOh+M0r9jJ3OisdVVhbm0VvwXKKyxfv9P\n9eqfXZ76m4ldJP1TrF6YcbKWbwrhN/hAfk4CXyGLduR9LXKofZVMqlI8ohiXD+j660Gv3Co833EL\noDbhhd00s7q+nm5/5T7uHuP/PtUVjI7PZjpRsWMUfh6P7ynjX0iiwuHpmtcpN7rxa/V/4P44247E\nzmifD1Lu5t2bHfPO8zRbR91nBRK7rUb3p9r50qt/dnnqbyR2YbHYK95ow2t+MSSvuSuKo0QyKfvv\n+XSfIpjgBjb48XveC1s2FV0Wl68jdb/FTqjQ/mRi+ck3Upa9p1e+qD9f+PsU6gW2fjKwgR6jdGyS\n524ct+9qjc9B1BpJw9yK7/njn46fXK3xzcuH4tfof//+OFM6vwpJWjI27krmb+N8cGP7+H0+3S98\n7tLPaAnPX1Mn1e4f+tr7v9M6X3r1zy5f/J1P7G5GBhYAACBGdjApEjsAAKCRHUyKxA4AAGhkB5Mi\nsQMAABrZwaRI7AAAgEZ2MCkSOwAAoJEdTIrEDgAAaGQHkyKxAwAAGtnBpEjsAACARnYwKRI7AACg\nXSc70F/XJddj2/fl9cjXv/S/q9V/fceBbQjfLyrX9b4rVn9dib9KX53y/mqdwleahHEsfmWXl4z1\nmq++UV/plD2j9ZUrI/3rfGVLaY62+ql163eeb8RzKB+bfv2h8TlT8/k7x29kfJrP3zO/nN3jg3OU\n5oYZn7/Lz8U/HQNt7/ocrJ/MxcL51ytvnb/epT72kU0xTnjGErAzmYE6YeLrfm4hg34sNymXs8wf\ncoUD0I7LQw6weGL1DkUnHHyV4gbdHtfeJY55e+33N74P117/evUNu+B2HMad+vZ5K55ffH+1/tj4\nnKf//Hb7B+ZnM7695+v7rZ1fxq7xwdmud758m5njZv5BDKzvprH6vfNv5Hy08zY7f1OXTuzMC7+P\nuIf6S5gLkzLZvE3nH9F7QtDkKi5odf/nzzLxdcDtoOn72I19qV/bNLJ+biD3P5WPRTbBwuHVOERt\nbEoz08enNWnr9MHux+S9cNzCysZDLay3rH+jB/cyvtV716ysr59v+5sGwPShvriz9nvV8fmQ0efX\n2m+V5udgfMvPH5tfrflx1PjgHNm+K+OTzYO21vki4nJ7ftj3mTGWuWFfc++zY29+XnN+hHvL6+F+\nyR7s18T7Hqptujw+365u6W9+6Thts/L8yIyfH9Up1ysX9j3t81dcOrGTQMWdfJmFFMfdDnY0MXX9\nEFhNv8+RRRjff1mUgdRLgm6Cm06EnyTQuv1B+fnryAQ4k41tNnGiyd+YWDbuhY7be5p4vvymtu7P\nsJpeSOsWZrl/fszflypXfV59MA/Vj9tQen6Ime9/FsNGfa82Pp/Sfn6//aI4foPjM9b/fKNO2ybX\nOeODc+h9V8ZnzTLQ9cP8CtL9x49zVC710+m04fywcyy6r/k5vKd9Pkp74vK8fX/buvOjrL0/hPGo\nnX+98qSNdh7U94/LJXZLUFSQ9W8j7yvuXDuwgTwnubfQi0zIM6OJL/WaC7PQxtLGUXz+SnLv85QO\nNdfud0wbE6t2cIbxDX1P7reKb1+y6NYszHL/0kSg/J5EZ3F1NeuX+mhe9TG0c+sl800v/qBcX4wl\nNudZlVitGL9MJb7955efPTI/jhgfnCMeG7nSPVifHeGK50/jfFFnheX3n0Cen0wXU776/NB1gt75\nWKpXavNF6bGLr2I8VnPrMbmXH7/Ru/f2h9CH8Az38zKHxsqX8Uzmn3LdT+yyhptFtXYS+oHRkucE\nAxNf6tUXpiz6PAPPFqZRfP5KMgFO4zeJtO2Vja8w8WsHp/Q7eb99Tu3gq/HtyO6/YmEW+5fXz9qr\n2fvUF1dXr34nPjbOG+r3E5tzDT+/1n/7euHQ0yrxbT9/fH715sfW8cE5kn03O182iM8XGcvwvwP1\nmjw/mVZbzo+kTszU18+PleqV2vxn7U3s+vtDtl+o9d8uHz9/xXUTO8N2NJrZslGWJ7WjP4q2A1PY\nwPVzHAlcutDtxmwCFyT1bNCj9uhF4ssLj688fx0Z1LPYuPc2PYlt5T3Vg1NNZHcPNWZ+gZQmrGtX\n7UD39dR8qd8nb7tud/KzbrtRu39Rr36h3Lazcv8sFivq634eqz5+QfH5K9pfHL8V41PrfxZTRddr\nxXHP+OAcEu9437Xxrw12Qe980eVh/gXJ8+182HB+mGfWzg55Xv1c6Z9vf9vo+VHf3/R+oH92YxrN\nEX3+9cpjtqx+Rl8nsbMNNYGQK5psIYiB/Tm8z15L51pl5gHljDd+T9wGuZ5PVyc83y82V24C/nLv\nD4spbOah/PmMN/eB568gdc+STUjFlof2vye4n/DvfulyT8U4f0xl4eixCVdyAx3jcmzr/VN9UG1I\nxzcv7+nVT+Jqr7T9aXnet3b9wfHZrbbx9Z/f639QG792fDvPH5pf7fmxb3xwqsHzpSWfP2r8kvPB\nXHJ+xPc+4fzQy6DZxt759ueNnB+1/U209werd/71yo1kHym242Kf2GGcDCoA4KpMokDShC8gO5gU\niR0AXJP+VAX4JLKDSZHYAQAAjexgUiR2AABAIzuYFIkdAADQyA4mRWIHAAA0soNJkdgBAACN7GBS\nJHYAAEAjO5gUiR0AANDIDiZFYgcAADSyg0mR2AHAKPV1Ud3/aLB/f+l99qvBVn4dXOkr6/70f7i4\nEd+vGflKsZbB+slcKHwXbK/8/dV09faRHUxKBh0AsIIciiaZ+HlWvlw9YQ7qAxMP+S7Y+Ev65dsp\n6l/a/xccG999/Pe8Rl/Oar89xLRvbITG6ofvAy59B6zolQs7bx7yPhK72yGxA4B15OC0h+bPM0+q\n1Jf4P3/SxEN/wf7apEwndvZ5rRO8IP2SeXM/lRjF5bb99n0muVAJQ/jKs6Q96lNF3b9wb3k93C9J\nLlT8sk/jOvFtWZ6XX2vHocwlZlk8JHb+x7aBxM7HtzrkvXJh32NiHv7pX9bIDiYlEwAAMO71CIeh\nTiokCYo/xVuSoiJzsK5NKHRiJ4nAmrxO1w+JQ5AmJnn7pX7yPN0Hk2jFxcX2+eTjfV/zc3jPyySa\n8e1sMvZ+/sr4ftzexE74Pr2vNPEK8Xj5pFr/mbVXnrTRjsNAYrc0houLi4vrkxc+QA7DKFORxOj9\noylLDnXh/2xbVHp/hz24ozFP6+ukIFzx4d1IHEpt9YlJ0E3s9Cdu5ioldsV+F+q6y7dxbXwVHbv4\nWjsOZfsTu/QTuvwTvNCH8Az38zKGY+VLPJPxV9hRJiUTAAAwJvz5MbnCwbs28aglOA1yML/rdA7m\nIT7xsEptVa/J8+uJnf5ErfB+Ue23qV+LlSjVW5HYnW9vYpfXt4lYVF//bPsfxbxdXkn8K+0jO5iU\nDCoAYEQ58Uj+NKsSLf2nzkQ1wamTgzs7+LPMqU7/qdMmHlF9Xa7bnzzfJg1RIqKTLF8+nti559Vj\nsjK+H5d/whbal/bIvy973b9f148DqBI58wZzjkc/98pjtoxP7G6HxA4ABvgkxX3KEQ7K+BMQ/5o9\nLMNr5no+3XtKyUcjwSmK7x3dLyQJI+x7wz3spQ72pJ/mkvbH99ZxeLk2hX7YRDMqfz7dzy43KX9i\npBO/ZhvXxPcrdB9LiVM9sXuXhfql96gY6Ph1yw2bMIb3FNtBYjctGVQAwOfpT9+uySQql0ma8Elk\nB5MisQOAT1Gf5lw8YdKf6uBvITuYFIkdAADQyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gU\niR0AANDIDiZFYgcAADSyg0mR2AEAAI3sYFIkdgAAQPsz2UH6dSnxVfkuthPEbdj7Xy2XexxOf+WL\nXPF3muwt735ly8hXurQl49z6SpfSd7UEST/S+dG9f8vu+PSev7f+AUbiW9Vrv/rKnsIzmv3rxt9o\ntn9v+/bPb2xUGnszP64hzAt1FoU2R+1M59fylWHJvI+vy/TxDlrrV5f5K9mDeut/b/niTyV2IZmS\n/yp32HOzL1buWL40eru4LVvJwB7ObiSN/u0sL35J8nvi+0OxWj5APT99Xj8pEGGDLBY37z9gV3yM\nzvP31t9nLL51vfHX5W6TS9ZRr3/N+Pfb345vr329/uFset+V+O/dh7XN50P4Htdkfpj2RvcrtV/m\nT5CdK3LPqBz7tNe/W+9LsV7vvfXf3x/az0/9yb/nSUCi+CRscMPmHh8C+guW/ZXcxx4cS1lt07h2\nYhf1T0+aveVKaWIncZH7de4Rs0lZPCB2zPKNVi+QN9/+2twYvX/Vrvisf/7e+ltV49vUG3+9cfr+\nrOnfYPxH25/Gt9e+/saNc8l4ZImPGufq/u/tOh9afFveiaH+2b7UPjeWcjMXzby6E7uWCvGVqxWT\nMzXXr58Py/j39rf1+0OrnMQuUgt8rPkb2c9PUib3Kz2nt0BHyIQ+nD34lv5lE2dvueUOQLco41ju\nT+zc88MnsH6hFP7UbttVGBi7eZjnvaTctq/0p5H+/at2xccYev7e+vvV4tu2dvx9++Pn9Po3FP9e\n+xvxTRTal9Tt1cfR9L4r4xMPT23+Bb1ysesTO2mMuac8Q+ag/Kjv5+a0v9Szk+RHleEoY+vfjUVc\nPrK/jewPY88nsQvCwlLk4/B4LJoL12fpy8BMlthp6iDM7Cr3B997Yo9M/L54c3u+ZDzyxKV2cIe6\noQ354hy7/7BV8XHGn7+3/nbtxKhmzfiX+yZW9a8S/7H219tQK7P3fb/m37M6TtgqnhtyJXOtt/8f\ncT60vO/vD24zT0TrfrY//n1CfnZ9MveIXr8DPXbxlYzjxwys/2S+uNda+9u6/aH1fBK7xe6FKwsy\nfa9MxtK4LAtwO5nQp7OJamOjOqQ8xKw/8deyC2XFwR02yvfzkvblavcftio+ue7zz25/RS2+baPj\n7w++gfuPxScvH25/Mb619uX9y+YbTpXsuzqpl7EsjPnnEzsj+stP+34y15byI84VrFBc/4Z9XecY\nvf1tw/5Qe75BYhfJ/48U+W8+yXvizUGCHL+3OLjOEQvw8MSuMEnsARcm1gnl6cR1Ezs+EJP6b/59\n2espe+9K/O19ywOTttGOb9rmoHX/ot3xSWXP31v/QNX4WrXx64//mjaPxKc8vyrtH4hvr336vu04\n4WgyPtnBGcW/t//vOh96ZH4V5sKS2EkSp56v6uj+ic2J5p9V2Z9W7K9uH9Bx7+9vzf1hxfPFn0rs\nXMBNsN6XCr4NXlxe2KTtYi3XT+//+H2a3+aWe8jCDGXxtW3hSd2jZfFRk2ZvuZ2ocXnWdx2jUmzq\niV16f13X10vuby59n2R80/Fv379vb3x6z99bf5/B+Ib3FcbPtLA+/mpc3pfeCN9lef/a8e+3vxnf\ngfZlzyjGAKeIx8fEPQhz0ert/zvPh6rkvv7wTl6T+5TOj+Wgz+b2+xpsA7z6/tTbXwP7vmxiiN75\n1t4fRp8v/uQndncgAwsAABAjO5gUiR0AANDIDiZFYgcAADSyg0mR2AEAAI3sYFIkdgAAQCM7mBSJ\nHQAA0MgOAAAAboLEDgAA4CZI7AAAAG6CxA4AAOAmSOwAAABugsQOAADgJkjsAAAAboLEDgAA4CZI\n7AAAAG6CxA4AAOAmSOwAAABugsQOAADgJkjsAAAAboLEDgAA4CZI7AAAAG6CxA4AAOAmSOwAAABu\ngsQOAADgJkjsAAAAboLEDgAA4CZI7AAAAG6CxA4AAOAmSOwAAABugsQOAADgJkjsAAAAboLEDgAA\n4Cbeid1///3HxcXFxcXFxcU10aXxid2kSoMJAAD+NrKDSZHYAQAAjexgUiR2AABAIzuYFIkdAADQ\nyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0AANDIDiZFYgcAADSyg0mR2AEAAI3sYFIk\ndgAAQCM7mBSJHQAA0MgOJkViBwAAtFtlB/9eD5vwyPV4/fOvHu/nee79R9w1sZPYhjH87/GyY/r8\n8YUHSZ5hrqPvv8frEbXt+fx9mhhs8/P7lHtsrp+K11ZyHXR/fF8y9/y1ZZ+TubJnf9xb/9uWtfL4\n1d1I9h698fx7/T5Cmb2eZhVHfp71uoNm3l/1fe11pc17p3iPLc9/v6dXyxe3zA5m3xhGyODejV24\n0UINE/3QtSub50WTET1vbTxUW18Ptdk3mY3gwL5m6+rCsbyqdeP3aXJwpO17PUJy4g6VkX21tf+O\n9H/2/bvafknMovUi6ztbT43N7ue5JIpb5tFd9tc98+Pa68/p9W+k/38msXO/jZqF8Yp+61G/UYXf\nWKVumPTxb03La/UNLvmtwkzwZWM09G9kOxaA1L+VzqYWpJ8qqEXaiW8yNpX3fFNzweq++SsJmXrP\n8+esxK5839766ZXf2sj4fV2e2Mk4xW1876OVaSqK83hF/4v1J1Jrv45lto4kRq0JsecTu969vdr+\nKm234/WMz890rnxqfy3Ft9u+kfmn36Pa3tq/4r7btsVjFe4TvxbeV1DqX6xXLv7UJ3ZuYPRgp5Pz\nHfxoMPR6qN2//BtR/FtWuiHa8o0TX9p4KybOvckq45e8R8Yqit9QfGXMN8b8E5KNtdDO+m+ccijH\n/Zefy/fYys3netus3vrpld/cvJ/YRcIBWJkDtf1R8IndEjO3nqJ4ZMlHJVbJ+wbnk3n23v1Vt7fY\nzw/sr7X4jrSvNf+Gzo/G/pU9T94bb27mF+342RLv0t5X61/QKxd/LrHTgZT3Jq8NLIDy/dVvX1q2\naMM1uDAVqXsrvbhL/Aqr4P3nidH4fmDjOYrMM72xVDemUvwO7usy7xtzvTeOA+vrzq6f2KXrpzVW\ndn4W3lPeH50/ndgZIWb2MmuoGQ+zVuLEIP6lr5QQNO3dXw1pe/es/MD+WovvSPuq8R49P5pxTH8x\nyn4pKjyjNI6t+SN65eKPJXb5b58jk0Er37+T2HXL15FJcSu9DaG78QzG9wMbz3HSjUJcI7Fr6K2f\ngfV1Z9dP7MbaJ3Mh+zTDa82Tv57YpbadGTb2jcSgaPf+6p7bPSs/sL/W4jvSvvr8Gzw/OvvXuw3m\nfemfy2VtpflH1l6v1r+gVy7+3Cd23U9ABg6e2v3lT7Hp6zKYy+LLy7e7XWJn9OKXJ+bpYhyKb3fj\ncc/c+knqdvnCt21VKz+JgWwe73ZK/bTN9l8NyPq6vX+leb96/XTX17fi/xn18buCfA5lZE7K+DTW\nUOvgGen/yMF1ZaPtl/UZL2/9p1C5T/7v0pkxesc+/t99e/dXaU/SnNJa/sD+WovvSPta82/o/Cj1\nOeFilsVSx8Wvo2x4jd786ZWLG2UHYcLoaxk4CfbT/AYSly8BKtdfAt+/v4g/LpdLD5wu3zrBpe4d\nNeMXDpVaudGKb15mrmxl7d94tinNL73RGnYzCuWqjUmZueQ/lyL/TDbabf2TzSS59/sK99m2frLw\nb2zfNFrj92Xp+igc8sXXc82DZ6D/IwfXlVXbn63PvPfJGFQSpHgtFm7RtHV/zZ4Z9SX09ez91f6i\nWrn/SPuszvzL+xDeM7p/+bYUCuI2yjoKuYh7a/n+pTbKfXrr457ZQUX+G8m8ZNAB4GpGDp6WvfW/\nbfb249pI7CIjvwnNhMQOwBWR2IVPZu7zQQKuYPlUj8TupkjsAACARnYwKRI7AACgkR1MisQOAABo\nZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyKxA4AAGhkB5MisQMAABrZwaRI7AAAgEZ2MCkS\nOwAAoJEdTIrEDgAAaGQHkyKxAwAA2js7kESBi4uLi4uLi4trnkvjY59JlQYTAAD8bWQHkyKxAwAA\nGtnBpEjsAACARnYwKRI7AACgkR1MisQOAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyK\nxA4AAGhkB5MisQMAABrZwaRI7AAAgEZ2MCkSOwAAoJEdTIrEDgAAaGQHk/pOYvfv9/X4zz7bXs8f\n//qon99nqGuvp3kldna59+/1+9hT/2s+FJ+fZ/Sex+/rn399dPyr9b8c36Rd/nr3QbfNX4+X6XWw\nt/0D9UMbV6+tI5zdv73lV3fm/qjL/DU8P/fWN5rrZwbfPL9G4i9H06Natib+JHYf9nroybCNDOpn\n+UXxnkhuoj6WU7tD15d5an5+T96zyxfy+uMhCygdC/v+gfrfsb//I/0LG0v0Nk/fvzz+9fpjzz+V\n3Rhr68/1Z2leJZ6b2z84fqYNrQ37PP350zbYv83lV6fbf/T+uHd+HjC/m+vn6s4en/3jp+Orx2NN\n/P9cYmeDZQJsLxP0lzng32OrMmI96HYg/OtLZh0FulX//SlResXjtobU/Sw9MaVLJgbxwm9yEzmJ\nicRreGHsre+FxRH+6V8uKdb/mrH+18tzWf9sTGpzcmD8m/VzH4+vb9/7aj3br9dWX9a1f3B+Gvb1\n0SAeZv38Se2dn3uf/21n74/K3vm5pb60R9ZNuKYZG/Ht80spxN+2J33BvKeeX7Ta/qcSO71h2kDG\nfyr6+UkOehmoOM5WCK4Javj5/Z6B+vN+YqflE7ltZONyi2+ZvHGs9tYX0T3sOJbGIr7HMWN1jF7/\nR+Ij6v0LG91L1oktj/+MqrnnZeupW/+L8VVj3tp43d5w9PzozU9H71OfMTp/avbOz73Pv5p8fbSt\n6//e+bmpvrQneq2ZuFzeGeMztr5FMf42vmHP9O1L8pPx+P+hxM4E3QShyWfRy8CUE7tkcGMD9e+R\n2PlJt2pR9xdGOlH1wttb37wSLya1SHJb+nimXv/78Unl/XPxWe5R3/zLsRmvLy4Q3+oc8G1rbvrr\n29+bn4F9X/PZZ3DtGZ8/Wq/+3vKZ+LFd1fY1/ff33zw/99b3pH3NPfSqBvqWcXVa4zO6vltlYQ+1\nucNL8onGL9eN+JPYvUm2nQZRgpzF3gQzGdy3sfrzJ3b+t5LmplDSWxh5uZ3kjfJ19X279fUuL7CJ\nemNhfdT6/qflBap/abyMYv/r4z9WP/Lt+NrnF9ajfb3wS522qv29+bkgsfNWPf8qztofI3vn51Hz\n25bPlth96/yKDMbf7gOt+Dbi/4cSOxeoZGD8INsAS5DMILzVgm8GM72HN1g//3f6ti2MbyR2dqJ2\nJ6Sb4PmE9q9HldPfcPzPuvz98976io69Ha90E6suzFNtjV+nfKR/+j02RsvP3fFv1f92fAvP1/Mn\ncP1U63K4/bXx888bmJ+118/VX1/OSfOzW3593fVhbY3fYt/83Fh/xfq5qrPHx/6sywsPK8Zfydq6\nMv5/KrETdnBMwMIVxz0E012P3+czDm75E5/x+p497MJ72oPbIvU/Kml3dCWdE7WFIXQMdf993VCe\n3WNvfccuiMJ7ktfttX18ttsTv3b5UP/UOL+Hd3T8a/WNb8c3XZ/mas2PbF6Ptr81fq35qcqK7zlb\nb36J8+bn2PMv6iP7o7Nvfm6vP7p+LukS55czFv987NbE/88ldnchAwsAABAjO5gUiR0AANDIDiZF\nYgcAADSyg0mR2AEAAI3sYFIkdgAAQCM7mBSJHQAA0MgOJkViBwAANLKDSZHYAQAAjexgUiR2AABA\nIzuYFIkdAADQyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0AANDIDiZFYgcAADSyg0mR\n2AEAAI3sYFIkdgAAQCM7mBSJHQAA0MgOJkViBwAANLKDSZHYAQAAjexgUiR2AABAIzuYFIkdAADQ\nyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0AANDIDiZFYgcAADSyg0mR2AEAAG2q7ODf\n62ETGrker3/+1eP9PM+9/xFI7AAAgDZldiAJ3rcSr9fj+fvj//c3kdgBAADtctmBfFoWPpX77/Ey\nidTjV+dw1cTu57nUNZd+z+uxvL58+rckas1PBP+9fh++LL6eS2VfvrT3fT/Tj6PJfQEAAGKXyg5s\nUvfOlEJitCax+0k+TZNELrqdE5K/kGyZn/V7qvc3mp/YmXvlyeQ5n/CR2AEAAO1C2cHP73Pwk61q\n4lX4VK2U2NWStmBzYmck5QPP2orEDgAAaDdK7Ex99emevO/TiV18/7M+rRMkdgAAQLtUdpD/v1El\nWcs/dSsmXvJpXZwY+k/vjk/souTR3Cv+d/QCm9Cd+GmdILEDAADa5bKD8H9wCNeSmLkkLy5z15JY\nSUK2vP74fT7dz+4e5fpr7m/ZZK5SFlQSviPJ8wEAAGJkB2c4+dM6QWIHAAA0soMDJf+pFnNlfwY+\nEIkdAADQyA4mRWIHAAA0soNPqPzHjZdr/b+PJ/UAAABiZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsA\nAKCRHUyKxA4AAGhkB5MisQMAABrZwaRI7AAAgEZ2cLD4K9H+4nfFJl/r9nj96gjUyytf6abu0bu/\n9f7Py+j/jIx+xtrybxtoX/jKu9Z/HTv5Wrzou4+Davy8zfE1Rtr3bdX4HDQ/GvEfmt+n2du/Xv29\n5fMbGt/a/Ehe99d7HenY+Us/o7X+mvcXE4xPq39dH5i/1faNjN+/9CtXG30ksTuQfPPEksy5gdo0\nvwbIwF6OnbTLZLbfxBEHoFmu4+Un8XD9hR2Hh2yg6cLS77c/vxdO/ry0/Nt67Rtb9OFgqRRbtfgF\ntXLbnp3t+7Z6fI6ZH834D87vc+zt3+D83Fx+AwPju2Z+pHr7p/9Z1p5cpQc075+391rjM9C/pv78\na/e/V7/XvsHxe//s3l/78IjE7kSySLPxO4hMjquxm1LcYfvJzrJR9MoT/lOh9O0D9cPm1NmkxLqF\n+W2ufclClj4W2qc3oDcbk8qhEfTi1yuP1OJXbd+3NeMzHv+q5v0H5/dp9vavV39v+fy649uZH6H8\nfbViU9g/g97+MHR/41r742Lb/rJ+/qX9H68/1L5s/PIPiux8qrSPxO40ZiBM0M8iC+9y7MYQ/nTg\nJnr+p4RGecROWn2odetHi8u+t3QougXiNi9dHpeVyr9p/8YRNoKXlNv+6dj34rc3vs62jfd87fiM\nx7+mG38b09b8PtPe/vXq7y2/gc74js2PZU2liUWquH961fU3dP/++v62bfvL6Pyr9X98/o60rzV+\njp8/lfuQ2J2h8dvSUWRiXZGbkG7iP18Sh3Rz6pU79Unbqp8sBrVJ5fwzooWXbmTthfN5+zeOELtw\njyRe+udC/PbGN9i28Z6vHZ/x+Ne07++E95Tm97n29q9Xf2/5PbTGd2R+JCQ+xXIXy9oaG15/1fsL\n/4wLjs+2/WXt/NP9H6/fb197/EZiT2J3MLcYz9+MZQO4OjuBGxtTtXwwMU7r60/b/NXaeOxzwljl\nC9OOZav+R+3fOLL+JP3vxW9vfBf9je072vEZj39N+/653vo51t7+9ervLb8fPb5r54crX79/Dq+/\n2v2DXvu+ZNv+smH+Jf0fr99tX3P8/D7c6R+J3YHsgJmB/AQ5VL/DTeDehms3qcbm0ip3Ze0DrXd/\nu6jiexQ2Ib2R6gW3bYPYqxZf/7pu35qNQ8fAxqiyMev4aRviG3wnrgOa8RmNf2N9rIh/d34fbm//\nevX3ls+iFp9UcXxb86Owvmrxcfeur11bL3mw0bt/oby2vr+t2L+3jfO32//x+dtun79vYfyKc6aC\nxO4oduDNgKlrZBC2kHt/R33jshP23fd8YvbKg9rEX1U/vC9qZ1pfLn0P37dC3c+px9f0QH1qFrdf\ntT1c+j72sFjKS/OzFr9gW3wH2/dtzfi04h+0xs9o3D+NX31+n2dv/3r195bPoB6fofFtzI9wsL+v\nyhyzz8kWdn/99e7fXt/fNrq/bJ+//f636o+2rzJ+al68r2ycHRK7ScmgAgAAxMgOJkViBwAANLKD\nSZHYAQAAjexgUiR2AABAIzuYFIkdAADQyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0A\nANDIDiZFYgcAADSyg4MlX8vyOO97Y6+X2OmvU1lisHxlivpald7X3mTl/a8c0vHPvzbGSL6eJf3+\nv3dZ5atavmvgK5dG2v/++rs1X4njdOM78vxW/L8paZe/dD92zY/e/Paq43O2gfnV1Ku/t/wmGvN/\naP/qzY9iuY6tv1pr2F5R+0bWx2Wt6H8jvvXxGbx/a/9oxndF+w0SuwPZ73gzgQ5kEjxOOrVkUK/F\nTbxlvvpD7P2C/tm9f4nP2vIl3u8I24WxLMbSd+6FhaleNgYP3a/p9X+8/VLv8ZA4xBvX3viOPb8e\n/wtQ/UvtnR86vnp+L8rjc7aB8W/q1d9bfg/N+T+wfwl5vTU/yuW9/dlZ0765jPVfVOPbHJ/e/f3P\n5j3V/aMZ3/H2CxK7E/2txE7xv/Us805PTBefZePulbuJnMRTFkK08dv3pzcwbYgWil046TNKahvq\nd/X7HzTbHzaPbBM5IL5e9fmD8f8a3773VYit2DY/evPbq47P2cbnV1mv/t7yG7BjWp//Q+vL3qMx\nP3rlQbY/G532hfLe+phCqf+iEb/R/c+q3d/o7Y9D8W3cX5DYHc5t4GFgziL3vzK7CFobi9/I6wek\nLh/Y+O3CCH868PWjPyWEg/QlC8uOUfnPgPdN7KJ7ZBvX/vgGteePxv9rVExsP1bFdw0fv+Q+rfE5\n28D4N/Xq7y2fX3f+2zFvra8oRva9en70yhel/Xmsff31MYNS/7vx647Ponx/p7p/rIhv6/6CxO5E\nMjDJRnWgayd2ftJXDz9fXt0USuXRogtkIah7uAlv6prr+ZLfapaFF8rCPWqL45iD+2hj/Re19if9\nzTau/fEN2s/vx/8yshg5++dHef63x+ds4/OrrFd/b/n8RuZ/eE9pffXmx/j88fNPzeGR9iWaz7iy\nVv/b8QsxKo3Ponz/YHj/qMa3fX9BYneqn9+n2ZjOIBPrspofE/tPNKuTsla+fuO3CyhaGHZRxu+3\n7RxPTL5rvP/l9kefJMfXu/7++Aa1+I3G/zJs+8b7N6Y2v3vjc7b145/q1d9bPr+18z9dX735sWL+\nVPbn1euzsj4ur9j/9euvtv/V4hsM7x+1+HbuL0jsDiQbUxzsfQdAm0y673AbcGvCu99q8gnpXq9P\nyHa5f25UaONbaUfxXnZBRBuVHByXS+xq8R3v/1D7bd/jMTogvl71+YPx/wrdNmNbfOvrozf/E9n4\nnG10/LfOz73ls6jFx1gx/7tzpTc/GuXu3oWyVvtWrI/vasTfq/Y/1olva3x697dx0xVXxHek/SR2\nh1JZvxmUs8j9v6O/cIoT1y6UKDbhCu/rlVv6t6p0ctvnVsre1HOW2/t+RWX2+vjG1Ypvq//j7U/i\nlJTvie/g86vx/76wWb+vpO2j8a2M39D8durjc7b2+Dtb56fYWz6DVnyMxvwf2r+M3vwYKq8tvEb7\n2uvjKjrxN5r9N2rxS17vjU92//7+MRrfXvsFid2kZOABAABiZAeTIrEDAAAa2cGkSOwAAIBGdjAp\nEjsAAKCRHUyKxA4AAGhkB5MisQMAABrZwaRI7AAAgEZ2MCkSOwAAoJEdTIrEDgAAaGQHkyKxAwAA\nGtnBKfxX43S+9mOPKRO70tcqlWKUvC/+/rzeVw6pr23pfa3L6vJvG/jKpRC7Vtvt9xJuqb/nK6F0\nmb+u+LVE1fgMxL/p6vE5s39ib/nkRve/vfOvWr+1v+l7+6s0/6r788UNxD/5Wq/a2qv2f2B8mvur\n14tva//2SOxOIIvn+TLBbw3eTjLo07ETtj4ZRVhYeej8phQV2O/Mey8+Xe4W2eO9KvaWf9tg/02b\n7RW9T5N6j4fEOR6LXv3B51fLXTyX4vz9V9GMT7V/PVePz9n921t+AwP7n5B+75l/Y/X1/jY2/+r7\n8wR68VflNr7D/R+c36ZueX91RuJbHt8Uid3BbNBlsUhW3RqdnWTwp5P8JmIuvSn58nLYegtHb0x+\nkRxW/m2u/0miKfEqtK+0Ib2FzUttYrFy/bGNa6R9lv+ts9bMr6nGZ2X/MlePz9n921t+A9IfM6bV\n/U/Y9+yYf9X6K/e30vyz97zgmh3l21+Lv41HGiATgyiGzf6Pz9/q/jwSX/ue0vimSOwOlEwMEruc\nmoxpYuDjZ35+yet28emPod3m5Mrkqk/ssNCKC8jaW/5pB2wc8T0aG0O9fiv+4+0Tdqyb4/cNrfis\n61/u6vE5u397y29A+hONqd7/khio947Fp1Vfc++t7W+l+dffny+uF39bHvrk4xP1sd3/kfFxavtr\nP77j40tidyA7YO9Dz13JQB9I7j09NTndZrLETG8u6UJsbUy+rLConL3l3xAt6mDLxhHi2dgYavXb\n8R9vX173GtrxWdO/kqvH5+z+7S2/IelfNMf2zr/R9R3uVY+tL1fzz91/aUPyvBkVYhT6KJf916l0\nYlft//j8be/P9fgmPzfHl8TuPHxi16c+6rYTN14Itjz9DSpeONn7Lf+pUjX2e8u/Ze/GoT9t89dw\n/V78x9vnxvVqf9LpxWdF/4quHp+z+7e3/IaS/W/v/Btd3wP7W2X+tffnCSXxz9l9cMf5VJu/zcSu\nev/x/VuQ2J1FBuXEnVkG9ZrcBM8mXGETsBO8OpENWRjRz3pB6J/twjBxqYV9b/lnVOIXXtf9Lyzs\n2saRsLEtb2q1+u34j7fPxbm+oV5CFp/R/u0fv+/E5+z+7S2fRSU+I/tfbPP88wrre3R/q86/zv58\nDcfEvxirZv/Hx8e+XhqENfEtjG+MxO4EduDMpJAryeAPJPe+psrCMsJieV+lRWUn7PKedP77e5fq\nq3rvK9xgb/nH1ONnGql+ayts/Em5uWobS1Y+Ur8Rf6vVvkV1Y7uIcnzESP+2jt/ie/E5u397y2dQ\nj8/Q/mfsm3+V+iv2t+b8a+7PV7A9/kncanOv2f/W+Ki9M1y6nQPxrc+PBYndpGRQAQAAYmQHkyKx\nAwAAGtnBpEjsAACARnYwKRI7AACgkR1MisQOAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCR\nHUyKxA4AAGhkB5MisQMAABrZwZHsd71FX/dx4lfiXDOx63zlTelrbeLvTNlbvvf5RvK1M5Wva/me\nga8UCn0sfRdN8J6nuv7YVxalcUy/y7AZv4H4f9XZ82O0/9XxOdvg+Ff16u8tn12vf+prp5K5oev6\nK5mDg/Frza+R/aOx/q9t7/wbjG8xPkeMX2t+pEjsjiQLprUgDiQDey1+0kX9t99pF09cO+Eri0Hs\nKj/++c3vTPy4Xv/GF73UezwkQUljofubxc8IiU3x9r349eL/bWfPj8H+18bnXAPrp2lwfm4un93a\n/rtDfvmucffzUr1yv4H4yev5/PL36+wfzfV/aWvjn8dvJL71+PTGb2379PxIkdgdicSuPfHtwWZe\nC5daFPvK3fOTiS7vX/F8uyjj8bO/2X7ycG0Z6J+nN6CEjYHpU/inf7mkNn61W3fj14n/13Xat3t+\njPQ/jEv4p3/5fOPzq6xXf2/57Fz/6vujPvhNDZlvtf7buddOsLL1K8K8Cv/0L8eq+4et037mdfXi\n3yvPZeVr4pONX2/+r5sfJHZH8oO1bN7nbcxy/+txk6/af7WZlBfG1vKBg2Ho/uGjc7/QL/OnhoH+\nebZfxd0luoeKxSIew7Q8bCQvub8tV7Hpxa8X/2/rta/Xv55u/0fG5yzj86usV39v+R109seEi0d5\nHZtSWYvF+vX1m8S4Mb9q+0d3/V9eL/4j47Njf4zk47d2/rfnB4ndmWRgKoHfSybO1aQHVXviWY3N\nxVpVvnZhGIX7uwXnFubzJYn6VTav8f41N+bQ315sw/hF9w+xCW3IN6flPUPx67bhy86eH+r+68bn\naBvWT6JXf2/5/Mb3R19W7XurbpDfY3R+tfeP9vq/sl78151ftfiOxKd0b/fa2PzPn62R2J3KZPcm\n+GeQCXQt+cS0E7sx+cwbzMHY2BhWla9ZGF7n/naht57/UeP9K2/M+rdRf3XjsyQu2Xiqcq0bv974\nftvZ8yO5/4bxOdSG9ZPo1d9bPru8f+X90c+DakJh2HmT/lmuKFmf4/OrvH8U2pvc/+p68R8dn4jq\n/3B8iuM3Ov8H5odBYneg4sTorr5tZFF+h5uApQmvN4Tk58Ikt+XhPnvLQ7v084frp+zYZYvvE2rx\n7fQvYl/vNVw2jTgpKcSnu1HZe6R1gix+K+P/cYfNj8r4re2/Hp/Tjc6vrfNzb/ksavHx/dH9i34e\n3XPc+9TcGFm/scb80u16W7H+v2d7/Js/j8R3MD7F8Qvt1s+P7j86PwSJ3cHs4Jjg28sMylnk/t9R\nXzjvsqj/8XvCxDyr3CwF9Vtpunh69e1CqtT9nFZ8W/1TsQ9X4T5JP6PytP9yFWJgN6vlPfEm04tf\nf/y+65j5UR+/0f7Xxud87fXjbJ2fYm/5DFrxUWs0fo9aV++rcIrb+VF7Palfjl95fg3uH431fw0b\n42+1y4fiOxCf2viZkvr8XzE/BIndpGRQAQAAYmQHkyKxAwAAGtnBpEjsAACARnYwKRI7AACgkR1M\nisQOAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyKxA4AAGhkB5MisQMAABrZwRmSr/84\n52txrpnYDXwlUIhN9lUouq6/kq91UV/50r1HJfbJ+BS+69B+51+j/tfsiW+k1b9W/SRu/iq9r3b/\n0frf0mzfyPxsWVH/a/NvcP1U9ervLZ/ZyPj3+598LV1t7jXmT7t+5/lXX79dg/Ort/42xXdk/Hvj\n4w3sDyR2B3MDUw/4UWTgr8UnXdFCt9+J956cY0nZ8rK+n/7Zvf8RZWX6O/jS5zth4WSPj0i9x+Mz\n4zhub3wX5f4N1Lcbez8m1fgN1v+aZvt687NnvP535l9vfvX06u8tn93a/a3QfzU/bXn0/kBeH1l/\nun7x58bz5zI+v6rx86rlzfgOrP+94xshsTuSZNJmonyCHL7Xsm7hlCZswv9WsrxNLwx5i5nchfsH\n2fPtwknvkQmL63KbmItvnMjatq6N70D/qvVtHVMWrlLsW/cfqf9Na9qXzc+VavVb8TvV+Pwq69Xf\nW34z2fi7/rf2T7vfxRPG3qO0xsrzZ6h+pLZ/Dq2Py+nH12rEzzoqvoX1v3d8YyR2RzLBfjzN9Z78\nhT/zHUTufz0u+VoWf3niNRMPz07yxsQ178gWqhO3Ia1v72kW8kueb8v1+ESHS2fhfN74wVeP71j/\nqvVVnXxj7Ny/W//LVrSvPz/byvW/Of/G51dZr/7e8nspj39n/7RzIuxZLl7pHhbF0L53bX1R3z/1\nPS+3frt651MnfofE1ymOf7d+r30LErsDucGKNifJuM3EP4M852rShe4nZiFBsO9rJnb1uo4vb24q\n+Xv0+OjFlfzcWTifFy3qQNpYiEEtvqP964+Pp+6xOn6Xi7FSbV9vfvaU6393/o3Pr7Je/b3ld1Ie\n/5H9080R87q5ni/5RGc5+EfmT6t+yj+/Ff/KM66qF99e/I6Lb3lsRav+yPMDErsD2cCbiRN7PWoL\nZx8Z+GtxkzXemEM8dPe7iUPhY+qF/62rVT+w91ELI25PUq5/m/NXa2P7qDy+dnEPx3e8f8OJXfKn\ngg3xS+pfUK199vXa/BxQrP/t+Tc+v8p69feW30hx/PP+1/bPwK7THesvrV9g29k4v2rr45J68e3F\n78D4Du4fe8aXxO5gMnmWATODYQJ/BhnU73ALpDShdEKgfw5qrwd2wRUWhHu9siAKm1A7kTPk4Kht\nXLbsG5tWLb7+dR3fgXEoavSvWL8Q39rzLX3/tfU/bUX7avPTqa+PoF3f+/j8G51fW+fn3vJZ1OKz\nqI2/7a/uf/RzrLkXis78yer39s+rr9+3evzXxLcXv9XxjdTGP7Z3fEnsDpdm1tWB2Unu/R2tjcuX\nhf4n71Flxfc4xQVnJ7KqK1f0PlsvKS9MfHWf0vgk9/n4xtWKr/6tLe7fyvhm5f36YbMplcVq8Rut\n/y2r+ldd2K3xc9r16/E7X2t+BVvnp9hbPoM946/WoLpHMi8asanNn179tDx/z9XXr9OKfzu+QS1+\nwdb4BvZ9hfFfVT+8r9IHErtJyaACAADEyA4mRWIHAAA0soNJkdgBAACN7GBSJHYAAEAjO5gUiR0A\nANDIDiZFYgcAADSyg0mR2AEAAI3sYFIkdgAAQCM7mBSJHQAA0MgOJkViBwAANLKDAyVfV+Kv5Eut\nDzRrYpd8LU3p61DCV35VvnKpW1/Y7zWU98Rfy6K/rshfyT3UV840vvbpOwa+cqkVv9LXsiXv6/W/\n8/zm/UfifxFJP+LvxzzoK682z8+z7e1fr/7e8tnt7//Q/ledv97m+T3/+Gw/P5x2/db+qWPnr/ge\n3f15PP4kdgd6PdJAyyQgsYvYibvEKP3OvIGkqll/Ia8/HrIA4/Fwi2J5u3+efv77Z/f+s8ZvPd0+\n3//3xrA+fqle/3vPN5r378X/GsLGnTdroP+Dts3Ps+3tX6/+3vLZHdD/gf2vPn+d7fP7BuOz6/ww\nmvV1fPT+ObC+m/unfp7/uRJ/ErsT6UTvSLI4Z2M3lXgi29+M8hjVFtxQ/bA4OovE1Y0Xml54/nmX\n2bjcRpAkmtLHQvtq8XOxMWXhSur2+j/w/Ob9lSz+F+DbX27TePyb7DO2zM+z7e1fr/7e8tm5/tUP\n5l65eUdv/7NzqjFnmuX723d1e8+Pdv2V50dpffvxeV+d2LbiT2J3Fhm4ZNSOJQM/HTtxw0f/fqMo\n/KnATthS7Lr13Wv2cCgszJhddI3y9/1PHMN1or4F0sfCwm7Hb+lze2PW/R94/or79+P/eWEjfkm7\n7eZamVtBJf510T1UrLTPx2dv/3r195bfgTv8l8Nbj2+n3M6Z+v7Xnr/98t3tu7pO/JI5aN+7Lv4p\nX145P4rrWz2zvH/GY1CPP4ndSWTgzswJZGBn5Ca0m5jPl/zWki+MamJitOoni6W0MN/ai+5dfqlD\nJdp0AuljoY2t+CWqMSr1f/z5b737n7lANghzK/QxmU9b+q8cNz/PsLd/vfp7y+eXHtT5GPfKRX//\nW2KYzt9++RHtu7p+/Nrrs1V/4WNTnbuDsau0wWk/g8TuJGf+GVbIxJqd3SgKE9e+PrBhpPX1b5P+\nKk380sfgb/4+l9uw3EIeOfhG4+fioONf6//489+K9zea8f8eu2nH/bHtDBv3hv4njpqfZ9nbv179\nveWzy/uXzrdeeU7vn+352ys/vn1Xt/n88HT8nYHzY3R92/cV9s9AjW+MxO4MEvCTd2WZdNfkNoDe\ngrebQmVyjyQmrfqWHAqVReHq5mXde35ELX7+9ahxNk6FOBfjV9gEdP12/zvPH7h/UIv/1+k+2DmU\nHnz9+I/Nf3fvdfPzXHv716u/t3wWtfj4/uj+rfg5VlyrzflrdMqPbN/31OMf6+71jfUpSvW79/SK\n67u3fxbK7X0q/SSxO4EE/Oz5LhPomjobm2m3u/Si8fXe5f6K7tOuv0jeV2uHHiC7kKN64Tp7IDOt\njUn/VhnHoB+/sPGUysb633p+5/6RYvyvQsUhbWa7/07/YNk0Pz9ib/969feWz6AVH7VGs/e0y5N5\nU4tNc/4azfJ97buGevyH4mfU1mez/orzo7a+e/tn+ny56n0gsZuUDCwAAECM7GBSJHYAAEAjO5gU\niR0AANDIDiZFYgcAADSyg0mR2AEAAI3sYFIkdgAAQCM7mBSJHQAA0MgOJkViBwAANLKDSZHYAQAA\njexgUiR2AABAIzs4WPK1Hyd+LdBVE7vka1Gyr3UZ/Mog+714pfJWfV3mr7gNpa99Scbo6l9pNNC+\n0MfW3CvGdyB+Rnt8vSTO6fcbWr3yb6u274D50Ryfb8+/vc/v1d9bPr9d+2MyL/01vH/pMn8lbbj/\n+Ozbv/bEZyT+3gH7D4ndgeykiRaaJHmt83UPGdjLsRNymWz6O/GKPxcmtrz+eMgCjCeu/w7Aan03\n6Zfi/P2mQnMxjLbvO3r9V9/jGPdbKcd3ffx0vETYOGuP75V/W6t9++ZHf3z23X+v3vzqGZyfm8tv\noLl+Bvrf3L969fvrO22Prn+D8dm5f+2Lz8D+atSfvy7+JHYHkkGJB0R+fpz0cYQM/tXYSZkG4PdR\n3YgqEzMsvmwTcxM7iae8p7ax+E+lkgVi72meGa5aXa+1cD5vvP+23fnO4FTjqxTi1x1fH9/ao7vl\n37ayfVvnR3N8Ip+dfyvXV6ZXf2/5/NrrZ6D/8rPsW+FKYrMyfqX9UUnnn7t/3P7Pzs/9du9fSik+\nu+LffP66+JPYHUp9VGqCfha5/+XYiRk+OvYTMftTWxwjnVhEi8PeKy5ft3DsItb3V/csL4xW+75p\nvP+2X8XdoRXfVD1+9fG1dUx7XvJ8G7907Hvl3zbWvv3zoz4+4lvzb936yvXq7y2/geb6Gei/rb/M\niT2JRXF9W635F5eVyi+uGX8fExOvbet/f/z7zx+PP4ndkcxAxvv1v9ezMDGOIQN7RW7Cuon3fMlv\nJaXFIfzCiiZ+MtnVJhbeP7Zw/L2rh6eXPSOWt++7xvtfSxza8Y3V49ca31AW2pg87/1zvfzb1rVv\n+/xoJ3bBp+ff+Pwq69XfW34PYY7l62dD/6X8PT/X1Hfvbc9B/56ofimR7M/ja6nHfykLMVy3/vfH\nv/f8NfEnsTvQz1MlMWZgk4E+kEyAq7MTsXVw24+jQ8z0byP+2rJwBv7MYOmP4rWkfd823v9y4tCL\nb2Qwfnp87UYU30/Fr1f+bavbt7H95fEp+Gh81hxMJb36e8vvJ10/G/pv58eG+rbe6P4Y5l9+/2y9\nTGbt/pXpxGdt/NvPXxd/ErsDSaCTwTIDmwz0geRQ/g43wXoL2k66ePIWFklrYtpFkSRd/rlRgO3C\nLNR3z47rGoXnJ/XXtu80tfiO99++rncNLYvvohg/JRtfoWNonxH93Cv/tlb7hudHf30Ux+fr8290\nfm2dn3vLZ1GLTypfP53+F+bH1vgV1/fA/LP30/ePfr6GrfE3dq3/nfEXrefbH8fjT2J3qM/9Owhy\n/++oLxw70Rp9T8vL7xHJ+5LnjMW3NuHDYi7fe7x952ptTK3++3pJubl641QrL8QvjU8lNnYzWt6X\n3aZX/m2N9o3Nj9r49cfn+/NvZH1tnZ9ib/kM6vHpr592/3v712j8xtZ3qb6aw8U58G174m/sWv/7\n4m8198fx+JPYTUoGFgAAIEZ2MCkSOwAAoJEdTIrEDgAAaGQHkyKxAwAAGtnBpEjsAACARnYwKRI7\nAACgkR1MisQOAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyKxA4AAGhkB5MisQMAABrZ\nwaRI7AAAgEZ2MCkSOwAAoJEdTIrEDgAAaGQHkyKxAwAAGtnBpEjsAACARnYwKRI7AACgkR1MisQO\nAABoZAeTIrEDAAAa2cGkSOwAAIBGdjApEjsAAKCRHUyKxA4AAGhkB5MisQMAABrZwaRI7AAAgEZ2\nMCkSOwAAoJEdTIrEDgAAaGQHkyKxAwAAGtnBpL6S2P087XOT6/njC0f8/D6T+k/zSmT0/sn7Hr+v\nf/51Ecpa7arV392/s3XiJ0b6/+/1+yjW79//3+uxlD9ev3Ho36rjM9D+U7Wer8v8lfTx3+/rEZWt\nnhud/n99/u0dn179veUXt3v8ev1vzT9d119qjQ6t3x37w1edHv/B/m+OX698fP8hsZuUDOzH2YVT\nmcwDfp7pZLQ/x5vLwP3DxpTP6bFJX69v7OzfuXz/qvEbX/RS7/GQOKR9XTs++v2iFd/u/U/Wfr7b\nVJdiHW/9s3v/I/mtoqU3fsZX599A+5p69feWT+DU/VHHR8+/3vw1BtavkNfz/WGC8flo/Ov9L8dv\n7/17458isVNsME3A7GWC+jIDFMfOBjeURwP3PtCeMrny8iCpb+79MM/YQup/nF04cfv3Lep04hq9\n+/vyaO4X6QX01qt/cP+O5RZ2spClvYU2VvsvbB/NvAz/9C+X6PGxczy+r/3NNLpHL75KNv4f1ny+\n/6176Ys+OH08hts/MH4+ft+Zf+Pzq6xXf2/5BA4ev3R+rpx/2fz1709fSNevsH0o7Q9ufOL6316/\nmVPjP9j/avxy6+6/bvxJ7CI2kFHkbOCiPyXVNp7AvX8ZTPk5fr/+OQzcFjJxP05N1uLE7nIT1C0+\nNfE79w8T+SWv2/rqz7CeHsegW/+Q/p1l/OCr9T+5h+rrojc+IWZ+I4piODY+jft/xNjzbV+a7cs3\n4raB8fvq/BufX2W9+nvLJ3DI+I2uDxev2vwrzl/bvvr6TcZA9cWJ29Zr3xecHv9e/3vxE3vuH2uP\nP4ndmwmqmQRV8ttNIYg/z/RgS95iBjfZqFYNXJvU/7rq5B3hJ2Zr4an7u81q2fxrh28tsRmt/7ar\nf0eLNo1A2leIX7v/vj/dvpXHJ8RQrudLfuNXid1wfAfG/1St5/uyyqa5re2uzsj4vX10/m1oX6JX\nf2/5hHaNn4tHuf+tMuHLG3tgff3W94c0Uao/4zIOjn+v/734pdbff5HX1Ujs3j6R2Cky+K1nNsjC\n/LrSR/lr2Prxb4yKur9dOPFkrtRvJjYD9d/29u9QbjGPHHzl/utfKvxVqP/WiY99zobxeeuVn632\nfPt6+mePhY9jubBhfPzebDs+Nf82tC/Rq7+3fEJ7x684PwfmX3P+LtL129sf8vHJ1vvVHBr/Xv/3\n7q+j8R3bf0jsIjLRk43FBzHEUP/7drbcBD6QgUjibTam+H5Zfdm4equvQibNRxU2GbsxZBPPTdDs\n9UL9ZOKO3F+/R+Kn6ghbrxTXVv3h/p2tEr/wetSvWvuq/Y/ZvkebXm98FFsWrQ1rZXxb9z/ciue7\nvuUHQrHPmY3j9/X5Nzq/ts7PveUXNzx+lfgNzM+x+Rfe105ouvfS+4Nh+6PHp9eYT/lA/Ff1f8P+\n2rv/6PgLEjvFDroJXriSINrBKZeHoL9fswPrfg7Jnb53b/G1SP1Pi/tor2zRiMrCMexEjeur/g/d\nP4qrXMv4+OdGZfbS96jWH+3f2erxM41XvxXG8Rvsv5GMQ1TeG5+0vDJ3G/Ht3f9so8+374sbLlS/\n3pd+3+bxu8L8a7fP2d6//eXXdur+ODz//H1qr7/r1mObvC9pp9pjPj4/284+n0b7X4vfrvuvGH9B\nYjcpGVQAAIAY2cGkSOwAAIBGdjApEjsAAKD9jeyg8O/Gpddc/y6HkHYDAADEyA4mRWIHAAA0soNJ\nkdgBAACN7GBSJHYAAEAjO5gUiR0AANDIDiZFYgcAADSyg0mR2AEAAI3sYFIkdgAAQCM7mBSJHQAA\n0MgOJkViBwAANLKDSZHYAQAA7Z0dSKLAxcXFxcXFxcU1z6XxsQ8AAMBNkNgBAADcBIkdAADATZDY\nAQAA3ASJHQAAwE2Q2AEAANwEiR0AAMBNkNgBAADcBIkdAADALfz+/g/B4F8AVL1l+wAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Here we just use Jupyter's tools to import an image of the Stata output\n",
    "from IPython.display import Image\n",
    "Image(filename='stata_reg.PNG')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pySpark ML Library\n",
    "\n",
    "Next we will compare the Stata output for this simple model to what we can do in pySpark using their ML library, short for Machine Learning.  The nature of distributed computing makes solving models like this more complicated than most researchers are used to: if your data is distributed across 100 computers and there's no possible way it can all be loaded in one place *how do you invert a matrix?*\n",
    "\n",
    "There is more than one approch to parameter estimation to address this, but the one used by pySpark 2.0 is called **L-BFGS**, which stands for Limited Memory Broyden–Fletcher–Goldfarb–Shanno.  Explaining how it works is beyond the scope of this simple tutorial, but suffice to say it's an algorithm designed to accomplish parameter estimation in a memory-efficient way that our distributed system can handle.\n",
    "\n",
    "# Configuring Independent Variables\n",
    "\n",
    "Before we run the analysis, we must first set up our data for ML to work with.  pySpark ML analysis tools, such as GLM or OLS, require a particular formatting of the independent variables. The tools they offer for this are powerful and flexible, but require the use of an excessive amount of obscure-looking code to accomplish what most social scientists will be used to achieving in one or two simple lines of Stata or SAS code.\n",
    "\n",
    "The standard formatting **requires exactly two columns** from a dataframe: *a dependent variable column* usually referred to as 'label', and *an independent variable column* usually referred to as 'features'. The dependent variable is simply a column of numerical data; the column for the independent variables must be a *sparse vector*.  \n",
    "\n",
    "*Optional explanation: A sparse vector is a memory-efficient version of a normal vector where, instead of filling in every zero, zero values are automatically filled in during calculation for any missing entries - this makes it more efficient the more zeroes there are in your independent variables.  A model that includes dummy variables is therefore ideal.*\n",
    "\n",
    "In order to skip over a lot of this tedious code, I've written a wrapper that handles it all behind the scenes for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Code copied and pasted from:\n",
    "\n",
    "    https://github.com/UrbanInstitute/pyspark-tutorials/blob/master/indep_vars/build_indep_vars.py\n",
    "    \n",
    "Then deleted from view for display purposes.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function, `build_indep_vars`, takes in as arguments a pySpark dataframe, a list of column names for the independent variables, and an optional list of any independent variables that are categorical.  It then handles in the background getting the data in the proper format and expanding the categorical variable columns into multiple columns of dummy variables.  When completed it returns the original dataframe with a new column added to it named 'indep_vars' that contains the properly formatted vector.\n",
    "\n",
    "It has an additional `summarizer` argument, which defaults to `True`.  This causes `build_indep_vars` to return two things: the modified dataframe and a \"summarizer\" object that knows some things about the data, and has the ability to create a more nicely-formatted summary output (which pySpark does not have natively).\n",
    "\n",
    "Let's demonstrate on some data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = spark.read.csv('s3://ui-spark-social-science-public/data/diamonds.csv', inferSchema=True, header=True, sep=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll trim the dataset to only the columns we're interested in, just for convenience, and add a column for the logged dependent variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = df[['carat', 'clarity', 'price']]\n",
    "\n",
    "from pyspark.sql.functions import log\n",
    "df = df.withColumn('lprice', log('price'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----+-------+-----+------------------+\n",
      "|carat|clarity|price|            lprice|\n",
      "+-----+-------+-----+------------------+\n",
      "| 0.23|    SI2|  326| 5.786897381366708|\n",
      "| 0.21|    SI1|  326| 5.786897381366708|\n",
      "| 0.23|    VS1|  327|5.7899601708972535|\n",
      "| 0.29|    VS2|  334| 5.811140992976701|\n",
      "| 0.31|    SI2|  335| 5.814130531825066|\n",
      "+-----+-------+-----+------------------+\n",
      "only showing top 5 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "df.show(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we call the `build_indep_vars` function discussed above.  The only arguments that *have* to be passed in are the first two; the dataframe and the columns that are independent variables.  All the others are left in here just for explanatory purposes.\n",
    "\n",
    "If you specify anything for `categorical_vars=[...]`, those columns must also be in the independent variables argument.  If you change `keep_intermediate=True`, the modified dataframe will include two extra columns per categorical variable that were intermediate steps to get to the final result.  And if you change `summarizer=False` then you'll need to change the assignment at the beginning to be only `df = build_indep_vars(...)` instead of `df, summarizer = build_indep_vars(...)`, because it will only return the modified dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df, summarizer = build_indep_vars(df, \n",
    "                                  ['carat', 'clarity'], \n",
    "                                  categorical_vars=['clarity'], \n",
    "                                  keep_intermediate=False, \n",
    "                                  summarizer=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----+-------+-----+------------------+--------------------+\n",
      "|carat|clarity|price|            lprice|          indep_vars|\n",
      "+-----+-------+-----+------------------+--------------------+\n",
      "| 0.23|    SI2|  326| 5.786897381366708|(8,[0,3],[0.23,1.0])|\n",
      "| 0.21|    SI1|  326| 5.786897381366708|(8,[0,1],[0.21,1.0])|\n",
      "| 0.23|    VS1|  327|5.7899601708972535|(8,[0,4],[0.23,1.0])|\n",
      "| 0.29|    VS2|  334| 5.811140992976701|(8,[0,2],[0.29,1.0])|\n",
      "| 0.31|    SI2|  335| 5.814130531825066|(8,[0,3],[0.31,1.0])|\n",
      "+-----+-------+-----+------------------+--------------------+\n",
      "only showing top 5 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "df.show(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll skip any more explanation of the contents of the `indep_vars` column, other than to note that it contains all the information necessary for the calculation of the right hand side of our equation.  \n",
    "\n",
    "# Running GLM\n",
    "\n",
    "We can now move on to analysis.  There are two steps to this: First a structure is specified for the GLM equation.  Note once again that we specify the names of the columns that have our **dependent and independent** variables using the common machine learning terms, **label and features**.  Then we fit that GLM equation to our data to get a model.  Most types of analysis in pySpark follow this general format.\n",
    "\n",
    "We will test it using a stanard gaussian distribution with an identity link function, making it equivalent to OLS.  From the pySpark ML documentation, other supported distributions and link functions are:\n",
    "\n",
    "    “gaussian” -> “identity”, “log”, “inverse”\n",
    "    “binomial” -> “logit”, “probit”, “cloglog”\n",
    "    “poisson” -> “log”, “identity”, “sqrt”\n",
    "    “gamma” -> “inverse”, “identity”, “log”"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyspark.ml.regression import GeneralizedLinearRegression\n",
    "\n",
    "glm = GeneralizedLinearRegression(family='gaussian', \n",
    "                                  link='identity', \n",
    "                                  labelCol='lprice', \n",
    "                                  featuresCol='indep_vars', \n",
    "                                  fitIntercept=True)\n",
    "\n",
    "model = glm.fit(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can query the results directly from the model, without using the `summarizer` that the `build_indep_vars` function generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DenseVector([2.0808, 0.7221, 0.8172, 0.5682, 0.8555, 0.9341, 0.9195, 0.9979])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.356165273724909"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[573.8145328923426,\n",
       " 51.24311447314843,\n",
       " 57.736801306869914,\n",
       " 40.069982524765486,\n",
       " 59.53629044537795,\n",
       " 63.12762399645197,\n",
       " 60.46906625364304,\n",
       " 60.69082330168219,\n",
       " 372.05231564539275]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.summary.tValues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.0036263332330212492,\n",
       " 0.014091572711412436,\n",
       " 0.014154515091987462,\n",
       " 0.014180739640881894,\n",
       " 0.014369486052126049,\n",
       " 0.01479661778476952,\n",
       " 0.015205659565030371,\n",
       " 0.01644307646451729,\n",
       " 0.014396269149497595]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.summary.coefficientStandardErrors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice however that pySpark does not relate model parameters to variable names.  They are returned in the same order as the independent variables we passed into `build_indep_vars`, but particularly when one or more columns are expanded into dummy variables, this relationship may not be clear.\n",
    "\n",
    "In the Stata output this is clearly done for our continuous independent variable, but not for our dummies.  If we are interested in the coefficients in Stata, we can get them easily because Stata stores the `encode` info in memory: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(filename='stata_list.PNG')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In pySpark, this is where the `summarizer` comes in.  We acess nicely formatted output by passing the `model` we just created into the `summarize` method, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------\n",
      "               |   Coef    Std Err    T Stat    P Val   \n",
      "---------------------------------------------------------\n",
      "      intercept|  5.3562    0.0036   573.8145    0.0    \n",
      "          carat|  2.0808    0.0141   51.2431     0.0    \n",
      "            SI1|  0.7221    0.0142   57.7368     0.0    \n",
      "            VS2|  0.8172    0.0142    40.07      0.0    \n",
      "            SI2|  0.5682    0.0144   59.5363     0.0    \n",
      "            VS1|  0.8555    0.0148   63.1276     0.0    \n",
      "           VVS2|  0.9341    0.0152   60.4691     0.0    \n",
      "           VVS1|  0.9195    0.0164   60.6908     0.0    \n",
      "             IF|  0.9979    0.0144   372.0523    0.0    \n",
      "---------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "summarizer.summarize(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The summarizer also keeps track of the relationships between the column names and their indexes used in the vector for calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 'carat',\n",
       " 1: u'SI1',\n",
       " 2: u'VS2',\n",
       " 3: u'SI2',\n",
       " 4: u'VS1',\n",
       " 5: u'VVS2',\n",
       " 6: u'VVS1',\n",
       " 7: u'IF'}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summarizer.param_crosswalk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, please note that the `build_indep_vars` and `summarizer` features are not official parts of pySpark, and were built at Urban to serve as a temporary work-around for the lack of an easy way to accomplish nice output.  In limited testing it has worked well for parametric models.  Given that SparkR has a cleaner interface that handles this already, it is expected that pySpark will catch up shortly.\n",
    "\n",
    "If you wish to examine the code behind `build_indep_vars` it is available in the Urban Institute public repository at https://github.com/UrbanInstitute/pyspark-tutorials/tree/master/indep_vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}