Skip to content

Instantly share code, notes, and snippets.

@kforeman
Last active May 4, 2017 00:02
Show Gist options
  • Save kforeman/1078e130759522e45963ff32f1493a71 to your computer and use it in GitHub Desktop.
Save kforeman/1078e130759522e45963ff32f1493a71 to your computer and use it in GitHub Desktop.
20170501 - Loess + Gaussian Process.ipynb
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd8lFXWx7/PpPc66ZUkQIAASegQuiiWWFdZXREbKu7q\nru+ur7vsrqtucXVf17bq6q6CZa1riQqCEkooUkKAJISQBpn0PimT6c/7R4opE0iZFJL7/XzyyeR5\n7sy9k2TOuffcc39HkmUZgUAgEEw8FKM9AIFAIBCMDsIBCAQCwQRFOACBQCCYoAgHIBAIBBMU4QAE\nAoFggiIcgEAgEExQhAMQCASCCYpwAAKBQDBBGbIDkCQpVJKk3ZIk5UiSlC1J0sMW2kiSJL0oSVK+\nJEmnJElKGGq/AoFAIBgatlZ4DSPwP7IsH5ckyQ1IlyTpW1mWT3dpsxaIaf+aD7za/v2C+Pr6yhER\nEVYYokAgEEwM0tPTa2RZVvan7ZAdgCzL5UB5++MmSZJygGCgqwO4FnhbbtOd+F6SJE9JkgLbn9sn\nERERHDt2bKhDFAgEggmDJEnn+9vWqnsAkiRFAPHA4R63ggFVl59L2q9Zeo2NkiQdkyTpWHV1tTWH\nJxAIBIIuWM0BSJLkCvwX+Lksy409b1t4ikUVOlmWX5dleY4sy3OUyn6tYgQCgUAwCKziACRJsqPN\n+L8ny/KnFpqUAKFdfg4ByqzRt0AgEAgGhzWygCTg30COLMvP9dEsBVjfng20AFBfLP4vEAgEguHF\nGllAi4HbgUxJkk60X/sNEAYgy/JrwDbgSiAf0AB3WqFfgUAgEAwBa2QB7cdyjL9rGxl4cKh9CQQC\ngcB6iJPAAoFAMEERDkAgEAgmKMIBCAQCwRhCazChNZhGpC/hAAQCgWAMkV3WiMls8ZiU1REOQCAQ\nCMYI52tbqG/Rj1h/wgEIBALBGKBRa6CgunlE+xQOQCAQCEYZk1kmq1SN2Tyy/QoHIBAIBKPM2com\nNLqR2fjtinAAAoFAMIpUN+korW8dlb6FA7iE2LJlC1u2bBntYQgEAiuhM5o4Xd5TPHnkEA5AIBAI\nRonTZY0YjCMc+O+CcAACgUAwCqjqNNQ2j1zKpyWEAxAIBIIRpllnJK+qabSHIRyAQCAQjCTmUUr5\ntIRwAAKBQDCC5Fc306w1jvYwAOEABAKBYMSoa9FTXKsZ7WF0Yq2awG9KklQlSVJWH/eXS5KkliTp\nRPvX763Rr0AgEFwqGExmTpeNXsqnJaxREhJgC/Ay8PYF2qTJsny1lfoTCASCS4rciqYRk3nuL1ZZ\nAciyvA+os8ZrCQQCwXijslFLhVo72sPoxUjuASyUJOmkJEnbJUmaPoL9CgQCwaihNZjIGcXTvhfC\nWiGgi3EcCJdluVmSpCuBz4EYSw0lSdoIbAQICwsboeGNXTqkHzZs2DCq4xAIBIMjp7wRo2lkCrwM\nlBFZAciy3CjLcnP7422AnSRJvn20fV2W5TmyLM9RKpUjMTyBQCAYFkobWkf9tO+FGBEHIElSgCRJ\nUvvjee391o5E3wKBQDAatOpNnK0c/dO+F8IqISBJkt4HlgO+kiSVAI8DdgCyLL8G3AQ8IEmSEWgF\n1smyPDbXRJcAIiwkEIxtZFkmu0yNaYyGfjqwigOQZfnHF7n/Mm1pogKBQDDuKa7T0KAxjPYwLoo4\nCXyJI2oECARji6ZRqO07WEYqC0gwCDIzM0lNTaWxsRGNRkNtbS2hoaGjPSyBQNAHZrNMdlnjkITe\nDCYzFWotUX6u1htYH4gVwBilw/jHxcWRnJzM2rVrycnJQaVSjfbQBAJBHxTWDE3ozWSWeX1fIbf9\n6zAtuuEXjBMrgDFGRzhHo9GQnJzM3r17AYiMjCQ+Pp7MzMxRWwVcaPNZbEwLJjoNGj3nhyD0Jssy\n73x/ngxVA79eOxUXh+E3z8IBjFGqq6u7HYTLzMzk5MmTHDt2DICYmJheIaKkpCTi4uJGa8gCwYTF\n2C70NpTcxs8yStmfX8NVcYH8ZEG49QZ3AYQDGKMolUqKi4sBUKlUqNXqzpl/XFwcu3fvprm5mbi4\nOHx8fFixYgUpKSmd9wUCwciRV9WMRj94obedpyvYllXBsslKrpsdZMWRXRixBzBGSUpKIiUlherq\nanJzc4mPj+f8+fNMmTIFpVKJg4MDLi4uKJVKFAoFkZGRJCcnk5aWNtpDFwgmFDXNOkrrWwf9/EMF\ntXx0rITEMC9umxdG+5nZEUE4gDFKXFwcK1euJDMzk2PHjnHq1CliY2M7VwFmsxmTqfuMIywsjOrq\n6tEYrkAwITGYzEMSejtV0sBbB4uYGuDGPUmRKBQjZ/xBhIDGNB1OAOCqq67q3BAGUCgU2NjYdGtf\nXFyM0E8SCEaOM+VN6AyDy/nMq2ritb2FhHo789MV0djZjPx8XDiAYcRamTExMTGkpKSg1+vx8fGh\ntrYWnU5HS0sL1dXV+Pj4UFRUREpKSqfDGO4xdUVsRgsmIhVqLZWNg9P4L6nX8FJqPl4udjy8MgZH\nO5uLP2kYEA7gEiA0NJTExESeffZZSkpKaGlpwcXFBVtbW/bt24eNjQ1arZaVK1eOuOHtel5BbEYL\nJgpag4kzFYML/dQ06/j7d3nY2yj4xerJuDvZWXl0/Uc4gGHC2rPiuLg4YmJiMBqNxMfHdxrbzZs3\nExsby6ZNm0Z8TABpaWm9ziskJyezfft24QAE45bTg9T4b2w18Ny3ZzGYzPzv5VPxdXUYhtH1H+EA\nBsHFwiiDnRVbMtBdycvLIz4+vjPO3/Vw2MXGN1wz9Z7nFUBsRgvGN6o6DXWD0Phv1Zt4flceDRoD\nj1w2mWAvp2EY3cAQWUDDQMeseCApmpakH1JTU7tJPzQ2NuLj49PteY888ghRUVHDMqb+0PW8Qgdi\nM1owXmnSGsirGrjGv8Fk5h978imtb+WB5VFEj4DOT38QDmAYGMysuC8D7ePj0zmTd3d3p7a2ex2d\n/hrb4Zqpdz2vYDabOzeje65eBIJLHbNZJqt04EJvZrPMG2mFnKlo4s7FEcQFe1ywvbL0O+x3/i9D\nOlbcT0QIaBgYzKy4PwY6JiaGjIyMzj2A/mb+DHZM/aEjfPTss8/S2Ng4apvRAsFwk1/dPGCBNlmW\neffweY4XN7BubigLJvlcsL1/8ddMP/wr5MBZoG8Bh+FdKQgHMED6s5HaMSvuSNvsj6Huj4HuOASW\nmZk5YGM7mDH1l67nFTZs2MCWLVtIT08XwnCCcUNdi57iQQi9fXmqnH15NVw5I4DVsf4XbBt47lOm\nHf0NDT4J2K/7BJdhNv5gvZKQbwJXA1WyLM+wcF8CXgCuBDTABlmWj1uj75Gkvxupg5kV99dAh4aG\ndjqCgRhYMVMXCAaHoV3obaDsPVtNyskyFkf5cH188AXbBhd8QGz676n1X8TJxa8w38FtsMMdENZa\nAWyhreTj233cXwvEtH/NB15t/35JMZCUx56z4gvRkbWzcuXKYTXQAxnTSCAkpAWXArkVTWgNAxN6\nyyiu593D54kL9uD2heEX1PcJPbuFKSf+THXgcjIXvYTZZuRSQ61VE3ifJEkRF2hyLfB2eyH47yVJ\n8pQkKVCW5XJr9D9SDHfKozUM9MVCVMLYCgT9p6pRS4V6YKd986qaeD2tkAgfF+5fOglbRd+5NhE5\nrxGd+RyVIZeTNf//kG3shzrkATFSewDBQNdSViXt13o5AEmSNgIbgV7GdiS40Kx0oBup1jK2RpOZ\nZp2Reo0evdGMySwjyzI55Y3YKCTsbBTY2yoozD3N4QP7LIao0tPTrTomgWC8ozOayKkYWMpnWUMr\nL6Xm4+1iz0Mro3HoS+JBlpmU/SKTTv+D8rBrOD3vr8iKkd+SHakeLa1/LOY4ybL8OvA6wJw5c4Y/\nD2oAWGMjtaeDsVT3NyQkhNpmHbUteupa9LTojMgylPWQnO0pQZvy1bfMTVpJdvr3tNRp+DBlBxq1\nhl279+DpPrwxReFYBOONnPImDMb+53zWa/Q8/10edjYKfr5qMm6OfUg8yDLRp54hIvfflEbeRE7i\nU6AY31pAJUDXOoYhQNkI9W01rL2R2nNTed6iJB77zWbO1WoInd/Q2U6WZWpb9JQbndHJNmjNNhiw\noeV4CQqFhK1CQpVzgrz9B5m1+ib0JjCZTbToTJglR04cSkfXqsGgbaGiVs2aVctJmD2r8/VHIxYv\nBOQEY5nShlZqmnT9bq/RG3n+uzw0BiOPrpmK0q2POL5sZkrGU4Tmv4cq+jZy438HUvcQkaez3YiJ\nw42UA0gBfipJ0ge0bf6qx2L8vz9GyZobqR2bytu/TUVV30qY2Z3AyKkUnc2htL6VEyUNnKlopLhW\nQ4veBPxQKUhCJi+7ElPnYRE/muxC+P37e7HznISbwoCnQoex8AjGikpmJ8zF3dOL4LiF/PujL8mr\nambp/AT83ByH9B4GgxCQE4xlWvUmzlb2P/RjMJl5eXc+FY1afr4qhjAfZ8sNzSZi039HcNEnnJty\nN/kzH4Uem8P+7o5MD3IfsboA1koDfR9YDvhKklQCPA7YAciy/BqwjbYU0Hza0kDvtEa/1mS4jVLH\nLLsr5RWVqBUeFFU3A1DbrEPlFMWphjwe/zIbgFAvJxLDvQjzdqb01EEcFUYcJRO2kszl168D2v4B\nv/r8U0pjfDAaC2hqUaN3C6S8tpGy40dwnXkZavt4AnUt+JncSViymgNpe/AJieaMTRPlai3eLoNT\nJBzM6kEIyAnGKrIsc7pcjamfQm9ms8y/9hdxtrKZjUmTiA10t9hOMhuZduR/CSz+ksJpD1I4/aFe\nxj/E24kp/m4jWhHMWllAP77IfRl40Bp9DRcDMUpDnfnLssy5mhaacObUmXwqjU7k6T355NNM9A0G\nXFyc+cmCcGaFeODp/ENWwI7TlrMR7GwUOCuMxAT5Mnn6LD588x8YW5qJcHHF0V3ixluvZ8+xbMqN\nLry6twBbyYx9Zj6RS9VMC3SnrllHfYuOUyUNRPq69B277Af9WUUJATnBWKW4TkN9i6FfbWVZ5v2j\nxaSfr+fmOSHMi/S22E4y6Zlx+H/wL9lBftwjnIu9v1ebCF+XUdEHEieB2xkpo6Q1mCitbyWksgmn\n0On86dV3aA6YjauHLUkBZgpPfUN0TBTLJitxtLPB09kON0dbnO1tyfdzxUYhYdM+Q1g9zR9ZljHL\nUBTghsksc92qhTSV5KIzmdEaTBzct5doZy0apyrMMkQuupI9x7I5YHbkhV15eDjZEWzyJspOzTtb\ntwJw+x13MEnpiqvDwP49+ruKGkg2lTgrIBgpmrQGCtpX4/1he1YFu3OruXyaP2umBVhsozDpiDv4\nM5Tle8id/RtUkzf0ahPt50qEr8tghz0khANoZ7hVLWVZprpZR02zjnqDPX/beZbcSgnH4Bn4Fu3B\nWVuD5+JlbHzwIRbPS8Df3bGXAba0MSRJEjZS2yrAzgaUbg54d9EYV8yJozgjDVmrwcbJHQ9jPX71\nWfzu7htocg5i277DHE4/zCFNA0o3B2ZFBVHVqKO6SUeghxOTlC793pDq7ypqOGUpBILBYDbLZJf1\nX+jtQH4Nn2aUMj/SmxsTQyy2URg1zDqwCZ/Kg+QkPklp1Lpu9yUJJvu7Eerdx57BCCAcQDvDaZTU\nrQYKq1vILy4jvaCcqiY9zu7ZXHnZCm68LZnUFC2ujrb89L578OlHgYies2FLqaQdchHhYWGd1cTU\n6kZ8HGHlyhV4h0bz9j+ew6YwjzXTZ1PhvIj82lZSzxxF+/F3rL9qKbIMlU1aJvm6oC4tsFiroOu1\nrKws7rvvvm5js7SK6m82lcgUEowUhTXNNGv7J/SWWapm66FzTAt0585FESgsxOxtDM3MTtuIZ+1x\nsuc9TXnEDd3vKySmB7nj5z7ySRhdEQ6gneHQyjGazORXN1NSp+FwUR0ZedU4RC9mjo8DNy+M5uTB\nXbRWBhLl54qjnU2/jH9PLIVdNm/eDPwgHmcpc8lgMvNKtYoZs2bi4OJBCLXMcLThuG0shw4e4Ize\nizXTArgyLoCdaUfJST/I5NhpBPn7sWLFCl577bXO1+7od/fu3ezcubPb+PpaRV0sm0pkCglGigaN\nnvP9FHo7V9PCq3sLCPFyZtPyKGwtFHK31TcSv+9u3OqzyJr/f1SGXdXtvpO9DXEhHrgPYa/NWggH\n0AVrpng2aPRklqopq2/lne/Pk15QQeDk2SwO0ONh28KShOkkxSjZvWsnzs79WwJaGpOlsEtHlbAO\nB2AJOxsFRq2G2dFhGMwyDRoDDRoDC/whUOeITZgXX2eWc6iwlvDKNK6/5mqyjh/CscVAREQELi5t\nMcuu1cnuuusu3nzzTebNmzfkVZTIFBKMBEaTmeyyxn5J71c36XghNQ93R1seXmW5kLudro74vXfh\n2phH5qIXqQ6+rNt9P3cHYgPdsbPgOEYD4QB6YI3NxqKaFgqrmzlaVMfWQ+cxyzJKfRlrAmLwcnXA\nz92VaUHumEwufPRBNY8//vjFX9RkBF0XRUKFLTi4Wdy89vHxobGxu3qhpbBReXk577//PoGBgcTE\nxBATHEJhSTkRIQFINUdwcnIkzy6G3ScKME1uJcxoQ2WjluPF9egMRmx75CqvWbOGL7/8clBy1T0R\nmUKCkeBsZTOt+osLvTXrjLyQmofJLPPQyhg8LBRyt2+tJmHvBpxaijm5+FVqA5d23lMo2uL9IV6j\nF++3hHAAVmLLli2YzGYSVl9PeYOWD4+p2Hu2mkhfFzYmTWJrtj0eNlpCvb06n9MrPCLLUJsPJceg\n6jRUn4G6ImipBm1D704lG5QnFBT/8TNWOTrRZKuEw3pWTPahYGpMn/HzjvDK4sWLKS8vJzg4mOzs\nbOrr66moqOCxxx7ju32HsGnS4WPMpMndhszcQnLcwpjtWENds55yta7Xacfi4mJmzJjRuaK5mDO9\n0H1RalIw3FQ36ShraL1oO4PJzCt78qlp0vHIZZMJ8uxdy9dBU0HC3jtwaK3kxJLXqfdf2HnPuT3k\nM5T06uFiXDqA0Ugd1BpMqOo0uFQ08fLufErqW7liegDXJwQRpXQjLz6OM9mZONvbdoZHNm/ezMwp\nkXDifcj/Fs7th+ZKAEzYoLbzxztmLrj4gYsvOHrQKatkNkBrPUneZ0k5lsNS3ypiXfMo+s9uXjum\nA0nBI/ODcJoeiX2QmpQvPwJZJm7mzG7hFS8vL/Ly8jh//jzff/89kyZNYv/+/bTU1xIdFExlo5aE\nqEBMtuc40aThiDmU+o/34aZuQtPSArb2TA4L7Bbu6RCeGwoiU0gwnOiNZnLKL67xb5Zl3jpwjrOV\nzdybFMlk/96aWo4tJSTsWY+9rp6Mpf9G7ZvYec/P3YFpge4W9wrGAuPSAYw0Nc06impaqNQ78MXX\nOZhlmZ+vimFupDczgjzwcLYjLCwUSWqbfTer69Hm7+dO/yxWyt/A5yZwDYDIpRCxBEIX8O7XB5Al\nGzbcvOGCfcetAjIz2zavCxpZkzQPW90p7lkaiq1qP766XJxOHCO53sz2J/9N3O03UJ1VTdj6W4Af\nNooNBgOurq78+Mc/7rWR7DEjBveQyZS/9hKGJgMnHYLwjprFjQnBHEnZyuFDh9C0/hDusYYDEAVs\nBMNJTnkj+n4IvX2WUcqRc3XcmBDM/Mje5Rydms6RuGc9NqZWji/fSqP3TKAtxTPGz61vWYgxgnAA\nQ0RVp+FsZRP5WjeO65T4u9sy25RLY2Yh81be322zZ7afzM1OVUS2HMde1tHk7MNp52XE3fxbCE7o\ndjRclr7v9xh6bl4/8cQThN3yW9555x2QZTYkLyUsbxfV//cy5KSgPF9N8e92sjpgCuecZ7M3t56I\niEnY2dl1FqTvupHs6mDLrVcsoTQ/m7pmPbFLk/lnWgHv5+uInZ3MbPs6Eq+4hUmhnt3GNdSV2Fgr\nYCMYH5Q1tFLdD6G3PblVbM+qYNlkJVdM733Qy0WdR8LeDUiyifTlb9PsGQuAg52CuODup/jHKhPa\nAfTHQF2oTX5VE0XVLXx1qpx0nR/+Ni385sp4cvbl4e/u2Gb8TUbI/Zq15c/jryvCINlzzmU2MT96\ngv+m5oIkEReS2O11h5r/3i1+LkngE0VxowLlCiPcfx9JiR+R8ukWVlZls9DlNB/n67BvCGDZzCQk\nuW1DzMfHh5KSkm7jMDbWEq4MICbAld9dNY3/HCnmYAHUmxxY2qzj6Lk6ZoV4XmBkAsHoojWYyO2H\n0NupkgbeO1LMzGAPbp0X1kufx7X+NAn77kSW7Ehf/i4tHtEAeLnYMyPYHQfb0ZF3HigT2gEMhfyq\nJgqrW/jgqIrUM1WE2zYy17kKV00Z2ccO8n2jGs3ZvSSZDhLnVIGzrQ+Hva8n33UBBoUj2kY3Unfv\n7mXkrZH/fsH4ucKGuMt+DAEzePaZZzA2lFKnPct93hqut/+SVtUeMt9I58P/ZJJ7thAPDw/i4+NZ\nu3YtmzdvJhaYH+lDTnkjdy6KoLXkDCd0vvx52xl+uiIasyzTqDVYLcdZzPwF1iS7rPGiQm/nalv4\n575CQr2c2bh0EjY9st3ca08Sv+9ujHauHF+2lVa3cKBNzydK6TKiYm5DRTiAQZBf1UxBVQtvHiji\ncFEdiW5NNB35jMMNVZw/nsZcPx33RmfTUFtNSrkf3PZHjte5ILfrfqtUKtRqtUUjb4389/7Ez+Pi\n4li5alXneD7JOU2TwhVl7VEO7XoTqdTE7y6fiTlmKvszc2lubu4MC9nbKpgV6klxrYYYRzUeCh3H\ndBH8aVsOD66IoqROQ5Cn0yWxBBZMHFR1Gupb9BdsU9Os46XUfFwdLOf6e1QfIz7tXvQO3hxf/jZa\nl2BsbSSmBbmPirT6UBEO4AJYCsW4BERSWN3MloPnOFxUR5KPBq+GfPwWJKA6uY/bQovIyMwm2yaI\nJQ9+QLJ9LNu/+QZn5x/2AvLy8njiiScsGvnB5r/3nCn3J37e9Xpm+0byyZMyN175AN57P+WukCJs\nda8yI2wKaV9/gI9fVLfzBWE+zvzm4Qc4VaLmioZWXtiVx/Pf5ZFo74pc34xZlvF2GbkC1wJBX7To\njORXXVjorUVn5IVdeeiNZv5n7dReuf5elQeZvf8BtM4BHF+2FZ1zAK6OtswM8cDZ/tI0pWMzN2kM\n0DUUk5yczNq1a/k4ZTvf7j/KO4fOc6iwlutmB+GnKeKBm6/kFrtUgipSWeyrJnDBjbzesAimXkVY\neHg3471hwwaioqL6NPKjlf/e4TAiIyPZ/Me/YRu9jH86/Ywsj1XMdSii+pv/Y1LeGyhdFGzZsqVz\nb8TT2Z55kd6E+7jw2BVTmaR04bA2ADlmOeUNWmqa+19VSSAYDmS5TejNZO479GMwmfnHnnyqmnQ8\nuCKqV66/T/keZqdtROMaRvqK99A5B+Dv7sjcCO9L1viDcAB90hGKUSqVKBQKbNz9iZ6znC2ffkNa\nfg1XxwVy85xQXKpOMnnHOqKaj6L3mULxtZ/SHLaKxqa22YYl430hI98Rv6+ursZsNnfG7zvE14Yb\nd3d3iouLiYmJ4VBmHjuMC3jVYRNS+EKOnzjFg+67mF/7CY6mHzbSHO1smBPhRYi3M79YPZm5EV58\ncryETJ03FQ3ai868BILhpLCmhcbWvjX+ZVlmy8G2XP+7FkUwNaB7URdlyU5mHXiQFo8Yji9/G4OT\nL9F+rsSFePTaH7jUsFZFsCuAFwAb4F+yLD/d4/4G4FmgtP3Sy7Is/8safQ+Wi2XadA3F1Gv0nK1s\n4sO00+QUlXDzZf7cOtuL+Sc3c+LcpxQHTCMz8HYcbSBlx24qKipwdXXt83DUhTZpRzv/PSYmhpSU\nFBwdHZk6dSoHDx4kPz+f5ORknCLn4eR6ipimA0Q3H4H9JliwCWztsbNREB/qyenyRu5dMglt5S4O\nHjjEWV0lOl0rl61cxjXLF1xSG2SCS58GjZ5zNS0XbPPZiVIOF9VxQ3ww8yd1z/X3L/6K6Yd/RaP3\nTE4kvQFOHswK9sB3EMKNY5EhOwBJkmyAfwCX0Vb8/agkSSmyLJ/u0fRDWZZ/OtT+LkZ/Uij7k2nT\nMUtv1hooV7fSXFjLiRoTvq4O3DNZw5xd96GoLyTphntJUYejbzQTHOxDYGAg7777Lm5ubmzfvt3i\n4aiLGfnRzH8PDQ3tlI9ubGzE3d2dm2++md/97nds2bKFQ8SS7bGCOXVfEPbd43DiPVj7DEStQNEu\ncVt49jShrYVMmxxNsfMVqJxNbN/xHUaTzPWrFvZyAqLoi2A4MJjMZJVeWOht79lqtmVWsDTGl7Uz\nuuf6Bxb9l2lHf0ODci4nlryGg6sHs0M9L+mQT0+s8U7mAfmyLBcCtBd+vxbo6QCGnf6mUPYn0yYp\nKYmPP/2Mggo1GudAPtj2PTYF+7l3Sh1zv7sZydEd1qcQF5n0w0ncxkbWrFlDcnIyoaGhnQatL6c0\nVg85dR2bJRrt/Ej1v5cNi4Jg+6PwznUw61a44i9ITp6U5GRw803X4/zNLtz0ak6ovdG7xpK6Zy9R\nU6YRF+wxYkWvBROX3IomtIa+hd4yS9W8d/g8M4LcuW1+eLeJSUj+e0w9/gS1/ks4ufgfeHt5MiNo\n7Eo6DBZrOIBgQNXl5xJgvoV2N0qStBQ4C/xClmWVhTZDor8plP3JtJkcO42M4nq+3PMKp+vzcLKT\nWWY6ijarhFcbo0na9HKb8efCM/YLOaVLnslr2uQr0v4Gac9B4R5Ifonq6mqS4mM59v0BJHUd0+Nm\n8tHR8xRl53O1uk1862JOQKwKBEOhXN1KhdpyDW1oq/37Wruu//3LorrF8sNy32TyyaepDlrJqYUv\nEubnSbSf67gMX1rDAVj6rfRcdH0JvC/Lsk6SpPuBrYDFKaYkSRuBjUAvI30x+ptCebFMG7NZJrNE\njWdwFA1xt+BTWcBtjW+zPqiSlpBVuFz/PClfbwMXn4vG5i/klDZt2tTn80bD8Fnq86LjsHOElb+F\nKWvhswfgvRtRVsRTnHc9vq4OKCSYNSOA5upSPjrhwBv7irh3aSQmWWZWiOclv4kmGHu06k2cqej7\ntG+9Rs8Us7/mAAAgAElEQVRLqXk429vw0Mrobrn+EadfJTrr71SGrCV7wbNMCfYZcxLO1sQa65kS\noGvlkRCgrGsDWZZrZVnuyAd8A+iufdC97euyLM+RZXnOQFMf+5tCebFMmzMVTdS16PnX/iIcZB13\n1b/C7UHnqIq6kWO+NxEZHUNycjJpaWkXHdN407VXqVSkpqby+eef88orr5CZmdl2IzgR7tsHix8m\nSXGclF9fjrE0A08nO1z1dShKjnPt5atIL67nX2lFVDfqyCiux2jqZxHWLnRNQxUIuiLLMlll6j5P\n+2oNJl7clYdGb+KhlTE/HFaUZaIynyM66++Uh1/L6UX/x4xw5bg2/mCdFcBRIEaSpEjasnzWAbd2\nbSBJUqAsy+XtPyYDOVbotxf9lRC+0Casqk5DWUMrKSfKqCg5xzaPZ3mloZKQJ/7DvvT6ztforxEf\nT7r2KpWKnJwc4uPjLe+x2DnCZU8SN+UqeOkn7Dv4BifNMWhbW/nRNVdg9grFJ7OCj9NLUEhF3L0k\nkhOqBkxmMzaK8RVbFYwOhTUtqDWWUz7NZpnX9xVS0tDKQytjfijGLsvEnPwL4We3UDLpZs7OfZK4\nEO9e9S7GI0N2ALIsGyVJ+imwg7Y00DdlWc6WJOlJ4JgsyynAQ5IkJQNGoA7YMNR+LTGQFEpLcfv6\nlrZ0zxOqBk5kneRLl7/ibWrEc+lGSl1nAns7n99fIz6auvZDCSNZem5eXh7x8fHdykBalKkIm0/c\nk0dwfymZBzWnwHMvzFxHjewKtGms//d4KZIEdy+ORFWnIczbZdBjFQjg4imfHxxTcapUzW3zw4gL\n9mi7KJuZcvxJQgv+Q3HMevITNjMr1GtQ9bkvRaySzyTL8jZgW49rv+/y+NfAr63R18UYbHaN1mAi\ns1RNTZOOXQcO8qnDU3gqjDTc9DGr9O4XNeJ99TXaef3WpLGxER+f7nnSfa6EnLzYrbyLqU37WVCY\nAq8twfemN5kVkgC0FT/7NKOUxrICDMf3YtJraG3VsHTp0kvydyMYXS5W2/e7nEpSz1Rx2TR/Vkzx\na7toNjHt2GaCzn3Kuan3UjTrV8wK9Zwwxh+EFhDQtjQ8VaJGozfyxZ5DvCk9ibudzPlrPiJmygJ8\n29sN1oiP5ZTPgeDu7k5tbW23lc8FV0KSxBn3JBb86CH4eANsvRqfNX9k5vQ7Aag6f5ZtO/cRM3UJ\ni/3MTIpfwq5du0bgnQguFfqbDXah2r4nVA18eExFfKgnP0oIAUAyG5h++FECVF9TMP1nFMf9jNmh\nXni5TCwBQxF4BU6XN9LYamDPkeP8tXkznrZGcta8Q+T0eZ1tOoz4ddddx6ZNmybkLDUmJoaMjIyB\ny1QEzuI9140UO8TCN4/hu2MTM/1s8Wku5Jqrr6HSNZpTBj9sPfyZtmBlvzbXBYIOapr7ru17vraF\n19MKCfd25p4lkSgUEpJJT9yhnxOg+pq8mb9CNfMhEsK8J5zxB7ECoKZJR4VaS0HBWe4/9wt8bDRk\nrnyHKdMX9Dr0cSnP3q1BR/nIzMzMAa+EDAonUv3uYkN0A6Q+hW/laZya5nDVzbdT/uFuzhq8+OxE\nKdfODORMUQkJ0ycP99sRjAMMpr5r+9a16HmxXdr5ZytjcLCzQWHUMvPQz/At30tu/G8pn3oH8WFe\nvZQ/JwoT2gFcdeOPOVXSQHNNKZcd24ifQs3JZVtQxszHw3li/kNcjNDQ0E5HMGCHKCnYUuBNoN/9\nXN70IWEF7+JfNJ1ZTnqMKNiWCbq6ClobNXz+9Q4wtA6qIppg4pBb0YTO0DuVWGsw8WJqHjqjiV9f\nEYuHkx0Ko4bZ++/Hq+owpxOfonrKj0kM88TNSsWLLkUmrANo0RnJKlODpoaZu+/An1r2z/8nHmHz\niFK6iJOow0i50xS4dR9JLTeQ+o9fsCQ+Cdn5cpwcHfjkvfeJDXAnLiqGqLAA1q5dO+CKaIKJQWWj\n1uJpX5NZ5rV9BZS1p3sGezlhY2gmPu0ePGpPkD3vr9RF30BiuBeuDhPWBAIT1AEYTGZOqhqQWhuY\nvON2lKZy/jv1eQIjFzEj2GNcHvkeDXpqIZWUlNDS0oKLi0vbzP66V1mpfJ607Z9g1mYRM/ce8j1c\nUAUsJsRNora5CRt3/wFXRBOMf3RGy6d9ZVnm/SPFZJU2cvuCcGYEe2CrayA+7W7c6nPIXPB31JOu\nIjHMC5cJbvxhAjqAjuIQWk0TU3fdhY/uPC/6PUX8rMuIVrqJfwor02H8fX19MRgMREREUFRURGRk\nJKn79rNy5aPE4cbG2o+QXD6mOWgOZVOiOFbZjC0yZyubmBqg7Ew17Y/aq2D8c6a8CYOxd+jn25xK\n9pyt5vLp/iybrMROW0fCvg24NBZwavHLNIdfRmL4+FL0HAoT5rfQEdJZcuVN1DY0MSNtE75NWfzW\n/lGSkq7Hy8WOMJ+2k4HCyFiPDi2kxx9/nISEBJRKJd7e3hQVFXXO7J1d51FvF0Ry60fEV3/OTc5h\n/NYmnsNaf+aWNFBZdh5ndy+eeuqpC59EFowr+voclqtbqW7qXWkuo7iej4+VkBjmxY0JIdi3VpGw\ndwNOLSWcWPJPWsOWkRjmhZO9jYXeJiYTxgEANGkNnKtqZMbhX+FffZBfGzcy9bJbcXKwYVpg28nA\n/kpKdzCR9goG8x47tJC6HiLz8fFBpVJ1HiILDw+nziEE1u9lqXYdqZ89x2+mzuYvjg/x4ucHmCsV\nsvby1Xz34b9I6M9JZMElT1+fQ53BRLNLUK/252paeGN/ERG+Lty1JAInTTmJe+/AXltNRtK/0Ics\nJDHcq1eR94nOuDwHsGHDhl7GSmc0UVqnYerxJwgo2c6fDLciJawn2MuJaKVr56ygZynIDiMjctMH\nR4cWUschMqDzMFnPQ2RbPvyCdP9bWHnTvVTknWTx9z/HreQg6cZgbH0iKKmqw8PTq9vrX8rCeoK+\n6etz+N9tbYWFulLbrOOl3fm4Odjy0xXReGpLmbP7Nux1tWQsfQtD6CISI4Txt8S4dAA9MZtlSupa\nucz4HSGFH/Cq6VoOB97GiilKvFzsfxCFYvypd442HVpIvr6+HD9+nNzcXNLT04mMjLR8iExSEHfP\nS8ReeR/PLDdzZOFeFgWaeX7XWRSOrmQXlWDuUtz7UhXWE1wYS59DG3clZRWV3a616k28mJqP3mjm\n4VUxBBpUJO6+DRtjC+nL3sYUMpfEcC8cbIXxt8SEcABnq5qYp93HKvM+PlOs4TWb27hzUQS2Ngqm\nBfYoAD2O1DvHAh0nqGtqaqisrCQ1NZWqqiqKiooueIjsvMtstgf9Ajt7B7Yq/sCs0o/IqTXzxccf\n8Z9P/sv58+cpKipi8+bNnSsLwfih5+dQazCx/0QOHl4/aFEZzWZe21tAhVrLA8uiiJFUJO75CZLZ\nSPryd5CC40kM98LedkKYuUEx7vcAqhq1GI+/zzXmHew2J/I/2vU8vDoSN0c7ov1ce20IjaZ653jF\nUonJriG6jsc9Nf7r7YOxu38vh/+aTPz5dwmMv46vtHfQem4vH3z0MXq9ntjY2M6DaYLxQ8/P4c7D\nmez/bgfx85cAHemeKrLLG1m/MJz5Tirid9+J2cae48u3YhcwlfhQz3FXwtHajGsHoDWYKMv4hplH\nf0MWMWzUP8zqaYFMD/LA09muW+ing/Gk3nkp0TPjo7a2ts2wu/hwTHkLS9f4E9+4nelGPe8GrWO6\nNh9vPyU+Hm6jPXTBMND1c1hRU8/UxCri5y8hcnIsADtPV7L3bDVrZwRwtXcp8XvuxmjnxvHlb+Pg\nF8VsYfz7xbh1ALIsk591hBn7H6TJNZL1Nb/ERWHk+vhgFAqYFuTe53PHi3rnWOJCv0dLGR+bN2/u\nvF9TV8/0X39A2rM38hN24Wdq4E2nH7H/8FGcFXpamptFuu4oMZxZcHFxcSxOWkZBTTNrrl3XeT39\nfD2fpJcwJ9yLu0LKiN97H3pHH9KXvY2zXwSzQ0Wp0f4ybl1k8fkCor+9E5ONEw8rfoNadmG+YyWp\nKR9x/NvPxUGQMYSljI/4+Hjy8vKAtnhwaWkpuQHXk2J7FWts0pmb9wIlJSp8wqK5+uprWLt2Lamp\nqT+UqBSMC8rUrchdznsV1jTz7/1FRPq68Gh0GQlp96Bz8ufYivdw8RfGf6CMSyvY2tSA9+c/wVbf\nyOuTXmbPKQfiHapxtzHgZG+Dj+vFZV/7M6MRB8asg6WMDx8fHxob21Qeu8aDDcoV/D3XyMGTn3Nn\nvIlzjpdR1awf1JmAiXSG41JEVadBo/tB47+mWcdLqfm4O9nyx2klzDn0c1rcJnF82RbcfQOZFeKJ\nQhj/AWEVByBJ0hXAC7SVhPyXLMtP97jvALxNWzH4WuAWWZbPWaPvXpiM2H56Fw7qs+xNfIkXvnci\nNsCFqOZ8JAUEeTpZRetnoAfGBH1jKfOqtrYWd/e2MF3PfRl7JxeqXadw56RSVMatvNd0K8W1mkGn\n61pyBMI5jC5ag4n86uZuP7+Umo/JLPO3aedYcPQxmjynkrH0TTx9/ZkZ7CGM/yAYsgOQJMkG+Adw\nGVACHJUkKUWW5dNdmt0N1MuyHC1J0jrgr8AtQ+3bIoYWJEMLOfF/4M9nQ1BIrdy5OJKjOzPwdXWw\n2mGQjrDF3r1tdYIvNgMVhqRvLGVe2dvb86tf/aqzTc9MopRtZv7ldyd31TzNvca3OFd4OZWOriJd\nd5yQU96Iqf3AlyzDWwfOUaZu5dW4fJZl/gG192wykt7A28eXOGH8B401VgDzgHxZlgsBJEn6ALgW\n6OoArgX+0P74E+BlSZIkWe6rgucQcPSg5cdf8NbXueRVlTDXsRJ1qSOZxw7gYNbj4eFOTEzMkLsR\nB8asx2Ayr+LjplFap+fhult5xvMD3D64ka3y9dy4bn1nGzGLv7To+Htdfv0t1Dbr2x+v46tTZaSf\nKOPZSSe4/Oyz1PvN5+TiV1H6+jA9yF2o9w4BaziAYEDV5ecSYH5fbWRZNkqSpAZ8gJqeLyZJ0kZg\nI9DLwPaXvBoNn58oJci2GeeaXE4cKWPRnARCAv2ora0lIyODzMzMIYVqxIEx69KfzKuu5wVCQ0OZ\nHZ/AY4//kZUnZ/Jj16PcHP4ujq0raNbNnPA67yONtZytwWQmt4vM8wlVA5+fKOP3/vv5Udkr1AQk\ncWrRPwhUehEb2Hcmn6B/WCMLyJL77Tmz70+btouy/Losy3NkWZ4zGGOqN5r57WdZONnbkOhQTWlx\nITdcdy1hwQEoFAqUSiXx8fFD1vbpCFsMuD6uwGrMnjWTG6+5gsvWJnNg9p/xVvoxc8+dVO58jlad\ncbSHJ7DAli1beh3460q5Wtup9VPW0Mq/9hfyv+47uUv9ClVBqzi5+FVC/L2F8bcS1pgmlQBdj2KG\nAGV9tCmRJMkW8ADqrNB3L/QmM1F+rqya6kdV+hn0rS0kxcfyn9PpnW061CiHgjgwZn0GM3u0s1EQ\n6e2Izqzltubf8hfpTVan/4mq2mxkwwIkOyfrD1RgVTocQr1GT1OrAQCN3sg/duezyeZzHtB/QEXo\nlWTPf5YwpQcx/uLwn7WwhgM4CsRIkhQJlALrgFt7tEkB7gAOATcBqcMS/wdcHWz58/VxfLTjADsO\n76e2vISn//JnJEnqlAzoUKMcKuLA2OhgSUbihVffwGCu55emX7JR+pQHzn3ENVIaacr13Z4r9gXG\nBj3/DkXnizl0PJOW5ma0ra2cNAaw1jaDh2w/oDz8WrLnPk2EnxvRfsL4W5MhO4D2mP5PgR20pYG+\nKctytiRJTwLHZFlOAf4NvCNJUj5tM/91fb/i0MnOyiTj8H5mzYxj6Zw4jEYjJ06coLy8nPLycvLz\n80lOTh7yPoBg7ODhZEe0l5l5s6fwzI4bOOUwgz8b/0Zy1QsY9rhgt/QRUAhFyLGISqXi8PFThEdP\nw93Ti3KXKKTtL3J91GGqEy7j9Ny/EK4Uxn84sMpOmSzL24BtPa79vstjLfAja/TVHw4dOMCy1VfQ\nfP4UNgoFiYmJ7N+/ny+++IJ58+Zx0003sX79epG3P87wcXVgcYwvWoOJ53fJXG3+My85vEbinj9i\nKtiFzY1vdLbteojP3b0tM0ysDkaHjMzThEZPw9Pbh2KDK941x3gy9gjvFbiTtOnvBHi7irDPMDEu\nUyVqaqpZFx/Lt6osoM3ABwUFYWNjw7XXXguIalLjlRh/NwwmGRn4vx1G7tA+wpvzy0nI+hOnHp/L\niYpYTlVLfP3114SGhhIdHU1GRgYfffQRERERzJs3b1jHNx6djEqlIi8vr5eQX3/eY5PWQGVdA5Fx\nXtSb7HHVl/KS/cuc8wxlV2sQV3u495JsF1iPcekAggP90aurul2rqqrC39+/2zWRtz8+6CnJsWTJ\nEpJiwjm6bxdprUH8LHsKdwT/laaMP/BwyDHe0gez+r5Hefq5FwFYsWIFCxYs4PvvvycnJ0foCQ0A\nlUrVq05zVyG/rvRcdUVHx5BZqsbJ2ZWaejU2jhpetXuRJo+pPF+ShJN7EXHBHiLPfxgZl2Jwy5Yu\n7ZWiWV9fj5+fX7d2Im//0qerJEdycjJr165l9+7dyHXniXQzkeRchlpr4JnPjuL6kzcp9F2Bov4c\ni7Mew8u2FQcHB5RKJUqlElmWrZIiPJHIy8sjvr1OsyUhvw56/p3i4uI4nHGK7MwsgsKiyM08xi+b\nn6NC9mZHzJPk5OYxb9Z0UcxlmBmXv92O7JzMzExSUlLYvn078+fPp6KiQuTtjzP6qh27f/9+Qjyd\niHQ18YvVk2lW1/FWlpZPuZxc98UUVTUTpT6EryYfSTZ3ag/5+Pj0uSq8WA77UNtfijQ2NuLj49Pt\nWlchvw66/p1KS0v5/thxKiqreO+ff8egb+KFqEN8XqBgw/5QDh45xvz4OKIiw0fyrUxIxmUICHqn\naG7ZsgWVSkVmZqZV8/bHUyz3UuRCkhzh4eGEeDqROMmHNQkxHK4sZ49bMDPC5vOiaQGNpi1E6M4y\nJed5UlUxxEybbbUU4YlAZmYm5eXlvP/++wQGBnZKrHQV8uug4++kUqk4lXkan/AYwqbN4eT+b0go\n3kpjhDemhfcyR3LniuWL2PrPlzideVKo7A4z49YBWCI0NLTzLIAw3OODi0lyKBQSccEeeNvqUZ7Z\nRkX4ZZzySCQ5sJmUChdcTUaQ8rkvtpIcKYRdGdX8+te/Jj093VJ3Q2I8yYd3vJfFixdTXl5OcHAw\n2dnZbN++nYyMDGJjY7u17/g7ncnNRRk5BRd3L6jK4lq7Q8yO9ueZoulMC3XH1FxL1rGDQmV3hJhQ\nDkAw/rhQDecOI65QSCROiwbgeE4quVkS//H0YeEVN+FiY6K4PIuvzhwk6PwbXBm3hrgZM3o5gKEa\n76eeeqrXZumlbNi6quF6eXmRl5dHXV0dW7duZfr06b3qNNfW1rJ582YKVeUsjInHufoE7rkfYevh\ny58LppOTk4OtowtmbTM/2/hiv1V2BUNDOADBJc2FJDm6GnFJanMCgUEhTG+CQ6ZJZOq1LHMuxT0g\niAiPNVzR9BFhrbvJfO4G9mQ40tDcikajwd/fn5KSkgvOSi+W3tl1sxQufcPWNfTWsbI2m82oVKpO\n49/1dxIaGkqdRs/xrDNkbt9KsncBZjdfjgesw6tFw0xHR55/5s9s2rQJtVrdrS+RrTd8jGsH0N8w\nz3jMzZ5IWJLk6KvI/OZfPMDZyiYO5tfwl69OsUcTzDKnUvLVBlI816OsTaHi8A5+Pi2A/CX3knj5\nWh599FHuuusuKisrgcEZb0ubpZeyYeuriE9f+ydNWgO2rr5sSArDvmAnUVGR/LVgCtclLuTbT99l\nYcIsIiMjWbhwISkpKd32dUS23vAxLrOABBObnimHkZGRpKSk8Pbbb/PKK6+gqzrHkhgly51KkYE9\nrcGoDXaoGnT8V+VJ8mP/ZoqHjuTKvxOpzyU4OJiCgoJufQzUeLu7u1NbW9vt2qVs2Hqq4VZXV5OR\nkWExq05rMFFa30qS6QCPKXfjHT6ZJwtncO5MJmWFuSxMmEVYWNuqITk5mUOHDolsvRFCOADBuKNn\nyuHZs2dZvXo1QUFBncXjWyuLiPayYblTKRJtTqDeYEdZdR1u8dfxZdAvabTzg/fXEU0x+T3y2gdq\nvGNiYsjIyLikDVvXtNaeqdaZmZnExsb2WhGZzTIni+u5zPAtyeYdfG2aR0rAT1m6aBHTZsZz7203\ndRp/AA8PDxISErqlcAuV3eFjXIeAutK1mEhXxlNmhqCNnTt34ujoiEKhIC8vjyeeeILdu3eTk5PT\nLXyjdHMGdEhSKXs0wextDULp6MHOI1k0mNzYHvAQt7t/T9RnH7C30J+wIH88lQHdNpr7S0dc3Npp\nyKNJzzKdlqhoaGHNwd8SYt7Hf4wr2CLdyEy7OmwUEvNmT+ebbV/12sBft25d5/6NCMsOLxPGAVhC\nFHYfn3SEW5RKJY2NjYSFhXXLTe96TkDp5oCMFmhzAiX+C/n6q69wtYO4SSEUJfyWs/vqWGfYScGB\nl9jOrEEb74mWhqxuaiK56V1CGs/wqvEa3uFq/nDLcmxtJU7v/RJXB3cSExMvuoEvGD4mnAPo+sEb\naGF3waVBR7glPj4eV1dXDhw40C03vWf4xs/NkfLSApTZ+8hpkNgJhCjqyc44gk6nZc29TxFnfz/6\nD9ZjkrJIVcwnPT1d/I90oefme0NtJfM4iSmshSd87+IL81KWOZeiUEhMDXCnuL1kp6ipMbpMOAfQ\nFVHYfXzSNdxSVlbGSy+9RGhoKMHBwd3CNx0G/KmnnqKyuIAF8TOItfEkrd6d0qx9LJ0SQOLaW4gI\n9wLHOL4OfIRVVW9wRcVLHPK5GdjQr/FYkp7uL5dShlrH+0yICeDWiL3Uq7XcfXoR9bYzWOVXhq0k\nM0npQpCnqNI2VpjQDkAUdh8/9DSQXcMtHWGGlJQUi+Gbrjn6pfWtLL18Gb/fqiUtP51rajVIksSc\ncC+u3/gYtN5H6QuXs6T2fdjuAWv+BDZ9f4xUKhVqtbpXmHE8FiNKS0vj8vgQQlOfxqyAh11+S2m0\nB8scK3FQuOHpbM8kpStwcYd2KTi88cCQsoAkSfKWJOlbSZLy2r979dHOJEnSifavlKH0aU1EYfeJ\nQUeY4brrrmPTpk29DG9UVBSPPPIIAEGejswO9WR1gB6TppG/7czlXE0Lx4vr0eiN4OTFd/73ke2+\nHA6/Bu/eAJq+y1vn5eVZFKsbj4qjldn7idp9Pxqc+LH+9xwxT2G5vwG5VY2Loy1Bno6jPURBD4a6\nAngM2CXL8tOSJD3W/vP/WmjXKsvy7CH2ZXVEYffxx2DCLV1XgpIkMT3IHVd9HbGeMo2SxN925vLL\nNVOQkJgT4YUs2XDU+3qmr7gZ0xcPoXkuES+/e6i3D+712h2b0F0Zd2FGWca8/wX8zn1Gtm8cvzQ+\nQIFZSZJTGbaNFXj6+RHq5SR0/ccgQ3UA1wLL2x9vBfZg2QGMWcQm1PihZ1ZXbW0tGRkZZGZmXvBv\n21NP6Ny5c5QXniFuyiQWXT6FZ3fk/uAEJDCYzNjZKCD+NrYfyWdF1b+5svx5DvjeSse+QMdYioqK\n+NOf/oQkSZ0hqUs9zNh1X8LGbGBh7Ycozh8ldu4y7sicRI3SxGK/MhwaSzh3Npu/PPE7sk6dGN1B\nCywyVAfgL8tyOYAsy+WSJPn10c5RkqRjgBF4Wpblz4fYr1URhn980DOrS6lUdhZ4udCqztJKcNq0\nWIKDQ4gNdONXPZzAuZoWInxdAKhxCOerwF+yvPpNlldvgW89yfS7ntQ9e4mLiyM6Ohqj0ciJEye6\nhRkHcoZgzNJUyeWVL+OnO0f+9Id4o+5KqlXH8Svey9mMInStGpTeHhz5/mCnFIdgbHFRByBJ0ndA\ngIVbluu+WSZMluUySZImAamSJGXKslxgqaEkSRuBjUCvpbNAcCEsZXX5+PigUqku+lxL9SMAZoV4\nYpbpdAJ/fOcb/FUHSNO1HRxsbKgnNDSUHQE/ZV7dp0w98DxpRV+Q/NgW9h5um/UmJiaSmZnJp59+\nilKpHBdhRm+dCt5YgZeuiq2KW9jfei2Hiir50WWLUB8oRqVoYdnCeYQE+hEYGMhbb72Fm5ubOGw5\nxrioA5BleXVf9yRJqpQkKbB99h8IVFlqJ8tyWfv3QkmS9gDxgEUHIMvy68DrAHPmzJEv+g4EgnYG\nKlDWH2xtFMwO9cRsllkXJfPy9ycpi1jJSn8DwTMWcOj1vwFtWUff+9zM1OU3U/3LTYR9cwfezj+i\nziGEuLg4rrnmGlJSUti0adOQ3uNo0XVvJdq+mp/4ZGGMDuZl23vYpZ9G1ulKVk714+q4QJ55v5Cl\nC+cSFhyASqWiqamJ1atXo1KpWLt2rThsOYYYqhZQCnBH++M7gC96NpAkyUuSJIf2x77AYuD0EPsV\nCHoxEIGygWBnoyAh3Iu6wlP8z723Ye/pT5ouhBYHH3zCppBzJveHxokbUK56kOLqZq4qf46pjftA\nli1WyRpNBlKussP4z54+mXsiS3CvTufnu+DO00mklLiSpfdhfqQ36+aGolBI2Jl0hAb6Az9kQcXE\nxNDc3GwxC2rDhg0iDDtKDHUP4GngI0mS7gaKgR8BSJI0B7hfluV7gFjgn5IkmWlzOE/LsiwcgMDq\n9Izlu7u7WxQoGwx2NgqcTC3Mj5vM2exP2KMJ5m87c5nvpERVdQStwYSjnQ0ASdeuJ2W7F3NqP2FB\nzScUvVhD5nF7YqbNGvI4BoKlQ2SDqVGclpZG8sIplL+9kSP51Uydv4pi3Sxyy9Vk1xQTPsmBu35y\nOwpJIjbQHV9vz35JcQhGnyE5AFmWa4FVFq4fA+5pf3wQEGs9wYjQH4GywRLg74cvTSgdzCynlEPm\naFIr7PDQmfn0qx3YmXWdMe6Va6/h2Wdy2FpjZrZtKg+FOdHk3TtNtC/GTDF5s5nqjG2E1Rzhi/MG\nbJqStB8AAB3GSURBVOfczlGn6dTauZJZksqk2QvwOr8PG8V6ov1cCfJ0GrAUh2D0mNAngQWCgZCU\nlMQ3277C2azFx86JO2Kc+dP2b2jAmSmT4gj1dWPhslWkfvcNK1euZOWqVcAqNqydR82/b8a3egt8\nrIYrngY3S3kVY4ctW7bgaqjhJtvdKMtSKY5dzndGM1McYjlvcONItUxYUADL/IwcyWkmzMe5MzOq\nv1IcgtFHOACBoJ90DTE1qNVM0esId9DQMO0GMhz9cTeVUW/jycIVl5OWloqzs3PbE/2n8XXgL5ih\n3kXima8h7ztY9ijMvx9s7a0ytqFoBvV6rsnANPVu4hu2gZ0DSet/x1v5ttTr95BdqSXT6IZtwX7u\nv389GQf34OPlwWR/t26vORApDsHoIRyAQDAAuoaYbrt9PavTdjMnQM8+LextDWJ+dQshnu7knS9l\nVuwPp5DvuPNu4G6oLYBvfg3f/g6Ov83/t3fvUVVd96LHvz8Qgc1TwQdPQUVExaAi8R3iO5omxrS5\nqWmOj9va1CRNk+ZlHT3n2J7b0Rw72tPbm+QmN01NT2xMmkRrjMZHjDGmNSpigooiPtmgoFseIiII\n8/7Bhm5wI+gG9gZ+nzEcshZr7/nbDFi/teac6zeZugKS7gcv96/NtHr1avpV5nJP7aekFR/B6p9E\n1A/fw6cylBjfA1z45AsObHqf3tGDmDwijqCgQHoHW3jyX3920/fVhy09lyYApZxozYnKx9uLmH69\n8bt2kXSLF7sqIlm19RgPJ/pS0zOI/OKrRIQ0qX8TNggeeQ9ytsCWFfDXRdAnCe56ri4R3KSw3K1o\n7UJH9cdVF+czvMcZxvQ/D0Pi+bTv9zntO4zRZUFcvFzBocoQqiY+RuKlI4Sd28upY1lk7+vLgnlz\nGDlyZJvErDqeJgDV5XTkVWZiYiJHjhyhb7xwl991vrrSmz+s3kVaQiS+fv5U1dQ0miHUYMgsGDwd\nDq+Dz/8T3l8CQZEw6nsw+lGgrpLo8ePHW30Srz+uX79+WK1WpwsdOcrLy6P0TBY/iMxmQtwxcst8\neePMEMalP8v7GzaRX5TNwVMXsPoOIPNKMAN6lJEa2xOvAZMI9OvBz59ehpeX1vfpzDQBKOWC+n7u\nY8dyyL9wibGjJ1IZl0h2SCI+VTaGmBK+OnWJ4ZHBhAf6Nn6xlzckfxuGPwDHNkHGati1CnatIq64\nP1l5QYwbOx2/funNrlbnbFW7559/niVLllBYWAg0XujIYrHgW1NO/JUDXNj3MffFVxDV28KH5ePZ\ndLEXR8+f4ZP/9RJDR49n9ORUjhDF/s+3kZ5+NwMpQgSC/H2IDvXXk38XoAlAKRfVD3gaYxh/z4NY\n3ljNvspyvrkWTmlNT6Zdu87BsyXE9LYwuG8g3k1PnF7ekPStun8lZ+HgO+z77e94Ij6f+Gs5XC4I\nI+jr/dw3aDCbP36P5KTEhsFjZ6vaRUVFceLECQID62rvc62cWJPPhX3rWJZQSHjVWQTDuiu+xC74\nHau2nOGb/BMkj7yDwis1BONNYdEFDvUYQmW4hYfmz6PnuW8QP396B/akf7CfVvbsIjQBKHUTtzK7\nRqRuucPoUF96lBZy5FoVh6vC+M3WYzx21yAAbFeuMSwimFBLM7N/QmMh/QW+/K9/EDtsIoXXsomo\nPE7Qsc3EXrFxYVcVVL0K4QnQK44Lu08Rm1jKsNKjGBH46hqDOU3utl38ZGIAIdVF8OunOXvpOn2K\nrkPCYA6GzuKsZSTnoo5xtu80Dh//V+KHDKfE+FNYWMSQsXeztyqC0hMHee7BOaREB/OX1z4jeehg\nwgN9W72Yi8c8y6CapQlAKbvWjh20dFxYoC89e3jhXVxCsFcVB4qjWPnREZZMjGNkdCj7TxcTEerH\noD6BN44N2AUHB3O2DK72mcLR4Cks+pd/4f/9+kVK+n4FE+6GomwoPkOfq7mc3fAr0nrZZxFtXseg\nC9V8fhJKBveH8H6UJNzHhsNlTP3XBXycfapubOHAMaxWK08981NOns4jfHAK5Rcvcr7kCgXnvQiN\nCye2Oo8xA3pRVHCWYQOjb+zCUp2eJgClmtHSTJqbLT4T5OfDwD7e+FyqYO7kYby26wT/e0cuM5L6\nMX90FOdKKiksqyQq1MKAMMsNicDxadqwsDBOnTnDrkNWkpKnw/SfNxw3eXIWG7Ztpba6gj69ezFx\nwgTW/MdL9LurD7+9UEbZiTJmhsQz9eG62Ndv3UV2djZDhicTNyIVv5Bwvlr5PJs/XEtNRDJViTMJ\nr75Ikm0XpYF+2M7nkXfwC+bMmkFGRkb7/9BVh9IEoJQTzgZXHQdhm1t85oMPPmiUFAYPTiAlNpSf\nzUnirxlWtmUX8o21hEfHD2Bo/2DyLlWQX1JB3yA/YnpZCLH4AI2fpq1foyApKemGmvpN6x+VVXuR\nOGxEo+Pq71hqaw3fHD5K//ghXMEC1YYpk6cxYup89u3fT3DcBMaE+zA5yoc/v/yfBAUHUXh4D3Nm\nzSA5OVkTQBekCUApJ5wNrtbPpElOTna6+Ez//v1Zu3YtaWlpjZKCuXSWEVGD8PXxIiU6lP/ec4bf\nbM1h0uBwHhgVRYi/D+dLKzlfWkmAbw8iQvy4dr2m0dO0jmsUNNXcWgb1LldWc760knOllRRcuETs\n8BAAKmu9eevvp8mNnI5/yGlirZ9y+XgJhRPuYvbcb3FH4kC8vISMjIzbenJXH/ryfJoAlHLC2eIy\njlUsnX2/qKiIqKiohkJnjiuSLVuWTLC/Dz29vRnUN4CN35xjy+Hz7D19iamJfZk9vD+Bfj24cu06\nuUXl5BaW4+vjTaBvDyy+3hw4+HWrH+zavv1TLpWW0tMvgD7RcUSn/XPRev+AQEpKSigKGEh2VW/M\nCRvj+xl6Tx9HgCUALy/hmceXsvnDtTe8t57Qux5NAEo54WxxGccqls6+X1RUxLhx4xrtc1yRLMTf\nhzsH9uboucs8ODqaSYPD2fB1AVsOn2dnThETBoUzJSGc6F51NYSuVddwrbqGo7kFHD1TSFDEQOKT\nw0gYM4l3/7aJs7YKEocNp9YYxs2ez5FDh/jTex9hiYinX2IvykqKOXH0EKdysokfkkR55XWuRaTw\n0den8BkcS0yvKzw0NprcfTsZdeck8nMPERHqT7/gJk8vqy5LE4BSTjRdKL5pFcum37fZbBQXFzNo\n0KCGB7AWLVrEqVOn2Lx5c8P7+nh7kRwdQnhpT3p4Cz+YPJA5yRFsyjrHrpwL7DhaRN/qQiqzM+hZ\naSM0wEL55VIe/9l/8M2+f1BRVYsE9SMx7W62fLoD375xDe+98/NdjJ08nW/2/QOA0N5hDBw6gi2f\nfobvBT/2nS6mqlcK/eN98c/7FO/jxZzucRdjJ0xm1uQ0IuelN7zXokWLyMrKahhb0KUcuyZNAEo5\n4WyheMcqls4Wn7nzzjvJyclpNmk4igjxp5elJ7lF5QD8YPJAHh5bzYbP9rBzZwZVsXfRI7gPpZfP\nYdv6Cp9lnYYaX4K8qjHG0DcimtJiW6P3LC220btfJMU1PbHV+OEbNZSsymLyv1pDeEARMT7lJPiW\nEhLkA3F34u0lLP3+EmJ7W+jh3bgYXUuD4LdKu488k0sJQES+A/w7dat+pdkXgnF23Gzg94A38IYx\n5teutKtUR2ipiqWzxWfqSx87SxpN+fl4MyIqhAFhFk5dvFJXZqEkl188uZB/fLWP89dtnO8ZSnHk\ncD5cv4HQid8FYPM7mVgqbVzJr8a6ORsvESqrazhx6iqfvrqFHqF1YxPBBWXE+V5lSMog4gNPc6kw\nn5NnTlJ19QrhvUK4Y0QSA/sEOo2tpUFw1TW4egdwCJgPvNbcASLiDbwMzACswD4R2aDLQqrO4Fav\nXG+n9HGQnw8jo0O5WlXD3yrLiIqJJWDfPxjUs4xBlDF4ZBTfZO4j9uo3ENwPS2gY2Xv2EzE0Bd8e\n3tTUGnoH9CQ0bTylJzIJ8r1ATO8AJoy7kz0799EnIY6t69dy8ZyVIQkJTJ50JxaLhczMTLKyspye\n0FsaBFddg6tLQmYDLdUFSQNyjTEn7ceuBe5HF4ZXqpF3//LflBRaGdCzgpgwC2UV1ZRfu46vnz/D\nklMoOnWAq1fKiZtwF9MemUf8kCS2rKubrTPrgYeBBE7lRPHumy+Tdagci6kgYeAALtsK6RtiYd7s\nhXh5eTUsz1g/Q8lZAmhpEFx1DR0xBhAF5DlsW4E7O6BdpTqdhIQENm78iKqqKiLDwrh48SIn8o7z\ngyefZu/+/VRdr2XqfQ9xvcbc8FpfHy/SxqRw2TqdAL8e/PD7/5NXXnmFef/jQQ7s2U2fPn3w8vJi\n1KhRZGVlkZ6e3jBDqamWBsFV19BiAhCR7YCzBUxXGGP+1oo2nN0e3Pjb+8/2lgJLgRtuQZXq6mJi\nYhqNIwQHBzMyeRgzJ6dRcKLupjk9sS81tYbqmlpO9AvES4S7h/ZtqDKa6VCzp74rJzg4GJvNRp8+\nfQgLC6OsrKxh25mWBsFV19BiAjDGTHexDSvg+Px6NFBwk/ZeB14HSE1NbTZRKNVVORtcbsrbS/D2\n8sa3h3fDdj3HsYf6rhzH2kK1tbWICJmZmSxfvrxVcegsnq6pI7qA9gEJIhIP5AMPAws6oF2l3KKj\nTpatWfaxvivHz8+PoUOH8ve//53c3Fzi4uJISkrSK/puztVpoA8AfwD6AB+LyEFjzCwRiaRuuucc\nY8x1EXkC2ELdNNA3jTGHXY5cKQ/SESd9xzZaO0/f2fMKDz300A1F5VrTpup6XJ0FtA5Y52R/ATDH\nYXsTsMmVtpTqyppezdtsNmJiYpo9Ad/KPP3WdCmp7kmfBFbKzZxdza9YseKmr7nVefpNV+nSK3sF\n4NXyIUqp9lR/NV8/TTM+Pp5Ro0Zx/PjxZl+j8/RVW9AEoJSbObuar5+q2Zz6wd0LFy5QW1vbME9/\n8uTJ7R2u6kK0C0gpN3N2NW+z2QgODm72NTpPX7UFTQBKuZmzp27ryzU449iPr/P0lSs6XQKorq7G\narVSWVnp7lDczs/Pj+joaHx8fNwdirpFjidxZ1fzztb/VaqtdboEYLVaCQoKIi4urqUidF2aMQab\nzYbVaiU+Pt7d4SgX1CeDm63r2xy98leu6HSDwJWVlYSFhXXrkz/UVWANCwvTOyGl1G3rdHcA0GL5\n6W5Dfw7up1fgqjPrdHcAnqCwsJAFCxYwcOBAxowZw/jx41m3bh07d+4kJCSEUaNGkZSUxMqVKwEa\n7R86dCjPPvusmz+BUkppArhlxhjmzZvHlClTOHnyJBkZGaxduxar1QrUzejIzMxk//79vP3222Rk\nZDTan5mZycaNG/nyyy/d+TGUUqpzdgG5044dO+jZsyePPfZYw74BAwbw5JNPsnPnzoZ9AQEBjBkz\nhhMnTtC3b9+G/f7+/qSkpJCfn9+RYSsP0pq6P9q1pDpCp04AKz86zJGC5p+WvB3DIoP5t28Nb/b7\nhw8fZvTo0S2+j81mY8+ePfz85z9vVJ+luLiY48ePM2XKlDaJV3UuTev+RERE8Kc//YmgoKBmSzor\n1V60C8hFjz/+OHfccQdjx44F6uq6jBo1ipkzZ/Liiy8yfPjwhv0jR46kf//+3HvvvfTv72yRNdXV\nOdb9yc/PJycnh+nTpxMZGck999zDjh07yMrKcneYqpvo1HcAN7tSby/Dhw/ngw8+aNh++eWXuXjx\nIqmpqUBdX//GjRtveF39/pycHCZNmsQDDzxASkpKh8WtPINj3Z/jx4+zcuVKPvvsM7Kzs29a0rmt\naNeScqR3ALdo6tSpVFZW8uqrrzbsq6ioaPXrhwwZwvLly3nppZfaIzzl4Rzr/pSVlREbG9uo7s/N\nSjor1dY0AdwiEWH9+vV8/vnnxMfHk5aWxsKFC2/phP7YY4+xa9cuTp061Y6RKk/kWMUzMDCQL7/8\nkszMTBISEgAt6aw6lqtLQn4H+HcgCUgzxuxv5rjTwGWgBrhujEl1pV13i4iIYO3atU6/l56e7nSf\n435/f3+dBdRNOdb9KSgo4A9/+AMxMTFERUU1lHTW1btUR3F1DOAQMB94rRXH3m2Muehie0p1eo5L\nNI4ZM4ZVq1axYcMGLemsOpyrawJng5YkUOp2OSYDHaBVHa2jxgAMsFVEMkRkaQe1qZRS6iZavAMQ\nke2As0nrK4wxf2tlOxONMQUi0hfYJiJHjTG7mmlvKbAUuGGZPKWUUm2nxQRgjJnuaiPGmAL7/0Ui\nsg5IA5wmAGPM68DrAKmpqcbVtpXyRNrdozxBuz8IJiIBgJcx5rL965nAL9q7XaW6oqZ1hLR0hHKF\nS2MAIvKAiFiB8cDHIrLFvj9SRDbZD+sH7BaRr4G9wMfGmE9cadfdAgMDne5//fXXGTp0KEOHDiUt\nLY3du3c3fC89PZ3ExERSUlJISUnh29/+NgDHjh0jPT2dlJQUkpKSWLq0boikoqKCRx55hOTkZEaM\nGMGkSZMoLy9v/w+nPJZjHaH77rtPS0col7k6C2gdsM7J/gJgjv3rk8AdrrTTGWzcuJHXXnuN3bt3\nEx4ezoEDB5g3bx579+5tqPuzZs2ahpIR9X784x/z9NNPc//99wM0/DH//ve/p1+/fg3bx44d07V/\nu6jWdgfV1xH6/PPPATqkdITq2vRJ4Dby0ksvsWrVKsLDwwEYPXo0Cxcu5OWXX77p686dO0d0dHTD\ndv0f8rlz54iKimrYn5iYiK+vbztErjoLxzpC9bR0hHJFpy4Gx+YX4Xwb3/72T4Z7fn3LLzt8+DBj\nxoxptC81NZW33nqrYfuRRx7B398fgBkzZrBq1Sqefvpppk6dyoQJE5g5cyaLFy8mNDSUJUuWMHPm\nTN5//32mTZvGwoULG8oFqO7JsY5QPS0doVyhdwDtyBjT6CG5NWvWcPDgQQ4ePMiqVasAWLx4MdnZ\n2XznO99h586djBs3jmvXrpGSksLJkyd57rnnuHTpEmPHjiU7O9tdH0V5AMc6QrW1tQ2lIyZPnuzu\n0FQn1bnvAG7jSr29DBs2jIyMjEZ1XA4cOMCwYcNafG1kZCRLlixhyZIljBgxgkOHDjFmzBgCAwOZ\nP38+8+fPx8vLi02bNpGUlNSeH0N5MMc6QmVlZVo6QrlM7wDayPPPP88LL7yAzWYD4ODBg6xevZpl\ny5bd9HWffPIJ1dXVAJw/fx6bzUZUVBRffvklxcXFAFRVVXHkyBEGDBjQvh9Cebz60hHz5s1j2bJl\nevJXLuncdwBuUlFR0Wjg9plnnuGZZ54hPz+fCRMmICIEBQXx9ttvExER0XCc4xhAeHg427dvZ+vW\nrTz11FP4+fkBdVd3/fv3Z+vWrfzoRz/CGENtbS1z587lwQcf7NgPqpTq0sQYz33YNjU11ezf37jC\ndHZ2tnaDONCfR/ezevVqQJ8mVs6JSEZrS+5rF5BSSnVTmgCUUqqb0gSglFLdlCYApZTqpjQBKKVU\nN9UtEsDq1asbZk4opZSq0y0SQFuzWq3cf//9JCQkMGjQIJ566imqqqpYvXo1TzzxhLvDU0qpVuny\nCaC+hvr69et55ZVXXK6dboxh/vz5zJs3j+PHj5OTk0N5eTkrVqxoo4iVUqpjdOkE0B4LaOzYsQM/\nPz8WL14MgLe3N7/73e948803qaioIC8vj9mzZ5OYmMjKlSsBuHLlCnPnzuWOO+5gxIgRvPvuu23y\n+ZRSyhUulYIQkVXAt4Aq4ASw2BhT4uS42cDvAW/gDWNMh1Rxa48FNJyVfQ4ODiY2Npbr16+zd+9e\nDh06hMViYezYscydO5czZ84QGRnJxx9/DEBpaalrH0wppdqAq3cA24ARxpiRQA6wvOkBIuINvAzc\nAwwDvisiLZfIbAPtsYBG0xLPTffPmDGDsLAw/P39mT9/Prt37yY5OZnt27fzwgsv8MUXXxASEnLb\n7SulVFtxKQEYY7YaY67bN/cA0U4OSwNyjTEnjTFVwFrgflfaba32WEBj+PDhNK1PVFZWRl5eHt7e\n3jckBxFhyJAhZGRkkJyczPLly/nFL35x2+0rpVRbacsxgCXAZif7o4A8h22rfV+7a48FNKZNm0ZF\nRQV//vOfAaipqeGnP/0pixYtwmKxsG3bNi5dusTVq1dZv349EydOpKCgAIvFwve+9z2effZZDhw4\n0FYfUSmlbluLYwAish3o7+RbK4wxf7MfswK4Dqxx9hZO9jVbglRElgJLgRu6b25VeyygISKsW7eO\nZcuW8ctf/pLa2lrmzJnDr371K9555x0mTZrEo48+Sm5uLgsWLCA1NZUtW7bw3HPP4eXlhY+PD6++\n+qpLn0sppdpCiwnAGDP9Zt8XkYXAvcA047y2tBWIcdiOBgpu0t7rwOtQVw66pfhaUr+ABrRd+dyY\nmBg++uijG/YvWrTIaRuzZs1i1qxZbdK2Ukq1FVdnAc0GXgDuMsZUNHPYPiBBROKBfOBhYIEr7d4q\nrZuuuhL9fVZtxdUxgP8DBAHbROSgiPxfABGJFJFNAPZB4ieALUA28J4x5rCL7SqllHKRS3cAxpjB\nzewvAOY4bG8CNrnSllJKqbbVKZ8E9uRlLDuS/hyUUq7odAnAz88Pm83W7U9+xhhsNlvDYvJKKXWr\nXOoCcofo6GisVqtLT/N2FX5+fkRHO3v2TimlWtbpEoCPjw/x8fHuDkMppTq9TtcFpJRSqm1oAlBK\nqW5KE4BSSnVT4smzaUTkAnDmNl8eDlxsw3DakqfG5qlxgcZ2uzw1Nk+NCzp/bAOMMa0qeezRCcAV\nIrLfGJPq7jic8dTYPDUu0Nhul6fG5qlxQfeKTbuAlFKqm9IEoJRS3VRXTgCvuzuAm/DU2Dw1LtDY\nbpenxuapcUE3iq3LjgEopZS6ua58B6CUUuomulwCEJHZInJMRHJF5EV3x1NPRN4UkSIROeTuWJoS\nkRgR+UxEskXksIg85e6Y6omIn4jsFZGv7bGtdHdMjkTEW0QyRWSju2NxJCKnRSTLvk7HfnfH40hE\nQkXkfRE5av+dG+/umABEJNH+86r/VyYiP3F3XAAi8rT99/+QiLwjIm1SBbJLdQGJiDeQA8ygbinK\nfcB3jTFH3BoYICJTgHLgz8aYEe6Ox5GIRAARxpgDIhIEZADzPOTnJkCAMaZcRHyA3cBTxpg9bg4N\nABF5BkgFgo0x97o7nnoichpINcZ43Hx2EXkL+MIY84aI9AQsxpgSd8flyH4uyQfuNMbc7rNIbRVL\nFHW/98OMMVdF5D1gkzFmtavv3dXuANKAXGPMSWNMFbAWuN/NMQFgjNkFXHJ3HM4YY84ZYw7Yv75M\n3cptUe6Nqo6pU27f9LH/84irFhGJBuYCb7g7ls5CRIKBKcAfAYwxVZ528rebBpxw98nfQQ/AX0R6\nABZusq76rehqCSAKyHPYtuIhJ7LOQkTigFHAV+6N5J/s3SwHgSJgmzHGU2L7L+B5oNbdgThhgK0i\nkiEiS90djIOBwAXgT/auszdEJMDdQTnxMPCOu4MAMMbkA78BzgLngFJjzNa2eO+ulgDEyT6PuFrs\nDEQkEPgA+Ikxpszd8dQzxtQYY1KAaCBNRNzehSYi9wJFxpgMd8fSjInGmNHAPcDj9i5IT9ADGA28\naowZBVwBPGasDsDeLXUf8Fd3xwIgIr2o68mIByKBABH5Xlu8d1dLAFYgxmE7mja6Verq7P3rHwBr\njDEfujseZ+xdBTuB2W4OBWAicJ+9r30tMFVE3nZvSP9kX5cbY0wRsI667lFPYAWsDndx71OXEDzJ\nPcABY0yhuwOxmw6cMsZcMMZUAx8CE9rijbtaAtgHJIhIvD2LPwxscHNMHs8+0PpHINsY81t3x+NI\nRPqISKj9a3/q/hiOujcqMMYsN8ZEG2PiqPs922GMaZOrMleJSIB9MB9798pMwCNmnxljzgN5IpJo\n3zUNcPtkgya+i4d0/9idBcaJiMX+tzqNunE6l3W6FcFuxhhzXUSeALYA3sCbxpjDbg4LABF5B0gH\nwkXECvybMeaP7o2qwUTgUSDL3tcO8DNjzCY3xlQvAnjLPivDC3jPGONRUy49UD9gXd25gh7AX4wx\nn7g3pEaeBNbYL9JOAovdHE8DEbFQN4vwh+6OpZ4x5isReR84AFwHMmmjJ4K71DRQpZRSrdfVuoCU\nUkq1kiYApZTqpjQBKKVUN6UJQCmluilNAEop1U1pAlBKqW5KE4BSSnVTmgCUUqqb+v9wsHWTSGIr\nNQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x113ea9dd0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"from scipy import linalg\n",
"# easiest way to install: conda install -c https://conda.binstar.org/pymc pymc\n",
"from pymc import gp \n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def loess(x_obs, y_obs, alpha=2./3., y_sd=None, x_pred=None):\n",
" \"\"\"\n",
" Local regression, with optional inverse variance weights\n",
"\n",
" Adapted from https://gist.github.com/agramfort/850437 with heavy modification\n",
" to allow for extrapolation, add inverse variance weights, and revert to\n",
" linear regression if alpha=1.\n",
" \n",
" Args:\n",
" x_obs (numpy.ndarray): one-dimensional array of independent variable\n",
" y_obs (numpy.ndarray): one-dimensional array of dependent variable\n",
" alpha (float64): proportion of window to include (0. < alpha <= 1.)\n",
" y_sd (numpy.ndarray or None): if given, standard deviation of y_obs\n",
" used to construct inverse variance weights\n",
" x_pred (numpy.ndarray or None): if given, a different set of x values\n",
" to output predictions over (otherwise uses x_obs)\n",
" \n",
" Returns:\n",
" numpy.ndarray of shape (len(x_pred),) of predicted values of y\n",
" \"\"\"\n",
" # setup predictions\n",
" n_obs = len(x_obs)\n",
" if x_pred is None:\n",
" x_pred = x_obs\n",
" n_pred = len(x_pred)\n",
" y_pred = np.zeros(n_pred)\n",
" \n",
" # generate tricubic weights\n",
" global_window = int(np.ceil(alpha * n_obs))\n",
" # just use equal weights (linear regression) if everything is inside the window\n",
" if global_window > n_obs - 1:\n",
" weight = np.ones((n_pred,n_obs))\n",
" else:\n",
" local_window = [np.sort(np.abs(x_obs - x_obs[i]))[global_window] for i in range(n_obs)]\n",
" weight = np.clip(np.abs((x_pred[:, None] - x_obs[None, :]) / local_window), 0., 1.)\n",
" weight = (1. - weight ** 3) ** 3\n",
" \n",
" # if extrapolating, keep constant weight in the tails\n",
" x_obs_min = np.min(x_obs)\n",
" for i in xrange(n_pred-1, -1, -1):\n",
" if x_pred[i] < x_obs_min:\n",
" weight[i, :] = weight[i+1, :]\n",
" x_obs_max = np.max(x_obs)\n",
" for i in xrange(0, n_pred):\n",
" if x_pred[i] > x_obs_max:\n",
" weight[i, :] = weight[i-1, :]\n",
"\n",
" # apply inverse variance weights if given\n",
" if y_sd is not None:\n",
" weight = weight / np.square(y_sd)\n",
" \n",
" # check to make sure there's still _some_ weight\n",
" assert np.all(weight.sum(axis=1) > 0.), \\\n",
" 'Some values of x_pred have 0 weight - increase alpha or decrease range of x_pred'\n",
" \n",
" # loop through values of x_pred to generate corresponding y_pred\n",
" for i in xrange(n_pred):\n",
" # subset weights to this point\n",
" w = weight[i, :]\n",
" # generate design and data matrices\n",
" b = np.array([np.sum(w * y_obs), np.sum(w * y_obs * x_obs)])\n",
" A = np.array([[np.sum(w), np.sum(w * x_obs)],\n",
" [np.sum(w * x_obs), np.sum(w * x_obs * x_obs)]])\n",
" # solve Ab\n",
" beta = linalg.solve(A, b)\n",
" # use to make predictions\n",
" y_pred[i] = beta[0] + beta[1] * x_pred[i]\n",
"\n",
" return y_pred\n",
"\n",
"def simple_gp(x_obs, y_obs, x_pred, y_mean_func, y_sd=None, scale=1., draws=100): \n",
" # create a simple interpolation of the mean function \n",
" def mean_func(x):\n",
" return np.interp(x, x_pred, y_mean_func)\n",
" \n",
" # build GP's mean and covariance functions\n",
" M = gp.Mean(mean_func)\n",
" C = gp.Covariance(eval_fun=gp.matern.euclidean, diff_degree=2, amp=1., scale=scale)\n",
" \n",
" # observe and make draws\n",
" gp.observe(M=M, C=C, obs_mesh=x_obs, obs_vals=y_obs, \n",
" obs_V=np.square(y_sd) if y_sd is not None else np.ones(x_pred.shape[0]))\n",
" y_draws = np.stack([gp.Realization(M,C)(x_pred) for i in range(draws)])\n",
" \n",
" return y_draws\n",
"\n",
"if __name__ == '__main__':\n",
" # make up some data\n",
" np.random.seed(1)\n",
" n = 100\n",
" x_obs = np.linspace(0, 2 * np.pi, n)\n",
" y_obs = np.sin(x_obs) + 0.3 * np.random.randn(n)\n",
" x_pred = np.linspace(0, 2.5 * np.pi, n*2)\n",
" y_sd = .1 + np.random.exponential(scale=.1, size=n)\n",
"\n",
" # use LOESS to get the mean function\n",
" y_loess = loess(x_obs=x_obs, y_obs=y_obs, alpha=0.5, y_sd=y_sd, x_pred=x_pred)\n",
" \n",
" # run GPR \n",
" y_gpr = simple_gp(x_obs, y_obs, x_pred, y_loess, y_sd, scale=10.)\n",
" y_gpr_mean = y_gpr.mean(axis=0)\n",
" y_gpr_ci = np.percentile(y_gpr, [2.5, 97.5], axis=0)\n",
"\n",
" # plot the results\n",
" f, ax = plt.subplots()\n",
" ax.errorbar(x_obs, y_obs, y_sd, fmt='ko', alpha=0.5, markerfacecolor='None', label='Obs')\n",
" ax.fill_between(x_pred, y_gpr_ci[0,:], y_gpr_ci[1,:], alpha=0.3)\n",
" ax.plot(x_pred, y_gpr_mean, label='GPR')\n",
" ax.plot(x_pred, y_loess, label='LOESS')\n",
" plt.legend()"
]
}
],
"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": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment