{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdftarget(x, normed=1):\n",
    "    return np.exp(0.4*(x - 0.4)*(x - 0.4) - 0.08*x*x*x*x)/normed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def propos(x, xmin, xmax, xhigh, yhigh):\n",
    "    m1 = yhigh/(xhigh - xmin)\n",
    "    c1 = -xmin*m1\n",
    "    m2 = -yhigh/(xmax - xhigh)\n",
    "    c2 = -xmax*m2\n",
    "    y = np.where(x < xhigh, m1*x + c1, m2*x + c2)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sampling_isosceles(nsamples, normalize, xmin, xmax, xhigh, yhigh):   \n",
    "    # start sampling\n",
    "    u = np.random.random(nsamples)\n",
    "    v = np.random.random(nsamples)\n",
    "    low = 0.5*xmin\n",
    "    width = 0.5*(xmax-xmin)\n",
    "    x = low + width*u\n",
    "    y = low + width*v\n",
    "    xsample = x + y # \"Fourier convolution\" to make an isosceles triangle\n",
    "    target_dist = pdftarget(xsample, normed=normalize)\n",
    "    ysample = propos(xsample, xmin, xmax, xhigh, yhigh)*np.random.random(nsamples)\n",
    "    \n",
    "    accept = ysample <= target_dist\n",
    "    \n",
    "    return xsample, ysample, accept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accepting rate:  0.302\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFX2//H3CQHCIiCIgsQgGJDFcUUEQY2AEkBEGBgW\nxx3FBVGHkcWvYlDHZWZk5CcuOC7IjAoKjIIoKGiAKPumQNghLCIo+05I7u+PSkjTZKkk1V1V3ef1\nPP3Yy031oYWT21WfuiXGGJRSSvlTjNsFKKWUKjlt4kop5WPaxJVSyse0iSullI9pE1dKKR/TJq6U\nUj5mq4mLSLKIrBaRtSIyOJ/X/yoiS0VkiYj8LCInRaSa8+UqpZQKJEXlxEUkBlgLtAV+ARYCvYwx\nqwsYfwvwuDGmncO1KqWUCmJnJt4cWGeMyTDGZALjgC6FjO8NfOJEcUoppQpnp4nXAbYGPN6W89wZ\nRKQCkAxMLH1pSimliuL0gc3OQJoxZp/D21VKKZWPWBtjtgMJAY/jc57LTy8K2ZUiIrpQi1JKlYAx\nRvJ73s5MfCGQKCJ1RaQcVqOeHDxIRKoCNwBfFFGI67dnn33W9Rq8ctPPQj8L/Sy8/1kUpsiZuDEm\nS0T6A99gNf33jDHpItLPetm8kzP0NmC6MeZoUdtUSinlDDu7UzDGTAMuDnpudNDjD4EPnStNKaVU\nUaLyjM2kpCS3S/AM/Szy6GeRRz+LPF7/LIo82cfRNxMx4Xw/pZSKBCKCKcWBTaWUUh6lTVwppXxM\nm7hSSvmYNnGllPIxbeJKKeVj2sSVUsrHtIkrpZSPaRNXSikf0yaulFI+pk1cKaV8TJu4Ukr5mDZx\npZTyMW3iSinlY9rElVLKx7SJK6WUj2kTV1EvKzuLiasmcuD4AbdLUarYtImrqDd05lD6TOpDj896\ncDL7pNvlKFUs2sRVVBuzbAyT0iex9YmtAAycPtDlipQqHm3iKmqlbUlj0LeDmNJ7CudWOpfx3ccz\nfcN0Ri8aXfQPK+URtq52r1Sk2bxvMz0+68HYrmNpXLMxANXiqjGl9xRaf9CahjUacmO9G12uUqmi\n6UxcRZ2Dxw/S+ZPODGk1hOTE5NNea1CjAZ/88RN6TezFut3rXKpQKftsNXERSRaR1SKyVkQGFzAm\nSUSWisgKEfne2TKVckZWdhZ9JvWhZXxLBlwzIN8xbeq1YXjScDp/0pl9x/aFuUKlikeMMYUPEIkB\n1gJtgV+AhUAvY8zqgDFVgR+Bm40x20XkHGPM7/lsyxT1fkqF0qBvB7Hwl4VM//N0ypUpV+jYAV8P\nYM3uNUztM5XYGN3zqNwjIhhjJL/X7MzEmwPrjDEZxphMYBzQJWhMH2CiMWY7QH4NXCm35SZRJvSY\nUGQDBxjRfgSgiRXlbXaaeB1ga8DjbTnPBWoIVBeR70VkoYjc4VSBSjkhMIlSo2INWz8TGxOriRXl\neU59R4wFrgTaAJWAuSIy1xizPnhgSkrKqftJSUkkJSU5VIJS+csviWKXJlaUG1JTU0lNTbU11s4+\n8RZAijEmOefxEMAYY14JGDMYiDPGDM95/C7wtTFmYtC2dJ+4CquDxw9y7fvX0veKvjzW4rESb+e7\nTd/Re2Jv0u5Jo0GNBg5WqFTRSrtPfCGQKCJ1RaQc0AuYHDTmC6C1iJQRkYrANUB6aYpWqrTsJFHs\n0sSK8qoim7gxJgvoD3wDrATGGWPSRaSfiDyQM2Y1MB34CZgHvGOMWRW6spUq2tCZQzl04hCjOo5C\nJN9JTLE82OxBbr7oZnpO6KlrrCjPKHJ3iqNvprtTVJiMWTaGF2a/wPy+820fyLTjZPZJOn3ciUY1\nGjGyw0jHtqtUYUq7O0UpXylJEsUuTawor9EzGFREKU0SxS5NrCgv0Zm4ihiFrYniNF1jRXmF7hNX\nlk2bICMj7/HVV0OlSu7VU0xZ2VncNv42aleuzehbRjtyINOOtxe9zWvzXmNe33lUi6sWlvdU0Uf3\niauibd4Mqal5tyNH3K2nmJxOotiliRXlNm3iyveKuyaK03SNFeUmbeLK10KZRLFLEyvKTZpOUb4V\njiSKXZpYUW7RmbjypXAmUezSxIpygzZx5TtOroniNF1jRYWbNnHlO24lUezSxIoKJ23iylfcTqLY\npYkVFS7axJVveCGJYpcmVlS4aDpF+YKXkih2aWJFhYPOxJXneTGJYpcmVlSoaRNXnublJIpdmlhR\noaRNXHma15ModmliRYWKNnHlWX5JotiliRUVCtrElSf5KYlilyZWVChoOkV5jh+TKHZpYkU5TWfi\nylP8nESxSxMryknaxJVnREISxS5NrCinaBNXnhEpSRS7NLGinGCriYtIsoisFpG1IjI4n9dvEJF9\nIrIk5/a086WqSBZpSRS7NLGiSqvIJi4iMcAooD3QFOgtIo3yGTrbGHNlzu0Fh+tUESwSkyh2aWJF\nlZadmXhzYJ0xJsMYkwmMA7rkMy7yv/8qx0VyEsWu3MTKsNRhfL/pe7fLUT5jp4nXAbYGPN6W81yw\nliKyTESmikgTR6pTES0akih2aWJFlZRTOfHFQIIx5oiIdAA+BxrmNzAlJeXU/aSkJJKSkhwqQflJ\nNCVR7ApMrMzrO49qcdXcLkm5JDU1ldTUVFtjxRhT+ACRFkCKMSY55/EQwBhjXinkZzYBVxlj9gQ9\nb4p6P+WS77+HWbPyHj/yCNSsGbK3G/TtIBb+spDpf54eVQcy7Rjw9QDW7F7D1D5TiY3R8/EUiAjG\nmHx3WdvZnbIQSBSRuiJSDugFTA56g/MC7jfH+uWwB6XyEa1JFLs0saKKo8gmbozJAvoD3wArgXHG\nmHQR6SciD+QM6y4iK0RkKfAa0DNkFStfi+Ykil2aWFHFYeu7mjFmGnBx0HOjA+6/AbzhbGkq0mgS\nxT5dY0XZpTvcoklaGhw6ZN2vXBlatw7bW2sSpfgCEytp96TRoEYDt0tSHqSn3UeT5cth3jzrtnx5\n2N5Wkyglp2usqKJoE1chF21rojhN11hRhdEmrkJKkyjO0MSKKog2cRUymkRxjiZWVEH0wKYKCU2i\nOE8TKyo/OhNXjtMkSujoGisqmDZx5ShNooSeJlZUIG3iylGaRAkPTayoXNrElWM0iRJemlhRoE1c\nOUSTKOGniRUFmk5RDtAkins0saJ0Jq5KRZMo7tPESnTTJq5KTJMo3qGJleilTVyVmCZRvEUTK9FJ\nm7gqEU2ieJMmVqKPNnFVbJpE8S5NrEQfTaeoYtEkivdpYiW66Exc2aZJFP/QxEr00CYeyfbuhY0b\n824nTpRsOydOkLV+HX3+cxstqzRlQMM7nK1ThYQmVqKD7k6JZMuXQ2pq6bezdy9D/3s3h/iFUVyL\nxP8MN9xQ+u2qkHuw2YOs+m0VPSf0ZGqfqcTG6D/5SKMzcVWkMavHMYl0JtCDcpRxuxxVTJpYiWza\nxFWh0rakMWje80yhNzWo6HY5qgQ0sRLZbDVxEUkWkdUislZEBhcy7moRyRSRbs6VqNxyKonS5nUa\nU9PtclQp5CZWhqUO4/tN37tdjnJQkU1cRGKAUUB7oCnQW0QaFTDuZWC600Wq8Dt44lBeEiWhjdvl\nKAdoYiUy2ZmJNwfWGWMyjDGZwDigSz7jHgUmALscrE+5IIts+nz7oK6JEoE0sRJ57DTxOsDWgMfb\ncp47RUTOB24zxrwF6CIaPjeUmRzKPKxrokQoXWMlsjiVN3oNCNxXXuC//JSUlFP3k5KSSEpKcqgE\n5YQxLGMS6cxPnqtrokSwEe1H0OnjTgycPpCRHUa6XY4KkpqaSqrNeLCdJr4dSAh4HJ/zXKBmwDix\npm3nAB1EJNMYMzl4Y4FNXHlLGlsYxLfM4m5qxFV3uxwVQrmJlRbvtmD0otH0a9bP7ZJUgOAJ7vDh\nwwsca6eJLwQSRaQusAPoBfQOHGCMqZ97X0Q+AKbk18CVhxw8CF98cerh5m0r6MFnjKWrJlGihK6x\nEhmKbOLGmCwR6Q98g7UP/T1jTLqI9LNeNu8E/0gI6lROO3YMli4F4CDH6cz7DKEVySS6XJgKp8DE\nSto9aTSo0cDtklQx2donboyZBlwc9Fy+Zw0YY+51oC4VJllk04dJtCSeAVzjdjnKBYGJlXl951Et\nrprbJali0IUUotxQZnKIE4yiI5JzPDozK4YZM8uxaB2UKQMtLi5LkhFiRL9kRSpdY8W/9LT7KJab\nRAlcEyVtSwKXvPUww/9ZkSNHrIUQB6acxSVvPsyszXVdrliFkq6x4k/axKNUbhIlcE2Uz1Y2odv4\nnvy93bfM+3ofL70E//gHLPl2Dy+2nUmvid0ZOU93uUQqXWPFn/Q7UxTazL4zkiipmy+k/9cdmXHn\nWC49byfQ7tR4Ebit0Wour/UrHT66nWMfH2SwrkQbkTSx4j86E48yVhLlk9OSKLuPVOCO/3XlP13/\nl9PA83dhtX3MuGMs73xZmzFjwlSwCjtdY8VftIlHkYKSKI981YkeTVZx80UbitxGnSoHmfLiCp58\nEpYsCWW1yk26xop/aBOPIvklUWZn1GXutnj+1mam7e00ufAIb7wB3bvDgQOhqla5TddY8Qdt4lEi\nvyRKthH+Mr09L7WdSYWyxfhHOmsWf1r7Am2rL+WvbZfA7t15r+3cCS+8kHdLS3P4T6LCSRMr3qdN\nPArkl0QBmLKmIQbofcnPxdtgdjacPMmrN01j+rr6fPNdwPFxY+DkybxbdrYzfwjlCk2seJ828QiX\nXxIFrF77Utp1DG2dRklXm61S/jijb/mShwZV5tgxhwpWnqNXBfI2beIR7ODJI2ckUXLNzqjLnqMV\n6NoovVTvkZy4nsuanuTVV0u1GeVxmljxLm3iESorO4s+6S8UuCbKGwub89g18ygTU/pT6V8dfph/\n/Qu2bi16rPIvTax4kzbxCDV05lAOZR09LYmSa9fhSnyz4SL+fOlPjrxXvbrZPPIIDBrkyOaUh2li\nxXu0iUegMcvGMCl9EhOappxKogT6cNlldG2cTtW4446956BBMGsWLP1ZTwKOdJpY8RZt4hEmbUsa\ng74dxJTeU6hRtuoZrxsD7y69kvuvdPZMnUqV4Omn4amXznJ0u8p7NLHiLdrEI8jmfZvp8VkPxnYd\nS+OajfMdM29bPDFiaBnv/A7svn1h7cYypG6+0PFtK2/RxIp3aBOPEAePH6TzJ50Z0moIyYnJBY4b\nv/ISejVdUeJYYWHKlYPnBx1i6My2GKeXHj96FNLT826//+7wG6ji0sSKN2gTjwBZ2Vn0mdSHlvEt\nGXDNgALHZRvhs1VN6HnJypDV0uu2YxzJLMvkNRcXPbg49u6F8ePzbqtXO7t9VSKaWHGfNvEIMHTm\nUA6dOMSojqOQQqbYaVsSOKfiERqdE7pZbEwMDE9K5bnZNzg/G1eepIkVd2kT97lTSZQeEyhXplyh\nY8etuISeTUM3C89168VrOJkdw1fr9KK70UITK+7RJu5jpyVRKtYodGy2ET5f3YgeTULfxGPE8Mz1\nsxk+K0ln41FCEyvu0SbuU3aSKIGW7KhNlfLHaVBjTxiqg26N0zmSWZbpC84Oy/sp92lixR3axH3I\nbhIl0JdrG9K54ZoQV5bHmo3PYviYujobjyK5iZXeE3uzfs96t8uJCraauIgki8hqEVkrIoPzef1W\nEVkuIktFZIGItHK+VAX2kyjBpqxtyC0N14awsjN1b7KKfYdimWn/ehMqArSp14aUpBRNrIRJkU1c\nRGKAUUB7oCnQW0QaBQ2bYYy5zBhzBXAf8K7jlSrAfhIl0PYDZ7Fp79m0SgjvClVlYgxP37mF4cPR\n2XiUebDZg9xU/yZNrISBnZl4c2CdMSbDGJMJjAO6BA4wxhwJeFgZ0CsBhEBxkiiBvlrXgOTE9cTG\nhP9/S88bd7Fzp7WuiooumlgJDzurFdUBAqdw27Aa+2lE5DbgJaAm0MmR6tQpuUmUWXfPOj2JMn06\n7N+f/w/99hsAX65rWPxo4bRpUL68db8UV3yIjYWnnoLnnoOkpBJvRvlQbmKlxbstGL1oNP2a9XO7\npIjk2JJzxpjPgc9FpDXwAnBTfuNSUlJO3U9KSiJJ/2UXqdAkyoYNsGtXgT+bmRVD6uYLee/WycV7\n0w0bSlBp/m6/3WriP/wArfRoSVTJTay0/qA1DWs05MZ6N7pdki+kpqaSmppqa6ydJr4dSAh4HJ/z\nXL6MMWkiUl9EqhtjzsizBTZxVbSSJFECzd8eT2L1PZxT8UjRg0OkbFkYMgSef96a4KvoEphYSbs3\njcTqiUX/UJQLnuAOHz68wLF29okvBBJFpK6IlAN6AadN60TkooD7VwLl8mvgqnhKmkQJNGNjfdrV\n2+hwZcV3112wahUsWOB2JcoNmlgJnSKbuDEmC+gPfAOsBMYZY9JFpJ+IPJAz7I8iskJElgCvA38K\nWcVRpCRJlGAzNtbnpovcb+Lly8PgwfDCC25XotyiiZXQsJUTN8ZMM8ZcbIxpYIx5Oee50caYd3Lu\n/90Yc4kx5kpjTCtjzNxQFh0NSppECXTgeHmW7zyPVhdscbi6krnvPli8GJYtc7sS5RZNrDhPz9j0\noOKsiVKYWZvrck2d7VQo641ZT1wc/PWvOhuPZrrGivO0iXtMcddEKcyMjfVpV9/9XSmB+vWDtDRY\nscLtSpRbdI0VZ2kT95DSJlGCzdjkvSZesSI88QT87W9uV6LcpGusOEebuEc4kUQJtONgZXYcPIsr\nau1woDpnPfwwzJwJa1YbOHmy4FtxztUP/LmsrNNfy84u+XZVyGhixRmOneyjSseJJEqg2Rl1ua5u\nBmVivNewzjoLHn0UXhx2lA+b/L3ggQMGQPXq9jY6ahTsy2kE558PDzyQ99qMGfDjj3mPH38cqlUr\nfuHKcQ82e5BVv62i54SeTO0zldgYbUnFpTNxD3AiiRJszpa6XJfgjVRKfh59FKbOiGPjXl1vPNpp\nYqV0tIm7zKkkSrA5WxK4LiHDse05rVo1ePiuw7w0p7XbpSiXaWKldLSJu8jJJEqgfcesGe6Vtb23\nPzzQ4/cfZtLqxmTsq+p2KcplmlgpOW3iLnE6iRLohy0X0LzOdsqW8faKwNXPNtx/5RJe+UFn40oT\nKyWlTdwFTidRgnl9f3iggS1/5NOVTdm0Vw80Kk2slIQ2cRc4nUQJ5vX94YFqVjpC/+YLeDZVlyhV\nFl1jpXi0iYdZKJIogY5mxrLs11q0iN/m+LZD5S8t5zJ9w0Ws2HWu26Uoj9DEin3axMMoVEmUQAu2\n1+GSc3dRqVxmSLYfClXKH2dIqzT+77s2bpeiPEITK/ZpEw+TUCVRgvlpf3igh65exNIdtZm7Nd7t\nUpRHaGLFHm3iYRDKJEowP+0PDxQXe5KUpFSGzmynZ8WrUzSxUjRt4iEW6iRKoJPZMczbFk9rH87E\nAe68bDm7Dldi6rqGbpeiPEQTK4XTJh5ioU6iBFr+63nEVzlAjYpHQ/o+oRIbk82rN09n4Dc3cyKr\njNvlKA/RxErBtImHUKiTKMH8uj88UIcG66l/9l7eWHC126Uoj9HESv60iYdIOJIowWZn1PXl/vBg\nI26ezotp1/Hb76H95qL8RRMr+dMmHgLhSqIEMgbStiRwXV1/z8QBGtf8nd6XrODZlyu4XYryGE2s\nnEmbuMPCmUQJtGb3OVQom0lC1f1he89QSklKZcKUcnpRZXUGTaycTpu4g8KZRAk2JyPB9/vDA1Wv\ncJSXnjnCAw+ceaEepTSxkkebuIPCmUQJZh3U9P/+8ED39DlBXBy89ZbblSgv0sSKxda1kEQkGXgN\nq+m/Z4x5Jej1PsDgnIcHgYeMMT87WajX5SZR5vedH5YkSrA5WxIY0jot7O8LwBdfQNmy1v0TJxzb\nbEwMjB4N110HXVvGUMexLUeQ+fNhzZq8x3/6E8TFuVdPmI1oP4JOH3di4PSBjOww0u1yXFHkTFxE\nYoBRQHugKdBbRBoFDdsIXG+MuQx4Afi304V6mRtJlEDbDlTh4PHyND7nt7C/NwBbt8LGjdZtm7ML\nbzVubF1YecAwvXBEvn7/Pe+z37gx6vY9aWLF3u6U5sA6Y0yGMSYTGAd0CRxgjJlnjMk9ojYPomfS\n5EYSJdicjASuq5tBmPfghM1TT8GK1WX5bGUTt0tRHhTtiRU7TbwOsDXg8TYKb9J9ga9LU5RfuJVE\nCRYJJ/kUJi4Oxo7cS/+vO7LjYGW3y1EeFM2JFVv7xO0SkRuBe4ACr7eVkpJy6n5SUhJJSUlOlhA2\nbiZRgs3ZksA9ly91tYZQu+bKTB68aiX3Te7C1D4fEaFfOlQpBCZW5t43l2px/r1aVGpqKqmpqbbG\n2mni24GEgMfxOc+dRkQuBd4Bko0xewvaWGAT9zM3kyiB9hwuT8a+alxR+1fXagiXp6+fzbXv38fb\ni5rx0E1uV6O86MFmD7Lqt1X0nNCTqX2mEhvj6Dw1bIInuMOHDy9wrJ3dKQuBRBGpKyLlgF7A5MAB\nIpIATATuMMZsKEHNvhLuNVEK88PGWlwTv43YGG9fFNkJZctk85+u/2NY6o0sX1/J7XKUR0XbGitF\nNnFjTBbQH/gGWAmMM8aki0g/EXkgZ9gzQHXgTRFZKiILQlaxy9xOogSbs/78iN4fHqzROb8zMnka\n3Z9pxP7IODlVOSzaEiu2vmsYY6YBFwc9Nzrg/v3A/c6W5j1eSKIEm7OhNi9et9ztMsKqzx9+5gfT\nknvvrcCECURsKkeVXG5ipfUHrWlYoyE31ovcC3HrGZs2eSWJEujIEfhpew2uiT/jEIV3nDgBhw/n\n3bKd2e0zov8mtm6Fv//dkc0V7MiR0+sPvGX65zqm0ShaEiv+3OsfZl5KogSaNw8uq7ObimU93EzS\n0qxbrkcfhRql3w1Vvpxh0iRo2RIuugi6l3qLBXj7bThwIP/X2re3ClCeFUmJlYLoTNwGryRRgs2e\nDdcn/uJ2Ga6Jj4fJk60zOudtruV2OcqjIn2NFW3iRfBSEiXYnDnR3cQBrrgC3n8fur7fmXW7q7td\njvKoSE6saBMvhNeSKIFOnIAFC6DVRZGfDy/KLbfA8x1/pN1/7mTT3sj7uqxKL5ITK9rEC+DFJEqg\nxYshMRGqVnBu1UA/69tiJYNb/UDbsXexZY9myNWZInWNFW3i+fBiEiXY7Nlw/fVuV+EtD1+9kAHX\nzOfGVzuzIeJPOVMlEYmJFW3iQbyaRAk2Z4428fw83mIeg9ov57rrrG8rSgWLtKsCaRMP4tUkSqCs\nLPjhB+tiCepM/a5P5803ITkZvvrK7WqUF0VSYkWbeAAvJ1EC/fwznHcenHuu25V41223WRccuv9+\neO6/9ck23vyFrNwTKYkVPdknR24SZdbds7yVRDl6FJYHnFZfpw6zZ1+gu1K2bLHOdsp1/PgZQ669\nFhYtgp4312De3D58eNv/qFnpSHjqO3gQVq7Me1y3LtSuHZ73VrbkJlZavNuC0YtG069ZP7dLKhFt\n4ng8iXL4MEyblvf4+uuZM+cCunQp+Eeiwtq11q0ItWvDzFcW8cyLFfjDWw/zeoev6BGG8ti79/T/\nb8nJ2sQ9KBLWWIn63Sl+SKIEys6GWbP0oGZxlI01vNxuBp/3Gsew1BvpdkdFNm1yuyrlFX5PrET1\nTNwvSZRAKzZVokoVSEgoeqw6XYv4bSztN5p/mIE0awb33mtdv/Pss92urHDGWNdD3rQJ9uyB/fut\nvWxly0LcylrU3nkBdavuo/ZZh3RWVkJ+XmMlqpu4H5Iowb5bdjZt27pdhX/FxZ7kmceP07d/BYYN\nsxbPuuceeOwxb/xiPHAAli614pFLlsBPP1kXsS9fHi68EM45B6pWhQoVrEUUj26qw46t55GxvypH\nM8tyzQ+xtLkJune3/mzKPr9eFcgfVYZAbhJlft/5nk6iBJu5pDp3POF2Ff5Xuzb8+9/wzDMwciRc\nfrkV2bzjDus0/ri40Newf7/VqBcvzmva27bBpZfCVVdBmzbwxBNWM65W0MRw6mJYuBCAXYcrMe/q\nR5mWGkurVlC/Pjz+OHTrBrFR+y+9eEa0H0GnjzsxcPpARnYY6XY5tkTlty8vr4lSmJPZMcz+uRo3\n+u/Yi2clJMCrr0JGBnTtCm+9ZTX4bt2sVWjX7Dq71PFEY2DH72X57jv45z+hd29o2BDq1IGnn7aC\nNsnJMHGi1djnzoVRo6zdPVddVUgDD3JupcPc2imLN9+E7dvhySfh9dehaVMrbmlMqf4YUcGPa6xE\n3e9nTydRirBw+/lceN5RatY8y+1SIs5ZZ8Hdd1u3nTvh22/hm2/g5S+7sedwHJfX+pWGNXaTUHU/\nF1TZT/UKR6l4dlUqGetg87Fj1u3AAfj1V9ixphq/zO3G2t01WPP7OcRViuHiptaqi8nJVvNu1AjK\nlAnNn6dMGeuXUteuVkhm4EDrF9R771m/PFTB/JZYiaom7rckSrDvNtWj7RV7AW3ioXTeefDnP1s3\nRrzHnp2ZLN1Ri417z2bL/qqkZlzIvmNxHN5Qi8PjrIYZF2fdKleGWrWgdvUs/nDRBvpfvYCLz9lN\n9a43QIsWrvx5kpOhbVt46SXrl8gbb0CPsOQs/SswsZJ2bxqJ1RPdLqlAUdPE/ZhECTZzU30GttsL\neOAIXBSpXuEobetvoi1BucTCruyz5SC8751rn5YtC8OGQadO1kHPxYvhb38L3TeBSOCXxErUNHE/\nJlECHc2MZcH2OlxvpsC/F+S9sGePe0WVxKef5h1lO1HEMrrjx9sfG+i336yjlrn27z/99U8+ydtu\nzZrWOfpO+Ogj65qccOYZpHPnWusl5KdWLejc2d577NgBX36Z93hfIQs4/fILTJ2a97hVK666qgkL\nF8Kf/mTt9x8/PjwHcf3KD4kV71UUAn5NogSas6Uul9X6lbMy94CHr4tcpJ07QzM2UGamdWTP6e0W\n5ddfrdPt87N//5m/THIVZzp8/Hjhf7bCxh4+DFgxxWnT4M47oWNH66DnWbqHrkBeT6xEfDrFr0mU\nYF+vS6RDov/OJlPeVK6c9cWhQQNo167ga0Er7ydWbDVxEUkWkdUislZEBufz+sUi8qOIHBORvzhf\nZsn4OYnAlYHpAAALc0lEQVQS7Kv1DejYYJ3bZagIUqaMFaO86iq49VbrLFCVPy9fFajIJi4iMcAo\noD3QFOgtIo2Chu0GHgX+4XiFJeT3JEqgDXvO5sDx8lxeS6+nqZwlYmXS69SxEiuZmW5X5F1eXWPF\nzky8ObDOGJNhjMkExgGnraFnjPndGLMY8MTq6pGQRAn09foGJCeuJ0b0bA3lvJgYGDPGuv/QQ3pS\nUGG8eFUgO028DrA14PG2nOc8y+9JlGBfrWtAx0TdlaJCp2xZGDfOih6++qrb1Xib164KFPZ0SkpK\nyqn7SUlJJCUlObr9SEiiBDp6FNK2JPDxHye6XYqKcJUrw5Qp1jlJiYnOJS8jUagTK6mpqaSmptoa\na6eJb+f0s0viKUXILbCJO82zV+cphZlzynJF7R1UizvmdikqCsTHw+efQ4cOVnKlaVO3K/KmUF8V\nKHiCO3z48IJrsbG9hUCiiNQFdgC9gN6FjHdl/0UkJVHYv9+6GjIw6X+xdGuU7nJBEezkycJPmMrO\nLvi1I0dO/9nCxpamhkClyQIePlzw+8TGQpUqADRrZu1S6dbNWiAx52kVxCtrrIixcRRDRJKBkVj7\n0N8zxrwsIv0AY4x5R0TOAxZhLeqRDRwCmhhjDgVtx9h5v+I6ePwg175/LX2v6MtjLR5zfPth99Zb\nsHMnmVkx1H71ryzpN5qEqgWcKKJUoCefhEqVrPubNsGHH9r7uXr14K67TnvqoYdg1y6YMMFKsaj8\nfbfpO/pM7BPSNVZEBGPyX07TVk7cGDPNGHOxMaaBMeblnOdGG2Peybm/0xhzgTGmmjGmujEmIbiB\nh0qkJVECzcq4kPpn79UGrlzx2muwdase6CyK24kV35+xGWlJlECT0hvzx8ar3C5DRany5a1Z+D//\nCTaPsUUtNxMrvm7iuUmUCT0mREQSJVC2Ef63uhHdGuv+cOWehAQYOxb69LHW3lIFG9F+BAADpw8M\n6/v6tolHypooBUnbkkDNikdoUMNnqxSqiHPzzdCvH/TqZR2DVflza40VXzbxiEqiFOA/yy/lz5f+\n5HYZSgHWlYji4qxrkqqCubHGiu+aeCStiVKQoyfKMDG9Cbf/QZu48oYyZeC//7VWPpwyxe1qvC3c\na6z4qolHchIl0Bc/1ePqOtupU6WAtamVckHNmtap+ffdZ6UXVcHCmVjxVROP5CRKoLHzG3Lnpd65\ntJdSua69FoYOtVY8DL54kTpduBIrtk72cezNSnGyz5hlY3hh9gvM7zs/Ig9k5tq8GZpdcpSMASOo\nVE7XBVXF1KqVtZoVWJduW7bM3s9VqwaXX573uGFDOP/8vMdpaaeOahoD3Z9tynlnn+DNJ9adOVad\ncjL7JJ0+7kSjGo1KtcZKYSf7+OLybJG4JkpB3noL7mqxRhu4KpkffijZz+3bd3oYvHLlM5v4MWv9\nHgHev3Yuzf79AB/H7qbPUzu0iRcg1GusgA92p0RDEiXX0aPw/vvw0HUr3S5FqUJVjTvOhB6f8ti0\nZFZt1CstFybUiRVPN/FoSKIE+vhjaN4cEs/VCx4q77us1k5eaTeD7oPrcygsi2z4VygTK55t4tGS\nRMmVlQX/+Ac88YTblShl371XLKXFJYe5776SL+IYLUKVWPFsE4+WJEquTz+FGjWgbVu3K1GqeN4Y\nvIUtW2DYMLcr8b5QJFY82cQjeU2U/GRlwfPPw7PP6pKfyn8qxBm++AI++QQ++MDtarzP6TVWPNfE\nI31NlPx89JGV8LrpJrcrUapkzj0Xpk6FIUNg5ky3q/E2p9dY8VQTj6YkSq6DB62TJ/71L52FK39r\n1MjaLdi7N8yd63Y13uZkYsUzTTzakii5nnsO2rWDa65xuxKlSu+GG6yla2+7DRYtcrsab3MqseKJ\nk32iLYmS68cfrUWFlusZ9sprvv/+9Ol0Mc6xT84YzbtdzqdTm7Z83ncqLev9mvdi27bQpImDhfpb\nYGJl7n1zqRZXrdjb8MRMPNqSKGCdIHfnnfDGG9b+RKU85fBh2L0771ac5TL27KHz+Yv58NZJdBnd\nkck/npO3nZyzPlWe0iZWXG/i0ZZEASuNcvvtkJxsXVFcqUiUnLieqX0+4sEvb+Ffc1sU6/dAtClN\nYsXVJh6NSRRjoH9/OHLEOpipVCS7us4v/Hjfe3yy4g90+7Qnew+UcbskTypNYsW1Jh6NSZSTJ+GR\nR2DJEvjii7zF5pSKZBdW28ece94nocp+Lu/dSC8qUYCSJlZsNXERSRaR1SKyVkQGFzDm/4nIOhFZ\nJiKX5zcmVzQmUXbtgltvhXXrYPp0qFLF7YqUCp/ysVmM7DCND57NYOBA6NIFfv7Z7aq8pySJlSKb\nuIjEAKOA9kBToLeINAoa0wG4yBjTAOgHvF3Q9ryQREkNXHIzxLKy4MMP4dJL4bLL4KuvrBN7vCJ1\n82a3S/AM/SzyhOqzaNP8ED/9BNdfb53c1r07fPedt9ddCWe/gOKvsWJnJt4cWGeMyTDGZALjgC5B\nY7oAYwGMMfOBqiJyXn4b80ISJRz/U3btgtdft9JU774LkyfDSy95bxeKNq48+lnkCeVnERcHAwfC\n+vWQlASPP25dV2LwYJgzBzI9tpR+uJs4FC+xYicnXgfYGvB4G1ZjL2zM9pzndgZvbFL6JOb3nR8x\nSZSsLNizBzZutHaV/PST9RcxPR06d4a337b+okZJclIp2ypXtg7yP/KIdWLQlCnw2GOwdq31zbVZ\nM6u516tn3WrWtL7Fem0iFCoj2o+g08edikyshP1knzqzpnDn9Bqn4kbGcNr94P+G4rWtW2HaNPs/\nF/xcVhYcOGBlvQ8dgqpVoX59SEyEpk2tGXfz5lCxYsk+I2rWhJgwHXOuXBlq1w7Pe3mdfhZ5ivNZ\nBP9Fr1ULTpywNxZrgnP11dbtueesf1tLlsDixbByJXz5pXXZwt27Ye9eqFDBauYVKkC5cnm3smWt\n/+b+08mdOJX0v7nWrLFqCb9YpMx43mvUotBRRV5jU0RaACnGmOScx0MAY4x5JWDM28D3xpjxOY9X\nAzcYY3YGbUuTokopVQKlucbmQiBRROoCO4BeQO+gMZOBR4DxOU1/X3ADL6wIpZRSJVNkEzfGZIlI\nf+AbrAOh7xlj0kWkn/WyeccY85WIdBSR9cBh4J7Qlq2UUgps7E5RSinlXa6vneImERkoItkiUt3t\nWtwiIn8XkfSck7QmikjUnYZk52S2aCAi8SLynYisFJGfRSR6lhQtgIjEiMgSEZnsdi0FidomLiLx\nwE1Ahtu1uOwboKkx5nJgHTDU5XrCys7JbFHkJPAXY0xToCXwSBR/FrkeA1a5XURhoraJA/8CnnS7\nCLcZY2YYY3LPl5sHxLtZjwvsnMwWFYwxvxpjluXcPwSkY53vEZVyJnodgXfdrqUwUdnEReRWYKsx\nRldvON29wNduFxFm+Z3MFrWNK5eIXAhcDsx3txJX5U70PH3g0BNX9gkFEfkWCDz1X7D+ZzwNPIW1\nKyXwtYhVyGfxf8aYKTlj/g/INMZ87EKJykNEpDIwAXgsZ0YedUSkE7DTGLNMRJLwcI+I2CZujMn3\n2vEicglwIbBcrMVb4oHFItLcGLMrjCWGTUGfRS4RuRvra2ObsBTkLduBhIDH8TnPRSURicVq4P8x\nxnzhdj0uagXcKiIdgQrAWSIy1hhzp8t1nSHqI4Yisgm40hiz1+1a3CAiycCrwPXGmN1u1xNuIlIG\nWAO0xTqZbQHQ2xiT7mphLhGRscDvxpi/uF2LV4jIDcBAY8ytbteSn6jcJx7E4OGvSmHwOlAZ+DYn\nSvWm2wWFkzEmC8g9mW0lMC6KG3gr4HagjYgszfn7EB0L/vtY1M/ElVLKz3QmrpRSPqZNXCmlfEyb\nuFJK+Zg2caWU8jFt4kop5WPaxJVSyse0iSullI9pE1dKKR/7/4+PajaZ7m91AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a4ff76c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin, xmax = -5, 5\n",
    "I = quad(pdftarget, -100, 100)\n",
    "nsamples = 1000\n",
    "xhigh = 0.0\n",
    "yhigh = 0.63\n",
    "xsample, ysample, accept = sampling_isosceles(nsamples, I[0], xmin, xmax, xhigh, yhigh)\n",
    "\n",
    "xaccept = xsample[accept]\n",
    "print \"Accepting rate: \", float(len(xaccept))/len(xsample)\n",
    "\n",
    "x = np.linspace(xmin, xmax, 1000)\n",
    "y_target = pdftarget(x, normed=I[0])\n",
    "y_propos = propos(x, xmin, xmax, xhigh, yhigh)\n",
    "plt.plot(x, y_target)\n",
    "plt.plot(x, y_propos)\n",
    "plt.hist(xaccept, normed=True, bins=50, histtype=\"stepfilled\", color=\"red\", alpha=0.5, linewidth=0)\n",
    "plt.xlim([xmin, xmax])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's try to make it more efficient\n",
    "\n",
    "NOTE: efficiency in accepting rate! not necessarily in computing time (need more testing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sampling_piecewise(nsamples, normalize, xmin, xmax, xhigh, yhigh):   \n",
    "    l1 = 0.5*(xhigh-xmin)*yhigh # left triangle area\n",
    "    l2 = 0.5*(xmax-xhigh)*yhigh # right triangle area\n",
    "    l = l1+l2 # total area\n",
    "    \n",
    "    # start sampling\n",
    "    # piecewise proposal distribution\n",
    "    xsample1 = xmin + (xhigh-xmin)*np.sqrt(np.random.random(int(nsamples*l1/l))) # left triangle\n",
    "    xsample2 = xmax - (xmax-xhigh)*np.sqrt(np.random.random(int(nsamples*l2/l))) # right triangle\n",
    "    xsample = np.append(xsample1, xsample2)\n",
    "    np.random.shuffle(xsample) # shuffle this array (for animation purpose)\n",
    "    \n",
    "    n = len(xsample) # may change a little\n",
    "    target_dist = pdftarget(xsample, normed=normalize)\n",
    "    ysample = propos(xsample, xmin, xmax, xhigh, yhigh)*np.random.random(n)\n",
    "    \n",
    "    accept = ysample <= target_dist\n",
    "    \n",
    "    return xsample, ysample, accept, l, n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.75675675676 0.382082506929\n"
     ]
    }
   ],
   "source": [
    "xmin, xmax = -5, 5\n",
    "I = quad(pdftarget, -100, 100)\n",
    "nsamples = 1000\n",
    "\n",
    "x = np.linspace(xmin, xmax, 1000)\n",
    "y_target = pdftarget(x, normed=I[0])\n",
    "high = np.argmax(y_target)\n",
    "xhigh = x[high]\n",
    "yhigh = y_target[high]\n",
    "print xhigh, yhigh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accepting rate:  0.518518518519\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VEUXwOHfJBCKIL13aQKhCUTpkV5EmnQLKCCGjiIg\nnwqCChZUIEEUpaqAggKSUCUgSgmdhCq9g0Do6fP9sVlTSFmS3b1bzvs89zG7O3vvyQqH2bkzc5TW\nGiGEEM7Jw+gAhBBCZJwkcSGEcGKSxIUQwolJEhdCCCcmSVwIIZyYJHEhhHBiFiVxpVQbpdQRpdQx\npdSYFF5vqpQKV0rtiT/+Z/1QhRBCJJclvQZKKQ9gJtAcuAiEKKVWaK2PJGu6RWv9vA1iFEIIkQpL\neuI+wHGt9RmtdTSwGOiYQjtl1ciEEEKky5IkXgI4l+jx+fjnkquvlNqnlFqtlKpqleiEEEKkKd3h\nFAvtBkprre8rpdoCvwGVrHRuIYQQqbAkiV8ASid6XDL+uf9ore8m+jlIKRWglMqvtb6RuJ1SSjZq\nEUKIDNBapzhkbclwSghQQSlVRinlBfQEViZuoJQqkuhnH0AlT+CJAjH8eP/99w2PwVEO+Szks5DP\nwvE/i7Sk2xPXWscqpYYA6zAl/e+01oeVUq+bXtbfAC8opd4AooEHQI/0ziuEECLzLBoT11qvASon\ne252op/9AX/rhiaEECI9brli09fX1+gQHIZ8Fgnks0ggn0UCR/8sVHrjLVa9mFLantcTQghXoJRC\nZ+LGphBCCAclSVwIIZyYJHEhhHBiksSFEMKJSRIXQggnJklcCCGcmLU2wBLuIC4Obt1KeKwU5M1r\nXDxCCEni4hHcvw9ffZXwOHduePNN4+IRQshwihBCODNJ4kII4cQkiQshhBOTJC6EEE5MkrgQQjgx\nSeJCCOHEJIkLIYQTkyQuhBBOTJK4EEI4MUniQgjhxCSJCyGEE5MkLoQQTkySuBBCODFJ4kII4cQk\niQshhBOTJC6EEE5MkrgQQjgxSeJCCOHEJIkLIYQTkyQuhBBOTAoli9TFxMDevQmPIyONi0UIkSJJ\n4iJ1UVGwerXRUQgh0iDDKUII4cQkiQshhBOTJC6EEE5MkrgQQjgxi5K4UqqNUuqIUuqYUmpMGu3q\nKaWilVJdrBeiEEKI1KSbxJVSHsBMoDVQDeillHoylXZTgLXWDlIIIUTKLOmJ+wDHtdZntNbRwGKg\nYwrthgK/AFetGJ8QQog0WJLESwDnEj0+H//cf5RSxYFOWutZgLJeeMJIb20ax0L2Gx2GECIN1lrs\n8yWQeKw81UQ+YcKE/3729fXF19fXSiEIawo6HsTy4ytZwW12coFptCYrnkaHJYRbCA4OJjg42KK2\nSmuddgOlngEmaK3bxD8eC2it9dREbU6afwQKAveAgVrrlcnOpdO7njBeREwE3gHezGj2KfV/2cFL\n/Eo4EfxMN4qSK6Fh7tzw5pvGBSqEm1BKobVOsXNsyXBKCFBBKVVGKeUF9ASSJGet9RPxRzlM4+J+\nyRO4cB5Tt06lZtGatH2iNXnJzgp60pxy1ONbtiUZWRNCGC3d4RStdaxSagiwDlPS/05rfVgp9brp\nZf1N8rfYIE5hJydunGDGzhnsfT1h4ysPFBPwpS7F6chiPuBZXqeO3PwQwgFYNCautV4DVE723OxU\n2r5qhbiEAbTWDA0ayugGoymVpxTcv5/k9eeoxF+8SmeWEMIF/HUPshsUqxDCRFZsiv/8duQ3Toef\nZmT9kam2qUgBttOfO0TR+IE/Z2+dtWOEQojkJIkLAO5F3WPE2hH4t/PHy9Mrzba58GIJL9A9S02e\nnvM0m05tslOUQojkJIkLACZvmUyj0o14ttyzgKn+w649Huy/XITYuIdHvxWK0V7PsqjzInot68W0\nbdOQmUdC2F+6UwytejGZYuiQDl87TOO5jTn4xkGK5S7GwoXw9ttQqEAcUZdvEBGThc9braVr1cNJ\n3xg/xfBM+Bm6LO1CpQKVmNNhDo95PWbMLyKEi8rsFEPhwrTWDA4czLtN3qVY7mJMmQITJ8Lvv8OB\nnREcGTKTRV2W8/aGlnywuWmK5yiTtwxb+20lm2c26n9XnxM3Ttj5txDCfUkSd3OLQxdz48ENBvsM\nZsUKCAiAP/+EOnUS2jQqfZa/Xv2en0K9CQipl+J5cmTNwdyOcxlUdxANvm9A0PEgO/0GQrg3SeJu\n7FbELd5a/xaz2s/idngWBg6EJUugWLGH2xbNdZdVvX5iQrAvOy+UeLgBpq98fvX8WN59Of1X9WfS\n5knE6Tgb/xZCuDcZE3djI9aM4E7kHb7r+B1DhkBcdCwBTZckNIiNhRNJh0Z+OVSV8X8048CgWWTL\nlzPpsvt792DFCgAuxtyk26WvKOiZmwVF3iCPZ86HA3j8cXjuOVv8akK4lLTGxKXavZvaf3k/Px78\nkTC/MM6fhx9/hGMHomHOsTTf90LVQyw6UIOpfzXivef2JH0xJgaOmd5fHNhED0ayBp+7Y/mVHlSl\nUNL2BQta8TcSwj3JcIobitNx+AX6MenZSRR6rBDTpkG/fpbn1Oltg/hqx9OcC8+dZjsvPPGnPeNo\nRFPm8QuHrBC9ECIxSeJuaP6++cTExdD/qf7cvAnz5sGoUZa/v3SeWwx4ag+T1j9jUfu+1GINfXiL\ndYxlA7HIOLkQ1iJJ3M3ceHCDcRvHEdAuAE8PT376CVq3hhIp36tM1dsN/2L5wYr8849l7etQnBAG\nEMJF2vID17mf/puEEOmSJO5m3tn4Dl2rdKVOcdMcwu+/h1czsGVZ/hwPGNZoDx98YPl7CvEYa3mR\nmhShLt+yN0a2tRUisySJu5GQCyGsOLqCyc0mA7B/P1y7Bs2aZex8wxrv5fff4fx5y9+TBQ8+pRVT\naUGr2wEs3L8wYxcXQgCSxN1GbFwsb6x+gynNp5AvRz4AfvgBXnwRPDNYdS1vjkheeQWmT3/093an\nGpseH8oHWz5gWNAwomOjMxaEEG5Okrib+Gb3N+TMmpOXa74MgNawfDm88ELmzjt8OHz3Hdy+/ejv\n9c5SnJABIZwKP0WzBc24fPdy5oIRwg1JEncDV+9d5f3g9/Fv549SpvUCBw6Y1vLUqpW5c5ctCy1b\nwty5GXt/3ux5WdFzBc3LNafet/XYdm5b5gISws1IEncDb69/m5dqvET1ItX/e275cujaFZQVaqz5\n+cHs2abefUZ4KA8m+E5gVvtZdFzckdm7Zsu2tkJYSJK4i9t6disbTm5ggu+EJM+vWAGdOlnnGo0b\nm/77598ZHFyP91yl5/jr1b+YsXMG/Vf2JyImwgrRCeHaJIm7sJi4GPxW+zGt9TRyZ0tYXXn5Mpw5\nA89YtlYnXUrBoEEw6/u0KwJZomKBimzvv507UXdoPLexlH8TIh2SxF3YjB0zKJKrCN2qdkvy/IYN\npmmFWay4c85LL0HQ+qxcvZf5ghC5vHKx5IUldK/aXcq/CZEOSeIu6sLtC3z454fMbDvzv5uZZuvW\nQatW1r1evnzQqX00C/bXtMr5lFKMbjiahZ0XSvk3IdIgSdxFvbnuTV6v8zqVC1ZO8rzWtkniAH17\nRzF/f80M3+BMSYsnWrCj/w5+OPgDvZf35l7UPeudXAgXIEncBW08uZHt57czvsn4h147eNBUGrNc\nOetft0nDWO5EZmPv5RSqSmSClH8TInWSxF1MZEwkgwMHM73tdHJmfbgQw4YNpnndtuDhAS/X3M/8\nfdYZUklMyr8JkTJJ4i5m2rZpVCxQkecrP5/i61u2QNOU6x1bxcs19/NTaHWiYjM33TAlUv5NiIdJ\nEnchp8NP89m2z5jeJuXNTOLiYOtWaNTIdjFUyH+DSgWuE3S8gs2u0bB0Q0IGhLDmxBo6L+nMrYhb\nNruWEI5OkrgLGbFmBCOeHkG5fCkPeB89aipr+ah7hz+qV2ruY/7+TK7nT0fx3MXZ9MomSuYuic8c\nHw5dk6pBwj1JEncRq4+tJuxaGKMbjk61zZ9/2rYXbta9WhgbT5XjxoMcNr2Ol6cX/u39GddoHE3n\nNWXZoWU2vZ4QjkgKJbuAB9EPGBo0lK+f+5rsWbKn2u7PP6FJk0RPREXBlClWjydP9khalT/BskNV\nGFBnT/pvyKS+tfpSvXB1ui7tSsjFED5s9iGeHtYfkxfCEUlP3AVM2TqFOsXr0Kp82pO/UxwPj4tL\nelhJL+9Qfgqtnn5DK6lTvA4hA0IIuRhC2x/acv3+dbtdWwgjSRJ3csevH8c/xJ8vWn+RZrsLF+DO\nHXjySfvE1a7icfZeLsrFO7nTb2wlhR4rxNoX11KzSE3qfluXvZf22u3aQhhFkrgT01ozNGgoYxqO\noeTjJdNsu3MnPP20dbaetUT2LDF0rHyUpWHV7HPBeFk8svBpq0+Z2mIqrRa1kvJvwuVJEndiyw8v\n59ztc4x4ZkS6bXfuhHr17BBUIr28D/JTqLd9Lxqve7XubHplk5R/Ey5PkriTuht1l5FrRxLQLoCs\nnlnTbR8SAj4+dggskeZPnOJ0eF5O3Mhn3wvH8y7sTciAEE7ePCnl34TLsiiJK6XaKKWOKKWOKaXG\npPD680qp/UqpvUqpnUqphtYPVST2weYPaFq2KU3Lpr/8Mi4Odu2yf088i0ccL1Q5ZNcbnMnlzZ6X\nlb1WSvk34bLSnWKolPIAZgLNgYtAiFJqhdb6SKJmG7TWK+PbVweWAlVsEK8Awq6GMXffXA6+cTDt\nhnFxsHAhxy89Tr6sLSi05jfTxt/Wcv8+zJ+f8Dgm5qEmvasfZODvHRjfeIvdxuOTM5d/q1u8Lh0X\nd2TSs5MYWGfgQ1v0CuGMLOmJ+wDHtdZntNbRwGKgY+IGWuv7iR7mAmRDCxvRWjM4cDDvNXmPormK\npv+GU6cI2ampV/isqZyPNcXGwqlTCce5cw81qV/qPHejvDh4tYh1r50B5vJv03dOl/JvwmVYksRL\nAIn/dp6Pfy4JpVQnpdRhYBXwqnXCE8n9ePBHbkfe5o16b1j8np0XSuBT4oINo0qdh9L0rBbKTweN\nucGZXMUCFdnRfwd3ou7QZG4Tzt16+B8eIZyJ1W5saq1/01pXAToBk611XpHgVsQtRq8fTUD7ALJ4\nWL7YNuRiceoVNyaJA/SqHsriMG+rFovIDHP5t25Vu+Ezx0fKvwmnZkkmuACUTvS4ZPxzKdJab1VK\nPaGUyq+1vpH89QkTJvz3s6+vL76+vhYH6+7e3fQu7Su255mSllc4jor15MCVItQpfsmGkaWtZpHL\nZM8Sw/bzJalf6rxhcSRmLv9Wu1htei3rxdsN32bkMyNlnFw4hODgYIKDgy1qa0kSDwEqKKXKAJeA\nnkCvxA2UUuW11ifif34K8EopgUPSJC4st/fSXpaELSHML+yR3nfoWiHK5g0nl1cURm2Vo1TCMnxH\nSeJm5vJvXZZ2IeRiCHM6zOExr8wXexYiM5J3cCdOnJhq23SHU7TWscAQYB0QBizWWh9WSr2ulBoY\n36yrUipUKbUHmAF0z3j4Irk4HYdfoB+Tn51MwZwFH+m9+y4XpXZR4+dH9/I+yNKwasTEOd7SBCn/\nJpyZRV0zrfUaoHKy52Yn+vkT4BPrhibM5u6di9aa15567ZHfu+9yUWqZk3hsLGxKNP4bG2ulCNNX\nscANSuW5xaZTZWlZ/qTdrmspc/m3Wbtm0eD7BszrOI+2FdsaHZYQ6ZKtaB3c9fvXeeePd1jTZw0e\n6tF7sfsuF+W5SsdMD7SGzZutHKHlensf5MfQ6g6ZxCGh/FvNIjXp/kt3BtUZxPgm4zP0uQthL/Kn\n08GN2ziO7lW7U7tY7Ud+r9amJF6ziPHDKQA9vMNYceRJImIcu++QuPxblyVdpPybcGiSxB3YjvM7\n+P3Y70xqNilD7z9zBh7ziqLQY/fTb2wHxXPfoVbRywQer2h0KOkyl38rkbuElH8TDk2SuIOKjYvF\nL9CPT1p+Qt7seTN0jn37SBgPdxC9qxu3s+GjkvJvwhlIEndQX+/6mtxeuelTvU+Gz7Fvv6KWgwyl\nmHWtcoh1J8pzOzKb0aFYrG+tvgT1CWLUulGM3TCW2Dj73RAWIj2SxB3QlbtXmLB5Av7t/DO1+MQR\ne+L5ckTgW/Y0vx62U4khK6lbvC67BuyS8m/C4UgSd0Cj14+mb82+VCucuao4+/Y7XhKHRPU3o6Ph\n5MmE47yNFwI9eJD0epcebRWrlH8Tjsixpwm4oS1ntrDp9CYODz6cqfPcvAnXr0P5/DetFJn1dKh0\nlEG/P8fVizEUXrAg4YXChcHPz3YXvnQJEl+vfPlH3prXXP6tXol6tFrUimmtpvFSTStu7yvEI5Ke\nuAOJjo3Gb7UfX7T+glxeuTJ1rv37oUYN0y6CjuYxr2jaVzrGz2FVjQ4lw6T8m3AUksQdyPQd0ynx\neAm6Vuma6XPt2wc1a1ghKBsxL/xxZonLvzVf0FzKvwlDSBJ3EOdvn+fjrR8zs+1Mq+ykFxoKNao7\nXi/crGX5kxz9tyCnwzM2fdJRmMu/NSvXTMq/CUNIEncQo9aO4o26b1CxgHUWwoSGgrcDT8f28ozl\nhaqHWOwkc8bTYi7/Nqv9LDou7sjsXbPRjrJ5unB5ksQdwLoT6wi5GMK4xuOscr64OAgLg2qZm9xi\nc728nWfhjyUSl38bsGqAlH8TdiFJ3GCRMZEMCRzCjLYzyJk1p1XOefYs5MkD+fJZ5XQ207jMWa7f\nz0nY1UJGh2I15vJvtyNvS/k3YReSxA322d+fUaVQFZ6r9JzVzunoQylmHkrTy/sgiw448B3YDJDy\nb8KeJIkb6NTNU0zbPo2v2nxl1fM6w1CK2Su19rPwQE1i41yrLJq5/NvCzgvptawX07ZNk3FyYROS\nxA00fM1wRj0zirJ5y1r1vM7SEwfwLnyVYrnvsPHUE0aHYhMtnmjB9v7b+eHgD/Re3pt7UfeMDkm4\nGEniBll1dBVHrx/lrQZvWf3czpTEAfrW3Me8fbWMDsNmyuYtK+XfhM1IEjfA/ej7DFszDP92/mTL\nYt3d/GJi4OhRqOpEiyF7eocSeLwi4fe9jA7FZszl316v8zoNvm9A0PEgo0MSLkKSuAE++vMjfEr4\n0OKJFlY/94kTULQoPOZEBdsL5HxAy/InWLqngtGh2JRSisE+g1nefTn9V/Vn0uZJxOk4o8MSTk6S\nuJ0du36Mr3d9zbRW02xy/rAw5xpKMetbcx/ztjvX9rQZJeXfhDXJLoZ2pLVmSOAQxjUaR4nHS1jn\npHPmmIppxgu9NMApk3jrCifov7oTR49C5crxT16/DsuXJzQqWBA6dzYkPmszl38buWYkPnN8+LXH\nr1Qt5ERjYMJhSE/cjn459AuX7l5i2NPDrHfSCxeSHM52U9Msi0ccL9Y7xvz5iZ6Mjk76+127Zlh8\ntiDl34Q1SBK3kzuRdxi1bhQB7QLI6pnVZtcJDXWeOeLJ9at/hHnzTLnbnSQu/zZuwzgp/yYeiSRx\nO5m4eSLNyjWjcZnGNrtGZIwnJ08mGo5wMlWL3aR8efj9d6MjsT9z+bedF3dK+TfxSCSJ20Ho1VDm\n75/PJy0+sel1jl0vQLlykD27TS9jU4MGwddfGx2FMaT8m8gISeI2prVmcOBgJjSdQJFcRWx6rdCr\nhZ1yPDyxrl1hzx7TVEl3ZC7/NqX5FFotasWiA4uMDkk4OEniNrbowCLuRt1lUN1BNr+WKyTx7Nnh\nlVfg22+NjsRYPbx7sOmVTUzcPFHKv4k0SRK3ofCIcN7e8Daz2s/C08PT5tcLvVbYaW9qAnD3Lqxb\nx8CqW5k7O5LILTsyd74bN2DduoRj1y7rxGknUv5NWEKSuA3974//8Xyl5/Ep4WOX6zl9T/z+ffj7\nbyqd3YB3vgv8ujSTvc/bt+HvvxOOQ4esE6cdSfk3kR5J4jay59Iefj70Mx81/8gu17sXlZWLd3JT\nwUVWrg+qs4tZu+oaHYZDMJd/C2gXIOXfxEMkidtAnI7Db7UfHzf/mAI5C9jlmof/LUTlAtfJ4iJr\ncDs9eYSTN/Ox91JRo0NxGB0qd5Dyb+IhksRt4Ls93+GhPOhbq6/drhl6tTDVCl+12/VsLatnHEPq\n7eSL7fWNDsWhSPk3kZwkcSv79/6/jP9jPAHtA/BQ9vt4Q68WxruQ6yRxgIF1drPqWCUu3clldCgO\nJXn5t+DTwUaHJAwkSdzKxm4YSy/vXtQqat8iB6FXC+PtQj1xgHw5IujtfRD/EPvcGHYmicu/9fyl\np5R/c2MuMoLqGLad20bg8UAODz5s92uHXStEtcLX4MiRhCdd4C/1sKd30Hjuq4xvvIUcRgdja+Hh\ncDnRNMJ8+aBI2gvEzOXfui7tSsjFEOZ0mMNjXk60mbzINIt64kqpNkqpI0qpY0qpMSm83lsptT/+\n2KqUqm79UB1bTFwMfoF+fNbqM/Jkz2PXa4dHZOfmgxyUzRsOixcnHEuW2DUOW6hc8Do+JS6wYH9N\no0OxvX/+Sfr/z8J57VL+zb2lm8SVUh7ATKA1UA3opZRKvnv/SaCJ1romMBlwu/V2s0JmkS97Pnp5\n97L7tQ9dK0SVQtfwUM7f807JmIZb+eTvhsTEKqNDcVhS/s19WdIT9wGOa63PaK2jgcVAx8QNtNbb\ntdbm8iTbAStVPHAOl+5cYuLmicxsNxOl7J9owq4Woloh19prO7HGZc5S8vHbLA4pb3QoDk3Kv7kn\nS5J4CSDxPKbzpJ2k+wNu1Q0YvX40r9V+zbDKLGHXClPNxWamJDe+8Z98tKY2cZKT0iXl39yLVWen\nKKWeBfoBD42bu6rg08FsObOFd5u+a1gMYdcKudzMlORaPnGCXNmi+fVXoyNxDubyb8VzF8dnjg+H\nrjnflgPCMpbMTrkAlE70uGT8c0kopWoA3wBttNY3UzvZhAkT/vvZ19cXX19fC0N1PNGx0QwOHMyX\nbb4kl5dxc5lNC31cdzgFQCn4X7u9TPiwNV26mB6LtHl5ehHQPoC5e+fSdF5Tvm7/NV2rdjU6LGGB\n4OBggoODLWprSRIPASoopcoAl4CeQJK7d0qp0sAy4CWtdZq3xhMncWf35fYvKZ2nNJ2fNK547437\n2bkX5UWpx13/K3OHGmd47y9YsQI6dTI6GufRr3Y/qhepTtelXdl1cReTm022y66aIuOSd3AnTpyY\natt0h1O01rHAEGAdEAYs1lofVkq9rpQaGN/sXSA/EKCU2quU2pnx8J3DuVvnmPrXVGa0nWHIzUyz\nsGuFqVromlv0TJWCjz6C8eMhVspQPhIp/+a6LBoT11qv0VpX1lpX1FpPiX9uttb6m/ifB2itC2it\nn9Ja19Zau/wSu5FrRzK43mAq5Dd220DTzBTXHg9PrG1bKFAAFiwwOhLnI+XfXJMsu8+Atf+sZe/l\nvYxtNNboUAi75nrL7dOiFEydCu+/DxGyid8jk/JvrkeS+COKiIlgSNAQZrSdQY6sxi8ED71ayOVv\naiZXvz489RT4+xsdifOS8m+uQ5L4I/rkr0/wLuxNu4rtjA4FgLCrrj9HPCUffWTqkf/7r9GROC8p\n/+YaJIk/gpM3T/LVjq/4svWXRocCwLVrEBXrSfHcd4wOxe6qVoXevU03OUXGSfk35ydJ3EJaa4YF\nDeOt+m9RJm8Zo8MBICwMvAtfdYuZKSmZMAFWrnS6+scOR8q/OTfZitZCK4+u5MTNEyzvsdzoUP4T\nFobbjYcnljfiMh+1OMSQblX5+83leBQvCl26GB1W6o4cgT/+SHhcpQo8+6xx8SRjLv/WaUknQi6G\nMLPdTLJnyW50WCId0hO3wL2oewxfMxz/dv54eXoZHc5/wsJwy/Hw/8TE8MoTf+IRE83XQWXgZqoL\nhR1DRARcvZpw3HG8YTBz+bdbkbek/JuTkCRugQ///JD6perTrFwzo0NJIjTUvXviAB5KM+f5lby3\n6VnOXJcybtaQyysXS19YKuXfnIQk8XQc+fcI3+z+hs9bfW50KEloLT1xs6qFrjGq/jYGzG/kCsWM\nHIKUf3MeksTToLVmSOAQxjceT/HcxY0OJ4krV0wLX4rkumd0KA5hdIO/uH43G3PmGB2JazGXf1t0\nYBG9l/fmXpT8eXM0ksTTsDRsKdfuX2Po00ONDuUhYWFQrZrs5meW1TOOBa9t5p134LD9S5y6tLJ5\ny/LXq3/h5ekl5d8ckCTxVNyOvM2b694koF0AWTwcbxKPOYmLBNVKhPPRR9CjBzx4YHQ0riVH1hzM\n6zhPyr85IMfLTg5iQvAEWpZvScPSDY0OJUWhoVCzJuBOKxbDw+H33xMe33v4q33//rBxI4wcCV8b\nv7WNfZ04kfRrSIUK8GTycrgZZy7/VqtoLbr/0p1BdQYxvsl4PJT0BY0kSTwFB68cZNGBRYT6hRod\nSqrCwqBPHyDY6Ejs6P79dFf2KAWzZ8PTT8M3pXMxMM3WLuby5aSfT44cVk3iZubyb91+7sbuS7uZ\n32k+ebLnsfp1hGXkn9BktNb4BfrxwbMfUPixwkaHkyKtTT1xb2+jI3FMefLAqlXw7rR8BJ8ua3Q4\nLknKvzkOSeLJLNi/gIiYCAY8NcDoUFJ15gzkzm3aV1ukrGJF+Gn6NXr+8gLHrssHZQvm8m9jG46l\n6bymLDu0zOiQ3JIk8URuPrjJmA1jmNV+lkOXrzpwAGrUMDoKx9esQQQfN99Ay4UvcSZcvu7bSr/a\n/QjqE8SodaMYt2EcsXFSdsmeJIknMv6P8XR+sjN1i9c1OpQ0SRK3XL/a+xj1zDZaLHyZS3dkRaet\nSPk340gSj7fr4i6WH17Oh80/NDqUdO3fHz8zRVhk+DM76FdrH03n9eN0eF6jw3FZUv7NGJLEgdi4\nWPxW+zGlxRTy58hvdDhpCw3lQEgkNbIfg4MHkXXmidy/b/pMzMfJk/+99E7jPxnis5PGc/sRetUx\nb1i7Ain/Zn8yxRCYs2cOXp5evFzzZaNDSdf9Jas4d+EtKu1ZDPvjjA7HsVy/DstSv7k27OkdFMx5\nn2bzX2Fep99oV14+P1vp4d2DaoWr0XlJZ0IuhPBZq8/I6pnV6LBcktv3xK/du8a7m94loH2AUyxa\nCL1amCc7wq89AAAWh0lEQVQL/ktWT0lAGdG7+kF+7bGYgas6MPHXGsQ50ccYEwNnz0JICPz9N2zd\nCocOwa1bRkeWMnP5txM3T0j5Nxty+574mA1j6FO9DzWKOMedwgNXClOjyBWjw3BqDUufI2TAN/Rc\n/RIbfeG770xTEh1JVKwn2w7mZ0eYKWnv2QPnzkHhwlC0KHh5mRY2Xb8O589DoUJQv0Ilns9ziQ6V\njvKYl2MUPjaXf/tg8wfU+7YeS19YSv1S9Y0Oy6W4dRL/6+xfrD2xlsODnWfHpAOXi0gSt4Jiue/y\nx7j1zLj5IvXrw+jRMHw4ZDeokI3W8M/ZbKzd6cPaE+XZcqYMlUtH0KAddOwIkybBE0+YkndycXFw\n7BhsnXeL+ctqMjiwHf1r72GUt6aI/X+Vh5jLv9UpVoeOizsy6dlJDKwzECW7t1mF448f2EhMXAx+\ngX583upzHs/2uNHhWGz/5SLULCJfS63B00MzYgRs3w7btkHlyrBgAUTbqRP7INKDwEAYMgTKlwff\nQZXZc6kYL1Y/wMlhX7Hzi7/48kt48UXT6vmUEjiAh4fp9f4drhDU5wf2DJzNvWgvvPs/zfTpEOsg\n07bN5d+m75zOgFUDiIiJMDokl+C2Sdx/pz8FcxakR7UeRodiMa1lOMUWKlSA336DH36A77839Xin\nTIFLl6x7Ha3h+PX8+O+sR/sfe1PkldZMnQqlSsGKFXB+9QG+77iCHt5hFMiZ8W0Yy+S9xcx2gWye\ntodffoGWLa3/u2SUlH+zPrdM4hfvXGTSlkn4t/N3qq90585B9iwxFHrsvtGhuIaYGNNdwfijUe17\nBAfDypVw9ChUrQq+vvDFF6Yx6Uft0d69EcX2DXfx/+wBPbpEUbxlVZoteIWQiyXoW3MfZ+esZ/Nm\nGDMGqle3/t7wVcvcY9Mm0+9Qt266e4fZjZR/sy63HBN/a91bDHhqAE8WtP4Ob7Z04ADULCq9cKs5\nc8aUoc0qVoQ+fahdG+bOhVmzYM0aCAqCb781/SNaubLpKFoU8uY17WETFweRkaa6xxcvwoULpinq\nly9moUq+f6lV9DLty5xh6otnKJMnPCFZ53rK5r+ipye8955phW/btrBoEbRubfPLpstc/q12sdr0\n/KUnYxqOYcQzI5yqU+Uo3C6J/3HqD/4+9zffdvjW6FAe2YEDUKOI1NS0l+zZoVMn0wFw86aph37s\nmKlYfXg4XLsGWbKYxqtz5YImTaBECShbFircD8VzxXJDfwezTp1MM1g6d4bFi6GZg9T8Npd/67Kk\nCyEXQ/i2w7c85vWY0WE5FbdK4lGxUQwOHMyXbb50yj8oe/ZA12JyU9Mo+fLBM8+YDoscsGk4j6xh\nQ/j5Z+jWzTQGX99BZvqZy78NWj2I+t/V59cev1I+f3mjw3IabjUm/sW2L3gi3xN0rNzR6FAyZPdu\neKq4JHGRcU2bwrx50KULnD5tdDQJpPxbxrlNEj976yyf/P0J09tMd8pxt+vXTUfF/LI7nMicdu1M\nN1M7doS7d42OJoG5/Nuy7svov6o/k7dMJk470ZJag7hNEh+xZgTDfIY57de0vXuhdm3TnGAhMmv4\ncNOMlf79HW8PtUalGxEyIISgf4LosqQLtyIcdF8BB+EWKSHoeBAHrhxgTKMxRoeSYbt3Q506Rkch\nXIVSMHOmae+VuXONjuZhycu/Hb7mPKuq7c3lb2w+iH7AkKAh+LfzJ3sWg9ZUW8Hu3aavv5xMt6nI\nqFOn4MsvU3+9fHno0CHh8bFjEBiYevuoKOvFZgM5cphmqjRtarrpWfnfv0wbtZg1bAj16hkWn7n8\n29y9c2kyrwmzn5tNlypdDIvHUbl8Ep/611RqFa1FmwptjA4lU3bvhokTkSRuSzExpnmDqbl3L+nj\n6Oi02zuBqlXhgw/gpZdg2+QIPBP/PpGRxgWWSL/a/ahepDpdl3Yl5EIIk5tNdujyifZm0XCKUqqN\nUuqIUuqYUuqhMQmlVGWl1N9KqQil1Cjrh5kxJ26cYObOmXzZOo3elRO4edM0L7lSJaMjEa5o0CDT\nHPcvl5c2OpRUSfm31KWbxJVSHsBMoDVQDeillEq+1PE6MBT41OoRZpDWmqFBQxndYDSl8pQyOpxM\n2bsXatUyrb4TwtqUMq1I/finMvxzw3ErW0n5t5RZ0hP3AY5rrc9oraOBxUCSidZa63+11ruBGBvE\nmCG/HfmN0+GnGVl/pNGhZJrc1BS2Vr48vNP7NANWdXC42SqJSfm3h1mSxEsAibcaOx//nMO6F3WP\nEWtHENA+AC/PVPbvdCI7dxp6f0m4ieGdz3ErIjs/hVY3OpR09fDuwaZXNjFx80SGBw0nOtYximAY\nwe43NidMmPDfz76+vvj6+lr9GpO2TKJR6Ub4lrX+uY2wbRtMnWp0FMLVeXrCjLaB9PilG89XPkou\nowNKh7n824vLX6T5guYs7baUormKGh2WVQQHBxMcHGxRW0uS+AUg8R2PkvHPZUjiJG4Lh68dZs6e\nORx846BNr2Mv586ZZqqVK2d0JIKLF+HXXxMeZ3ZmyunTSc93PZ2bdf/8AwcT/bmuUMG0h21qjh6F\n27cTHleubJqOkoaGpc/RrNwpJm9pwpTSoaYdvlJTvbopBgMlL//2c7efeaakpZvbOK7kHdyJEyem\n2taSJB4CVFBKlQEuAT2BXmm0N2xNu9aawYGDea/pexTLXcyoMKxq+3bTRkVOuFOA67l9G/bvt975\nbtwwHZb699+k18+ZM+0kfvWq6TDLkyfdJA4wtcV6qs/y49Xae6l0OY3ft1gxw5M4JC3/9vxPz7td\n+bd0x8S11rHAEGAdEAYs1lofVkq9rpQaCKCUKqKUOgeMBMYrpc4qpez+bWxx6GJuRtzEr56fvS9t\nM9u2PcKueUJYQbHcdxnbaCuj1jrAxuOPoEPlDmx9davblX+zaJ641nqN1rqy1rqi1npK/HOztdbf\nxP98RWtdSmudV2udX2tdWmtt1611bkXc4q31bxHQLoAsHq6zhsncExfCnob67ODQtUJsOlXW6FAe\nSaUCldyu/JvL7J3yfvD7tCnfhvqlXCfjRUaavj3LzBRhb9myxPJR842MXt+KOO1cwxLm8m8vVH3B\nLcq/uUQS3395Pz8e/JEpLaYYHYpV7d1rWqX5mPPVrxAuoHu1MJTSLAmtZnQoj0wpxdsN32Zh54X0\n/KUnX2z7Au3IE+AzwemTeJyOwy/Qj8nNJlPosUJGh2NVMpQijOShNJ+2XM87fzQnMsY5lwuby78t\nPLCQPsv7cC/qXvpvcjJOn8Tn75tPTFwM/Z/qb3QoVvfnn9CggdFRCId15YppOa/5OH8+7faXLiVt\nf/FiupfwLXsa78JX8Q/xsVLQ9mcu/5bVMyv1v6vPiRsnjA7Jqpz6DuCNBzcYu3Esgb0D8VBO/+9R\nEnFxsHkzfPWV0ZEIh3XypOmw1PHjpsPM1xeKF0/3bVNbbKDpvL70rbWP/DkePHqcDsBc/i0gJIAG\n3zdgXsd5tK3Y1uiwrMKpM987G9/hhSovUKe4620sEhZmKsxbsqTRkQh3V7XQNTo/eYSP/2xkdCiZ\n4qrl35w2iYdcCGHF0RVMbjbZ6FBsIjjY1FESwhFM8A3m+321OROex+hQMs3Vyr85ZRKPjYvljdVv\nMLXFVPLlyGd0ODYhSVw4kuK57+BXN4R3NzUzOhSrcKXyb06ZxL/Z/Q05s+bkpRovGR2KTZjHw5s2\nNToSIRKMbvg3606UZ99l19hkylz+bWzDsTSZ14Tlh5cbHVKGOF0Sv3rvKu8Fv4d/O3+X3RshNBTy\n5pXxcOFYHs8WybtNNjNmQwujQ7GqfrX7EdQniJFrRzJuwzhi42KNDumROF0Sf3v927xc42WqF3H8\nPY8zau1aaNUq/sG//yY9XHTBgkhBRETS//fJa3waYGCd3Zy8mY/1J54wOhSrSlz+rd2P7Zyq/JtT\nTTHcenYrG05u4PBg5x2/skRQEIw0FySaPdtUkFe4n4MHk2496wCyesbxcfONvL2hJbvfOON8vcA0\nmMu/jdswjnrf1mNZ92XULlbb6LDS5TT/D2LiYvBb7ce01tPInS230eHYzJ07EBICzVzj/pFwQV2r\nHCKbZyw/ritodChWZy7/9nHzj52m/JvTJPEZO2ZQJFcRulXtZnQoNrVxo2nrWdkvRTgqpeDTluv4\n3+xSRLjobq/OVP7NKZL4hdsX+PDPD5nZdqbL3sw0CwqCtq6xkEy4sMZlzlKz4j38/Y2OxHbM5d/+\nufkPzRc05/Ldy0aHlCKnSOJvrnuT1+u8TuWClY0Oxaa0hsBASeLCOXz8xjmmToWbN42OxHbyZs/L\nql6raFauGfW+rcf289uNDukhDp/EN57cyI4LOxjfZLzRodjczp2mYZQnnzQ6EiHSV7XcAzp1go8/\nNjoS2zKXfwtoF8DzPz3P7F2zHWpbW4dO4pExkQwOHMxXbb4iZ9acRodjcz//DN26ST1N4TwmTIDv\nvoOzZ42OxPYctfybsue/KEop/SjX++jPj9h2fhureq2yYVSOQWtTRfuV3RdRo2CiLULv3zcuKOHa\nfH2T7u2wcaNp/2NLZc0KWbMyYX1DDl0pwNK596CO621Gl9zdqLv0W9GPM+FnWNZ9GaXylLL5NZVS\naJ1yiSWH7YmfDj/N59s+Z3qb6UaHYhchIZAtG1TPe86UuM2HEI4qOhru32eMzyZ2nSvCur/cY0qV\no5V/c9gkPmLNCEY8PYJy+coZHYpdLF0qQynCOeXIGsP0tkEM+awMkZFGR2MfjlT+zSGT+Opjqzl0\n7RCjG442OhS7iI6GRYvgJdfcz0u4gecqHaNK2Qg++8zoSOzLEcq/OVwSfxD9gKFBQ5nZbibZs2Q3\nOhy7WL0aKlaEyq49g1K4uC9HnuGLL+Cff4yOxL4Sl39r8H0Du5d/c7gk/vHWj6lTvA6tyrdKv7GL\n+O47eO01o6MQInPKFY/kf/+Dvn0h1rk2Asw0c/m3gU8NpMH3DQg6HmS3aztUEj9+/TgBIQF80foL\no0Oxm5MnYds203i4EM5u2DDw8IDp7jEfIQmjyr85TBLXWjM0aChjG42l5OPus5H2l19C//6yV4pw\nAXv24PHzEuZ2/Z0P34/kyNozRkdkCHP5t8DjgXYp/+YwSXz54eWcv32e4U8PNzoUu7l503RDc+hQ\noyMRwgquXIHDhyl/cxcfNllHjyGF3HaWbPHcxQnuG2yX8m8OkcTvRt1l5NqR+LfzJ6tnVqPDsZtp\n06BjRyhRwuhIhLCugXV2U73CA4YMMToS49ir/JtDJPEPNn9A07JNaVrWfYpKXrkCAQHw/vtGRyKE\n9SkFX489zY4dphv37ixx+bd3Nr5j9fJvhifxsKthzN03l09bfmp0KHb1wQfw4otQtqzRkQhhG7ly\nxrFsGbzzjmlFvzszl3/bcWGH1cu/GZrEtdYMDhzM+03fp2gu16igbYkdO2D5cumFC9f35JOm1ci9\nesH+/UZHYyxz+bcahWtQ79t67L201yrnNTSJ/3jwR25H3uaNum8YGYZdRUTAgAHw+eeQP7/R0Qhh\ne02bgr+/aZ/8AweMjsZYtij/Zlih5PCIcEavH83yHsvx9PA0Kgzr27MHoqISHj/1FHh5/fdwxAjT\nysxeveKfCAlJujIiJsY+cQpx/jxsT1Tk4MIF657/1Kn//mx3KwVx7z9Jq1Z5WbUK6tWz7qWcTQ/v\nHlQrXI3OSzoTciGEz1p9luFJHYZtRTssaBgPoh/w7fPf2u36djFtGty+nfD4zTcht6mw8+zZph74\nrl3w+OPxr0+dCg8e2D9OIeytQwdWXazDq6/CjBnQs6fRARkvPCKcPsv7cCfyDku7LU11WNnhtqLd\ne2kvS8KW8HELFy8JksgPP5huZgYGJkrgQriZDh1g/XrTzc7hw2W3ZWuUf7MoiSul2iiljiiljiml\nxqTSZrpS6rhSap9SqlZq54rTcfgF+vFhsw8pmLPgIwfsbOLiYMoUGDsW1qyBChWMjkgIY9WqZfo2\neu0aVK8O69YZHZGxkpd/+2b3N4/2/vQaKKU8gJlAa6Aa0Esp9WSyNm2B8lrrisDrwNepnW/u3rkA\nvFr71UcK1JqCg4Ptcp1/buSndecc/PabaX+U6tXtctlHEnz6tNEhOAz5LBLY+rPInx9+/BG++goG\nDzYVGNq40VThytHYK1+Yy799teMrBqy0vPybJT1xH+C41vqM1joaWAx0TNamI7AAQGu9A8ijlCqS\n0sne+eMdAtoF4KGMmxhjy/8pWsOO8yV4dUVHnpnTnxa+sWzdCiUddDsYSVwJ5LNIYK/P4rnn4PBh\n6NfPNLxSsaJp2HH3btO3WEdgryQOUKlAJba/tp3wyHCazG3CuVvn0n2PJbNTSgCJz3QeU2JPq82F\n+OeuJD9Z96rdqV2stgWXdXwxMXD5Mhw/DseOmW70b1rxGl4qhldq7uP40OnkGzkYsmQzOlQhHFaW\nLPDKK/Dyy7BzJ/z0E/TpAzduQN26puEXb28oVcp0FCtmKmXoqnJny83SF5by6d+f4jPHh5+6/pRm\ne7tPMTz69STa+yd8bcrIfzPzXjBV5l67NmPvjYyE8HC4dct0U6ZQIVPvoWJF8PGB0VW2UOXxi/Fl\n1vKZ9uVMS5EiGFrTKlcu098KIZ9FYrb4LHLmTPNlpeDpp00HmP6e7tljWiT022+mGZHnz8OlS6a/\nVo8/bpr4lSuXqWZzliymw9Mz4WcPj9RLHqZVCjHxa8eOmWYCP8p7Mk8Bb1Mm91O0uZn2NJ50pxgq\npZ4BJmit28Q/HgtorfXURG2+BjZprZfEPz4CNNVaX0l2Lgcc8RJCCMeX2hRDS3riIUAFpVQZ4BLQ\nE+iVrM1KYDCwJD7phydP4GkFIYQQImPSTeJa61il1BBgHaYbod9prQ8rpV43vay/0VoHKqXaKaX+\nAe4B/WwbthBCCLDzik0hhBDWZfhWtEZSSr2plIpTSrntVlRKqU+UUofjF2ktU0q53XpSSxazuQOl\nVEml1B9KqTCl1EGl1DCjYzKaUspDKbVHKbXS6FhS47ZJXClVEmgJuGchwATrgGpa61rAcWCcwfHY\nlSWL2dxIDDBKa10NqA8MduPPwmw4cMjoINLitkkc+AIYbXQQRtNab9D6v5Lc2wEHXZZkM5YsZnML\nWuvLWut98T/fBQ5jWu/hluI7eu2AOUbHkha3TOJKqeeBc1rrg0bH4mBeBYKMDsLOUlrM5raJy0wp\nVRaoBewwNhJDmTt6Dn3j0LD9xG1NKbUeSLz0X2H6n/E/4B1MQymJX3NZaXwW47XWq+LbjAeitdY/\nGhCicCBKqVzAL8Dw+B6521FKtQeuaK33KaV8ceAc4bJJXGvdMqXnlVLeQFlgv1JKYRo+2K2U8tFa\nX7VjiHaT2mdhppTqi+lrYzO7BORYLgClEz0uGf+cW1JKZcGUwBdqrVcYHY+BGgLPK6XaATmA3Eqp\nBVrrlw2O6yFuP8VQKXUKeEprfdPoWIyglGoDfA400Vpbr3qrk1BKeQJHgeaYFrPtBHpprQ8bGphB\nlFILgH+11qOMjsVRKKWaAm9qrZ83OpaUuOWYeDIaB/6qZAczgFzA+vipVAFGB2RPWutYwLyYLQxY\n7MYJvCHQB2imlNob/+ehjdFxibS5fU9cCCGcmfTEhRDCiUkSF0IIJyZJXAghnJgkcSGEcGKSxIUQ\nwolJEhdCCCcmSVwIIZyYJHEhhHBi/weivJ8/P/EtDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a211d8350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xsample, ysample, accept, l, n = sampling_piecewise(nsamples, I[0], xmin, xmax, xhigh, 1.05*yhigh)\n",
    "xaccept = xsample[accept]\n",
    "print \"Accepting rate: \", float(len(xaccept))/len(xsample)\n",
    "\n",
    "y_propos = propos(x, xmin, xmax, xhigh, 1.05*yhigh)\n",
    "plt.plot(x, y_target)\n",
    "plt.plot(x, y_propos)\n",
    "plt.hist(xaccept, normed=True, bins=50, histtype=\"stepfilled\", color=\"red\", alpha=0.5, linewidth=0)\n",
    "plt.xlim([xmin, xmax])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Make an animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_samples(xsample, ysample, accept, normalize, xhigh, yhigh, xmin=-5, xmax=5, nbins=50, write=False, filename='plot_rejection_sampling.png', trace=False):\n",
    "    nsamples = len(xsample)\n",
    "    ofile = '/home/ridlo/project/stats/rejection_sampling/'+filename\n",
    "    \n",
    "    xaccept = xsample[accept]\n",
    "    n_accept = len(xaccept)\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(xmin, xmax)\n",
    "    \n",
    "    x = np.linspace(xmin, xmax, 1000) # draw target dist line\n",
    "    target_dist = pdftarget(x, normed=normalize)\n",
    "    propos_dist = propos(x, xmin, xmax, xhigh, yhigh)\n",
    "    \n",
    "    ymax = 1.1*np.amax(propos_dist)\n",
    "    ax.set_ylim(0, ymax)\n",
    "    \n",
    "    ax.plot(x, target_dist, 'k')\n",
    "    ax.plot(x, propos_dist, 'g')\n",
    "    \n",
    "    if nsamples > 0:\n",
    "        if n_accept > 0:\n",
    "            ax.hist(xaccept, normed=True, bins=nbins, histtype=\"stepfilled\", color=\"blue\", alpha=0.5, linewidth=0)\n",
    "    \n",
    "            if trace:\n",
    "                sample_x = xsample[-1] # last sample \n",
    "                sample_y = ysample[-1]\n",
    "                sample_accept = accept[-1]\n",
    "                ax.axvline(x=sample_x, ymin=0, ymax=(sample_y)/(ymax), c='k')\n",
    "                ax.axhline(y=sample_y, xmin=xmin, xmax=(sample_x-xmin)/(xmax-xmin), c='k')\n",
    "                if sample_accept:\n",
    "                    ax.plot(sample_x, sample_y, 'k.')\n",
    "                else:\n",
    "                    ax.plot(sample_x, sample_y, 'r.')\n",
    "\n",
    "            rate_accept = float(n_accept)/float(nsamples)\n",
    "            text = r'$n_{sample} = '+'{0:d}$'.format(nsamples)+'\\n'\n",
    "            text += r'$r_{accept} = '+'{0:0.2f}$'.format(rate_accept)\n",
    "            ax.annotate(text, xy=(0.7, 0.97), xycoords='axes fraction', ha='left', va='top') \n",
    "    \n",
    "    if write:\n",
    "        plt.savefig(ofile, bbox_inches='tight', dpi=400); plt.close()\n",
    "    else:\n",
    "        plt.show(); plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD7CAYAAACRxdTpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYVEffxvHvLAgWVMQC2DvGErtRQUXsXWM3do2aJxrf\nJyaPPWLXxFSNJiax995iiwUVjEHFLsUSG2IviIJS5v0DRESUBXb37C7zua69snt25swNwR+H2XPm\nCCkliqIoinXRaR1AURRFMTxV3BVFUayQKu6KoihWSBV3RVEUK6SKu6IoihVSxV1RFMUK2Wod4CUh\nhDonU1EUJY2klCKl7WZ15C6l1PQxYcIEzTOYy+Nt34tZsyQTJrx6BAVpn1X9XKjvRWb9XryLWRV3\nRVEUxTBUcVcURbFCZjPnbg48PT21jmA29P1e+PtDcPCr12XLQrlyxsmkFfVz8Yr6Xrxi7t8Lkdq8\njakIIaS5ZFHe7ttv4cmTt7/foAE0bGi6PIqSmQkhkJbwgaqiKIpiGKq4K4qiWCE1564oitWIi4tj\nxowZlChRgoiICD7++GO9t1kbdeSuKIrVWLlyJUWLFqV79+5cvHiRa9eu6bXt+vXrWkc3OFXcFUWx\nGn5+fhQuXBiAYsWKcejQIb23WRs1LaMoSpr4+fmxdu1aPD09kVJy7tw5xo0bp3UsAHLmzElMTAwQ\nf8V7aGio3tusjSruyls9fw5+fm9uy4gXL8DX9/VtdetC1qwZ269ieoUKFaJmzZps3rzZIPsLDw/n\n888/5/fff0/x/fPnz/PXX38hxJtn/vXp04fcuXPTs2dPDh06ROPGjTl9+jRly5bVe5u1UcVdeasX\nL+DgQcPuMzr6zX1Wq6aKuyVxd3dn+vTp1KxZk/DwcLJnz26Q/WbPnp0SJUq89f3y5ctTvnz5d+6j\nUqVK3L9/nx07dlCoUCEqVqyo9zZro4q7oihpEhkZmVjQt2/fTsuWLfHz86Nu3bqsXLmSrFmzUq1a\nNY4cOUJAQABt27Zlz549dOzYkbVr19KyZUvOnz9P//79OXbsGFu2bKFFixbs37+fmjVrEhQUxJEj\nR3B2dqZFixaJ4748ck9OCEHv3r1xdHRk9+7d3Lhxg/79+7Nz504aNWqk9zZro4q7oihpcu7cOerX\nrw+Ag4MDV69epUKFCgghOH/+PPny5cPX15dvv/2WkJAQSpYsiYuLCzlz5qRMmTLodDqyZcsGQNGi\nRXFyciJPnjzY29tTrVo1vvzyS4YNG4aDg8Nr4+pz5F6mTBkCAwOZN28eXbp0wdbWVu9t1kav5QeE\nEM2BH4g/u+YPKeXMt7SrCRwGukopN6Sxr1p+wMw8eRK/3EBapLb8wNOn8M03r2/7v/8DR8e051PM\ny9q1aylbtiy3bt3i0qVL1KhRg+3bt2NnZ0fnzp05duwYDRo0YPv27bi4uNC6dWtWrVrFs2fPiIuL\nIzIykg8++ICQkBAqVqxIXFwc1apV0/rLMmvvWn4g1eIuhNABIUAj4CZwFOgmpQxKod1fQCSwQEq5\nQd++Cf1VcTczqrgrxvT8+XMGDhzIjz/+iJOTk9ZxLFJG15apBVyQUl6VUkYDq4B2KbQbBqwD7qSj\nr6IomYy9vT1Lly5Vhd1I9CnuhYCkl2/dSNiWSAhREGgvpZwHiLT0VSzPXQI5znwk6i8tRTFXhvoU\n4QdgZEZ34u3tnfjc09PT7NdLzowkkh0M5RanuM1pmvMjOmy0jqUomYKPjw8+Pj56tdWnuIcCRZO8\nLpywLakawCoRf3VBPqCFECJGz76JkhZ3xTxd5i8ec52hBLOOLqyjKx+yDFvUieqKYmzJD3onTpz4\n1rb6TMscBUoLIYoJIeyAbsCWpA2klCUTHiWIn3f/j5Ryiz59FcshiWMPo/BiKtnJSw+2o8OWZTQj\nikdax1MUJYlUi7uUMhYYCuwGzgGrpJSBQojBQohBKXVJra9Bkismd4416LChPJ0AsMWejqzAhWos\npB7h3NA4oaIoL+k15y6l3Am4Jdv261va9k+tr2J5YnnBPsbRml8RST4zF+hoxnf8zbcswJ1KT3cA\n777QRFEMafLkyVSuXJmzZ88yZsyYFNucOXOGSpUqcenSJQoXLoy9vb1e/SyZWvJX0UsAv5OHkpTk\nzcu0BYK6fIEXUxlxuiF+1/xS2IOiGN7evXsBaNu2LdHR0fgmX5UugaenJwULFmTTpk3Y29vr3c+S\nqeKupOoFERxkMo2Z8c5279OTUW5L6bC6A5uCNpkonZKZ+fn5UbVqVQCqVq3Kvn37Umw3e/Zsbt68\nyYgRI9LUz5JZ34IKisEd4QeK0QBXUr8UvKZTU3bU2UGblW24HXGbwTUGmyChYkq+vr5s27aNR48e\n8fjxYz799FM8PDw0yXLnzh1y5MgBxK9zc+vWrRTbHTt2DEdHRwIDAxkxYoTe/SyZKu7KOz3jHkf4\ngYEc0av9nTtQoEB1/vA4xJADzTl/4yY/tPNOcQ1uxTLlz58fBwcHvLy8aNCgAfb29gbZr7+/P7Vq\n1QL0W7sd4u+ZamMTf51FbGxs4vPkvv32W4QQXLlyhV27dundz5Kp4q680yGmUYEuOFEaKSUBAb/h\n7z+bJ0/CKFKkLg0aTKBgweqJ7QMD4x9Qii74sT6oFc90N5nXeh62OvXjZg3c3Nw4duwYI0eOJEuW\nLAbb7/bt2xOLuz4rQAI4Ozvz9OlTIP5mH/nz53+jzaJFi4iNjWXAgAFkzZqV06dP4+Likmo/S6f+\ntSlvdS38KqdYzH84h5RxbNs2hLCw47RqNQ8np9IEBW1i+fIWNGv2He+/3/ON/jkowJiC+1kX3okO\nqzuwutNqwDA3dlC0I6XkxYsXiYX9ypUriWu3T5gwgQcPHnD48GGePHlCs2bNEp97eHgkrtOeP39+\ntmzZQsuWLXn06BFSSoQQhIeHkytXLr3Wbgfw8PDg2LFjtGjRAn9//8R12a9evUqxYsUAyJcvX+Iv\njStXriReBHT06NE3+lkTVdyVt5p2eAI1+AQHXPD1m8nt26fp2/cAdnbx62zXqDGEYsUasGSJF1mz\n5qFs2VZv7COrzoGt3bcycOtAGi1pxOp224C8Jv5KFEO6du0a1au/+mvtxx9/5LvvvuPChQtky5aN\nr7/+mp9++gkhBMOGDUt83q9fv8R12h0dHXFycsLR0ZHg4GBsbGzo27dv4jy4vkfuXl5e7Nixg3Xr\n1iGEoGnTpjx69IgePXrgl3CPyFatWjF79mxy5cpF4cKF8fLyQkrJ9u3bX+tnbfRaz90U1JK/5uXs\nnbM0XOTFgGcXCL9zncWLvRg06Di5cxd5o+3163+zalU7Bg8OIFeuwq+9V7o09OwZf7Q3Zu8Y1gdu\npOX9nThSPLGNWvLXss2dOzdx7fYuXbrwzz//UKVKFeLi4jh9+nTi88DAQCpWrEhsbCwXLlxIXMe9\nT58+LF68mPfee49atWoZdKrH2mVoPXdTUcXdvLRd2ZY6rp5E7f8vy5Y1pWzZtnzwwbC3tj9wYDI3\nbvzNRx9tf237y+L+0qxDs5m0dyY9+BMXKgOquGc2ah13w8noeu5KJuN7zZdTt08xsPJ/uHbtEI8e\nXaFGjSHv7OPhMZKHDy9x4cKOd7b7pNowmvE9S2nCv+w3ZGzFQqh13E1DFXflNVJKRu0ZxUTPiWS1\nzYqf39fUrfslNjbv/lPZxsaOpk2/Zffuz4mLi3ln2wp0pjNrWE83zrHGkPEVRUmgirvymm0h23gU\n9Yhe7/fi8uVLhIb6U7lyb736linTiuzZ83P27OpU2xbHk178xW5G8MuJHzMaW1GUZFRxVxLFxsUy\neu9opjWaho3OhhUrFlOp0kfY2uq3VrsQgvr1x3Ho0FSkjEu1vTPv0w9fFp75hZF/jSROjz6KouhH\nFXcl0bLTy3DM6kibsm2Ii4tj1aolVKnSN037KFmyCXZ2DgTpubaMI8XY0dmXQ9cO0XdTX6Jjo9OR\nXFGU5FRxVwCIioniK5+vmNF4BkIIDhw4QO7cjri4VE7TfoQQ1KnzOf7+s/Xu45QtL3t67+Hx88e0\nXtmaJ8+fpDW+oijJqOKuADDv6Dzed34fj6LxC0CtW7eOTp26p2tf7733IffuBXPnzjm9+2TPkp31\nXdZTLHcxGi5uyO2I2+kaW1GUeHoVdyFEcyFEkBAiRAjxxo2whRBthRCnhBAnhBD+Qgj3JO9dSfqe\nIcMrhvE46jEz/GYwzWsaEH/GzJYtW2jVql269mdjY0f16oM4evTnNPWz1dnya+tfaVO2De4L3Ln4\n4GK6xlcURY/iLoTQAXOAZkAFoLsQolyyZnuklJWllFWBAcDvSd6LAzyllFWllLUMlFsxoFmHZ9G8\ndHMqOVcCICAggGzZslGmTPpvoFW9+iDOnl1FZGTapliEEEzwnMD/3P9H/YX1OXbzWLozKEpmps+R\ney3ggpTyqpQyGlgFvHZIJ6V8luSlA/EF/SWh5ziKBm5F3GLusblM8pyUuG3z5s20a9cuQ8v05sxZ\nkGLF6nH06Lp09R9UfRDzWs2j5fKW7Lq4K905FCWz0qfoFgKuJ3l9I2Hba4QQ7YUQgcBWIOl9VCXw\nlxDiqBDi44yEVQxv8oHJ9H6/N8UciyVu2759O23atMnwvitX7sPBg4u5e5fEx/37+vdvV64dm7pt\nos+mPiw9tTTDeRQlMzHYqpBSyk3AJiGEBzAFaJLwlruUMkwIkZ/4Ih8opbS+GxZaoIsPLrL63GqC\nhgYlbnvw4AEhISHUrl2b588ztv8yZVqxbdtgpkz5lzx5SqRrH3WL1GV/n/20WN6CsIgwvqz7pbrx\nh6LoQZ/iHgoUTfK6cMK2FEkpfYUQJYUQTlLKB1LKsITtd4UQG4mf5kmxuHt7eyc+9/T0TFx3WTGO\n8fvHM/yD4eTLni9x2/79+/Hw8MDOzi7Dxd3W1p4KFbpy+vQyGjQYn+79vJf/Pfz6+9FieQtCw0P5\nvvn36ISa6VMyHx8fH3x8fPRqm+qqkEIIGyAYaASEAf5AdyllYJI2paSUlxKeVwM2SymLCCGyAzop\nZYQQIgewG5gopdydwjhqVUgTCggLoNWKVlwYdgGHhPXZAT755BNKly7NiBEjePIEvv02Y+OEhh5l\n/fruDBt24a1H3PquCvko6hHtV7XH2cGZJe2XYG9rmNu7KYqlytCqkFLKWGAo8YX5HLBKShkohBgs\nhBiU0KyjEOKsECIAmA10SdjuDPgKIU4AR4CtKRV2xfRG7x3NuHrjXivsAHv37jXoXWkKFqyBTmfL\nzZtHM7wvx6yO7Oy5kzgZR/PlzXkc9dgACRXFOqn13DOhff/u4+OtHxP4aSB2NnaJ20NDQ6lcuTJ3\n7txBp9MZ5MgdYP/+r4iOfkbTprNSfD+t67nHxsXyfzv/j4PXDrLjox0UzFkw4yEVxQKp9dyVRC+X\n9J3ScMprhR3g8OHD1K1bF53OsD8WFSp04dy5NXotJqYPG50NP7X4ie4Vu1P3j7oE3QtKvZOiZDKq\nuGcyGwI3EBMXQ9eKXd9472VxN7T8+StgZ+fAjRv/GGyfQghGecSvO++5yJO/r/9tsH0rijVQxT0T\niYmLYey+sUxvND3Fs02MVdyFEFSo0JVz5wx/Y44+VfqwsN1C2q5qy9bgrQbfv6JYKlXcM5GFJxZS\nMGdBmpZ6807vkZGRnD17lho1ahhl7AoVOnP+/FqDTc0k1aJMC7b32M6gbYP4PeD31DsoSiaginsm\n8Sz6GRMPTExc0je5Y8eOUaFCBbJnz26U8fPnL0+2bHm4fv2wUfZfs1BNDvY9yHTf6Uw6MAn14byS\n2aninknM/mc2tQvXplahlNduM9aUTFLvvdeJwMANRtt/mbxlONz/MJuDN/PJn58QGxdrtLEUxdyp\n4p4JPIx8yKy/ZzHVa+pb2xw+fJg6deoYNUe5cu0JCtpk1KNqZwdnfPr4cPnhZTqu6UhkdKTRxlIU\nc6aKeyYww3cGHcp1wC3f25fwPXbsGLVqGXdFZmfn9wHJnTtnjDpOTvucbOuxDQc7BxovbcyDyAdG\nHU9RzJEq7lbuRvgNfj/xOxMaTHhrm1u3bhEZGUnx4sWNmkUIQblyHfS+v2pG2NnYsaTDEuoWrovH\nAg+uPb5m9DEVxZyo4m7lJvpM5ONqH1Mo1xurNCcKCAigWrVqJllt8eXUjCnohI5vmn7Dx9U+xn2B\nO2duG/cvBkUxJwZb8lcxP0H3gtgUvImQoSHvbPeyuJtCkSJ1CQ+/zqNHV3B0LG6SMf9b57+45nSl\n8dLGrOm0hgbFG5hkXEXRkiruVmzsvrF8WfdL8mTLA8CKFXAxyW1JPT2hfv344t6165tXrBqDTmdL\n2bJtCAraTO3awwH46afX23TuDO+9Z9hxu1XsRoEcBei8tjNzW82lU/lOhh1AUcyMmpaxUv/c+Af/\nUH+G1RqWuE1KiIt79Xh50oopj9whfmomOPjV1EzSTElzGZpXCS9299rN8J3DmeM/xziDKIqZUMXd\nCkkpGbV3FBMaTCBblmzvbHv//n0ePHhAqVKlTJQOSpZsQlhYAM+epeGeewZSxaUKvv18me0/m7F7\nx6qLnRSrpYq7Fdp1aRe3Im7Rt0rfVNsGBARQtWpVg68E+S5ZsmSjZMnGhIRsM9mYSZXIUwK//n7s\n/Xcv/bf0Jzo2WpMcimJMqrhbmTgZx6g9o5jqNRVbXeofqQQEBFC9enUTJHudm1t7goI2mnzcl/Jl\nz8fe3nu5+/Qu7Va14+mLp5plURRjUMXdyqw6u4qstlnpUK6DXu2PHz9u0vn2l8qWbcW//+4jOvqZ\nycd+KYddDjZ124SrgysNFzfk7tO7mmVRFEPTq7gLIZoLIYKEECFCiJEpvN9WCHFKCHFCCOEvhHDX\nt69iOC9iXzB+//i3Lg6WklOnTlGlShUjJ3tTtmxOFCxYncuX95h87KRsdbb83vZ3mpVqhvsCdy4/\nvKxpHkUxlFSLuxBCB8wBmgEVgO5CiHLJmu2RUlaWUlYFBgC/p6GvYiDzj8+nbN6yeBb31Kv98+eR\nXLt2jbJlyxo32Fu4ubUjKGizJmMnJYRgstdk/lv7v9RbWI+AsACtIylKhulznnst4IKU8iqAEGIV\n0A5IvLeZlDLp39YOQJy+fRXDiHgRwdRDU9nx0Q69+2zZEoijY2nmzIm/3V7JktC+vbESvsnNrR2+\nvtOJi4tFp7Mx3cBv8UnNT3B2cKb5suYs/3A5TUo10TqSoqSbPtMyhYDrSV7fSNj2GiFEeyFEILAV\n6J+WvkrGfff3d3iV8KKKi/5TLDdunCFfvkqEh0N4OESaeAHFPHlK4ODgwo0bR0w78Dt8+N6HrO+y\nnp4be7L89HKt4yhKuhnsClUp5SZgkxDCA5gCpPmwx9vbO/G5p6cnnp6ehopn1e4+vctP//yE/8f+\naep3585ZChSoZKRU+nFza0dw8GaKFnVPvbGJ1CtWj32999FieQtuRdxiRN0RWkdSFAB8fHzw8fHR\nq60+xT0UKJrkdeGEbSmSUvoKIUoKIZzS2jdpcVf0N/XQVLpX7E7JPCXT1O/OnTPUqjXUSKn04+bW\njg0betCkydea5kiuQoEK+PX3o/ny5oQ+CWVW01kp3ndWUUwp+UHvxIkT39pWn5/Wo0BpIUQxIYQd\n0A3YkrSBEKJUkufVADsp5QN9+ioZc+XRFZaeXsq4+uPS3PfOnTMUKFDRCKn05+pajejoZ9y7Z34f\nwxTJXQTffr4cu3mMnht68jzmudaRFEVvqRZ3KWUsMBTYDZwDVkkpA4UQg4UQgxKadRRCnBVCBACz\ngS7v6muEryPT+mr/VwytORRnB+c09YuMfMDz50/InbuYkZLpRwhB2bJtTbYMcFrlyZaHXT13ERUT\nRcsVLQl/Hq51JEXRi15/Z0opd0op3aSUZaSUMxK2/SqlnJ/w/GspZUUpZTUppbuU8u939VUM4/Tt\n0+y+tDtdc8Lx8+0VTbKGe2rKlYufd9fMkyfw99/x/01BtizZWNt5LeXylqPBogaEPQkzcUBFSTs1\niWjBxuwdw2iP0eSyz5Xmvrdvn9H8w9SXihf35N69ICIibpl+8CdPoF69+LWP69V7a4G30dkwp+Uc\nOr3XCfcF7oTcf/ca+YqiNWEuq+IJIcwjiJKp1AYOAHYAWbLAwYNQu/Y7+yw4sYAxe8ewudtmPij8\ngQlSKkrKhBBIKVP889usbtZhLr9ozJ2UEvcF7nxS4xN6Ve6ld7/ly+HChfjnCxfWw9NzEiVKNEx8\n380Nund/1f7JE/j2W0OlfrezZ1dx+vQyevTYRpcuUL68acZ9eeT+/NQp7MuXhwoVUu3Sv2p/nHM4\n02ZlGxa2W0irsq1MEFRR0kZNy1igLcFbiHgRQY9KPdK9j3v3gsmXz3xWgihdugVXrx7kxYsI0w6c\nMyccOkQDgEOH4l/roVXZVmztvpUBWwaw4MQCo0ZUlPRQxd3CxMbFMmbfGKY3mo5NOi/Zj4x8SExM\nFA4OLgZOl35Zs+amcOHaXLq02/SD58zJPwn/TYsPCn/Agb4HmHJwClMPTlV/eSpmRRV3C7Pk1BLy\nZstLyzIt072P+/eDyZfPzSzOlEnq5dWqlsQtnxt+/f1YF7iOoduHEhsXq3UkRQFUcbcoUTFRTPCZ\nwMzGMzNUmO/dCyZvXjcDJjMMN7e2hIT8SUxMjNZR0sQ1pysH+h4g6H4QXdZ1ISomSutIiqKKuyX5\n2f9nqrlWo06ROhnaz/375lncc+cugqNjMQICfLWOkma57HOxvcd27GzsaLq0KQ8jH2odScnkVHG3\nEI+jHjPTbybTGk3L8L7u3w8hb15t1nBPjZtbe/bts6ypmZfsbe1Z/uFyahSsQb2F9bj++HrqnRTF\nSFRxtxBf+31N67KtKZ8/4+cIvpxzN0flyrVj//7NFvvhpE7o+K7Zd/Sr0g/3Be6cu3NO60hKJqWK\nuwUIexLGL8d/YaLn21eA01dcXCwPHlzCyamMAZIZXoEClZBScvbsWa2jZMiIuiOY3mg6Xku88L1m\nedNMiuVTxd0CTDowiX5V+lEkd5EM7+vx42tkz54PO7scBkhmeEIIGjZsx+bNljk1k9RH73/Esg7L\n+HD1h2wM3Kh1HCWTUcXdzF24f4G159cy2mO0QfYX/2Gqec63v+TlZR3FHaBJqSbs7LmToTuGMu/o\nPK3jKJmIKu5mbtz+cXxe53PyZs9rkP3Ff5hqnvPtL1WvXo/Lly9z48YNraMYRDXXahzqd4jvjnzH\n+H3jLfbzBMWyqOJuxo7fPI7vNV+GfzDcYPuMX3bAvIu7ra0tLVu2ZMsW67mvS8k8JfHr78euS7sY\nuGUgMXGWdS6/YnlUcTdjo/aOYnz98eQw4Py4uZ7jnly7dtYzNfNSgRwF2NdnH2ERYbRf1Z6nL55q\nHUmxYnoVdyFEcyFEkBAiRAgxMoX3ewghTiU8fIUQ7yd570rC9hNCiLTdwTkT23N5D1cfXWVA1QEG\n3a85nwaZVLNmzfj77795/Pix1lEMysHOgc3dNpMvez4aLWnEvWf3tI6kWKlUi7sQQgfMAZoBFYDu\nQojkywleBupLKSsDU4D5Sd6LAzyllFWllLUME9u6xck4Ru0ZxRSvKWSxyWKw/UZFPeXZs3vkypXy\nWTdhYbBp06vHzp0GGzpN/P1h796clC1bj8mTd3LOyk4Vz2KThYXtFuJVwgv3Be5ceXRF60iKFdJn\nPfdawAUp5VUAIcQqoB2QeEdjKeWRJO2PAIWSvBao6Z80WXd+HQCdyncy6H5v376Ak1NpdG9ZTTI8\nHE6eNOiQ6XLlSvx/XV3bsWvXZpo376rPMusWRQjBtEbTKJizIO4L3Pmzx59UcamidSzFiuhTdAsB\nSa+jvsHrxTu5gcCOJK8l8JcQ4qgQ4uO0R8xcomOjGbtvLDMaz0AnDPs7MSzMMubbX3Jza8fFizt4\n8cJ6F+IaWmsoPzb/kaZLm7Lv331ax1GsiEHvxCSEaAj0AzySbHaXUoYJIfITX+QDpZQpXrLn7e2d\n+NzT0xNPT09DxrMIf5z4g+KOxWlcsrHB921pxd3BwRkXl6ocPbqTli3bax3HaDqV70T+7Pnpsq4L\nPzb/kW4Vu2kdSTFTPj4++Pj46NU21XuoCiFqA95SyuYJr0cBUko5M1m794H1QHMp5aW37GsC8ERK\n+V0K78nMfv7v0xdPKTO7DFu7b6V6weoG37+7+0fkz9+UKlX6GHzfxnLs2K88ferD/v0rjTpOwr0o\njTpGas7cPkPLFS0ZUWcE/1f7/zTNoliGd91DVZ+/+48CpYUQxYQQdkA34LUTkIUQRYkv7L2SFnYh\nRHYhhEPC8xxAU8CyFw0xoh//+ZF6xeoZpbADhIWFWMSZMkm9996HHD26g2fPnmkdxegqOVfCr78f\n84/P58vdXxIn47SOpFiwVIu7lDIWGArsBs4Bq6SUgUKIwUKIQQnNxgNOwNxkpzw6A75CiBPEf9C6\nVUqpwX3UzN/9Z/f57u/vmNJwilH2L6W0uGkZgBw58uPmVovt27drHcUkiuYuim9/Xw7fOEzvjb15\nEftC60iKhUp1WsZUMvu0zBe7v+Dpi6fMa22c9UfCwsIoV64y//3vHaPs35jCw3/n+vVdrF271mhj\nmMO0TFKR0ZF0X9+dZ9HPWN9lPTnt03Z/VyVzyOi0jGJk1x5fY+HJhXzV4CujjREcHIyLi2Udtb/k\n4fEhu3fvJiIiQusoJpMtSzbWdVlHCccSeC725HbEba0jKRZGFXcz4O3jzZDqQ3DN6Wq0MYKDg3F1\nNe/VIN8mVy4n3N3d2bZtm9ZRTMpWZ8svrX+hnVs76i6oy4X7F7SOpFgQVdw1dv7uebaFbON/7v8z\n6jghISG4ulrmkTtAly5dWL16tdYxTE4IwVcNvmKU+ygaLGrA0dCjWkdSLIQq7hobs3cMI91Hkjtr\nbqOOE3/kbrnFvX379uzbt8/q1prR18fVP+aX1r/QckVLdl7UaF0IxaKo4q6hw9cPExAWwKe1PjX6\nWJZe3B0dHfHy8mL9+vVaR9FMW7e2bO62mb6b+rL45GKt4yhmThV3jUgpGbVnFBM9J5LVNqtRx3rx\n4gXXr19UNZ+DAAAgAElEQVSnQIGSRh3H2Hr37s2SJUu0jqGpukXq4tPXhwk+E5jhO8OszvBRzIsq\n7hrZfmE79yPv07tyb6OPdenSJYoUKYKtrZ3RxzKmli1bcvbsWa68XFkskyqXrxyHBxxmxZkVDN85\nnNi4WK0jKWZIFXcNxMbFMnrvaKZ5TcPmLSs0GlJISAhubpY7JfOSvb09Xbt2ZdmyZVpH0VzBnAU5\n2O8gZ+6codv6bkTFWO/iakr6qOKugRVnVpDTPidt3dqaZLzg4GCrKO7wampGTUeAY1ZHdn60E4Gg\n+bLmPIp6pHUkxYyo4m5iz2Oe85XPV8xoNAMhUrywzOCsqbjXqlULIQT//POP1lHMgr2tPas6raKy\nc2XqL6xPaHio1pEUM6GKu4n9cuwXKuSvQL1i9Uw2ZnBwMGXLWuYFTMkJIejTpw+LF6uzRV7SCR0/\nNP+Bnu/3xH2BO4F3A7WOpJgBVdxNKPx5ONN8pzG90XSTjmstc+4v9ezZkzVr1vD8+XOto5gNIQT/\nc/8fkxpOouHihhy+fljrSIrGVHE3oW8Pf0uzUs2o5FzJZGM+fPiQqKgoXFxcTDamsRUtWpSqVauy\nceNGraOYnd6Ve7Oo/SLar2rPluAtqXdQrJYq7iZyO+I2c47OYVLDSSYd9+V8u6nm901l0KBB/Prr\nr1rHMEvNSzdn+0fbGbJtCPOPz0+9g2KVVHE3kSkHp9Dr/V4Udyxu0nGtab49qfbt2xMYGEhwcLDW\nUcxSjYI1ONjvIF/7fY23j7c6uygTUsXdBC4/vMzKsysZW2+syce2pjNlkrKzs6Nv377Mn6+OTN+m\ntFNp/Pr7sS1kG4O3DSYmLkbrSIoJ6VXchRDNhRBBQogQIcTIFN7vIYQ4lfDwTbifql59M4Px+8fz\n2QefkT9HfpOPbW0fpgI8eQKnT0O9eh+zcOESAgLUBTxv4+zgzP4++7n6+Cod13TkWbT1365QiZdq\ncRdC6IA5QDOgAtBdCFEuWbPLQH0pZWVgCjA/DX2t2slbJ9n37z4+r/O5JuNb45H7nTuwYQMcPVqK\nPHmqMnNm5l1MTB857XOytftWctnnovGSxtx/dl/rSIoJ6HPkXgu4IKW8KqWMBlYB7ZI2kFIekVK+\nXIv1CFBI377WbvTe0YytNxYHOweTjx0bG8vFixcpU6aMycc2lerVB+Pv/4vWMcyenY0di9svpl7R\nengs9ODqo6taR1KMTJ/iXgi4nuT1DV4V75QMBHaks69V8bniQ8j9EAZVH5R6YyO4du0a+fLlI0eO\nHJqMbwpubm15+PAKAQEBWkcxezqhY2aTmQypPgSPhR6cvn1a60iKERn0A1UhREOgH5Ap59aTklIy\ncs9IJjecjJ2NNqsxWuOUTHI2NlmoU2cY33//vdZRLMbw2sOZ1WQWjZc0xueKj9ZxFCOx1aNNKFA0\nyevCCdtek/Ah6nyguZTyYVr6vuTt7Z343NPTE09PTz3imaeNQRt5EfuCbhW7aZbBGj9MTUmtWh/z\n008luXnzJgULFtQ6jkXoWrEr+XPkp8vaLvzc8mc6V+isdSRFDz4+Pvj4+OjVVqR2/qsQwgYIBhoB\nYYA/0F1KGZikTVFgL9BLSnkkLX2TtJXWci5uTFwMFedW5IfmP9C8dHPNcnz66ae4ubnx2WefAbB8\nOVywwHsse3hA48avXl+6BEuXvnqdNSvcuDGU3LlzM3Xq1HSNIYTIlOeCn7p1ilYrWjHSfSTDPhim\ndRwljRJ+blO8QjHVaRkpZSwwFNgNnANWSSkDhRCDhRAvJ5PHA07AXCHECSGE/7v6ZvgrMnOLTi7C\nNacrzUo10zSHtV7AlJLhw4czf/58nj1Tp/qlRWWXyvj29+Xnoz8zes/oTPkLzlrpMy2DlHIn4JZs\n269Jnn8MfKxvX2sWGR3JxAMTWdd5nUkv+Y+OhmvXXt927lwwWbO6celS/GtrrntlypShTp06LF68\nmE8++UTrOEYTFQWhySY2S5XK2D6LOxbHt78vbVa2oe/mvvze5ney2GTJ2E4VzelV3BX9zfafTa1C\ntfig8AcmHffZs9enKl68eMq9e/fYt68oukxyHfLIkSPp2bMnAwcOJEsW6yxO9+69/v9ZCJgwIeP7\nzZc9H3t776Xruq60WdmGdV3WaXL6rmI4meSfvWk8jHzIN4e/YapX+uZ9DenBgws4OZVGZ4Lb+JkL\nd3d3ihcvzooVK7SOYpGyZ8nOxq4bKZyrMA0XN+TO0ztaR1IyQBV3A5rpN5P2bu0pl0/7i3Dv3Qsm\nb97MMd+e1Pjx45k2bRqxseqm0elhq7Pltza/0aJ0C9wXuHPpwSWtIynppIq7gYSGh/JbwG94e3pr\nHQWA+/eDyZs303zUkahhw4bky5ePNWvWaB3FYgkhmNRwEiPqjKDewnocv3lc60hKOqjibiATD0xk\nYNWBFMplHhfgZtbiLoRg/PjxTJkyRR29Z9CQGkOY22ouLZa3YPel3VrHUdJIFXcDCL4XzMagjYzy\nGKV1lET374eQL1/mK+4AzZo1w8nJiaVJP3lU0qV9ufZs6LqBXht7sez0Mq3jKGmgirsBjN03li/q\nfEGebHm0jgLEL32QWefcIf7ofebMmXz11VdERkZqHcfieRT1YF/vfYzdN5ZZh2epc+EthCruGeQf\n6s+RG0fM6uq+iIhb2Nraky2bk9ZRNFO3bl2qV6/Ozz//bLQxHj+G27dfPZ48MdpQmqtQoAJ+/f1Y\nfGoxn+/6nDgZp3UkJRWquGeAlJJRe0YxocEEsmfJrnWcRJl1vj25adOmMXPmTB4+fJh643TYuxfm\nzXv18PMzyjBmo3Cuwhzse5DjYcfpsb4Hz2Oeax1JeQdV3DNg96Xd3Hxyk35V+2kd5TXxUzKquL/3\n3nt07NiRCYa4ykcBIE+2POzutZuYuBhaLG/B46jHqXdSNKGKezrFyThG7R3FVK+p2OrM60Lf+/eD\nM+2HqclNnTqV1atXc/LkSa2jWI2stllZ3Wk15fOXp8GiBoQ9CdM6kpICVdzTafXZ1djZ2PHhex9q\nHeUN8cVd+wupzEHevHmZMmUKQ4cOJS5OzRMbio3OhtktZtOlQhfqLqhL8L1grSMpyajing4vYl8w\nbv84ZjSaYdLFwfR1716QKu5JDBgwgBcvXrBkyRKto1gVIQRj6o3hq/pf0WBRA47cOJJ6J8VkVHFP\nh9+O/0YZpzI0LNFQ6yhviImJIjw8FEfHElpHMRs6nY5ffvmFkSNHcvPmTa3jWJ1+VfuxoN0C2qxs\nw7aQbVrHURKo4p5GES8imHJoCtMbTdc6SooePLhInjwlsFFLtr6mWrVqDBkyhEGDBqnztI2gZZmW\n/NnjTz7e+jF/BPyhdRwFVdzT7Pu/v6dh8YZUda2qdZQUWeOUTHQ0RES8ekRFvf6+lK+/HxEByafX\no6Jg+PCxXL8eyq+/LuZ5CmfxJe3/9Knxvh5Te/789a8t+dceG/vm9y89v/9qFarFgb4HmHpoKlMO\nTlG/RDVmXqd5mLm7T+/y4z8/8s/Af7SO8lb37gWRN691Ffd//ol/vM3z5zBr1uvbPvsMnJJcw7V1\nK5w7Z0edOosYMaIx9vYe9OtX+rU+SfdhZwdjxhggvBnYtQsCAl69rl4d2rR59fr2bZg///U+6T17\ntGzeshwecJgWy1sQGh7KnJZzsMlEy06bE72O3IUQzYUQQUKIECHEyBTedxNCHBZCRAkhPk/23hUh\nxKmkt9+zVNMOTaNbxW6UcsrgrW+MyBqP3A3JxaUyDRp4M3VqZ7U0gZG4OLhwoO8BLjy4QKe1nYiM\nVt9nLaRa3IUQOmAO0AyoAHQXQiSvHveBYcA3KewiDvCUUlaVUtbKYF7NXH10lSWnlzC+/nito7yT\nOsc9dTVr/odChcoyfPhwraNYrVz2udj+0Xay2WajydImPIh8oHWkTEefI/dawAUp5VUpZTSwCmiX\ntIGU8p6U8jgQk0J/oec4Zu0rn6/4tOanODs4ax3lreIXDAtSV6emQgjBf//7OwcOHOC3337TOo7V\nsrOxY9mHy6hduDb1Ftbj+uPrWkfKVPSZcy8EJP2/coP4gq8vCfwlhIgF5kspLe5f05nbZ9h5cScX\nhl3QOso7PXlykyxZcpDNTFanNGfZs+dk69at1K9fn6JFi2odx2rphI5ZTWfh6uCK+wJ3tn+0nYoF\nKmodK1MwxQeq7lLKMCFEfuKLfKCU0jelht7e3onPPT098fT0NEG81I3ZN4bRHqPJZZ9L6yjvpK5M\nTZuyZcuyfv16OnTooHUUqzei7ghcc7rSaEkj1nZeS/1i9bWOZJF8fHzw8fHRq60+xT0USHpoUzhh\nm16klGEJ/70rhNhI/FF/qsXdXPhe8+XM7TOs67xO6yipUlMyaefu7s7PP/9Mly5duHPnHAUKVNA6\nktXqUakHBXIUoNOaTsxrNY+O5TtqHcniJD/onThx4lvb6jMXfhQoLYQoJoSwA7oBW97RPvF6fCFE\ndiGEQ8LzHEBT4KweY5oFKSUj94xkUsNJ2Nvaax0nVepMmfTp3LkzAEuXNuHOHYv58bRIjUs2ZlfP\nXXy28zPmHp2rdRyrluqRu5QyVggxFNhN/C+DP6SUgUKIwfFvy/lCCGfgGJATiBNCDAfKA/mBjUII\nmTDWcimlxdyMcWvIVsKfh/NRpY+0jqKXe/eCKF26hdYxLFbTpt+yZEljunbdQKlSdU0y5ty5cOfO\nq9ft20OVKiYZWjNVXatyqN8hmi9rTmh4KFO8ppjlGk2WTq85dynlTsAt2bZfkzy/DRRJoWsEYJE/\nqrFxsYzZO4YZjWdYzEUYas49YypV6k7WrLlZtaodrVv/BHTXOlKqIiMf4ud3nqCgIEJDQ3n8+DFx\ncXHkzJkTV1dXKlSowPPnNYFsWkd9Tck8JfHr70erFa0Iiwjj19a/kkUtmWFQ6grVt1h6eil5suWh\nVZlWWkfRy9OnT3n69C65c6szPzKiTJmW9O69j1Wr2vDFF8eZOnUq9vbmMyUXEfGE8+d38e+/+7l2\n7SCPHl3Fx6c85cqVo0iRIri6uqLT6Xjy5AnHjx9n0aJFnDkTSJEinlSvPpjSpZsD5nGwkj9Hfvb3\n2U/ntZ1pv7o9azqtIYddDq1jWQ1V3FMQFRPFBJ8JrOy40mL+XLxwIRgnp9LoLOSvDHPm7FyJTz89\nzrlzA6lduzbLly+nfPnymuW5f/8+W7duZcOGDezf74Ozcx1KlmxClSp9cXWtyqRJ7/5nvGLFQ9as\n2cyBA97s2fM/dLpvadOmuYnSv1sOuxxs7raZQdsG4bXEi23dt5E/R36tY1kFi7+4yBjmHp1LFZcq\n1C1imnlXQwgKUmd6GFL27HnZsGEDn3zyCQ0aNOCLL74gPDzcZOOHhYUxd+5cGjduTMmSJdm6dStd\nu3blyJFr9Oy5i7p1v6BQoZrY2KR+fObgkIcqVfoycKA/jRpNZ/78YfTo0YPHj83jFnlZbLKwoO0C\nmpRsgvsCd/59+K/WkayCKu7JPI56zEy/mUzzmqZ1lDQJCjpH/vyquBuSEIJBgwZx9uxZHj58iJub\nG9OnTzfaDbcfPbrKsmXf4+HhQfny5Tl8+DD/+c9/CAsLY/369Xz00Ufkzu2Y7v0LIXBza8tPP53C\n0dGRatWqERgYaMCvIP2EEEzxmsLwD4bjsdCDE2EntI5k8dS0TDLfHP6GlmVaUsHCjoIDA89RoMAA\nrWNYJWdnZ/744w/Onj3LN998Q6lSpejYsSOFC/cgLq5+uqfCoqKiCAry5fjx3Vy+/Bfh4Tdo0qQt\nY8aMoVGjRkab67e3z87cuXNZunQpnp6ezJu3FjCPi4o+rfUpLg4uNFvWjBUdV9C4ZGOtI1ksVdyT\nCHsSxrxj8zgx2PKOGgIDz9K2rWX9QrI0FStWZPHixYSGhrJixQrmzRvBrVvXKFasHkWKeKDTvUep\nUmVwcXEhR44c6HQ6YmJiiIyM5NatW1y9epXLly9z4sQJjh07xvnz53FxqUKRIk1o0WIOhQt/wIcf\n2prsVMhevXrh6upK164dadduE0WLuptm4FR0LN+R/Dny03ltZ35o9gPdK5n/WUvmSBX3JCYfnEzf\nyn0pagFnnISEwObN8c+fP48gLOw2efKU1DaUGfntN9AlmXRMfoOPo0fh1Km07fOPP+DBA4hfbulL\n+vX7kocPb3D16iGuXz/Mpk27WLjwAnfv3iUyMhIbGxtiYmLIli0bzs7OFC9enOLFi1OlShV69epF\n5cqVWbw4x2vnuScXEAB79756HRubtsypady4MXPmLGfgwA707LkLVzO5CU39YvXZ02sPLVe0JCwi\njM/rfJ56J+U1qrgnuPjgImvOrSF4qGXcxT029tXdgkJDA8mb102dKZNEaku1R0fHP9K6z+R3aMqV\nqzCVKnWnUqXu1K4NzRNOQomLiyM6Oho7O7sMnXEVE2P8u0I1aNCUli1/ZvXq9gwc6I+Dmax8Wsm5\nEof7H6bZsmbcfHKTr5t8jU6ojwn1pb5TCcbtG8d/a/+XvNnzah0lze7eVR+mmhudToe9vb3FnEpb\noUJnKlfuy5o1HxITk8I9CDVSJHcRfPv7cuTGEXpt7MWL2BdaR7IYqrgDx28e5+DVg/xf7f/TOkq6\n3LlzVhV3JcM8PSeQPXt+9uwZpXWU1zhlc+KvXn/xLPoZrVa0Ivy56U5JtWSquAOj945mfP3xFnt1\n3N275yig1shWMkgIHW3b/kFg4Dp27zavJaCyZcnGus7rKJWnFJ6LPLkVcUvrSGYv0xf3vZf38u+j\nfxlYbaDWUdJNLVWrGEr27Hlp124RAwb05969e1rHeY2NzoZ5rebRoVwH6v5Rl5D7IVpHMmuZurhL\nKRm1dxRTGk6x2EWLnj8PJzLyAY6OxbWOoliJkiUb0alTZ0aMGKF1lDcIIRjfYDxj6o2hwaIG+If6\nax3JbGXq4r7u/DriZBydK3TWOkq63b17nnz5yiHUWQSKAU2ePBkfHx/27dundZQUDaw2kPmt59N6\nRWt2XNihdRyzlGlPhYyOjWbsvrH83PJniz696s6ds2pKxkyFhcHGja9e29jA4MHv7rNvHxw+/Op1\naqd0GouDgwNz5sxhyJAhnD59mqxZs2oT5B3auLVhS/cttF/VnhmNZ9C3Sl+tI5kVy61qGbTgxAKK\n5i5Kk1JNtI6SIfFnyqgPU81RdHT8jThePu7eTb1PePjrfZ48MX7Ot2nTpg2VKlVi1qxZ2oVIRe3C\ntfHp68PEAxOZdmgaUkqtI5kNvYq7EKK5ECJICBEihBiZwvtuQojDQogoIcTnaemrhWfRz5h0cBIz\nGs/QOkqG3bp1wmyuKlSszzfffMP333/PzZs3tY7yVuXylcOvvx+rz61m2I5hxMYZ+DJeC5VqcRfx\nk7lzgGZABaC7ECL57X7uA8OAb9LR1+R+PPIj7kXcqVGwhtZRMkRKya1bp3B2rqx1FMVKlSxZkoED\nBzJu3Dito7xTwZwFOdj3IOfvnqfruq5ExUSl3snK6XPkXgu4IKW8KqWMBlYB7ZI2kFLek1IeB2LS\n2tfUHkQ+4Lsj3zHFa4qWMQzi0aMr2Nk5kEPd3EAxojFjxrB9+3ZOnDDvBfVyZ83Njo92YKOzodmy\nZjyKeqR1JE3pU9wLAdeTvL6RsE0fGelrFNMPTafjex0pm7esljEM4tatk7i4WOQtahULkjt3bry9\nvfnyyy+1jpIqe1t7VnZcSVWXqtRbWI8b4Te0jqQZszpbxtvbO/G5p6cnnp6eBt3/9cfXWXByAWc+\nOWPQ/WpFFXfFVAYMGMA333zD/v37adiwodZx3kkndHzf7HtmHZ6F+wJ3dny0g/L5tbtNoiH5+Pjg\n4+OjV1t9insokHQN3MIJ2/SRpr5Ji7sxePt4M7j6YArmLGjUcUzl9u2TvP9+b61jKJlAlixZ8Pb2\nZvz48Rw6dMjsF0QTQvCl+5e4OLjQcHFDNnTZgLuZrFefEckPeidOnPjWtvoU96NAaSFEMSAM6Aa8\na/X8pP/X09rXaM7fPc/WkK2EDLOeS5Zv3TpJ06bfaR1DSXD+PK+tzZ58DXljkBKWLHl9W+vW4ORk\n+LF69OjB9OnT2bVrF46OzQlJ8k+pdGmoa4a3HO5VuRfODs50WN2B+W3m075ce60jmUyqxV1KGSuE\nGArsJn6O/g8pZaAQYnD823K+EMIZOAbkBOKEEMOB8lLKiJT6Gu2reYex+8byP/f/4Zg1/fegNCeP\nHj0gMvIhefKU0DqKkiA8PP5hapcvv/76uZFW7LWxsWHSpEmMGzeOCROacfnyq+O43LmNM6YhNC3V\nlO0fbaftyrbcirjFkBpDtI5kEnrNuUspdwJuybb9muT5baCIvn1N7e/rf3P85nFWdlypZQyDCgo6\niYtLZbXsgGJSH374IVOnTuWff7Zga6vpiW9pUqNgDQ71O0SzZc0IexKGt6e32U8tZZTVV4aXi4N5\ne3qT1db8LqFOr6Cgkzg7qw9TFdPS6XSMHz+e1aunWNzVoKWcSnF4wGH+vPAng7YOIiYu+Znb1sXq\ni/uOizu4+/QuvStb1weP584dw9W1mtYxlEyoffv2PH/+jMuX/9I6SpoVyFEAn74+XA+/TofVHXgW\n/UzrSEZj1cU9TsYxeu9opjWahq3OrM76zLDTp/0pXPgDrWMomZBOp6NLlzEcOjRV6yjp4mDnwNbu\nW8mTNQ+NljTi3jPzWrfeUKy6uK84s4IcWXLQzs1y5gb1cf/+fR48uEPevJp+lKFkYvXqdSU8PJRr\n13y1jpIuWWyysLj9YhoUa4DHAg+uPLqidSSDs9ri/jzmOeP3j2dG4xlW98HJsWPHqFixBjqdjdZR\nlEzKxsYWd/eRFnv0DvHnws9oPIP/1PwPHgs8OHXrlNaRDMq65iqS+PX4r5TPX576xeprHcXg/P39\nqVixptYxlDSKjYX161/f9vix4cfZuxeyZXv1+kayK/D//ff1HCmtGb9hw7vHCA2FypV7c/DgJG7e\nPE7VqtXTH1hjn33wGS4OLjRZ2oTVnVbTsIR5X4GrL6ss7k+eP2HaoWns7mVeN/k1FH9/fxo37svD\nh1onUdJCSjhjgpUvLl589/sPHsQ/3kWfnLa29tSp8wW+vtNo1Wp96h3MWJcKXcifPT9d13VlTss5\ndKnQRetIGWaV0zLf/v0tTUo14X3n97WOYnBSSvz9/alUqZbWURSF6tU/5to1P65dO6d1lAxrWKIh\ne3rv4fNdn/PjkR+1jpNhVlfc7zy9w2z/2UzynKR1FKO4du0aQghcXAprHUVRyJIlOx98MJz166dr\nHcUg3nd+H7/+fsw7No+Rf40kTsZpHSndrK64Tzk4hZ6VelLCSi/L9/X1pW7dulb3IbFiuWrW/A8B\nATu5dOmS1lEMophjMfz6+3Ho2iH6bupLdGy01pHSxaqK++WHl1l+Zjlj64/VOorRHDp0iHr16mkd\nQ1ESZc2am+bNP2HmzJlaRzGYvNnzsqf3Hh5FPaL1ytY8ea7hzWzTyaqK+1f7v+KzWp9RIEcBraMY\njSruijlq02Y469ev50byU3MsWPYs2dnQdQNFcxWl4eKG3I64rXWkNLGa4n7q1in2XN7D53U+T72x\nhbp//z7Xr1+nShW1poxiXnLlyke/fv2YNWuW1lEMylZny/w282ldtjXuC9y5+CCVU5HMiNWcCjl6\n72jG1htLTvucWkcxipAQWLHClzJlarN3ry3372udyLrFxMCuXa9vi4jQJoslCA2FatVGMHhwBerU\nGUPhwgVwt/x7YwDxFzt5e3rj6uBK/YX12dJ9CzUK1tA6VqqsorgfuHKAoHtBbOq2SesoRnPtGuze\nfYhcuerx999ap7F+cXGo73Ma3LkDd+64Uq5cN2bP/p6PPppuNcX9pcE1BuPi4ELL5S1Z2mEpzUo3\n0zrSO1n8tIyUkpF7RjK54WTsbOy0jmNU164dolgxNd+umC939/8REDCfZ8+s8wq7duXasbHrRnpv\n6s3SU0u1jvNOehV3IURzIUSQECJECDHyLW1+EkJcEEKcFEJUTbL9ihDilBDihBDC31DBX9oUtImo\nmCi6V9Lk7n0mExHxiLt3z1O4cG2toyjKWzk6Fqds2Tb4+MzROorRuBd1Z3+f/YzbP46v/b4223Xt\nUy3uIv5WP3OAZkAFoLsQolyyNi2AUlLKMsBgYF6St+MATyllVSmlQS+rjImLYcy+MUxvNB2dld+R\n6MSJvRQt6oGtFd1wRLFOHh6j8fGZTYQVf0hRPn95Dvc/zNLTS/m/nf9nlhc76VMRawEXpJRXpZTR\nwCog+Rq67YAlAFLKf4DcCfdVhfgbZhul8i4+uRjnHM40L93cGLs3K8eP76ZkyaZax1CUVOXL50bZ\nsp78+uuvqTe2YIVyFeJQv0OcvH2S7uu78zzGSDevTSd9im4h4HqS1zcStr2rTWiSNhL4SwhxVAjx\ncXqDJhcZHYn3AW+rXNI3OSklx47tolQpVdwVy9Cs2Ri+/fZboqKitI5iVI5ZHdnVcxdxMo7my5vz\nOMoIy3ymkynmMtyllNWAlsCnQggPQ+x0jv8cahasSe1MMAd94cIFYmNjyJ+/vNZRFEUvRYpUoVq1\navz+++9aRzG6rLZZWdVxFRXyV6D+ovrcfHJT60iAfqdChgJFk7wunLAteZsiKbWRUoYl/PeuEGIj\n8dM8Kd6+xdvbO/G5p6cnnp6eKQZ6FPWIbw5/w4G+B/SIb/l27dpF9epNrf4vFMW6TJo0idatW9Ov\nXz9y5MjxxvunTsGTJFf1lywJBQuaMKAB2ehsmN1iNjN8Z1D3j7rs7LmTcvnKpd4xjXx8fPDx8dGr\nrUjtk14hhA0QDDQCwgB/oLuUMjBJm5bAp1LKVkKI2sAPUsraQojsgE5KGSGEyAHsBiZKKd9YaF0I\nIfX91Hn0ntHcfXaX39ta/1EBQNOmTaldexA6XSeto1itiRMFEyaY51kPlsjZGT75BLp06ULVqlUZ\nPTfekX8AAAuVSURBVHr0G23++AOuJ5nMbdECPrCC2wIvOrmIUXtGsbHrRuoUqWPUsYQQSClTPOpL\ndVpGShkLDCW+MJ8DVkkpA4UQg4UQgxLabAf+FUJcBH4F/pPQ3RnwFUKcAI4AW1Mq7Glx88lN5gfM\nx9vTOyO7sRgPHz7kyJEj1Kxp/R8aK9Zn8uTJfPfddzzMRHeW6VulLwvbLaTtqrZsDd6qWQ69rlCV\nUu4E3JJt+zXZ66Ep9PsXMOhCKBN9JjKg6gAK58oc65lv27aNhg0bki2bg9ZRFCXN3NzcaN++PV9/\n/TXTp1vHmu/6aFGmBX/2+JN2q9ox+elkBlYbaPIMFnVyePC9YDYEbWCUxyito5jMxo0b+fDDD7WO\noSjpNmHCBObPn29VK0bqo1ahWhzse5DpvtOZdGCSyS92sqjiPm7/OEbUGYFTNieto5jE48eP2bdv\nH23atNE6iqKkW+HChRkyZAgjR6Z4cbtVK5O3DH79/dgUtIlP/vyE2LhYk41tMcX9aOhR/r7+N599\n8JnWUUxm3bp1eHl54eSUOX6ZKdZr9OjRHDhwAD8/P62jmJyLgwsH+h7g0sNLdFzTkcjoSJOMaxHF\nXUrJqL2j+KrBV2TPkl3rOCazePFi+vTpo3UMRckwBwcHvv76az777DNiY0139Gouctrn5M8ef5LD\nLgeNlzbmQeQDo49pEUv+/nX5L26E36B/1f5aR0mXhw/hypVXr+3toXyy65FOn4akP/M63WUCAwNp\n0aKFSTIqiiFFRsKJE69vK1euO3Fxc/nqq9/o1GkIT59qk00rdjZ2LO2wlJF/jcRjgQc7e+6kaO6i\nqXdMJ7Mv7nEyjlF7RjHVayq2OrOPm6LQUNi8+dXrPHneLO7bt0PSK7VDQ+fRp08f7OysexljxTqF\nh7/+Mx9PUKfOL/zwQ0OePWtF7txFUupq1XRCxzdNv8E1pyvuC9zZ3mM7lZwrGWcso+zVgNacW4Ot\nzpaO73XUOorJvHjxlPXrF/Lpp59qHUVRDKpAgYrUqjWMP/8cYrZL5ZrC53U+5+vGX9N4aWMOXDHO\nlfZmXdxfxL5g3L5xmWJxsKROn15G9eoelChRQusoimJwHh6jCA+/wcmTi7SOoqnulbqz4sMVdF7b\nmXXn1xl8/2Zd3H8P+J1STqXwKuGldRSTiY2Nxs9vJgMHfqF1FEUxChsbOzp0WMaePf/j7t3A1DtY\nsUYlG7G7126G7xzOHH/D3uDEbIt7xIsIphycwoxGM7SOYlInTy7Eyak0NWoYZPFMRTFLzs6V8PKa\nxrp1XYiOfqZ1HE1VcamCbz9fZvvPZszeMQabrjLb4v7DkR9oULwBVV2rpt7YSrx48ZSDB6fg6TlR\n6yiKYnTVqg3E2bkyW7YMRJrhnYxMqUSeEvj282Xvv3vpt7kf0bH/3979x0Zd33Ecf74O21laIJiU\nllmG2hIZPxYYSpgFFBxSMRRiDNDQkMk0RhzToBJAwCmNy5Y1M2NhZPxYYGJRN6osDihrMZRuuLEe\nDAsMmBQ2sLXMFLPKgHLv/XEHraxc74De97j7PJKmvev3c/e6T9v3ffr9fj7f74Xrfsy4LO6nvzjN\n67tfZ9m4ZV5HiamdO0vo338M/br4THKOEw8kMXnyKpqb66mqWuJ1HM9lpmdSNauKpi+amLJxCi3n\nr2+uaFwW99eqX2P64Onk3ZbndZSYaWjYj9+/moceKgWC0yePHGn7SKKT6jlJJCUljRkz3uPAgbd5\n550v73P+/PMv/w0cP+5RyBhKT03n3envkp2Rzbh142hqabrmx4q7iePHm4+zbt866ubUeR0lZlpa\nWigrm86ECT8hIyMbgKoqj0M5Toykp2dSXFxBWdmDZGf/lxdeCE4mOHYMysvbtsvMhGSYHZzSLYU1\nhWtYsmMJ+Wvz2Vq8lbt63xX148TdyP3lD15mzj1zyA4VuUQXCASYPXs2OTkjGTbMnWrASU69e9/J\nypU7Wb16NfPnz6e1tdXrSJ6SRMn4Ep4b9RxjfjWG2k9qo36MuCruH336EVuObuHF/Be9jhITgUCA\nuXPn0tDQwNSpv/A6juN4qk+fHKqrq6mtraWgoIDPPvvU60iem3PvHJY/vJyJb0xk+z+2R9U2ouIu\nqUDSIUmHJXV43k5JP5N0RNJeScOiaXvJospFLMhfQM+v9IzqRdyMzp49y8yZM/H7/WzevJmUlDSv\nIzmO5zIzM9m2bRujRo3isceG4vevTeqVrACPfv1RNk3bRHF5MRv+tiHidp0Wd0k+4OfARGAwUCRp\n4BXbPAzkmtkA4ClgZaRt29vXuI+n73064vA3WqQXnr1ex479iREjRuDz+aisrKRXr14xed5o1Nd/\n4HWEuOH6ok0s+qJbt26UlJSwYsVW9uxZyapV93Dw4CYCgfiaLhmregEwpv8YKmdVsrByIaV/LI2o\nTSQj95HAETM7bmYXgI3AlCu2mQKsBzCzD4FekrIibHvZqw+8yq233BpR8K7QlT+sQOAiH39cSVlZ\nIevWTWPp0qVs2LCBtLT4HLG7gtbG9UWbWPbFwIHDeeKJ3Ywdu4Rdu37IK6/ksnjxYvx+f1wU+lgW\nd4AhfYZQM7uGtXvXMm/bPAKdrA2IZLbM7UC7a5TzL4JFu7Ntbo+w7WXF3yiOIE78a21t5eTJkxw9\nepRDhw7x/vs1VFfvoEePrzJ8+Hd58sm3mTHDuzcxx7lZSD4GDpzK3XdP4dy5vZw9+wZFRUU0NTUx\nevRohg4dyqBBg8jNzSUrK4usrKy4HTDdCP169WPX47so3FhI8abw9bKrpkJe01m+CicXAlzexxbu\ncyTbRNvmxIkTVFRUXNPznD9/nubmZs6cOUNLSwt9+/YlLy+PAQMGMHbst7nvvmVkZuYC0LODQwrZ\n2XDuXLQ91nUyMqBvX69TxNbVXm8y9sXVdGVfdL/iOjxpae2fS/TuPZxp04ZTWlrKqVOnqKmpoa6u\njvLycurr62lsbKSxsZHU1FTS09Pp3r375Y/U1FR8Ph8+nw9Jl79uf1tSVCcoPHz4MHv27Llhrz8a\nGcqgKjf8fGl1drBC0ijgB2ZWELq9ADAz+1G7bVYCO8zsrdDtQ8D9wJ2dtW33GMl91MRxHOcamFmH\n70iRjNz/AuRJ6g98AswAiq7YZjPwDPBW6M2g2cwaJZ2OoG3YgI7jOE70Oi3uZnZR0veACoIHYNeY\n2UFJTwW/bb80s99LmiTpKNACPB6ubZe9GsdxHAeIYLeM4ziOc/OJqxWq8UTS85ICkm7zOotXJP1Y\n0sHQwrTfSkr81WXtRLMAL5FJypFUJalO0n5J3/c6k9ck+STVStrsdZarccW9A5JygAlAEpyHLqwK\nYLCZDQOOAAs9zhMz0S7AS3CtwDwzGwx8C3gmifvikmeBA16HCMcV9479FEiOE9yEYWZ/sLarKOwG\ncrzME2NRLcBLZGbWYGZ7Q1//BzhIcA1LUgoN/iYBq73OEo4r7leQVAj808z2e50lzswGtngdIoau\ntjAvqUm6AxgGfOhtEk9dGvzF9QHLuDufeyxI2g5ktb+L4A9qMbCI4C6Z9t9LWGH64iUz+11om5eA\nC2b2pgcRnTghKQP4DfBsaASfdCQ9AjSa2V5JDxDH9SEpi7uZTejofklDgDuAfQouVcsB/ipppJkl\n5PlHr9YXl0j6DsF/QcfHJFD8OAl8rd3tnNB9SUnSLQQL+6/N7D2v83goHyiUNAlIA3pIWm9mszzO\n9X/cVMgwJB0DvmlmSXmRO0kFQCkw1sz+7XWeWJLUDfg78CDBBXh/BoqSdZ2GpPXAaTOb53WWeCHp\nfuB5Myv0OktH3D738Iw4/rcrBpYDGcD20LSvFV4HihUzuwhcWoBXB2xM4sKeD8wExkvyh34XCrzO\n5YTnRu6O4zgJyI3cHcdxEpAr7o7jOAnIFXfHcZwE5Iq74zhOAnLF3XEcJwG54u44jpOAXHF3HMdJ\nQK64O47jJKD/AY8w82ahirVPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a1f30ca10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# test\n",
    "xmin, xmax = -5, 5\n",
    "I = quad(pdftarget, -100, 100)\n",
    "nsamp = 1000\n",
    "\n",
    "x = np.linspace(xmin, xmax, 1000)\n",
    "y_target = pdftarget(x, normed=I[0])\n",
    "high = np.argmax(y_target)\n",
    "xhigh = x[high]\n",
    "yhigh = y_target[high]\n",
    "\n",
    "xsample, ysample, accept, l, n = sampling_piecewise(nsamples, I[0], xmin, xmax, xhigh, 1.05*yhigh)\n",
    "plot_samples(xsample, ysample, accept, I[0], xhigh, 1.05*yhigh, trace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .\n"
     ]
    }
   ],
   "source": [
    "# animation plot\n",
    "# make an animation\n",
    "divisor = 10\n",
    "for i in range(n):\n",
    "    plot_samples(xsample[0:i], \n",
    "                 ysample[0:i], \n",
    "                 accept[0:i], \n",
    "                 I[0], \n",
    "                 xhigh, \n",
    "                 1.05*yhigh,\n",
    "                 write=True,\n",
    "                 filename='plotsample_nsamp_{0:04d}.png'.format(i),\n",
    "                 trace=True)\n",
    "    \n",
    "    if (((i+1) % divisor) == 0):\n",
    "        print \".\","
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/Pk2pSMDIXww\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f1a1f30c190>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "YouTubeVideo(\"Pk2pSMDIXww\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "- SEAYAC Workshop, *MC methods in astronomy*, Tri L. Astraatmadja (MPIA Heidelberg)\n",
    "    Krabi, 4 December 2015"
   ]
  },
  {
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}