Skip to content

Instantly share code, notes, and snippets.

@tanutarou

tanutarou/gibbs_sampler.ipynb

Created September 5, 2017 09:30
Show Gist options
  • Save tanutarou/aa2fd52ba301508f91eb9be1054a3b32 to your computer and use it in GitHub Desktop.
Save tanutarou/aa2fd52ba301508f91eb9be1054a3b32 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
gibbs_sampler.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from __future__ import print_function, division\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import numpy as np\n",
"import math\n",
"import itertools\n",
"import copy\n",
"import seaborn as sns\n",
"sns.set(context='paper', style='darkgrid', rc={'figure.facecolor':'white'}, font_scale=1.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2変量正規分布(連続変数)\n",
"$$\n",
"p(x_1, x_2) = \\frac{1}{Z} \\exp \\left(- \\frac{x_1^2 - 2b x_1 x_2 + x_2^2}{2} \\right)\n",
"$$\n",
"b : x1とx2の相関の強さ\n",
"\n",
"条件付き確率:\n",
"$$\n",
"p(x_1 | x_2) = \\frac{1}{\\sqrt{2\\pi}} \\exp \\left( - \\frac{(x_1 - bx_2)^2}{2} \\right)\n",
"$$\n",
"$$\n",
"p(x_2 | x_1) = \\frac{1}{\\sqrt{2\\pi}} \\exp \\left( - \\frac{(x_2 - bx_1)^2}{2} \\right)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mcmc_gibbs_normal(init, max_iter, burn_in, b=0.8):\n",
" \"\"\" \n",
" ギブスサンプリングによるMCMC(2変量正規分布)\n",
" init : 初期状態\n",
" max_iter : 繰り返し回数\n",
" burn_in : 捨てる初期サンプル数\n",
" b : 相関の強さ\n",
" \"\"\"\n",
" samples = []\n",
" x = init\n",
" for i in range(max_iter):\n",
" # 交互に条件付き確率からのサンプリング\n",
" if i % 2 == 0:\n",
" x[0] = np.random.normal(b*x[1], 1)\n",
" else:\n",
" x[1] = np.random.normal(b*x[0], 1)\n",
" \n",
" # サンプル列に追加\n",
" samples.append(copy.deepcopy(x))\n",
" return samples[burn_in:]"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"b = 0.0\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD7CAYAAABNPKDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGwpJREFUeJzt3X1QU2e+B/DvCfUdL4G2SbgdhfKipBrstFe5viwFKiId\nLVbxhTqMnY4y2/pStL27s3JlR8ftTne0bmf8o1Pv3Lm9vTuOy7ajs86U1kEQdVrpdreE3oIiGui1\nJKkGkAiCwrl/0KSEJBA4OeSc5Pv5y+Q8hOeU8uU8z3nO8xNEURRBRDRBmlB3gIjUjSFCRJIwRIhI\nEoYIEUnCECEiSRgiRCSJ5BDZs2cPFi9ejL1797rfM5vNWLNmDXJzc3H8+HGfX9fW1ob169cjNzcX\n5eXl4J1mInWSHCJbt27FO++84/HeoUOHcOzYMVRWVqK2thZXr171+rojR46gtLQU586dQ2dnJ2pq\naqR2hYhCQHKIZGRkYNasWe7XNpsNoigiNTUVUVFRWLt2rVdAiKKI+vp6ZGZmAgDWrVuH6upqqV0h\nohAI+pyI3W6HXq93vzYYDLDZbB5tOjo6oNVqR21DROrwSKg7EKiHDwcQrtMmUVECBgbC7+R4Xuoz\nZUrUuL8m6CGi0+k8riqsVit0Op1Hm9jYWHR2do7aZiRRBDo7e4LbWYXQameG5bnxvNTn8cdnj/tr\ngj6ccQ1lmpubMTAwgLNnzyI7O9ujjSAIMJlMqK2tBQCcPn3aqw0RqYMg9SnekpISmM1m9Pb2IiYm\nBu+//z76+/tRVlaGvr4+FBQUYPfu3QCAsrIybNmyBSaTCRaLBfv27cPdu3exdOlSHDx4EBqN/0x7\n8GAgbNM/XP+y8bzUZyJXIpJDZLIwRNSH56U+ihjOEFFkYYgQkSQMESKSRDXrRCj0RFFEU2sHLLZu\nJOpnIy0hFoIghLpbFGIMEQqIKIr4sLIJl8ztGBQBjQCsSI/HK/nGUHeNQozDGQpIY2uHO0AAYFAE\nLpnb0WhxhLZjFHIMEQpIq63bHSAugyJgsXWHpkOkGAwRCkiifjY0I6Y/NMLQ+xTZGCIUkLSEWKxI\nj3cHiWtOxJgYF9qOUchxYpUCIggCXsk3IsOod9+dYYAQwBChcTImxjE8yAOHM0QkCUOEiCRhiBCR\nJAwRIpKEIUJEkjBEiEgShggRSSLbOpH6+nqUl5e7Xzc3N+Pjjz+G0fjzU585OTmIjo6GIAjQ6XQ4\nceKEXN0hIpnIFiKLFi3CmTNnAAC3bt1CcXGxR4C4VFRUYNq0aXJ1g4hkNinDmcrKSuTl5U3GtyKi\nSTZpIZKfn+/zWFFREQoLC1FZWTkZXSGiIJP92Zlbt27B4XAgPT3d69jJkyeh1+ths9mwbds2GI1G\nJCQk+PycqCgBWu1MubsbElFRmrA8N55XZJA9RD777DO/QxlXtTy9Xo/ly5ejsbHRb4gMDIhhW+sj\nXOuY8LzUR5F1Z/wNZXp6euB0OgEA3d3dqKurQ3JystzdIaIgk/VK5IcffoDD4YDJZHK/t2PHDhw+\nfBj9/f3YuXMngKFNgLdu3YrU1FQ5u0NEMmAZTQUI18tjnpf6KHI4Q0ThjSFCRJIwRIhIEoYIEUnC\nECEiSRgiRCQJQ4SIJGGIEJEkDBEikoQhQkSSMESISBKGCBFJwhAhIkkYIkQkCUOEiCRhiBCRJAwR\nIpKEIUJEkjBEiEgSWTdqXrBgAVJSUgAACxcuxO9+9zuP42azGfv370dfXx8KCgqwa9cuObtDRDKQ\nNUS0Wq27Hq8vhw4dwrFjx5CUlISioiLk5uZi/vz5cnaJiIIsZMMZm80GURSRmpqKqKgorF27FjU1\nNaHqDhFNkKxXIl1dXVi/fj2mTZuG0tJSZGRkuI/Z7XZ3BTwAMBgM+OKLL/x+Fstoqg/PKzLIGiJV\nVVXQ6/W4fv06SkpKcObMGcyePf66FgDLaKoRz0t9FFd3xnWlkZKSgnnz5sFisbiP6XQ62Gw292ur\n1QqdTidnd4hIBrKFSFdXF/r7+wEMzX9cu3YNc+bMcR93BUxzczMGBgZw9uxZZGdny9UdIpKJbMOZ\nlpYWlJeXQ6PRQKPRYP/+/dBqte5avHq9HgcOHEBpaan7Fi/vzBCpD2vxKkC4jrF5XuqjuDkRIgp/\nDBEikoQhQkSSMESISBKGCBFJwhAhIklkXfZO4UsURTS1dsBi60aifjbSEmIhCIJXm0aLY9Q2pH4M\nERo3URTxYWUTLpnbMSgCGgFYkR6PV/KNHm3e/6QBVX9r89uGwgOHMzRuja0d7gABgEERuGRuR6PF\n4dHm/Ndto7ah8MAQoXFrtXW7w8FlUAQstm6PNgODo7eh8MAQoXFL1M+GZsTUhkYYen94myjN6G0o\nPDBEaNzSEmKxIj3eHSSu+Q5jYpxHm5xn547ahsIDH8BTALU+0DX8zouvcNBqZ+KLb/5v1DZqpNaf\nVyAm8gAe787QhBkT48YMhkDakLoxRCjoXGtIrF0/wBAz3e8akrHWmZA6MEQoqAJdQzJWG1IPTqxS\nUAW6hmSsNqQeDBEKqkDXkIzVhtRDtuHMjRs3sH//fjidTkyZMgW/+c1vsGTJEo82Y5XZJPVxrSEZ\nHhK+1pCM1YbUQ7YQmTZtGt5++20kJSWhpaUFr732Gj7//HOPNmOV2ST1ca0hGTnfMXINyVhtSD1k\nC5EnnnjC/e+kpCQ4nU6IosgZeIUYfnckQRcNAGi1Oyd8p2T452UY9ViSpoO16z76ex/A2tGDY6f+\nAUEQoI2eiiVGPbatTkOGUe+xhoR3bNRpUu7OVFVV4amnnvL6H2K0MpsjsYxm8LiesD3/dZvX8y1R\nGiDn2bl4bUP6hD8vSgNkPzsHgiDgXF2bV/tasxW5i4e+x9JRPmO8/ZgsLKPpSfYVq7du3cKrr76K\nDz74AAkJCR7HbDZbwGU2uWI1eL6zOPDuqW+8JjddNALw5uanRx1eDL9qEETgLxdaPD5PEICx/s96\na/MiCILg/oyKmhYM/5JA+hEKXLHqSdYrEafTiddffx0HDhzwChDAd5lNk8kkZ5cIvu+ODOe6U+Lv\nl3fkOg/fbcbux8e1N9Bq9d+XsfpByiDbLd6BgQG88cYb2Lx5M1asWOF1fKwymyQfX0/hDicIgN3R\ng0aLA74uVEeu8/D3GWO52T56mPGOjTrIdiVSW1uLL7/8Erdv38apU6cAAB999BGKi4tx5swZv2U2\nSX4j746MJIrAhfp21Na34xeLhlaSuoYvN9q78HXTj2MGyApTPJz3+/GPa3f8thntakUjAMtNBgDA\np1daOdGqYHyKVwFCNcYe/hSuKIr4n3PXYHX0erQRBOAXJgM6nP1ouBH4ilKddjrsnffH3acnDbMx\nVx+Nf5n/OL66+qMil8ZzTsQTn52JYGkJsQDgntgcGSDA0NVCrdk67s+eSIAAwE1rN1pt3bhz9z6+\nvdnhfn9QBGrr26HXzsDqf03gFYmCMEQiVCCTo6EyKMIjQIaruHADts5eRVyR0BA+OxOhvrM4UFuv\nvAAJRG09H9ZTEoZIBBJFEZ/U3gh1NyThw3rKwRCJQI2tHbBYR/8l3PRcEjIXxY96KziUBDVeQoUp\nhkgEarV1j7kYTNQIeCXfiDc3Pw1TkvIWe3111R7qLtBPOLEagXw9ij+cRgAEcWh9RoIuGlDgH/2b\nVif+98ZtfG+/h6a2DqTNjcWqJXOg0fDv4mTjOhEFmOx1B647Mxfr233mgyFuBuwdvYqfdI2e8Qic\nvQ/drw1xM/B2ydJRviI4uE7EE2M7AgnC0FBlY1ay9zEMrRdReoAA8AgQYKjflV+2hqg3kYshEsFE\nH5OmKsiOUTW1+V5fQvJhiEQwXw/iqX0haNrc2FB3IeJwYjVCuB6gu2m9Cw0EDEJEon42VpjicanB\n8/mUNls3LFZnqLvsl0YAli3U4/qtux5L9Q1xQ0vih+NuafJjiEQAfxOprtDYt2kRvmoaumX67LzH\n8LUCb5/qtNNgnBuHDmcf0uZoMdcwG/GPRcN+pwcOZx+Mc2N9Bgjr28iPIRIBXPt/jJzvcD3UZm65\njU7nAwBDWwAokb2zD/bOob6Zhz1NLABYmBQHURh6Knn4lYa/+jYZRj03OgoihkgEGGsnM1eAqJEI\noOGGAw03HF5XGhbrXb/1bRgiwcOJ1Qgw1k5m4WJkJT0NfJ+0oPZbUArDEIkArp3MIiBHPCrpDfq5\nYe3r1jZNHIczEcC1uMxV58V+pwcXzMqc+wgG176sTxr+iZX2JoGsVyLV1dXIy8vDqlWrUFFR4XXc\nbDZjzZo1yM3NxfHjx+XsCgEwJsYhPyMBxavnQx83I9TdkZ3rCsw1lGOlPXnIdiXy8OFDvPPOO/jo\no48wa9YsbNiwAStXrkRs7M+LgQ4dOoRjx44hKSkJRUVFyM3Nxfz58+XqEv2kqa0Tdh9bIYaL4ROn\nGUa9exi3OE2Hp558NHQdC1OyhYjZbMa8efOg0+kAAFlZWbh8+TLWrFkDYKhMhCiKSE1NBQCsXbsW\nNTU1DJFJ0GrrVv3y9tEI4tAakf/6tBEXh+0POyiKDBEZyBYidrvdXZwKAAwGA2w226jHv/jiC7+f\nxzKawbMg+TH8pbolbINk2owpaL19zyNAAOCi2YrnlyQgPeVxSZ/PMpqeVDOxOjAghu3j15P9aPmc\nR2diYVLcuEpAqElf7wPU/O17n8eqv2rD3MdmSfp8bgXgSbaJVZ1O53HlYbVa3UObQI5TcImiiEaL\nA59eaUVTawfeKExHugJ3LAuUIACJhmifx0QBEXE7WylkC5H09HRcvXoVdrsd9+7dQ3V1tUc5TddQ\nprm5GQMDAzh79iyys7Pl6k5Ecz1DcvTUN6iobsHRU9/gvz+7ij2F6TDEqvMujS52BjZmpXgtonPd\nwl2c5vsP0hKj3uf7NHGyDWceeeQR/OpXv0JxcTEGBwexfft2xMbGYseOHTh8+DD0ej0OHDiA0tJS\n9PX1oaCggJOqMvH3DEmbtRvWDnXepbF39EIURY9yoMNv4YqiiMxF8bhoboco/lTJj7d3ZcHtERVA\n7jH2p1daUVHdItvnh8rG7GTkZyR4lAMdGRKjHZsozol4Us3EKk3cWBszq5Xd0eN+ctdfQBgT43j1\nITM+OxMBfK3cfNKg/qXfF+rbcfTUN/iwsinUXYlovBKJACOfnUnUz8ZN613cHKOAlRpwj5DQY4hE\nkJGX9uEyxOEeIaHF4UyEmj9XC51Kb+/6wj1CQodXIhGqqa0TNpXe3vXl3E8rVFkFb/Lxv3aECqQe\nr5p0OPvx55oW/Pt/XAl1VyIOQyRCDF/23mhxIEEXrfoaM76wCt7k43AmAoiiiA8/bcLFhp9Xb65Y\naMAvTPGoDcMdzpraOrzKR5B8GCIR4DuLwyMsRBG42GDFv215GhlP6VHXOPQgpF47A3++cCNU3Qwa\nVsGbXAyRCOAqTDVSXaMN2/KN7lujoijC2tnrft5EjXxVwSN5MUQigL+pj47uPoii6C725FqUZoid\niT/XqONZG40AJBpmI3rGFKT5qIJH8mOIRIDFaTqfle0abjrwYWWTV1lJpZZUeC49HrpHZ7p3aw/2\ng3U0MQyRCGBMjEPmonjUjggS8acl40vSdBAEwf1LmaDzvdlPqOkenYn8jJ+vNBgeysAQiQCuYYoA\n71q7gyLw8YUWtNqc7j05lpviETt7Kjq6+0PTYT+4KlWZuE4kgiwx6r12AhME4KbV6blhUUO74gIE\nUO4wK9IxRCKIry0BEn1sCTCeOzPpSXEozErCrBnyX9Sycp0ycTgTQXxtCQAAR0994/E0ryAEHiR5\nS+ZCBHCv92HwOzxMehI3F1IqWUKko6MDpaWl+PHHHxEVFYXXX38d+fn5Xu1ycnIQHR0NQRCg0+lw\n4sQJObpDIwzfEsDfPqUAPPZlHU2rTfq+JPrYGX4fCDQlxaF009OSvwfJQ5YQ0Wg0ePPNN5Geno47\nd+7gpZdeQnZ2NqZPn+7VtqKiAtOmTZOjGxQAX1cnw0tQWmzdEET4XTfi+hope5NkLorHttVpeK+i\nHuafauEIPw21Cp9L5hWIwskSIjExMUhPTwcAPProo9Bqtejq6vIZIqQMIzcsGr5/91x9NExPxqLh\nZofH12gEIEEXDVEUkaCPxk2rc9zfd1NWsnuBWOmmp2XZWJnkJfucyHfffYfBwUGPkpnDFRUVQaPR\nYPv27Vi9erXc3aEAuOrUeAxxTPEwDauaN3Qr2IC6JnvAwx6f32vEHRdurKw+Ew4RV2HukT755BNM\nnToVAHD37l38+te/xqFDh3y2PXnyJPR6PWw2G7Zt2waj0YiEBN/LllmLd/LUX/8Rlxs869Rc/rYd\nB17NwIYcoOVWF5KfiIEI4PB/Xhk1QKI0wNa8NPxw5x6qvvreY8I2SjNUF1hJ5x4Ipf28Qm3CIXL2\n7NlRjz948AC7d+/Gtm3b8Mwzz/hs47o60ev1WL58ORobG/2GCGvxTp7vWm5jYNDzvYFB4H9bbiM/\nI8Fdy/bTK61e7YaL0gwtXMta9M8QRRH9fQ89rm6Wm+Ix97FZijr3QCjt5xVMiqo7c/DgQSxcuBCF\nhYU+j/f09GBwcBDR0dHo7u5GXV0dtmzZIld3aBx8TZS6ylMG0q7wuWSIwtBVhitwRpvAJXWTJUSu\nXbuGiooKzJ8/H5cuXQIAHD16FCkpKe4ymv39/di5cyeAoTH41q1bkZqaKkd3aJxci9J8lacMpJ1r\notTXX2zOeYQfltFUAKVeHgd6p8RfO6Wel1Thel6AwoYzpH6BXjXw6iKy8dkZIpKEIUJEkjBEiEgS\nhggRScIQISJJGCJEJAlDhIgkYYgQkSQMESKShCFCRJIwRIhIEoYIEUnCECEiSRgiRCQJQ4SIJGGI\nEJEkDBEikkS2nc0CKZHZ1taG0tJSdHd3Y+nSpTh48CAEgaXfidRE1u0RxyqReeTIEZSWliIzMxN7\n9uxBTU0NsrOz5ewSEQVZyIYzoiiivr4emZmZAIB169ahuro6VN0hogmSNUSKiopQWFiIyspKr2Md\nHR3QarXu1waDATabTc7uEJEMZCujOZ4SmYFgGU314XlFBtnKaI5VIjM2NhadnZ3u11arFTqdzu/n\nsYym+vC81GcidWdkGc709PTA6XQCgLtEZnJyskcbQRBgMplQW1sLADh9+jQnVYlUSJYQuXPnDl5+\n+WW8+OKLePnllz1KZJaVlaGhoQEA8NZbb+GPf/wjVq5ciZiYGGRlZcnRHSKSEctoKkC4Xh7zvNRH\nMcMZIoocDBEikoQhQkSSMESISBKGCBFJwhAhIkkYIkQkCUOEiCRhiBCRJAwRIpKEIUJEkjBEiEgS\nhggRScIQISJJGCJEJAlDhIgkYYgQkSQMESKSRJYKeO3t7fjlL3/pfm2xWHD06FGsXLnSo11xcTFu\n376NqVOnAgDOnDkjR3eISEayhEh8fLw7EHp7e5GTk4Nly5b5bHv8+HGvneCJSD1kH85cuHABixcv\nxsyZLPZDFI5kD5HKykqsXr3a7/G9e/fipZdewp/+9Ce5u0JEMphwyYixymgCwP3795GTk4OqqirM\nmDHDq63NZoNer0dXVxdKSkqwb98+ZGRk+PzcwcFBDAyoorrFuEVFaTAwMBjqbgQdz0t9pkyJGvfX\nyFZGE/h5KOMrQICfS23GxMQgLy8P3377rd8QYRlN9eF5qY/i6s5UVlYiPz/f57GHDx/C4XAAAPr7\n+3Hx4kWkpKTI2R0ikoEsd2eAoaFMXV0dfv/733u8X1ZWhi1btiA5ORnbt2/HgwcPIIoiVq9ejeee\ne06u7hCRTFhGUwHC9fKY56U+ihvOEFH4Y4gQkSQMESKShCFCRJIwRIhIEoYIEUnCECEiSRgiRCQJ\nQ4SIJGGIEJEkDBEikoQhQkSSMESISBKGCBFJwhAhIkkYIkQkCUOEiCRhiBCRJAwRIpJEcogcPnwY\ny5Ytw6ZNmzzeb2trw/r165Gbm4vy8nL42srV4XCguLgYq1atwq5du9DX1ye1O0Q0ySSHyAsvvIAP\nPvjA6/0jR46gtLQU586dQ2dnJ2pqarzanDhxAi+88AI+//xzzJkzBxUVFVK7Q0STTHKIPPPMM9Bq\ntR7viaKI+vp6ZGZmAgDWrVuH6upqr6+trq7G2rVrR21DRMomS92Zjo4Oj2AxGAyw2Wxe7e7du4fo\n6OhR27hMmRI1oe3s1SJcz43nFf7GDJFAau4SUeQaM0QCqbk7UmxsLDo7O92vrVYrdDqdV7uZM2fC\n6XQiOjrabxsiUjZZbvEKggCTyYTa2loAwOnTp5Gdne3VLisrC3/9619HbUNEyia5jGZ5eTnOnz+P\nzs5OxMXF4be//S2ef/55WCwW7Nu3D3fv3sXSpUtx8OBBaDQavPfee1i4cCGef/55OBwO7NmzBzab\nDampqXj33Xcxffr0YJ0bEU0C1dTiJSJlUsWKVSkL2tQkJycHL774IgoKCrBjx45Qd0eS6upq5OXl\nYdWqVWG1/mfBggUoKChAQUEBysrKQt0dSfbs2YPFixdj79697vfMZjPWrFmD3NxcHD9+PLAPElXg\n66+/FhsaGsSNGzd6vL97927xwoUL7n+fP38+FN0LmuzsbPH+/fuh7oZkDx48EPPy8kSbzSY6nU4x\nLy9PdDgcoe5WUCxbtizUXQiaL7/8UqyqqhJLS0vd723YsEG8du2a+PDhQ3Hjxo1iU1PTmJ+jiisR\nKQvaaPKZzWbMmzcPOp0Os2bNQlZWFi5fvhzqbtEIGRkZmDVrlvu1zWaDKIpITU1FVFQU1q5d63Ol\n+UiqCBFfAl3QpjZFRUUoLCxEZWVlqLsyYXa7HXq93v06XH42ANDV1YX169ejqKgIV65cCXV3gmqi\nPzdZVqyOV6QsaBvrPE+ePAm9Xg+bzYZt27bBaDQiISFhkntJo6mqqoJer8f169dRUlKCM2fOYPbs\nyF69qogQkXNBm5KMdZ6uvwJ6vR7Lly9HY2OjKkNEp9N5/AWzWq1YsGBBCHsUPK6fUUpKCubNmweL\nxQKTyRTiXgWHr59bIL9Tqh3OBLqgTS16enrgdDoBAN3d3airq0NycnKIezUx6enpuHr1Kux2O+7d\nu4fq6mqsWLEi1N2SrKurC/39/QCG5g+uXbuGOXPmhLhXweMKyObmZgwMDODs2bMB/U6pYp3IeBe0\nqdH333+PnTt3AhiaNN66dSu2bNkS4l5NXFVVFf7whz9gcHAQ27dvx+bNm0PdJcn+/ve/o7y8HBqN\nBhqNBrt27cLKlStD3a0JKykpgdlsRm9vL2JiYvD++++jv78fZWVl6OvrQ0FBAXbv3j3m56giRIhI\nudT5Z5uIFIMhQkSSMESISBKGCBFJwhAhIkkYIkQkCUOEiCRhiBCRJP8PH5WLVwZVLywAAAAASUVO\nRK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc1a92b99e8>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"b = 0.25\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD7CAYAAABNPKDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHNRJREFUeJzt3XtQVFe+L/DvbuITCHTHdDeTERBEabExJ3OR+AgKEdGU\nBmM0SryM56YiNRMfg2ZqpkaOTGk5mckcHSdV/pGK9585uXM9GU9SWmPVIcklIJqKcmfmhsYTfERt\nSSm90TQgiILS+/5hukM/aVi96df380/SvTfN2iF8WWvttddPUhRFARHRGGnC3QAiim4MESISwhAh\nIiEMESISwhAhIiEMESISIhwiO3bsQEFBAXbu3Ol6z2KxYNWqVSgtLcXhw4d9fl17ezvWrl2L0tJS\n1NbWgneaiaKTcIhs2rQJb7/9ttt7+/btw6FDh1BXV4empiZcvHjR6+sOHDiA6upqfPrpp+ju7kZj\nY6NoU4goDIRDpLCwEImJia7XsixDURTk5OQgISEBq1ev9goIRVHQ0tKCoqIiAMCaNWvQ0NAg2hQi\nCoOQz4l0dnbCYDC4XhuNRsiy7HZOV1cXUlNTA55DRNHhsXA3IFgPHw4hVqdNEhIkDA3F3sXxuqLP\nhAkJo/6akIeIXq9361XYbDbo9Xq3c7RaLbq7uwOe40lRgO7u/tA2NkKkpk6NyWvjdUWfJ59MHvXX\nhHw44xzKXL58GUNDQzh58iSKi4vdzpEkCWazGU1NTQCA48ePe51DRNFBEn2Kt6qqChaLBffu3UNK\nSgreffddDA4OoqamBgMDAygvL8f27dsBADU1Ndi4cSPMZjOsVit27dqFO3fuYMGCBdi7dy80Gv+Z\n9uDBUMymf6z+ZeN1RZ+x9ESEQ2S8MESiD68r+kTEcIaI4gtDhIiEMESISAhDhIiEMESISAhDhIiE\nMESISAhDhIiEMESISAhDhIiEMESISAhDhIiEMESISAhDhIiEMESISAhDhIiEMESISAhDhIiEqFYy\noqWlBbW1ta7Xly9fxocffgiTyeR6r6SkBElJSZAkCXq9HkeOHFGrOUSkEtVCZN68eThx4gQA4MaN\nG6isrHQLEKdjx45h0qRJajWDiFQ2LsOZuro6lJWVjce3IqJxNm4hsnLlSp/HKioqsG7dOtTV1Y1H\nU4goxFQvo3njxg3Y7Xbk5+d7HTt69CgMBgNkWcbmzZthMpmQkZHh83MSEiSkpk5Vu7lhkZCgiclr\n43XFB9VD5OOPP/Y7lHFWyzMYDFi0aBHa2tr8hsjQkBKztT5itY4Jryv6RGTdGX9Dmf7+fvT19QEA\nent70dzcjOzsbLWbQ0QhpmpP5ObNm7Db7TCbza73tmzZgv3792NwcBBbt24FACiKgk2bNiEnJ0fN\n5hCRClhGMwLEaveY1xV9InI4Q0SxjSFCREIYIkQkhCFCREIYIkQkhCFCREIYIkQkhCFCREIYIkQk\nhCFCREIYIkQkhCFCREIYIkQkhCFCREIYIkQkhCFCREIYIkQkhCFCREIYIkQkRNWNmvPy8jBz5kwA\nwNy5c/Gb3/zG7bjFYsHu3bsxMDCA8vJybNu2Tc3mEJEKVA2R1NRUVz1eX/bt24dDhw4hKysLFRUV\nKC0txezZs9VsEhGFWNiGM7IsQ1EU5OTkICEhAatXr0ZjY2O4mkNEY6RqT6Snpwdr167FpEmTUF1d\njcLCQtexzs5OVwU8ADAajfjiiy/8fhbLaEYfXld8UDVE6uvrYTAY8PXXX6OqqgonTpxAcvLo61oA\nLKMZjXhd0Sfi6s44exozZ87ErFmzYLVaXcf0ej1kWXa9ttls0Ov1ajaHiFSgWoj09PRgcHAQwKP5\nj0uXLmH69Omu486AuXz5MoaGhnDy5EkUFxer1RwiUolqw5krV66gtrYWGo0GGo0Gu3fvRmpqqqsW\nr8FgwJ49e1BdXe26xcs7M0TRh7V4I0CsjrF5XdEn4uZEiCj2MUSISAhDhIiEMESISAhDhIiEMESI\nSAhDhIiEMESISIiqD+ARBUNRFFy43gWr3IsMfRIA4HpnHzINycjN0EKSpDC3kAJhiFBYKYqCP9Vd\nwBlLBxwea6c1ErA4Pw3/vNIUnsZRUDicobBqu97lM0AAwKEAZywdaLPax79hFDSGCAVFURS0We34\nz3PX0Wa1Q/SRK+fnfdLc7jNAnBwKYJV7hb4XqYvDGRqR55BDdJgRaAjjSSMBmYaxbWRF44M9ERqR\n55BDdJgRaAgznDOsTJm6MX0fGh/sidCIrsu9Xr/wzmHGWH7BfX0eAORn6VA2Px3Ao8/ONCQzQKIA\nQ4RGlGlIhkaC2y++yDDD3+eVzU93hQbDI3pwOEMjys3QYnF+GjTfLdeQAOTN0MEq97omWUcz8Zqb\nocVicxqcyz+kIIctoZ7cpdBgT4S8DF/85Vzw9c8rTSg0GXDNdgeX2rtx/podrVft0EjAIrMRkiS5\n5jkkAHOzdHh6th6GxycBcF88BgCQAGcGfP9P7+/rXGgW6sldCh3VQuTq1avYvXs3+vr6MGHCBPzq\nV7/C/Pnz3c4ZqcwmjT9fv6yLzEY8O8cIq9wLDSS0XrXD2QdwKMBpiw3S8FAA0Hr1UcgM5/zFn28y\n4Iylw+3Y6ZYOdPcOuH22OUuH6vXzIEmS38ndQpOBQ58wUy1EJk2ahLfeegtZWVm4cuUKfvrTn+KT\nTz5xO2ekMps0/nz9sp622HCm1YZAo4dgRhbOX/zu3gGviVUFgMUjdFqv2vHOsRZUv/J0yCd3KXRU\nC5GnnnrK9e9ZWVno6+uDoih8DiLC+btzEqrpB4fiHRaBWK7a0Wa1h3xyl0JnXOZE6uvrMWfOHK8A\nCVRm0xPLaI6PvOxp+OjUFQw5wt2S79l67mNNUTaev2LHZ39vx5ADSNAAJT9Kx4Knfzju7Ymkn1ck\nUL1kxI0bN/Daa6/hvffeQ0ZGhtsxWZaDLrPJkhHjYzSrSceDRgLe3PC0a8jSZrUHtYYk0CStqEj6\neYXaWEpGqNoT6evrwxtvvIE9e/Z4BQjgu8ym2WxWs0k0AkmS3O7E/K2tE1a5LySfvSQ/Dadb/YeT\nBCAzLRnXOh49K+NrxaopUxfUrWDeyRk/qq0TGRoaws9+9jNs2LABixcv9jo+UplNCi9Tpg6ZxsdD\nFiDP5RtRYNIjw5Dk87hGAp6bl4Z/+fF/wytLs5GfpcO6JdnYvCJ31N8r1Mv0KTDVeiJNTU04e/Ys\nbt++jQ8++AAA8P7776OyshInTpzwW2aTIofVdkf4M8wztMhN1+Jvl27hTGuLa4JWlzwRc2foYNQl\nwgEFM4yPIzdD69aDOH/NDpu9H4VzDKPasIh3csYXy2hGgEgcYyuKgneOtYzqToqn5KkTMO3xSbhm\n89+bkSQgw5iM+bP1uPltHz5vlRHM/5CBhihtVjsOfvCl152c4XMrIiLx5xUqETcnQtGr7XoXzl/z\nDpCJj0kYfBjc353e/gfo7X8Q8BxFAawdvbB2jG7PEIcCNLV0wKidirLCdLceiXOZvuecCHsh6mCI\nxKmR7l74Wy8SbICMl780XoGtq9+tRzJ8cphPA6uPIRKHgrl74WtxV6RqaumABGC+yeAWhsHcySFx\nfIo3DrVd78LpFo+l7S3udy88n9yNdKdaOnDwgy/xp7oLAc/jk8Chx55IHLLa7nhNXip4dPciN0Pr\nGuYUzH4Sliu30d0XeF4jUoz0UB7Xj6iDIRKHNPDdvZBv92H3kbOQ7ffGuUWh41CAj5vbAcBrnodP\nAquDw5k45PBzE7WpVY7qAHGyXLX7HNoEWj9CY8cQiUMzjI9HzVzHWPlapeqcLB6OTwKLY4jEodwM\nLRaZjeFuRkhIElzL5D159jI8J4u5fiQ0OCcShyRJQkGuHqcttnA3RVhK0kSUFaYjw5iM89fsAfcb\n4foRdTBE4pCiKPjf/+dyuJsREt29g3jnWAtmpacib4YO57/bXjFQL8Pf+hGHw4FPmr/BhfYu5KZr\nsXz+dGg07KyPhCESh9qud8EWAxOoTpardliu2iFJj/ZlnZ2hHXUvw+Fw4F/+5znXfxfLVTuaLDfx\nVtUCtZodMxizMczfwqrrMXo3Qvlu68VAAeLvv8nHzd94BavNfg91Z6+r3u5ox55IjPJcWCVJQKYx\nGeuWZCP9ycRwN09VzW3yqBebXWzv8vlZF9q7sOJZ7w216HsMkRjlubBKUYBrHb3413//EkbdlPA2\nLkwCLTbLTdf63PYgN107zq2MPhzOxCh/T+ECiJn5EIN2ss/355sMPoctgRabLZ8/3Stcjbop7IUE\ngT2RGBVNT+GOldx13+u9onlpXjukOYcthSaD37ITGo0Gb1UtQN3Z6667MwyQ4DBEYpRzYdVpS0fI\nasZEKkkCir6rrGfK1OErq93nsGV+rn7EzYpWPJvB8BglVUOkoaEBv/vd76AoCrZs2YL169e7HbdY\nLNi9ezcGBgZQXl6Obdu2qdmcuOJcWDU/V48PT11xbVGokQC9dgpk+72gtiGMBooC6HVTXWHgb9hy\nvbOPi81UoFqIPHz4EG+//Tbef/99JCYm4uWXX8ayZcug1X4/UbVv3z4cOnQIWVlZqKioQGlpKWbP\nnq1Wk+LSnBlPYM6MJ7zqtThfSwqgSIDkUPCXU1fD3dwxG74ydaRqedysKLRUm1i1WCyYNWsW9Ho9\nEhMTsXTpUnz++eeu47IsQ1EU5OTkICEhAatXr0ZjY6NazYl7pkwdVhZmIDdD6xYoZYXpKCuYjq/8\n3OKMBjOMSW6hwGdkxpdqPZHOzk5XcSoAMBqNkGU54PEvvvjC7+exjKY4RVHw7ketrlKUGgl4etaT\nsNnv4uat6Ny9XJKA5/7ph2houYnsp1Jgzp4GSZJQXfEjlBTcwpUbPch+KgX5M58M2fdkGU13UTOx\nOjSkxOw2/aEsQeBrA2YAuHC9C+faZDS1dLjOdSjAPy7eCsn3DRd96hT8r7oLPncqS5+WiPRpjxbW\nhfL/HZaMcKdaiOj1ereeh81mQ15eXsDjer1erebEBV8rMheZ09Ddex+t16J3uBKI3PX9mhfuVBYe\nqs2J5Ofn4+LFi+js7MTdu3fR0NDgVk7TOZS5fPkyhoaGcPLkSRQXF6vVnLjga0XmaUtHzAaIL9yp\nbPypFiKPPfYYfvGLX6CyshJr1qzBa6+9Bq1Wiy1btrh6IHv27EF1dTXKysqwePFi3pkRFGiVarzg\nTmXjj2U0I0Coxti+ykf6M8OQiOTESV4b+USz8dq9nXMi7qJmYpVG5lk+MpACkxErns3AV9e+xX+c\nugJrgHq5kcyom4LK5bO5eCyMGCIxZPj2f81tMk4NuxPjyYFHD6ida5OjJkASJz+G7B88Dm3SRHT1\nDbo938LwCB8OZyKAGt1jRVGw/9/+hmt+CmVnGpOiJjw8Fc0Lb8EpDmfccSuAGCVJEtYtyfZ5LEFC\n1AYI8Kj27vBSEE4skRkeHM7EsNwMLYrmec+RDMXA75ZV7nUbwrBEZviwJxLDnHMkb254GkXz0sLd\nnJDyvI3rb9cyXz0WCi2GSBwwZepg0MXWsx7Xbb1B71pG6uJwJk5E6wIso26Kz+0c/9J4BUBwu5aR\nutgTiRO5GVo8lx9dpTM1ElCU/4OAdYOdwxZFUfj4f5iwJxInJEnC/3hhDp6dY0RdcztafexsHmkW\n56ehrDAdtq7+gAvouGtZeDFE4owpUwer3BvRIWLQTcGPl892hcDwcJAU4D9OXeGuZRGEw5k4FOnF\nq4rm/cDr9q1TuiEJi8xGDlsiCHsicei6HLkLzSSPyVB/6z/e3PA0hy0Rgj2ROOJc0Xmm1fczNamJ\nj4W9Op55hs61Gxvgf/0HAKwszGCARACGSJxw/kU/+MGXfivgLS/IwFtVC7AkP3wL085fs+NPdRdc\nr7n+I/IxROKE5190T8NLRv54ZS6K5qUhwJ3VoCVNGd2I2XOlqbP8w3Bc/xFZGCIxxt9DaP52PTPq\npuCVpdl4q2qB2/uFJgPmZokNFSQJ+NEs37usBxo2De9pOPdIGZ4jeR5DHgovTqzGkEAPofkr6FQ5\n7Faqr88QMXligs8naZ2LyDxv1Q4/7uxpSJKEzSty0d07AMt3t6X/67shDx+uiwyqhEhXVxeqq6tx\n69YtJCQk4I033sDKlSu9zispKUFSUhIkSYJer8eRI0fUaE7c8DcJWWgyeO165uvWqKIo+Phcu1tZ\nCRH3BobQZLF5vf+kdgp++ORUzJ2hcwWDk7Ndw4tsScqjuRIn7uoeWVQJEY1GgzfffBP5+fn49ttv\n8dJLL6G4uBiTJ0/2OvfYsWOYNGmSGs2IO4EmIU2ZOq8Vnc5f1Gu2O5AU4P9eugWrn02MQkm238Mf\njrW6vZdpTEaBSe9q10i9oeHXReGlSoikpKQgPz8fAPDEE08gNTUVPT09PkOEQmekGrTA9ys6ncOW\n0y0dEVHYu13uxfql2TBl6vCV1T7icIqTq5FD9YnVr776Cg6Hw61k5nAVFRVYt24d6urq1G5KzBtN\nDVrn0CcSAgRwn0wdqfQFV6lGljH3RFatWuXz/Y8++ggTJ04EANy5cwe//OUvsW/fPp/nHj16FAaD\nAbIsY/PmzTCZTMjIyPB5LmvxBifYGrRyz80RJ04lAJUrc5H1VApOnL6K/6diyc0EDZCXPQ2pqVOR\nlz0NH526giGH+/FNZblQgJDX1h11W1mL182YQ+TkyZMBjz948ADbt2/H5s2b8cwzz/g8x9k7MRgM\nWLRoEdra2vyGCGvxBi+YGrTGlMleQx9P5iwdluSnQZIkLPunp/DlpVsYzbalEuDq6TyWIOGhn30Z\nneU+06cloru7H9OfmIpF5jSvcqBL5/3A9TXh/H+BGzW7U+0W7969ezF37lysW7fO5/H+/n44HA4k\nJSWht7cXzc3N2Lhxo1rNIQ/OoU+gOzHnh91Kzc3QItOY7Hf3eE9F89JQXDAd/3byK1yz9fkNkPws\nHcrmp7sNTYaXvuDzMZFPlRC5dOkSjh07htmzZ+PMmTMAgIMHD2LmzJnYsmUL9u/fj8HBQWzduhXA\no1uLmzZtQk5OjhrNIR+cv6jzc/X4sOmqz3BwKI92Vjdqp6KsMB3rlmTjX//9y4Cfm6abgv/+3dqT\n67fvBnzYTyPBK0CG42P90YF1ZyJAJHSP26x2fNzc7rVuw6loXho2r8gNeOvVnKXDzleedr1uaLmJ\n9//zgveJiO7d2CPh56WWiBrOUHRx/sX3V5vXubjLOcw41ybjtKXDbY7kv67Z0Wa1uz4r+6kUr3kX\nSQKK8tMwnwvFYgafnSEX5zyJL8NvwTp3j/fsw3o+XWvOnuZ1y/m5/DRsXmligMQQ9kTIxTlPYtRO\nde2m7uS5uCuYhW2cII0P7ImQl7LCdBTNC7xobTQL20yZOm4gFMM4sRoBInWizvkAXKAeRKBzIvW6\nRMXqdQGcWKUQC+YWK2/DEoczRCSEIUJEQhgiRCSEIUJEQhgiRCSEIUJEQhgiRCSEIUJEQhgiRCSE\nIUJEQhgiRCSEIUJEQlR7AC+YEpnt7e2orq5Gb28vFixYgL1790KSQlGLnojGi6pP8Y5UIvPAgQOo\nrq5GUVERduzYgcbGRhQXF6vZJCIKsbANZxRFQUtLC4qKigAAa9asQUNDQ7iaQ0RjpGqIBCqR2dXV\nhdTUVNdro9EIWZbVbA4RqUC1MpqjKZEZDJbRjD68rvigWhnNkUpkarVadHd3u17bbDbo9Xq/n8cy\nmtGH1xV9xrI9oirDmf7+fvT1Pap85iyRmZ2d7XaOJEkwm81oamoCABw/fpyTqkRRSJUQ+fbbb/Hq\nq6/ixRdfxKuvvupWIrOmpgatra0AgJ///Of44x//iGXLliElJQVLly5VozlEpCLu9h4BYrV7zOuK\nPhEznCGi+MEQISIhDBEiEsIQISIhDBEiEsIQISIhDBEiEsIQISIhDBEiEsIQISIhDBEiEsIQISIh\nDBEiEsIQISIhDBEiEsIQISIhDBEiEsIQISIhqlTA6+jowE9+8hPXa6vVioMHD2LZsmVu51VWVuL2\n7duYOHEiAODEiRNqNIeIVKRKiKSlpbkC4d69eygpKcHChQt9nnv48GGvneCJKHqoPpw5deoUCgoK\nMHUqi/0QxSLVQ6Surg4rVqzwe3znzp146aWX8Oc//1ntphCRCsZcMmKkMpoAcP/+fZSUlKC+vh5T\npkzxOleWZRgMBvT09KCqqgq7du1CYWGhz891OBwYGoqK6hajlpCgwdCQI9zNCDleV/SZMCFh1F+j\nWhlN4PuhjK8AAb4vtZmSkoKysjKcP3/eb4iwjGb04XVFn4irO1NXV4eVK1f6PPbw4UPY7XYAwODg\nIE6fPo2ZM2eq2RwiUoEqd2eAR0OZ5uZm/Pa3v3V7v6amBhs3bkR2djZef/11PHjwAIqiYMWKFViy\nZIlazSEilbCMZgSI1e4xryv6RNxwhohiH0OEiIQwRIhICEOEiIQwRIhICEOEiIQwRIhICEOEiIQw\nRIhICEOEiIQwRIhICEOEiIQwRIhICEOEiIQwRIhICEOEiIQwRIhICEOEiIQwRIhIiHCI7N+/HwsX\nLsQrr7zi9n57ezvWrl2L0tJS1NbWwtdWrna7HZWVlVi+fDm2bduGgYEB0eYQ0TgTDpEXXngB7733\nntf7Bw4cQHV1NT799FN0d3ejsbHR65wjR47ghRdewCeffILp06fj2LFjos0honEmHCLPPPMMUlNT\n3d5TFAUtLS0oKioCAKxZswYNDQ1eX9vQ0IDVq1cHPIeIIpsqdWe6urrcgsVoNEKWZa/z7t69i6Sk\npIDnOE2YkDCm7eyjRaxeG68r9o0YIsHU3CWi+DViiARTc9eTVqtFd3e367XNZoNer/c6b+rUqejr\n60NSUpLfc4gosqlyi1eSJJjNZjQ1NQEAjh8/juLiYq/zli5dir/+9a8BzyGiyCZcRrO2thafffYZ\nuru7odPp8Otf/xrPP/88rFYrdu3ahTt37mDBggXYu3cvNBoN3nnnHcydOxfPP/887HY7duzYAVmW\nkZOTgz/84Q+YPHlyqK6NiMZB1NTiJaLIFBUrVkUWtEWTkpISvPjiiygvL8eWLVvC3RwhDQ0NKCsr\nw/Lly2Nq/U9eXh7Ky8tRXl6OmpqacDdHyI4dO1BQUICdO3e63rNYLFi1ahVKS0tx+PDh4D5IiQJ/\n//vfldbWVmX9+vVu72/fvl05deqU698/++yzcDQvZIqLi5X79++HuxnCHjx4oJSVlSmyLCt9fX1K\nWVmZYrfbw92skFi4cGG4mxAyZ8+eVerr65Xq6mrXey+//LJy6dIl5eHDh8r69euVCxcujPg5UdET\nEVnQRuPPYrFg1qxZ0Ov1SExMxNKlS/H555+Hu1nkobCwEImJia7XsixDURTk5OQgISEBq1ev9rnS\n3FNUhIgvwS5oizYVFRVYt24d6urqwt2UMevs7ITBYHC9jpWfDQD09PRg7dq1qKiowLlz58LdnJAa\n689NlRWroxUvC9pGus6jR4/CYDBAlmVs3rwZJpMJGRkZ49xKCqS+vh4GgwFff/01qqqqcOLECSQn\nx/fq1YgIETUXtEWSka7T+VfAYDBg0aJFaGtri8oQ0ev1bn/BbDYb8vLywtii0HH+jGbOnIlZs2bB\narXCbDaHuVWh4evnFszvVNQOZ4Jd0BYt+vv70dfXBwDo7e1Fc3MzsrOzw9yqscnPz8fFixfR2dmJ\nu3fvoqGhAYsXLw53s4T19PRgcHAQwKP5g0uXLmH69OlhblXoOAPy8uXLGBoawsmTJ4P6nYqKdSKj\nXdAWjb755hts3boVwKNJ402bNmHjxo1hbtXY1dfX4/e//z0cDgdef/11bNiwIdxNEvaPf/wDtbW1\n0Gg00Gg02LZtG5YtWxbuZo1ZVVUVLBYL7t27h5SUFLz77rsYHBxETU0NBgYGUF5eju3bt4/4OVER\nIkQUuaLzzzYRRQyGCBEJYYgQkRCGCBEJYYgQkRCGCBEJYYgQkRCGCBEJ+f/c/Gkku9ExvgAAAABJ\nRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc1a93d2588>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"b = 0.5\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD7CAYAAABNPKDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG7xJREFUeJzt3W1QVPe9B/Dv2fUJgcCSZHe5uQryIBAF27TK+BAiRIJm\nNBijicThktur3jY+FG2nM5VKR8e2SWtqM+OLNL7KzWSc1CZXW+eGJCUgxptImtzLkoqID4jXuLsq\ny8MKgu6e+0J3w7Jnl4XDYc/Z/X5eyZ7D8j8CX87/6fwEURRFEBGNkS7cDSAibWOIEJEsDBEikoUh\nQkSyMESISBaGCBHJIjtEtm/fjvnz52PHjh3e1ywWC1auXIni4mIcPHhQ8vM6OjqwZs0aFBcXo7q6\nGpxpJtIm2SGyYcMGvPrqqz6v7d27FwcOHEBNTQ0aGhrQ2trq93n79+9HZWUlPv74Y3R1daG+vl5u\nU4goDGSHSH5+PmJjY70f22w2iKKIzMxM6PV6rFq1yi8gRFFEU1MTCgoKAACrV69GXV2d3KYQURiM\n+5iI3W6HyWTyfmw2m2Gz2XzOcTgcSExMDHoOEWnDpHA3IFR377oQqcMmer0AlyvyLo7XpT2TJ+tH\n/TnjHiJGo9HnrsJqtcJoNPqcYzAY0NXVFfSc4UQR6OrqG9/GqkRi4vSIvDZel/Y8/HD8qD9n3Lsz\nnq5MW1sbXC4Xjh8/jsLCQp9zBEFAbm4uGhoaAABHjx71O4eItEGQu4t38+bNsFgs6O/vR0JCAt54\n4w0MDg6iqqoKAwMDKC0txbZt2wAAVVVVWL9+PXJzc9He3o6dO3eip6cHCxcuxJ49e6DTBc60O3dc\nEZv+kfqXjdelPWO5E5EdIhOFIaI9vC7tUUV3hoiiC0OEiGRhiBCRLAwRIpKFIUJEsjBEiEgWhggR\nycIQISJZGCJEJAtDhIhkYYgQkSwMESKShSFCRLIwRIhIFoYIEcnCECEiWRgiRCQLQ4SIZFGsZERT\nUxOqq6u9H7e1teG9995DTk6O97WioiLExcVBEAQYjUYcOnRIqeYQkUIUC5F58+bh2LFjAICrV6+i\nvLzcJ0A8jhw5gqlTpyrVDCJS2IR0Z2pqalBSUjIRX4qIJtiEhciKFSskj5WVlWHt2rWoqamZiKYQ\n0ThTvIzm1atX0dnZiby8PL9jhw8fhslkgs1mQ0VFBXJycpCSkiL5Pnq9gMTE6Uo3Nyz0el1EXhuv\nKzooHiIffvhhwK6Mp1qeyWTC4sWL0dLSEjBEXC4xYmt9RGodE16X9qiy7kygrkxfXx+cTicAoLe3\nF42NjUhPT1e6OUQ0zhS9E/nmm2/Q2dmJ3Nxc72ubNm3Cvn37MDg4iC1btgAARFHEhg0bkJmZqWRz\niEgBLKOpApF6e8zr0h5VdmeIKLIxRIhIFoYIEcnCECEiWRgiRCQLQ4SIZGGIEJEsDBEikoUhQkSy\nMESISBaGCBHJwhAhIlkYIkQkC0OEiGRhiBCRLAwRIpKFIUJEsjBEiEgWhggRyaLog5rnzJmDjIwM\nAMDcuXPxq1/9yue4xWLBrl27MDAwgNLSUmzdulXJ5hCRAhQNkcTERG89Xil79+7FgQMHkJaWhrKy\nMhQXFyMrK0vJJhHROAtbd8Zms0EURWRmZkKv12PVqlWor68PV3OIaIwUvRPp7u7GmjVrMHXqVFRW\nViI/P997zG63eyvgAYDZbMZnn30W8L1YRlN7eF3RQdEQqa2thclkwvnz57F582YcO3YM8fGjr2sB\nsIymFvG6tEd1dWc8dxoZGRmYPXs22tvbvceMRiNsNpv3Y6vVCqPRqGRziEgBioVId3c3BgcHAdwb\n/zh37hxmzJjhPe4JmLa2NrhcLhw/fhyFhYVKNYeIFKJYd+bChQuorq6GTqeDTqfDrl27kJiY6K3F\nazKZsHv3blRWVnqneDkzQ6Q9rMWrApHax+Z1aY/qxkSIKPIxRIhIFoYIEcnCECEiWRRdbEYUClEU\ncfayA+22XqSa4pGdYoAgCCMeI3VgiFBYiaKIt2rO4lPLNbhFQCcAS/KS8dKKnKDHSD3YnaGwarns\n8IYEALhF4FPLNbS0dwY9RurBEKGwumzr9YaEh1sE2m29QY+RerA7Q2GVaoqHToBPWOiEe697/h3o\nGKkD70QorLJTDFiSlwzd/bFSz7hHTmpS0GOkHlz2rgKRuox6NNfV0t7pnYEZHhLBjoVDpH6/gLEt\ne2d3hlQhJzUpYEAEO0bhx+4MEcnCOxFSjCiKPl0ROQvFuOhMvRgipAhRFPHG+82o/XtHyAvFPEFx\nydoDHQS4IWKW+QFkzUzEf3zYykVnKsUQIUW0XHbgky87/BaKmQ3TId6fpvXcTYiiiDOXbuKdv7XB\n2tnv8z46AZg7KwlfX+r0e6/8HBOyUwyy7lB4hyMfQ4QUcdnWC5fb9zW3CPyp/gKAb+8mKpZn460P\nzqLBck3yfdwiYLnov0LVLQKXrD043WIb8x0Kl9WPDw6s0rjwjH98cPoyWto7kWKMgz7IT5fnbuLD\n0x042SwdICO53tkva1k8l9WPD8XuRC5evIhdu3bB6XRi8uTJ+PnPf44FCxb4nDNSmU3SBqm/6Itz\nzSj63kzvmIgUtwh8cdaOsa5UcjgHJJfF1zR24LK11zumEqiLEmxZPaeUQ6dYiEydOhW//vWvkZaW\nhgsXLuBHP/oRPvroI59zRiqzSdog9Rf9VLMV1f+Wj++kJ6Hd1gvBLeJIw0W/wGi3jm0fjE4Asmca\nfMZKPJovdqL5fhcoWBclxRgn+d5cVj86inVnHnnkEaSlpQEA0tLS4HQ6oZHFsTRKgf6iX7jajZzU\nJCxfMBPXHH2SdxyBfiIS4yYH/HrC/WAoyZ/psyxeSrAuCn8ax8eEDKzW1tbi0Ucf9bulDFZmcziW\n0VQXURTRfOEGLlztRsy0ydDr4DOQqtcBs2cYkJg4HU3nr+O/v7aG9L4xU/VYnJcMQRDwceOVAF8b\nmDplMgyGWFSWfQ9F86/j2MmL+J/W65Lnu0XA2n0bC4f9H9t7vpE8X+rcobT4/VKS4iFy9epV/O53\nv8Obb77pd2w0ZTZZRlM9ho+BCABMSTGwO/qHjIkkY07ag+jq6sOZCzf8ZmoC6R9w4W9f/N+I5/3t\niw58J/3ecviZD8Vi2XcfQdO56wHHX8wJ0/z+j80J0yR3CUudO5TWvl+jobqSEU6nEy+//DJ2796N\nlJQUv+PBymyS+nhmYN6qOYuGpm/HQEQAdkc/nitIQ0GeGXNnJUEA0HTeDlEUA449xMWM/W+YCOD0\nGau3i+zZ8Tsa3CU8PhS7E3G5XPjxj3+MF154AUuWLPE73t3djZiYGEyZMkWyzCapy/C7j+HcIvBF\nqx3tVqf3tRNN15A7y4De23cl39PZL/16qBosVkAQ8NKKHAiCgIrl2Th3pctvwRogPeMi3P/c/ByT\nqnYJa41iIdLQ0IDPP/8cN27cwLvvvgsAePvtt1FeXo5jx44FLLNJ6jR8BkbK0ADxaL7kULBVQEPT\nt6tgBRGwOfwDBAg+48JdwvLweSIqoIU+9genL+NI3YWwfO1UUxzabf4BFapZ5jjsfmnByCeGSAvf\nr7FS3ZgIaZ9nHMTeGZ5fmrmpiQg6hxuCtUszxqk1JIV7ZyigkcZBJsLVG33ocg6G54tTSHgnQpJE\nUcSHpzt8ZmGAewu90szSsy1KcDgHZS8K49PhlcUQIT+eOxDPjlvfY8A/G+Ogld3yfDq88hgi5Kfl\nsgMnm6R31uoEoLN3YMyb5iYa130ojyFCftqtPQG7EEZDDL5WeNp2vJgN0/hskAnAECE/QpC7DKmF\nXGqVNcPgfb6JRlYyaBJnZ8iHKIr4otUe7mbINmWSDieb+cSyicA7EfLRctkhufJUawbvuvnEsgnC\nECEflyN0OpSFwJXDEIkSw5+BGmiMINCOW60QBCA3LclvkSunepXDMZEoEOyp5p5waWyxwdE7gJ5b\n2l4dKorAFbsTj6Yk4h/tXRDBLf5KY4hEqKH1VIT7YwLDxwgWZBvR2GIPWK5Bq7qcg+hyDmJWcjzi\nYyYje6YBJfkzw92siMUQiUCh7Hlxi0DjWfuYyzVowaVr98ZAvr7UCaujj7MzCuGYSAQK5dkfANCl\noZWncnB2RlkMkQgk9fR1KVKV5SIVZ2eUwxCJQKmmeLmP4Ig4nJ1RDkMkAg1/AHG0EvgA5gmh6MBq\nXV0dXnnlFYiiiE2bNmHdunU+xy0WC3bt2oWBgQGUlpZi69atSjYnagx9APHpFhsaAuzIjXQFeckw\nJk3nA5gVpliI3L17F6+++irefvttxMbG4rnnnsOyZctgMBi85+zduxcHDhxAWloaysrKUFxcjKys\nLKWaFHVyUpOidhxAJwALckwMjwmgWHfGYrFg9uzZMBqNiI2NxdKlS3Hq1CnvcZvNBlEUkZmZCb1e\nj1WrVqG+vl6p5kStaBkfMcRPYf2YMFHsTsRut3uLUwGA2WyGzWYLevyzzz4L+H4soxm6oSUu0/7p\nARR9fwZq/34lIqdzH8t6GKsen4V5GUZYzl/HhavdSH8kAXkZDyv2NVlG05dmFpuxjGZoAi1xfzzX\nfK/YU4T5qvU64qZNQspDcZj5UCxmPhQLAIr+rLBkhC/FujNGo9HnzsNqtcJoNIZ8nMZm+EIztwic\nbLqGK/Zb4W2YghqaruHMpZvhbkbUUixE8vLy0NraCrvdjlu3bqGurs6nnKanK9PW1gaXy4Xjx4+j\nsLBQqeZEDamFZiKAS9bIHmB9r+FiuJsQtRQLkUmTJuFnP/sZysvLsXr1avzgBz+AwWDApk2bvHcg\nu3fvRmVlJUpKSrBkyRLOzIyDaBlIHa7d2stl7WHCMpoqoPSYyKOpBs08XDmQuJhJIxYAX1eYjhX5\nKYq3hWMivjQzsEqhGV7pPsUYh9NnbJLnhvKLqQYF85KxINuI/e82BTxH4LL2sOGy9wiVk5p076+y\nIOBks/SsjBYCxJQUg5dW5CAnNQkF85IDnpdq5qrUcGGIRDitPzPV7uhHS3un9w7r+aXpkuetfUL6\ndVIeQySCiaIYtIaMFojDtvCX5M9Ewbxkn9WpBfO4OjWcOCYSoTwDrIHKYWrJ0LGO4WM+3FwXfgyR\nCOVZdKbxGxHkpSVJhkROqvTrNPHYnYlQoT7dTM2MhhhUPv+dcDeDRsAQiVBarx+TFD8Fr/z7wnA3\ng0LA7kyE0upNyCS9gMVzzPiXFdnhbgqFiCESobQ6tXvXJeJk8zWIAnxKPAyto5Nqikd2igGCEIXr\n+1WIIRKB3G43vjx7PdzNGDNPiYf8+08mC1bBj8KPYyIRRhRFvP5ni6p37ZoNMSOeM7TEg9TjDVhH\nRj0YIhGm5bIDzSqvJ2N19Id0nmehnNRME+vIqAdDJMJodSxEinh/yEPq8QasI6MeDJEIk2qKRyQM\nNw4NieF1dPggZnXhwGqEyU4x4PF5yZqvNTNnVhKyU+6VF+FSd3XjQ4lUQImH3Lz1Xy04YVF/kJiT\nYlD+VBYuWXvwZet1XLp2rzum5hkYPpTIF7szEerhpJFnQNTAfn+QNdX8AC4PmVHiDIx2KNKdcTgc\nqKysxPXr16HX6/Hyyy9jxYoVfucVFRUhLi4OgiDAaDTi0KFDSjQnKs0yPwCdANXvnxk6yxJoBoZd\nF3VTJER0Oh1+8pOfIC8vDzdv3sSzzz6LwsJCTJs2ze/cI0eOYOrUqUo0I6p5BiOHL9AyG6bjT/UX\nwt08r6EDqMNDjzMw2qBIiCQkJCAvLw8A8OCDDyIxMRHd3d2SIULKCDQYKYoirJ19aBjH8ZJZ5njJ\ntRyCgKBV94bOsoiiKBl6vAtRP8VnZ86cOQO32+1TMnOosrIy6HQ6bNy4EcuXL1e6OVFn+HM3BEHA\ngkdN4xYiqeZ4/KLi+5LL0ovmz8Q/Ltzw3k14wmzovz1t4wyMdo15dmblypWSr7///vuYMmUKAKCn\npwcbNmzAnj178Nhjj/mda7PZYDKZYLPZUFFRgT/+8Y9ISZF+5L/b7YbLpfIO/hjp9Tq4XO4J+3r/\neeI83v7grN/ry+bPAACcPmNF7607Qd9DAPDdrIfxi3/N9742vBbuRF/XRInU6wKAyZP1o/6cMd+J\nHD9+POjxO3fuYNu2baioqJAMEODbKngmkwmLFy9GS0tLwBBhLd7xY06YJjnomhQ7BSX5M/Hd9Aex\n/93/9euKPJGXjPk5Rly2O5FivDcgfvjDFu+u2qG1cB2OW7hys897JxJJu245xetLse7Mnj17MHfu\nXKxdu1byeF9fH9xuN+Li4tDb24vGxkasX79eqeYQvt1Of8nag7mzktB8sdPnuSN/PnEBVkcfKpZn\n4/G8ZJxs8n28omeLfsXy7KC7aj27bk81X4PLHXzNB7f4a58iIXLu3DkcOXIEWVlZ+PTTTwEAr732\nGjIyMrBp0ybs27cPg4OD2LJlC4B7P0gbNmxAZmamEs0hSFfGSzHHo11ibUZ+jgkvrcjxm8nxHDcb\npkvuqvVs3Q+069ZzPFib1LrAjAJTJERmz56N1tZWyWND14L85S9/UeLLkwSpX+x2iccFDF2bIUrc\nELhFoKXDEXRNR7Bdt0NDJNSwIXXjitUoEeqDm4euzQhUHLy3707QXbWh7rrlFv/IwBCJEoECYShh\n2NqM7BQD5s7yvyO4bO3F3FlJAXfVeha66XXSx4O1iQvMtIe7eKPE8BWsUgryklExZDxCEARkpRhg\nGfaQIxFAVooBJQtmSq7p8Kz5GLpORKp7EmhVLbsy2sIQiRJDF3M1ttjQYLnmM4WrE4AFOf4LAj13\nC1LL0UcqIJWX8bB3ynekNnGBmXYxRKJMTuq953SIQEh3ABNxt8BqdtrG54moQLgWL7W0d4Z8BzCa\ncz0idVFWpF4XoLLFZqR+o7kD4N0CBcLZGSKShSFCRLIwRIhIFoYIEcnCECEiWTg7Q5Lb8QFwiz6F\nhCES5aS24y/ONUMQBG7Rp5AwRKKc1Hb8kxarz0OWuUWfguGYSJQL9IiA4euYuUWfAmGIRLlAjwgY\nPvzBLfoUCEMkynk22A19NsjjeWY8Puw1btGnQBQbEwmlRGZHRwcqKyvR29uLhQsXYs+ePZwBmGDB\ntuNziz6FQtGB1ZFKZO7fvx+VlZUoKCjA9u3bUV9fj8LCQiWbRAFIbbDjpjsKRdi6M6IooqmpCQUF\nBQCA1atXo66uLlzNIaIxUjREysrKsHbtWtTU1PgdczgcSExM9H5sNpths9mUbA4RKWDM3ZmRymge\nPnzYp0RmTk5OwOp2odDrBSQmTh/z56uZXq+LyGvjdUUHxcpojlQi02AwoKury/ux1WqF0WgM+H4s\no6k9vC7tGcuTzRTpzvT19cHpdAKAt0Rmenq6zzmCICA3NxcNDQ0AgKNHj3JQlUiDFAmRmzdv4sUX\nX8QzzzyDF1980adEZlVVFZqbmwEAP/3pT/GHP/wBy5YtQ0JCApYuXapEc4hIQXxQswpE6u0xr0t7\nVNOdIaLowRAhIlkYIkQkC0OEiGRhiBCRLAwRIpKFIUJEsjBEiEgWhggRycIQISJZGCJEJAtDhIhk\nYYgQkSwMESKShSFCRLIwRIhIFoYIEcnCECEiWRSpgHft2jX88Ic/9H7c3t6O1157DcuWLfM5r7y8\nHDdu3MCUKVMAAMeOHVOiOUSkIEVCJDk52RsI/f39KCoqwqJFiyTPPXjwoN+T4IlIOxTvzpw4cQLz\n58/H9Oks9kMUiRQPkZqaGixfvjzg8R07duDZZ5/FO++8o3RTiEgBYy4ZMVIZTQC4ffs2ioqKUFtb\ni5iYGL9zbTYbTCYTuru7sXnzZuzcuRP5+fmS7+t2u+FyaaK6xajp9Tq4XO5wN2Pc8bq0Z/Jk/ag/\nR7EymsC3XRmpAAG+LbWZkJCAkpISfP311wFDhGU0tYfXpT2qqztTU1ODFStWSB67e/cuOjs7AQCD\ng4M4efIkMjIylGwOESlAkdkZ4F5XprGxEb/5zW98Xq+qqsL69euRnp6OjRs34s6dOxBFEcuXL8cT\nTzyhVHOISCEso6kCkXp7zOvSHtV1Z4go8jFEiEgWhggRycIQISJZGCJEJAtDhIhkYYgQkSwMESKS\nhSFCRLIwRIhIFoYIEcnCECEiWRgiRCQLQ4SIZGGIEJEsDBEikoUhQkSyMESISBaGCBHJIjtE9u3b\nh0WLFuH555/3eb2jowNr1qxBcXExqqurIfUo187OTpSXl+Opp57C1q1bMTAwILc5RDTBZIfI008/\njTfffNPv9f3796OyshIff/wxurq6UF9f73fOoUOH8PTTT+Ojjz7CjBkzcOTIEbnNIaIJJjtEHnvs\nMSQmJvq8JooimpqaUFBQAABYvXo16urq/D63rq4Oq1atCnoOEambInVnHA6HT7CYzWbYbDa/827d\nuoW4uLig53hMnqwf0+PstSJSr43XFflGDJFQau4SUfQaMURCqbk7nMFgQFdXl/djq9UKo9Hod970\n6dPhdDoRFxcX8BwiUjdFpngFQUBubi4aGhoAAEePHkVhYaHfeUuXLsVf//rXoOcQkbrJLqNZXV2N\nTz75BF1dXUhKSsIvf/lLPPnkk2hvb8fOnTvR09ODhQsXYs+ePdDpdHj99dcxd+5cPPnkk+js7MT2\n7dths9mQmZmJ3//+95g2bdp4XRsRTQDN1OIlInXSxIpVOQvatKSoqAjPPPMMSktLsWnTpnA3R5a6\nujqUlJTgqaeeiqj1P3PmzEFpaSlKS0tRVVUV7ubIsn37dsyfPx87duzwvmaxWLBy5UoUFxfj4MGD\nob2RqAFffvml2NzcLK5bt87n9W3btoknTpzw/vuTTz4JR/PGTWFhoXj79u1wN0O2O3fuiCUlJaLN\nZhOdTqdYUlIidnZ2hrtZ42LRokXhbsK4+fzzz8Xa2lqxsrLS+9pzzz0nnjt3Trx79664bt068ezZ\nsyO+jybuROQsaKOJZ7FYMHv2bBiNRsTGxmLp0qU4depUuJtFw+Tn5yM2Ntb7sc1mgyiKyMzMhF6v\nx6pVqyRXmg+niRCREuqCNq0pKyvD2rVrUVNTE+6mjJndbofJZPJ+HCnfGwDo7u7GmjVrUFZWhtOn\nT4e7OeNqrN83RVasjla0LGgb6ToPHz4Mk8kEm82GiooK5OTkICUlZYJbScHU1tbCZDLh/Pnz2Lx5\nM44dO4b4+OhevaqKEFFyQZuajHSdnr8CJpMJixcvRktLiyZDxGg0+vwFs1qtmDNnThhbNH4836OM\njAzMnj0b7e3tyM3NDXOrxofU9y2U3ynNdmdCXdCmFX19fXA6nQCA3t5eNDY2Ij09PcytGpu8vDy0\ntrbCbrfj1q1bqKurw5IlS8LdLNm6u7sxODgI4N74wblz5zBjxowwt2r8eAKyra0NLpcLx48fD+l3\nShPrREa7oE2Lrly5gi1btgC4N2i8YcMGrF+/PsytGrva2lr89re/hdvtxsaNG/HCCy+Eu0myffXV\nV6iuroZOp4NOp8PWrVuxbNmycDdrzDZv3gyLxYL+/n4kJCTgjTfewODgIKqqqjAwMIDS0lJs27Zt\nxPfRRIgQkXpp8882EakGQ4SIZGGIEJEsDBEikoUhQkSyMESISBaGCBHJwhAhIln+H7H1psSETxCx\nAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc1a921fa58>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"b = 0.75\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD7CAYAAABNPKDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH3BJREFUeJzt3W1wU9e9LvBny7xjasskkmgCNn4BG7DJpAccAiEYYgwM\nxIRCAmG49HQSpg2BOqSnncINHRjaTnqhaWf4kElu555MpkNTTtIwZU4cOMbYkEnilDQWJMLmxcKE\nWHJAsrHAYLDW/eBIWNKWLGlrW3vLz+8T3tqW1sbjx3utvdb6S0IIASKiOBmS3QAi0jeGCBEpwhAh\nIkUYIkSkCEOEiBRhiBCRIopDZOvWrZg1axZeeukl/zGr1Yrly5ejvLwc+/fvl/2+1tZWrFq1CuXl\n5di5cyf4pJlInxSHyPr16/Hqq68GHNu9ezdee+01VFdXo76+Hk1NTSHft3fvXlRVVeHo0aPo6OjA\n8ePHlTaFiJJAcYiUlpZi7Nix/q+dTieEECgoKEBaWhpWrFgREhBCCDQ2NmL+/PkAgJUrV6K2tlZp\nU4goCRI+JtLe3g6z2ez/2mKxwOl0BpzjdruRmZkZ8Rwi0odhyW5AtO7e7UWqDpukpUno7U29i+N1\n6c/w4Wkxf0/CQ8RkMgXcVTgcDphMpoBzjEYjOjo6Ip4TTAigo+NmYhurEZmZY1Ly2nhd+nP//eNi\n/p6Ed2d8XZlz586ht7cXhw8fRllZWcA5kiShuLgY9fX1AID3338/5Bwi0gdJ6SreTZs2wWq1oru7\nGxkZGXj99dfR09ODHTt24Pbt26isrMSWLVsAADt27MDatWtRXFwMu92Obdu24fr165gzZw527doF\ngyF8pt2505uy6Z+qf9l4XfoTz52I4hAZLAwR/eF16Y8mujNENLQwRIhIEYYIESnCECEiRRgiRKQI\nQ4SIFGGIEJEiDBEiUoQhQkSKMESISBGGCBEpwhAhIkUYIkSkCEOEiBRhiBCRIgwRIlKEIUJEijBE\niEgR1UpGNDY2YufOnf6vz507h3fffRdFRUX+YwsXLkR6ejokSYLJZMKbb76pVnOISCWqhcjMmTNx\n6NAhAMCVK1ewYcOGgADxOXjwIEaOHKlWM4hIZYPSnamurkZFRcVgfBQRDbJBC5GlS5fKvrZu3Tqs\nXr0a1dXVg9EUIkow1ctoXrlyBS6XCyUlJSGvHThwAGazGU6nExs3bkRRURGys7Nl3yctTUJm5hi1\nm5sUaWmGlLw2XtfQoHqIfPjhh2G7Mr5qeWazGXPnzoXNZgsbIr29ImVrfaRqHRNel/5osu5MuK7M\nzZs34fF4AABdXV1oaGhAXl6e2s0hogRT9U7km2++gcvlQnFxsf/Y888/jz179qCnpwebN28GAAgh\nsH79ehQUFKjZHCJSActoakCq3h7zuvRHk90ZIkptDBEiUoQhQkSKMESISBGGCBEpwhAhIkUYIkSk\nCEOEiBRhiBCRIgwRIlKEIUJEijBEiEgRhggRKcIQISJFGCJEpAhDhIgUYYgQkSIMESJShCFCRIqo\nulHz9OnTkZ+fDwCYMWMGfvOb3wS8brVasX37dty+fRuVlZV48cUX1WwOpSAhBM5ecsPu7EKOeRwK\ns42QJClh59PAVA2RzMxMfz1eObt378Zrr72G3NxcrFu3DuXl5Zg6daqaTaIUIoTAW9VncdLaBq8A\nDBIwr2QCfrQ0tOZzPOdTdJLWnXE6nRBCoKCgAGlpaVixYgWOHz+erOaQDtkuuf2BAABeAZy0tsFm\ndyXkfIqOqncinZ2dWLVqFUaOHImqqiqUlpb6X2tvb/dXwAMAi8WCjz/+OOx7sYym/qh9Xc7Ob/yB\n4OMVgKPzFubIfG6s54eTqj+veKkaIjU1NTCbzTh//jw2bdqEQ4cOYdy42OtaACyjqUdqX5clYxQM\nEgKCwSD1HZf73FjPDydVf16ABuvO+O408vPzMWXKFNjtdv9rJpMJTqfT/7XD4YDJZFKzOZRiCrON\nmFcyAYbvxkV9YxxFOVkJOZ+io9qdSGdnJ0aPHo0RI0bA6XSiubkZEydO9L/uC5hz584hNzcXhw8f\nxu7du9VqDqUgSZLwo6VFKC0y+5+2yAVC/ycypUVmzC404VK7J+z5Sg21J0CqhciFCxewc+dOGAwG\nGAwGbN++HZmZmf5avGazGa+88gqqqqr8j3j5ZIbiUZSTFTYMhBD4zw9sOGF1+I89VmLBvy+bpkpb\nBvvztIC1eDUgVfvYWriuL1uuYd87jSHH/2PtQ3HfhUS6LjU+bzBpbkyEKNk+O9sue7zB5pQ9rrfP\n0wKGCKW0wR6JSN2Rj/AYIpTSZhXKP/GbXWSWPa63z9MChgiltKKcLMyfOQG+hyOSBMyfqd5j3cH+\nPC3gwKoGaGEAUg1aui6b3RXxMXAsormuRH7eYIpnYFXVGatEsfB6vTjScBlnW90onGTE4tkTYTAk\n5ma5MNsIALA7uwK+Vms+R6THzqmGIUKa4PV68b//76dwuLoBANaLLtRbv8FvN81R/N5yq3fnFlsg\nSRJX9CYAx0RIEz5suOwPEB+HqxvVn1xS/N5yq3dPWB04wRW9CcEQIU1oanXLHj8rc1wIAZvdhQ8+\nvQSb3YVww3q+8440tIas3u17PfBrr7jX3aHosTtDmlA4yQjrxdC7gMJJxoC1KNmmdHxqc+Kj046I\n3ZDgLkw0DBKQY45vlflQxhAhTSif9SCOnroMd1eP/5glazQqSicFhIEkBd5B+LohpUVm/0CmEAIf\nftqK+sa2qD+fK3rjxxChpBNC4O0jzej03AuQHHM61pTl463qswFhINdz8Yq+aeW+Jy7B3xOJhL4F\ncqXTLAyQODFEKOmCBz4BwO70YO87X8iGhpx6axsE+maGnrRGfwciAFxu9+BHyxgg8WKIUNJdcnZF\nNfAZifiuWyMBUY+B+LQ4PLDZXbwTiRNDhAZF8EY9Uydloqm1A3ZnFyTR161QOnXaFx7BWyBGo8Hm\nHDKbCCUaQ4RU43vE2uK4jubWDpxpcfmfqJiMo9Hu7vZ/bc6693W8DFLfAjgBxPRUBgDq+o2hpPom\nQonGECFVCCHw+nunUfPP0DkaXoGAiWVeAbS7u7H68TwICXBevYH60w7Eal7JBEybPB7TJo9HaZEZ\nHza0yj427i/4aQ/QNxHtEQ60Ro0hQoqE20/UdsmNY6fkJ3nJ8QqgV3iRBgM+P381rrZYjGMghIAk\nSf4A8N399Pf0gjxkW8bB7uzC2UtunJYJmgabkyESJdVC5OLFi9i+fTs8Hg+GDx+OX/3qV5g9e3bA\nOQOV2SRti1RR7pKzC73e2N6v9vMrcPWbJxKrvx2/AIf7pn/imW939/5dm8mWdGRb+sKuKCcL7S5t\nrDLWM9VCZOTIkfjtb3+L3NxcXLhwAT/96U9x5MiRgHMGKrNJ2hauolxpkRk55nFIMyDqIEkfPUxR\ngPictLZhdqEJkiT5d3f/tyn34S815+F0daPF4cHev36Bx2b2hd2sQlPAeIhPKm8ilGiqhcgDDzzg\n/3dubi48Ho//VpNSg9yjWd/6kyWzJ2HhDybh6GetUb2Xp/tuQtrkFcC79RdxydHXNglA5rgRATNh\nfQOvvlmu82dOwAlrG8R3M2If48zVmAzKmEhNTQ2mTZsWEiCRymwGYxlN7Zmedx/eq7sQcLeRZug7\nbjSOxeY1MwFJ4GjD5UFrkwTA7ujyD5YKICBAfLwC+NeFqxgzdgSyv5+BnO9nwCsE8h/MREn+/RE/\nQ68/L7WovrPZlStX8OMf/xhvvPEGsrOzA15zOp1Rl9nkzmbaE2lMRAiBA8fOo+afl2OaNKbUZEs6\nWhyemL8vlv1E9PrziobmdjbzeDx44YUX8Morr4QECCBfZrO4uFjNJlECRapA95Xdhf/5bPDuQHzS\nx4yIa7KZ3EI+io5qIdLb24uf/exneOaZZzBv3ryQ1wcqs0n64dsKsP/kslNnv01KW85cdKE4N0v2\n0e5AfOM5DJHYqBYi9fX1+OSTT3D16lW88847AIC3334bGzZswKFDh8KW2SR98nVtTjS2KZ6+rqgd\nAKZMykThxEz84+NL6O7pjfp7uZ9IfLjbuwbouY/tm2z2qc3pf8KRbJas0SFbLQ6EYyJ9NDcmQqkt\nnt3DBkOsAfL4zAmYzbGQuHGPVYqb3D4gemTKGsMAUYAhQnELtw+InnAcRDmGCIUVaVd1IQQknQSI\nMX142NemT87yb6tI8eGYCMmKNJHM6/XiT/9llV39qiXD0iQ8NW8ymi53wO2Rb+uXLS68VX2WRasU\nYIiQrHCL62YXmnDks8uaDxAAuNsrcLDuYsRzOMlMOXZnSFa4xXUNZ9sH3OhHb1i0ShmGCMnKMY+D\nIWjBtUFCbLsn6wQHV5VhiJAs34Y+viDxjYm4Pcr3/NASFq1SjmMiJEtucZ0AsO+vXyS7aX7po4bB\nlDkKE03puD9rDCZbvgchBN6tu+BfyWuQgLnFFnR4evzjOAbp3lOZ/osGKT4MEYrIt7gOAP77E3tS\n18UE89y6C4/Dg4sOD+bPnIBlj+QAAKZNHg+b3RWysljuGCnHEKEB+eaLxFJZbrCdCHrC0j/8fOSO\nkXIMkSGu/27t2aZ0AMCldo//33ZnF5paOzT/SFcIoMVxnSGRBAyRIUyrC+ji1dzagWWPJLsVQw+f\nzgwB4aavp8oCOp8zLS7Y7Nq+Y0pFvBNJYb7weLf+on/z4uDaMKkSIAB3JksWhkiK8nVV6oNqqgTX\nholnP1KtkiT4x3KChavUR8oxRFKUr6sip39tmOAKcZEMHybhzl3tJo74blr+tMnjg46HX0xIyqk6\nJlJbW4uKigosXrwYBw8eDHndarVi+fLlKC8vx/79+9VsypATqavi+4vtm1D28jMPYfWCXDxeMiHi\ne2o5QHxOWttCxkXCLSbk+EliqBYid+/exauvvoq3334bf//73/HnP/8Zbrc74Jzdu3fjtddeQ3V1\nNerr69HU1KRWc4YcubUvPr6/2D6F2Ua0u7tRp+F5INGSW0wXqVIfKadaiFitVkyZMgUmkwljx47F\nggUL8NFHH/lfdzqdEEKgoKAAaWlpWLFiBY4fP65Wc4ac4LUvwU5a2/Dlxauo/uQS9rz1T5yQqUer\ndcb04bLXJwkEbKAUbjEhF90lhmpjIu3t7f7iVABgsVjgdDojvv7xxx+HfT+W0Yxd1bofYOGsb3Ho\nxEX8qymwDoxXAG8ctqHr5p2Ef+5gmfxABsZ/bwyOnWoNKOX5bv0FuG/cwU9/WAIAeGTmaCy64PKf\nl2YAFv5gEuY89GBcn8symoF0M7Da2ytSdpt+NUsQTBw/BgUTvhcSIgB0HSAA8EXzVbz8zEMwjh2O\nvx2/4D/e6wVq/tmKh/LuTXNftygfD+VlBaydiff/nCUjAqnWnTGZTAF3Hg6HAyaTKerXSTkhBP7z\nA1vALxjQdytvTB+RpFYljlcADTYnvDLLAuXGPIpysrC0NJvzSBJMtRApKSlBU1MT2tvbcePGDdTW\n1gaU0/R1Zc6dO4fe3l4cPnwYZWVlajVnSPrK7sIJqyPk+IzJWSiZnBqbE9c1tqG5tYNjHkmkWndm\n2LBh+MUvfoENGzbA6/Xiueeeg9FoxPPPP489e/bAbDbjlVdeQVVVFW7fvo3KykpMnTpVreYMSZ/1\newLTX6ptb3imxYXpk7Pw5Xf1d7nR0OBiGU0NUKuP/dYHNtTp8KlLPNaU5SHHPG5Q9gvhmEgg3Qys\n0sCCl/Xfnzky2U1KOEkK3ebV13XhfiHJwRBJEam2rF9OSW4Wyv/tQXzW9G3IFHaGR/IwRFJEqi3r\nl3OmxYXMcSND9n5lgCQXQyRFpNKyfkkCHiuZAItxTMDj6eBCUwwPbeCmRCki0loZvckxj0NpkTnq\n+R+UXAyRFDHQWhk9aXF0Yd87X3D+h06wO6Mz4TZWzjGPw8YlhSgtMuPDhlbdzwXxCs7/0AuGiI5E\negITvNHOme9+8fRismUcWhyB3RSv6LvDWjJ7EgdRNYzdGR2J9ASm/0Y7hdlGzCuOvMGQlhjTh2P1\ngryQroskAe2uvkldS2ZPYoBoFENERwZ6AqPXQcfMcaNkx3SE6Fsbs++dL/BW9dnkNZAiYndGR6LZ\nWDnHPA5f2V2o19EuZZccXTh7ye2f//GpzRmwwbTvLstiHAPx3cAqN1rWDoaIjkydlAmTcTQcru6I\n531mc0Z8XWv6l3ooysmSvZvyCvjnjHCjZW1hd0ZHzrZ2oN0dOUBaHNfRdLljkFqUGMGPbQea88KN\nlrWFIaIjA42JGCTAAAlO963Ba1QCzC22BAyaRjPnRa/jP6mI3RkdiTQm4rvFFzocJgge2/CVsvCt\nj5EE8F91FwKum5POtIMhoiO+v9D9V7BOn5yFwmyjfw6Fze6Kuqrd04/n4uipr+H29Kjf+Aj6r4fp\nzzdGIoSAw32TK3c1ipsSaUCsm9zY7K6wk6+EEPh///0VTp4eeHD1sWIzJmSNxd/qLsbc5kRbU5aH\npaXZEc+JdN2DiZsSBeKdiA7I1ZGN9EvUeSO6XdxPRBE0iSY3MzXarglX7mqTKiHidrtRVVWFb7/9\nFmlpaXjhhRewdOnSkPMWLlyI9PS+co4mkwlvvvmmGs3RDbmwABAw1V0CMCP3XhcmeL6E7ZIbpzWy\nbqbgwQycv9IZsBOZ3dGFktws/7R8dk30T5UQMRgMePnll1FSUoJr167hqaeeQllZGUaNGhVy7sGD\nBzFyZOpt4xercEWnZxWaAiZeCQCnL7pw+qJLdr7EJY08sSjOzcJDU00493VnwHEBYMqkTBROMuJs\nqxuFk4yoKJ2UnEZSQqgSIhkZGSgp6as+Nn78eGRmZqKzs1M2RKiPXNHpE41taI4w5yN4kx4guU8s\nJAmY/13wSZKEf128FrInqgSgubXDfydypsUFh/smJ47pmOrzRL766it4vd6Akpn9rVu3DqtXr0Z1\ndbXaTdE0uTkgAhhwdmrwfInCbCOKk1BTxvDdbmT/a0khGs62Y987X+B/Gi4HBIhB6uuK9V9hzIlj\n+hf3ncjy5ctlj7/33nsYMaKvutr169fxy1/+Ert375Y998CBAzCbzXA6ndi4cSOKioqQnS0/Qp/q\ntXin592H9+ouBNSUlQCZvb2CvtcATM+7z/9/I4SAeXw6Tre4VWuvnNwHMlC17gdoPP8tPjoduNJY\nArBo1kTMm/l9XLjSGTJm4xWAo/MW5ujk58tavIHiDpHDhw9HfP3OnTvYsmULNm7ciIcfflj2HN/d\nidlsxty5c2Gz2cKGSKrX4p04fgzmFofOAflSZl8QXxfBIAFziydg0n1j/f83X9ldOHbq8qBfw/mv\nO/HxF1/D7uwKCEKgLwiN6SMw6b6xuOG5HTKPxSABloxRuvn58hFvINUe8e7atQszZszA6tWrZV+/\nefMmvF4v0tPT0dXVhYaGBqxdu1at5mhe8CxN35MXucHWSDudJ3PDZl+b5ELCN1YjN2GOT2f0TZUQ\naW5uxsGDBzF16lScPHkSALBv3z7k5+f7y2j29PRg8+bNAPpuwdevX4+CggI1mqMrwXMhwpVHCPdL\nF812AWrwBcVAISEXlgwQfeOMVQ1I5O1xrEWsTFmj0T7A4K2ccaOHoav7LgD5pfk2uwuOzluwZIxK\nuZBgdyYQZ6ymmOC/9O3XbqIuzAZFBgkompgZV4gseSQ7Yu3bopwszEnhXza6h1sBpKiinCwsLc3G\n/Vmjw54zr2QCZheZYy4z0b/27dLS7JS706DYMERS3GTL92RD4ukFefjR0qKQvTsMEpA7IfwtLQdC\nKRi7Myku3EDnkkf6HqVH+1QoeMsBIh8OrGrAYAzUxbOMXunS+1QdgEzV6wI4sEoRxLOMnkvvKRoM\nkRQXqewmyy5QIjBENCTcfiLBxwb6xfe9T4vjesCK2f5iLbsg1zYGEAEMEc2Q209kbrEFkiSFDIpG\n+sWPdrKZ3DYCsbSNdV/Ih494NUJ2PxGrAydCjrXhq5ZrUb9PJNGWXbBdcuNEY+heJ1y+TwBDRDPC\nLZwLfnYmBPBu3YWY30dOtHub2h3XQ7YkEGDdF+rDENGIcFXf5EYdWhwevPWBDTa7C/2f0AshIMUQ\nINFOGjPItgJRfxalNo6JaITcpLC5xRZ8/e0NtLSF/sWva2zDCWubf2yi/7hFf8ETxQDEPPfDG2Zr\nJD0WyqLEY4hoRLgl8ja7C//nr1/Ifk//wVEByI6FrH48zz871SfWuR++qfOsQEdy2J3RmOBFbYXZ\nRsyfGb4urW9wNOyYSgLuFuTW13D9DPnwTkTj+t+hNNicqLe2hWx+7LsjUOtugRsJUSQMEZ0oyukb\n1+jfbel/RyCEUH3bQU6DJzlcgKcBiazFq5V6tUDqLlRL1esCNLYAL5oSma2traiqqkJXVxfmzJmD\nXbt2cSp1FCLdEfBugQabqt2ZgUpk7t27F1VVVZg/fz62bt2K48ePo6ysTM0mEVGCJe3pjBACjY2N\nmD9/PgBg5cqVqK2tTVZziChOqoZIpBKZbrcbmZmZ/q8tFgucTqeazSEiFahWRjOWEpnRSPUymql4\nbbyuoUG1MpoDlcg0Go3o6LhX8d7hcMBkMoV9v1Qvo5mK18br0p94ns6o0p25efMmPB4PAPhLZObl\n5QWcI0kSiouLUV9fDwB4//33OahKpEOqhMi1a9fw7LPP4sknn8Szzz4bUCJzx44dOH36NADg5z//\nOf74xz/iiSeeQEZGBhYsWKBGc4hIRZxspgGpenvM69IfzXRniGjoYIgQkSIMESJShCFCRIowRIhI\nEYYIESnCECEiRRgiRKQIQ4SIFGGIEJEiDBEiUoQhQkSKMESISBGGCBEpwhAhIkUYIkSkCEOEiBRh\niBCRIqpUwGtra8NPfvIT/9d2ux379u3DE088EXDehg0bcPXqVYwYMQIAcOjQITWaQ0QqUiVEJkyY\n4A+E7u5uLFy4EI8++qjsufv37w/ZCZ6I9EP17kxdXR1mzZqFMWNY7IcoFakeItXV1ViyZEnY1196\n6SU89dRT+Mtf/qJ2U4hIBXGXjBiojCYA3Lp1CwsXLkRNTQ1Gjx4dcq7T6YTZbEZnZyc2bdqEbdu2\nobS0VPZ9vV4vent1Ud0iZmlpBvT2epPdjITjdenP8OFpMX+PamU0gXtdGbkAAe6V2szIyEBFRQXO\nnDkTNkRYRlN/eF36o7m6M9XV1Vi6dKnsa3fv3oXL5QIA9PT04MSJE8jPz1ezOUSkAlWezgB9XZmG\nhgb87ne/Czi+Y8cOrF27Fnl5eXjuuedw584dCCGwZMkSPP7442o1h4hUwjKaGpCqt8e8Lv3RXHeG\niFIfQ4SIFGGIEJEiDBEiUoQhQkSKMESISBGGCBEpwhAhIkUYIkSkCEOEiBRhiBCRIgwRIlKEIUJE\nijBEiEgRhggRKcIQISJFGCJEpAhDhIgUYYgQkSKKQ2TPnj149NFH8fTTTwccb21txapVq1BeXo6d\nO3dCbitXl8uFDRs2YPHixXjxxRdx+/Ztpc0hokGmOESWLVuGN954I+T43r17UVVVhaNHj6KjowPH\njx8POefNN9/EsmXLcOTIEUycOBEHDx5U2hwiGmSKQ+Thhx9GZmZmwDEhBBobGzF//nwAwMqVK1Fb\nWxvyvbW1tVixYkXEc4hI21SpO+N2uwOCxWKxwOl0hpx348YNpKenRzzHZ/jwtLi2s9eLVL02Xlfq\nGzBEoqm5S0RD14AhEk3N3WBGoxEdHR3+rx0OB0wmU8h5Y8aMgcfjQXp6ethziEjbVHnEK0kSiouL\nUV9fDwB4//33UVZWFnLeggUL8I9//CPiOUSkbYrLaO7cuRPHjh1DR0cHsrKy8Otf/xqLFi2C3W7H\ntm3bcP36dcyZMwe7du2CwWDAn/70J8yYMQOLFi2Cy+XC1q1b4XQ6UVBQgD/84Q8YNWpUoq6NiAaB\nbmrxEpE26WLGqpIJbXqycOFCPPnkk6isrMTzzz+f7OYoUltbi4qKCixevDil5v9Mnz4dlZWVqKys\nxI4dO5LdHEW2bt2KWbNm4aWXXvIfs1qtWL58OcrLy7F///7o3kjowKlTp8Tp06fFmjVrAo5v2bJF\n1NXV+f997NixZDQvYcrKysStW7eS3QzF7ty5IyoqKoTT6RQej0dUVFQIl8uV7GYlxKOPPprsJiTM\nJ598ImpqakRVVZX/2A9/+EPR3Nws7t69K9asWSPOnj074Pvo4k5EyYQ2GnxWqxVTpkyByWTC2LFj\nsWDBAnz00UfJbhYFKS0txdixY/1fO51OCCFQUFCAtLQ0rFixQnameTBdhIicaCe06c26deuwevVq\nVFdXJ7spcWtvb4fZbPZ/nSo/GwDo7OzEqlWrsG7dOnz66afJbk5CxftzU2XGaqyGyoS2ga7zwIED\nMJvNcDqd2LhxI4qKipCdnT3IraRIampqYDabcf78eWzatAmHDh3CuHFDe/aqJkJEzQltWjLQdfr+\nCpjNZsydOxc2m02XIWIymQL+gjkcDkyfPj2JLUoc388oPz8fU6ZMgd1uR3FxcZJblRhyP7dofqd0\n252JdkKbXty8eRMejwcA0NXVhYaGBuTl5SW5VfEpKSlBU1MT2tvbcePGDdTW1mLevHnJbpZinZ2d\n6OnpAdA3ftDc3IyJEycmuVWJ4wvIc+fOobe3F4cPH47qd0oX80RindCmR5cvX8bmzZsB9A0ar1+/\nHmvXrk1yq+JXU1OD3//+9/B6vXjuuefwzDPPJLtJin3++efYuXMnDAYDDAYDXnzxRTzxxBPJblbc\nNm3aBKvViu7ubmRkZOD1119HT08PduzYgdu3b6OyshJbtmwZ8H10ESJEpF36/LNNRJrBECEiRRgi\nRKQIQ4SIFGGIEJEiDBEiUoQhQkSKMESISJH/D+aCa+Km25EkAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc1a91bba90>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"b = 1.0\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD7CAYAAABNPKDeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wU3W+N/D3SUWElqUpmoSrtKWl0Agtji6t/BBppRQY\nsC7ij8pw2esAs/4AK7t3d1YG9oGH3Rl3YF2f4Znx0Wd2rtc74yhXLzzyrEVvKRS9Qv0xNmXpD340\njYskAZKWFgqV5Nw/yolNck6S9uQ0Oen79RfN+ebkeyj58P39EURRFEFENEyGRFeAiPSNQYSIVGEQ\nISJVGESISBUGESJShUGEiFRRHUQ2b96MOXPm4OWXXw68ZrPZsGLFClRUVGDv3r2y73M4HFi1ahUq\nKiqwfft2cKaZSJ9UB5E1a9bg1VdfDXpt586deO2111BbW4uGhga0tbWFvW/37t2oqanBp59+iq6u\nLhw5ckRtVYgoAVQHkdLSUqSnpwd+drlcEEURBQUFSEtLw8qVK8MChCiKaGpqwsKFCwEAjz32GOrr\n69VWhYgSIO5jIm63G2azOfCzxWKBy+UKKuP1epGZmRmxDBHpw22JrkCsbt70IVWHTdLSBPh8qfdw\nfC79GTMmbcjviXsQMZlMQa0Kp9MJk8kUVMZoNKKrqytimVCiCHR1XYtvZZNEZub4lHw2Ppf+3HXX\nhCG/J+7dGakrc/r0afh8Phw8eBBlZWVBZQRBQFFRERoaGgAA+/fvDytDRPogqN3Fu3HjRthsNvT1\n9WHixIl444030N/fj61bt+LGjRuoqqrCpk2bAABbt27F008/jaKiItjtdmzZsgVXrlzB3LlzsWPH\nDhgMyjHthx98KRv9U/V/Nj6X/gynJaI6iIwUBhH94XPpT1J0Z4hodGEQISJVGESISBUGESJShUGE\niFRhECEiVRhEiEgVBhEiUoVBhIhUYRAhIlUYRIhIFQYRIlKFQYSIVGEQISJVGESISBUGESJShUGE\niFRhECEiVTRLGdHU1ITt27cHfj59+jQ++OADWK3WwGvl5eXIyMiAIAgwmUx46623tKoOEWlEsyAy\ne/ZsHDhwAABw/vx5rF27NiiASPbt24exY8dqVQ0i0tiIdGdqa2tRWVk5Eh9FRCNsxILIsmXLZK9V\nV1dj9erVqK2tHYmqEFGcaZ5G8/z58/B4PCguLg679u6778JsNsPlcmHdunWwWq3IycmRvU9amoDM\nzPFaVzch0tIMKflsfK7RQfMgcujQIcWujJQtz2w2Y/78+WhpaVEMIj6fmLK5PlI1jwmfS3+SMu+M\nUlfm2rVr6O3tBQD09PSgsbER+fn5WleHiOJM05bI999/D4/Hg6KiosBrGzZswK5du9Df348XXngB\nACCKItasWYOCggItq0NEGmAazSSQqs1jPpf+JGV3hohSG4MIEanCIEJEqjCIEJEqDCJEpAqDCBGp\nwiBCRKowiBCRKgwiRKQKgwgRqcIgQkSqMIgQkSoMIkSkCoMIEanCIEJEqjCIEJEqDCJEpAqDCBGp\nwiBCRKpoelDzzJkzMW3aNADArFmz8Pvf/z7ous1mwyuvvIIbN26gqqoKL774opbVISINaBpEMjMz\nA/l45ezcuROvvfYa8vLyUF1djYqKCsyYMUPLKhFRnCWsO+NyuSCKIgoKCpCWloaVK1fiyJEjiaoO\nEQ2Tpi2R7u5urFq1CmPHjkVNTQ1KS0sD19xudyADHgBYLBZ88cUXivdiGk394XONDpoGkbq6OpjN\nZpw5cwYbN27EgQMHMGHC0PNaAEyjqUd8Lv1JurwzUktj2rRpmD59Oux2e+CayWSCy+UK/Ox0OmEy\nmbSsDhFpQLMg0t3djf7+fgAD4x/t7e2YMmVK4LoUYE6fPg2fz4eDBw+irKxMq+oQkUY0686cPXsW\n27dvh8FggMFgwCuvvILMzMxALl6z2Yxt27ahpqYmMMXLmRki/WEu3iSQqn1sPpf+JN2YCBGlPgYR\nIlKFQYSIVGEQISJVGESISBUGESJShUGEiFRhECEiVRhEiEgVBhEiUoVBhIhUYRAhIlUYRIhIFQYR\nIlKFQYSIVGEQISJVGESISBUGESJSRbMzVs+dO4dXXnkFvb29GDNmDH7729+ipKQkqEy0NJtElPw0\nCyJjx47FH/7wB+Tl5eHs2bN47rnn8MknnwSViZZmk4iSn2ZB5O677w78OS8vD729vRBFEYIgaPWR\nRJQAmmbAk9TV1eHee+8NCyCR0myGYhpN/eFzjQ6ap4w4f/48nn32Wbz55pvIyckJuuZyuWJOs8mU\nEfrD59KfpEsZ0dvbi+effx7btm0LCyBA5DSbRKQPmgURn8+Hl156CU899RQWLFgQdj1amk0i0gfN\nxkQaGhpw/PhxXLp0Ce+99x4A4J133sHatWtx4MABxTSbRKQvTKOZBFK1j83n0p+kGxMhotTHIEJE\nqozIOhGiSERRRGunF3ZXD3LNE1CYY4y4KHGo5UlbDCKUUKIo4u3aVnxmuwC/CBgEYEHxZPx8mTWo\njBQ0ckwZOH7Kic+aXYHrDxVb8E/L701E9QkMIpRgLZ3eQAABAL8IfGa7gFKrGdbcrLAgIwhA6FTA\nMZsTD95rgTU3a+QfgDgmQpGJoogWuwcfn+hEi92DeE/mdbp6AgFE4hcBu6sHQHiQUfr4xhaX/AXS\nHFsipCiWroZaOaYM2ddzzQNTjXJBhpILWyKkSKmr0WL3xO0zosWHXPMEGGIYMy2xmuNSHxo6BhFS\nZHdeke1qdDivxO0zHO5e+c++1Z0pzDFiQfHkQCAxCIAlaxykyRhBABbOnszxkARid4ZkiaKINkeX\n7LU2RxeWPxifz5FaGoODlUH4sTsjCAJ+vsyKUqs5MKVrzc1Ci90T9DMlDlsiJKul04uT5+S7LX/r\n8MStSyPX0phfZAGAwGCu3+8Peo8oirDmZmFZaQ4DSBJgS4Rkdbp6FMcrpNmTWL7ASgvDBr9eajWj\npNCETncvckwZaGx1Y8973wYGc03GcXB7+zQb3CV1GERIllw3QzK4uxGJ0uzOuqWFirM+p+yesMFc\np6cvcM/QdSSUeOzOkCypmxE6MSJ94WP5AtvOXgoLCA1NF/D6viY0NAW/fqxpYNYnlindwetIKPHY\nEiFZoQOaggiIt1ogsbYAzp3vlg0INpmxFhEIdHmUWkCSWFtCNDIYRCgia27WsLsN+XdPjBoQBhPE\nH1tAg7s6cmMi7MokDwYR0sysvEmYOTULzQqzPKFEgVO6esQgQpoQRRH/5z9O4m8dsU8FD+6ihLaA\n1LSISFuaDqzW19ejsrISS5Yswb59+8Ku22w2rFixAhUVFdi7d6+WVaER1tLpxeGvHWFdGR77kXo0\na4ncvHkTr776Kt555x2kp6fj8ccfx+LFi2E0GgNldu7ciddeew15eXmorq5GRUUFZsyYoVWVKAah\nZ3cAAwOeBgjwQ8RUy09iOgSo09UDnz/89YeKLOjq7ZcdXI117QklF82CiM1mw/Tp02EymQAAixYt\nwueff44VK1YAGEgTIYoiCgoKAAArV67EkSNHGEQSKHRdh5xYF3vlmicgzYCwQOLt7cfSkmyc7PAo\nLnUnfdEsiLjd7kByKgCwWCxwuVwRr3/xxReK92MaTe01nbmIz5uVAwgwMNPyefMFlM/JRvG0u8Ku\ni6KI5rOX4Oy6jry7J+L0d91B10+e88A8aTzum27Ct+1u+EUgzQCUP5CNuffdE+9H0kSy/L6ShW4G\nVn0+MWWP6U+WFASnzl6S7YKE8vmBv529hOw704NeDzuFTOa9IoD//PLvMAjArKlZmJFjDMy4JMPf\nQSyS5felhaRKGWEymYJaHk6nM9C1ieU6jbxYz+5Q6nqEnUIW4R5+ETjZ4eGUbQrQLIgUFxejra0N\nbrcbV69eRX19fVA6Takrc/r0afh8Phw8eBBlZWVaVYcGUTryMHRHrRKlxV5DPYXMLwKHGh2aHLtI\nI0ez7sxtt92GX//611i7di38fj/Wr18Po9GIDRs2YNeuXTCbzdi2bRtqampw48YNVFVVcVBVQ9Ks\nS4fzCtodXYGBzcEDpYMXep1ocaGh6ULYfW6/zYBSqxmiKAZmaKR7uz3XZA9SjsR2zoOTHZ6gwVqm\nhNAXptFMAlr3saPNuhgE4JdP3RfUuvj/X3Tgg6MdivdcOHvgSx/LjE4spDoU5hg1P9dVLY6JBOMu\n3lEgdKwilNyu2Eve6xHvKZ21qnTv4rws2R3ATy7KR3FeeFdIqsNInOtK8cUgMgpEG6uQGyj19N6I\neE/pS6907+nZmVg8JzvoxLIFxZOx9MEcVJZkh427SHWIlkKCko9upnhp+KIdMCQ3UGrNNkbcODc4\n8MidkTrV8hPMve8e3JefFbZxTm6n7uA6RDpzlZIPg8goIPelnTk1C4WD1miEWlIyBQ2274NOFZMM\n/tKLohgxIMhtnFPaqatUV279T24cWE0CIzFQJ4oiDp1woNXhRWG2EZWl2VFnPHw+H/7wr1/ju0tX\nkTFuDBbf/w+AwYCplp+Efamlrfo5pgwIggC7qwcz8+/ElEnjhzWzksxb/zmwGoxBJAmM9OxMLDMe\nPp8Pm18/hr7+4CWs0qxMLJ+TZgDmFyXXzEo8MIgEY3dmFJCb8WhougCLcTwq5tyDT7/8O1ocXmRl\n3A6/KMJgMMBz5XpYAAEG3jf4kOTBazqEWzMp0uf4/DxUeTRgEBkFlGZQ3j9yFh/9lx19/b4h3U/a\nsi+KIv7l4xYcszkVyw4lvQTpE6d4R4FIe2KGGkCk+wHAKbsnYgABOLMyGjCIpAil/TDArRmPoslx\n+Zxxt6cFPu/LVnfEsmkGzqyMBuzOpICYBk7jtPWkr9+HPe99K5uTZrDivCysKi8IOy6AUg9bIikg\n2lJx6Xq8SPc3Z45TLFOYbZQ9tIhSD4NICoi2VFxpYDXPkiF7P0vWOCwstkT8TL8I+AUgV+EeIjfd\njhoMIilAbuB08ICm0vWfFpplB1ydnj54eyLvnZGWtj+xaJrsNQ6mjh4MIikg9DCh0KXicocNzZqa\nhSUlUwbGNmQCSXOHV/HzBt+/MMeIh0JaLfOLLBxMHUU4sJoCIu1Fka6vW1oIb8+NwKa6kx0e/Ouh\ntoHDiAAclTmASE5xXhYqS7LD7i8dRiTcymJHowdbIinEmpuFZaU5KMwxhk33tnR6g7LRDR58LbHK\nd2tCY4FBQFgAkQZtpRllked/jDpsiehILMcGKk333pV5h+Lga+WcKUE5cwUAC4otEAQh6m7aSIO6\nc+P9F0BJSZMg4vV6UVNTg4sXLyItLQ3PP/88li1bFlauvLwcGRkDuz5NJhPeeustLaqTEmLdRKe0\nTyZjXPivWhAA1+Wr+F//bgtqpeROngDLpHTkmiegpNCETnev4m5aubNKOLA6umgSRAwGA375y1+i\nuLgYly9fxs9+9jOUlZXhjjvuCCu7b98+jB07VotqpBSltSChm9uUpnN7+26GvSaKQIPMsvWOCz3o\nuNAT025fnv9BmgSRiRMnori4GAAwadIkZGZmoru7WzaIUDi5botSt6HDeQUAAmVzTBmKp5gNlVKg\nGizaoC6lPs3HRE6dOgW/3x+UMnOw6upqGAwGrF+/HkuXLtW6OklPqdtSemvwc3BwEAC0O7rw4dFz\nP5YtmhzWMsg2T4DdObwzSmPdhSt3ghmNDsMOIlJi7lAffvghbr/9dgDAlStX8Jvf/AY7d+6ULfvu\nu+/CbDbD5XJh3bp1sFqtyMnJkS07WnLxhubDlXLfls2Zgkd+mo3DXzvg8w9sbps9zYSmM+7gsicv\nYNuzpSj76RR83nwBEIF5xRb83wOn8P2lq0OvmwGYmX/nsP7uUzVnbao+13ANO4gcPHgw4vUffvgB\nmzZtwrp163D//ffLlpFaJ2azGfPnz0dLS4tiEBktuXjl8uH6/MCps5dR/ci0oIOP7a4efNPuDit7\n8sxFuL19gdbI4a++w4LiyVgwy4KvWl045+yNuW7ziyYj+870Yf3dp+oJYKn6XECS5Z3ZsWMHZs2a\nhdWrV8tev3btGnp7B/4x9/T0oLGxEfn5+VpVRzeiLWGX1oJYc7MUyxogyA7C5lgmYOu6ObJ5X5SU\nWuW7oUQSTYJIe3s79u3bh2PHjqGqqgpVVVU4c+YMAGDDhg1wuVy4fPkynnnmGTz66KN45plnsGbN\nGhQUFGhRHV0JXaIuYOBkdrurR/6cEJnl7qLMwKo0tiEIAmbkGGOuD/O9UDQ8qDkJyDWPW+yeiHlz\nQ8sOnhlpsXuw571vw9ZurH44H6IACCLw70fPxjSD889P3zfsAdNUbfan6nMBPKg5pVhzsyACgZkX\nYKA1ccx2ASWFJtw7dVJQ2cIcI1o7vfj4RCey70qHyTguKGfM2DFpgcAhADBnjYPb2xcxkBTnccaF\nomMQSWJya0NEEfjg6NmgIBJLUu3BZ6mKANzevkDLxO25JrsBbyjdHhq9GESSkLTYzO25BgEDX/rB\nOpy9ONVxOZAkShCBY00XwspF4hcHDg5aWpKNQyccYde5dJ1ixSCSZGJpVQDAv33aHrU7Ek2OKQNv\n17biWEgrhEvXaSgYRJJM6B4ZJXI5cofqyxYXjtmcYS2Y1Q/nY+mD8ut1iEIxiCQBKd2D3dUDt+da\nXPa9xOKoQs4Yno9KQ8EgkmCiKOKND5tR95VjYOZExRf4yUX5aHV4YTs3/AOBOBZCQ8WTzRKspdOL\nw187Aq0PNat2ciwTUFmSHbaKVcDACe7RcCyEhoMtkQTrdPWE7ZUZLrurB0tLsmXP9ygpNGH3e01h\n7xEE4IlbU73cxk/DwSCSYLnmCUgzIOZAIghA0dQsNHd4glotUjdE7nyPwhwj3v64VfZ+DxVP5iAq\nqcLuTIIV5hhR/sCPXZBIYyIGYeBL/9ITs/FQhBQRQPBGvZZOLz5rDl9M9uSi/IinlhHFgi2RBBME\nAc89XhzY4q+0ejQ0VcNQThNTOjKRszAUDwwiSUI6GazF7sGxkHUicqkaBr8nGh6mTFpidybJRMtm\nlyz3JJLwKIAkoHQUQLwPPtbinpGk6pb5VH0ugEcBpBQtDj7mYcqkBQYRHYiU+S6WrHhEWmIQSXKR\nMt/FmhWPSEuaBZFYUmQ6HA7U1NSgp6cHc+fOxY4dO/i/aIhIme9EIKaseERa0rQlEi1F5u7du1FT\nU4OFCxdi8+bNOHLkCMrKyrSsku5ESpgt/VnuGoMIjZSETfGKooimpiYsXLgQAPDYY4+hvr4+UdVJ\nWpFSSERLL0E0EjQNItXV1Vi9ejVqa2vDrnm9XmRmZgZ+tlgscLlcWlZHlyKt8eD6D0oGmqXRHEqK\nzFiMljSacmqqH0D5nIs4e74b+XdPRPG0u2K6lmipmm4yVZ9ruDRLoxktRabRaERXV1fgZ6fTCZPJ\npHi/0ZJGU0n2nenIvjMdAMLKRrqWSKm6KCtVnwtIojSasaTIFAQBRUVFaGhoAADs37+fg6pEOqRJ\nEImUInPr1q1obm4GAPzqV7/Cn//8ZyxevBgTJ07EokWLtKjOiJLOS/34RGdQ2kul14n0jntn4khp\n8de6pYURF4WlavOYz6U/SdOdGa2UFoYdOuGQfb3FPvwDlYmSBYNIHCktDGt1eCMuGCPSMwaROFJa\n/FWYbeSiMEpZDCJxpLT4q7I0m4vCKGVxYFUDSof/KL2eqgN1fC794aFESULp8B8eCkSpiN0ZIlKF\nQYSIVGF3ZohCjyOckZ2JNkdX1OMJox1xOHi8hEcckp4wiAyB3IpUk3Ec3N6+iMcTRjvi8I0Pm1H3\nlYNHHJIusTszBHIrUp2evqgrUZVWsrbYPWjp9OLw1w6uZiXdYhAZAqV0lIPJrUSNdMRhp6snLJk3\nV7OSnrA7MwRy6ShDya1EjZbGMs2AoEDC1aykJ2yJDIHcilRL1rioK1GjHXFY/kA2V7OSbnHFqgJp\nNqXDeQUGCPBDxFTLT1CYYwybZTl0woFWhxeF2UYsKZmiOFsTacXqF9/+fURTXI6EVF3ZmarPBQxv\nxSqDiIzQ2RRJ6MxJaDkBgDkr+mxNqFT9R8nn0h+eJxInobMpktCZk9ByImKbrSFKJQwiMiLNwgye\nORnubA1RKtFkdubChQv4xS9+EfjZbrdjz549WLx4cVC5tWvX4tKlS7j99tsBAAcOHNCiOkMWaRZm\n8MzJcGdriFKJJkFk8uTJgYDQ19eH8vJyzJs3T7bs3r17w06CTzRpNkVpTEQa+AwtpzQmkioDpURy\nNF8ncvToUcyZMwfjx498sh+5/SoAFGddpFkUQRDw82VWlBSa8GWrG96eGzBOGIs5hSbcO3VS4P5S\nuVKrOWhmRWkWhigVaR5EamtrsXTpUsXrL7/8MtLS0rB69WqsWbNmyPdX2tgmt19l3iwzzpy/Aqen\nL+gecrMofr8fn3z5HWznBgZFBQEQRQQFEUnoOSE8N4RGk2FP8UZLowkA169fR3l5Oerq6jBu3Liw\nsi6XC2azGd3d3di4cSO2bNmC0tJS2fv6/X74fMFVlTavHf7aAZ9/YOVn+QPZeO7xYjSduYhdfzkR\ntqRcSZoB2PZsKYqn3QVRFPH7vzTim9MXw8r9j/WlcU9VmZZmgC/WiuoIn0t/xoxJG/J7NEujCfzY\nlZELIMCPqTYnTpyIyspKnDx5UjGIyKXRPGX3BHa/AgNLx+u+cuC+/CzYZfakROLzA387ewnZd6bj\nlN0jG0AAoP5LRyBlZbyk6roDPpf+JN06kdraWixbtkz22s2bN+HxDHQV+vv7cezYMUybNm1I94+0\nsU3u5PVIBs+idHJKlihmmgWR69evo7GxMSw1ppRGs7+/H+vXr8fKlSuxatUqPPDAA3j44YeH9BlK\nKRqksZHB+1UiCZ1FiTQlW2I1D6mORKlO18veRVHEv3zcgmM2Z+C1oqlGFOZmIceUAQD4stUNALh/\nxp14Y/8p9PX7AmVvMwA/ezgfOaYMCIIQdFrZ27WtQfcFgIWztTksKFWbx3wu/RmVp70LghCYOQGA\n5g4vmju8QWUMwsCS9P+95WG8/dcWtP29CzPuycS65VbFU8f+afm9ePBeCxpbXAAGWiCccSEKp+sg\nIu1didaWkvawlFrN+MdlhWixe/Blqxtvf9wCU+Y4HGu6AFGmLKdqiaLTdRCJZe+KxC8CHc4rOH7K\nGdZNkStrd/XAmps17IOZiUYLXQeRWPauSAwCYIAQNYBIZXPNE4Z9MDPRaKLrXbyxzsBIX3Y/Ymu2\nSDM1wz2YmWg00XVLZPDelQ7nFbQ7utDc4YEoDixTL5qahRk5xkC348/7mqLe8+HZk7HuVstiKFv9\nOXZCo5WuWyKiKOJUx2U0trjQ5uiC7ZwnMMgqihjY9+IXA92S5nORWwwCgteBxLJgjVv9abTTTUvE\nduYipkwaH5Q17u2PW9FguxDxffuOnov5M8xZ44JaFKFb/ZXGRNgKodFMN0Hkf/7lBOYX/Zg17tAJ\nR9QAMlROTx9qj3eisjT71voTbvUnikY3QcTnHxjELCk0obHVjYam+AYQyftHzsLpvRY048Kt/kTK\ndDUm4heBxlY3PotzCyQUZ1yIYqerIGIQBgY/Y11gNlw8XJkodroJImmGgUHMEqs54ozJk4vysbDY\nouqzOONCFDvdBJFtz5bi58usigvMDMLALtslJVPQ9l1X1PtZssbhn5++D6sX5aEoL4tpLImGSTcD\nq5LQGRNBBMRbLQdrbhY+Pt4Jl/d61Pu4vQPnrC5/MBfLH1ROcUlEkekmiPy/Y2dRNvsf4Lh4NbAC\nVe7L3ubwyrw7nLQhT7oHZ1yIhkc3QeSbtkv4pu0SgMgb3wqzjYET2qNpd3Rh+YNxrSbRqKObMZHB\nIm18W1IyBZas4IOhMzPGhL0GACc7PJzKJVJJNy2RUH4RONTogHS6Y6e7N3Ak4oLiyWjv7ILTexUG\ngwFubx+6en+QvQc3zxGpozqI7Nq1C3/9619xzz334P333w+87nA4UFNTg56eHsydOxc7duwIO7zH\n4/HgpZdegsvlwvTp07Fnzx6MHTs25s+2nfPE3HVRwqlcInVUd2eWL1+ON998M+z13bt3o6amBp9+\n+im6urpw5MiRsDJvvfUWli9fjk8++QRTpkzBvn371FaHiEaY6iBy//33IzMzM+g1URTR1NSEhQsX\nAgAee+wx1NfXh723vr4eK1eujFhGa1yZSqSOJmMiXq83KLBYLBa4XK6wclevXkVGRkbEMlpKMwCz\nZ5iHdUx+vCVDHbTA50p9UYNILDl3R8JHe6pG7LOIKHZRg0gsOXdDGY1GdHX9uPTc6XTCZDKFlRs/\nfjx6e3uRkZGhWIaIkpsm60QEQUBRUREaGhoAAPv370dZWVlYuUWLFuGjjz6KWIaIkpvqNJrbt2/H\n4cOH0dXVhaysLPzud7/DI488Arvdji1btuDKlSuBKV6DwYDXX38ds2bNwiOPPAKPx4PNmzfD5XKh\noKAAf/rTn3DHHXfE69mIaAToJhcvESUnXSx737VrF+bNm4cnn3wy6HWHw4FVq1ahoqIC27dvh97j\nYXl5OR599FFUVVVhw4YNia6OKvX19aisrMSSJUtSav3PzJkzUVVVhaqqKmzdujXR1VFl8+bNmDNn\nDl5++eXAazabDStWrEBFRQX27t0b241EHfj666/F5uZm8Yknngh6fdOmTeLRo0cDfz58+HAiqhc3\nZWVl4vXr1xNdDdV++OEHsbKyUnS5XGJvb69YWVkpejyeRFcrLubNm5foKsTN8ePHxbq6OrGmpibw\n2uOPPy62t7eLN2/eFJ944gmxtbU16n100RJRs6CNRp7NZsP06dNhMpmQnp6ORYsW4fPPP090tShE\naWkp0tPTAz+7XC6IooiCggKkpaVh5cqVsivNQ+kiiMiJdUGb3lRXV2P16tWora1NdFWGze12w2z+\nMQlYqvxuAKC7uxurVq1CdXU1Tpw4kejqxNVwf29JsYs3WRa0aS3ac7777rswm81wuVxYt24drFYr\ncnJyRriWFEldXR3MZjPOnDmDjRs34sCBA5gwYXSvXk2KIKLlgrZkEu05pf8FzGYz5s+fj5aWFl0G\nEZPJFPStd9vUAAABLUlEQVQ/mNPpxMyZMxNYo/iRfkfTpk3D9OnTYbfbUVRUlOBaxYfc7y2W75Ru\nuzOxLmjTi2vXrqG3txcA0NPTg8bGRuTn5ye4VsNTXFyMtrY2uN1uXL16FfX19ViwYEGiq6Vad3c3\n+vv7AQyMH7S3t2PKlCkJrlX8SAHy9OnT8Pl8OHjwYEzfKV2sExnqgjY9+u677/DCCy8AGBg0XrNm\nDZ5++ukE12r46urq8Mc//hF+vx/r16/HU089legqqfbNN99g+/btMBgMMBgMePHFF7F48eJEV2vY\nNm7cCJvNhr6+PkycOBFvvPEG+vv7sXXrVty4cQNVVVXYtGlT1PvoIogQUfLS53/bRJQ0GESISBUG\nESJShUGEiFRhECEiVRhEiEgVBhEiUoVBhIhU+W+yXVhzgnXTqQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc1a92c4eb8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# サンプリングの散布図の描画\n",
"init = [5, 9]\n",
"bs = [0.0, 0.25, 0.5, 0.75, 1.0]\n",
"for b in bs:\n",
" print(\"b = {}\".format(b))\n",
" samples = np.array(mcmc_gibbs_normal(init, 1000, 0, b))\n",
" plt.axes().set_aspect('equal')\n",
" plt.scatter(samples[:, 0], samples[:, 1])\n",
" plt.xlim([-10, 10])\n",
" plt.ylim([-10, 10])\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mcmc_gibbs_ising(init, max_iter, burn_in, theta=0.2):\n",
" \"\"\" \n",
" ギブスサンプリングによるMCMC(イジングモデル)\n",
" init : 初期状態\n",
" max_iter : 繰り返し回数\n",
" burn_in : 捨てる初期サンプル数\n",
" theta : 相関の強さ\n",
" \"\"\"\n",
" samples = []\n",
" x = init\n",
" L = int(np.sqrt(len(x)))\n",
" assert int(L*L) == len(x)\n",
" \n",
" for i in range(max_iter):\n",
" idx = np.random.randint(len(x))\n",
" adjs = [idx+1, idx-1, idx+L, idx-L]\n",
" sum_adj = 0\n",
" for a in adjs:\n",
" if a >= 0 and a < len(x):\n",
" sum_adj += x[a]\n",
" \n",
" p_1 = (1 + np.tanh(theta * sum_adj)) / 2\n",
" if np.random.rand() < p_1:\n",
" x[idx] = 1\n",
" else:\n",
" x[idx] = -1\n",
" \n",
" # サンプル列に追加\n",
" samples.append(copy.deepcopy(x))\n",
" return samples[burn_in:]"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"/* Put everything inside the global mpl namespace */\n",
"window.mpl = {};\n",
"\n",
"\n",
"mpl.get_websocket_type = function() {\n",
" if (typeof(WebSocket) !== 'undefined') {\n",
" return WebSocket;\n",
" } else if (typeof(MozWebSocket) !== 'undefined') {\n",
" return MozWebSocket;\n",
" } else {\n",
" alert('Your browser does not have WebSocket support.' +\n",
" 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
" 'Firefox 4 and 5 are also supported but you ' +\n",
" 'have to enable WebSockets in about:config.');\n",
" };\n",
"}\n",
"\n",
"mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
" this.id = figure_id;\n",
"\n",
" this.ws = websocket;\n",
"\n",
" this.supports_binary = (this.ws.binaryType != undefined);\n",
"\n",
" if (!this.supports_binary) {\n",
" var warnings = document.getElementById(\"mpl-warnings\");\n",
" if (warnings) {\n",
" warnings.style.display = 'block';\n",
" warnings.textContent = (\n",
" \"This browser does not support binary websocket messages. \" +\n",
" \"Performance may be slow.\");\n",
" }\n",
" }\n",
"\n",
" this.imageObj = new Image();\n",
"\n",
" this.context = undefined;\n",
" this.message = undefined;\n",
" this.canvas = undefined;\n",
" this.rubberband_canvas = undefined;\n",
" this.rubberband_context = undefined;\n",
" this.format_dropdown = undefined;\n",
"\n",
" this.image_mode = 'full';\n",
"\n",
" this.root = $('<div/>');\n",
" this._root_extra_style(this.root)\n",
" this.root.attr('style', 'display: inline-block');\n",
"\n",
" $(parent_element).append(this.root);\n",
"\n",
" this._init_header(this);\n",
" this._init_canvas(this);\n",
" this._init_toolbar(this);\n",
"\n",
" var fig = this;\n",
"\n",
" this.waiting = false;\n",
"\n",
" this.ws.onopen = function () {\n",
" fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
" fig.send_message(\"send_image_mode\", {});\n",
" if (mpl.ratio != 1) {\n",
" fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
" }\n",
" fig.send_message(\"refresh\", {});\n",
" }\n",
"\n",
" this.imageObj.onload = function() {\n",
" if (fig.image_mode == 'full') {\n",
" // Full images could contain transparency (where diff images\n",
" // almost always do), so we need to clear the canvas so that\n",
" // there is no ghosting.\n",
" fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
" }\n",
" fig.context.drawImage(fig.imageObj, 0, 0);\n",
" };\n",
"\n",
" this.imageObj.onunload = function() {\n",
" this.ws.close();\n",
" }\n",
"\n",
" this.ws.onmessage = this._make_on_message_function(this);\n",
"\n",
" this.ondownload = ondownload;\n",
"}\n",
"\n",
"mpl.figure.prototype._init_header = function() {\n",
" var titlebar = $(\n",
" '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
" 'ui-helper-clearfix\"/>');\n",
" var titletext = $(\n",
" '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
" 'text-align: center; padding: 3px;\"/>');\n",
" titlebar.append(titletext)\n",
" this.root.append(titlebar);\n",
" this.header = titletext[0];\n",
"}\n",
"\n",
"\n",
"\n",
"mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
"\n",
"}\n",
"\n",
"\n",
"mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
"\n",
"}\n",
"\n",
"mpl.figure.prototype._init_canvas = function() {\n",
" var fig = this;\n",
"\n",
" var canvas_div = $('<div/>');\n",
"\n",
" canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
"\n",
" function canvas_keyboard_event(event) {\n",
" return fig.key_event(event, event['data']);\n",
" }\n",
"\n",
" canvas_div.keydown('key_press', canvas_keyboard_event);\n",
" canvas_div.keyup('key_release', canvas_keyboard_event);\n",
" this.canvas_div = canvas_div\n",
" this._canvas_extra_style(canvas_div)\n",
" this.root.append(canvas_div);\n",
"\n",
" var canvas = $('<canvas/>');\n",
" canvas.addClass('mpl-canvas');\n",
" canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
"\n",
" this.canvas = canvas[0];\n",
" this.context = canvas[0].getContext(\"2d\");\n",
"\n",
" var backingStore = this.context.backingStorePixelRatio ||\n",
"\tthis.context.webkitBackingStorePixelRatio ||\n",
"\tthis.context.mozBackingStorePixelRatio ||\n",
"\tthis.context.msBackingStorePixelRatio ||\n",
"\tthis.context.oBackingStorePixelRatio ||\n",
"\tthis.context.backingStorePixelRatio || 1;\n",
"\n",
" mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
"\n",
" var rubberband = $('<canvas/>');\n",
" rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
"\n",
" var pass_mouse_events = true;\n",
"\n",
" canvas_div.resizable({\n",
" start: function(event, ui) {\n",
" pass_mouse_events = false;\n",
" },\n",
" resize: function(event, ui) {\n",
" fig.request_resize(ui.size.width, ui.size.height);\n",
" },\n",
" stop: function(event, ui) {\n",
" pass_mouse_events = true;\n",
" fig.request_resize(ui.size.width, ui.size.height);\n",
" },\n",
" });\n",
"\n",
" function mouse_event_fn(event) {\n",
" if (pass_mouse_events)\n",
" return fig.mouse_event(event, event['data']);\n",
" }\n",
"\n",
" rubberband.mousedown('button_press', mouse_event_fn);\n",
" rubberband.mouseup('button_release', mouse_event_fn);\n",
" // Throttle sequential mouse events to 1 every 20ms.\n",
" rubberband.mousemove('motion_notify', mouse_event_fn);\n",
"\n",
" rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
" rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
"\n",
" canvas_div.on(\"wheel\", function (event) {\n",
" event = event.originalEvent;\n",
" event['data'] = 'scroll'\n",
" if (event.deltaY < 0) {\n",
" event.step = 1;\n",
" } else {\n",
" event.step = -1;\n",
" }\n",
" mouse_event_fn(event);\n",
" });\n",
"\n",
" canvas_div.append(canvas);\n",
" canvas_div.append(rubberband);\n",
"\n",
" this.rubberband = rubberband;\n",
" this.rubberband_canvas = rubberband[0];\n",
" this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
" this.rubberband_context.strokeStyle = \"#000000\";\n",
"\n",
" this._resize_canvas = function(width, height) {\n",
" // Keep the size of the canvas, canvas container, and rubber band\n",
" // canvas in synch.\n",
" canvas_div.css('width', width)\n",
" canvas_div.css('height', height)\n",
"\n",
" canvas.attr('width', width * mpl.ratio);\n",
" canvas.attr('height', height * mpl.ratio);\n",
" canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
"\n",
" rubberband.attr('width', width);\n",
" rubberband.attr('height', height);\n",
" }\n",
"\n",
" // Set the figure to an initial 600x600px, this will subsequently be updated\n",
" // upon first draw.\n",
" this._resize_canvas(600, 600);\n",
"\n",
" // Disable right mouse context menu.\n",
" $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
" return false;\n",
" });\n",
"\n",
" function set_focus () {\n",
" canvas.focus();\n",
" canvas_div.focus();\n",
" }\n",
"\n",
" window.setTimeout(set_focus, 100);\n",
"}\n",
"\n",
"mpl.figure.prototype._init_toolbar = function() {\n",
" var fig = this;\n",
"\n",
" var nav_element = $('<div/>')\n",
" nav_element.attr('style', 'width: 100%');\n",
" this.root.append(nav_element);\n",
"\n",
" // Define a callback function for later on.\n",
" function toolbar_event(event) {\n",
" return fig.toolbar_button_onclick(event['data']);\n",
" }\n",
" function toolbar_mouse_event(event) {\n",
" return fig.toolbar_button_onmouseover(event['data']);\n",
" }\n",
"\n",
" for(var toolbar_ind in mpl.toolbar_items) {\n",
" var name = mpl.toolbar_items[toolbar_ind][0];\n",
" var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
" var image = mpl.toolbar_items[toolbar_ind][2];\n",
" var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
"\n",
" if (!name) {\n",
" // put a spacer in here.\n",
" continue;\n",
" }\n",
" var button = $('<button/>');\n",
" button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
" 'ui-button-icon-only');\n",
" button.attr('role', 'button');\n",
" button.attr('aria-disabled', 'false');\n",
" button.click(method_name, toolbar_event);\n",
" button.mouseover(tooltip, toolbar_mouse_event);\n",
"\n",
" var icon_img = $('<span/>');\n",
" icon_img.addClass('ui-button-icon-primary ui-icon');\n",
" icon_img.addClass(image);\n",
" icon_img.addClass('ui-corner-all');\n",
"\n",
" var tooltip_span = $('<span/>');\n",
" tooltip_span.addClass('ui-button-text');\n",
" tooltip_span.html(tooltip);\n",
"\n",
" button.append(icon_img);\n",
" button.append(tooltip_span);\n",
"\n",
" nav_element.append(button);\n",
" }\n",
"\n",
" var fmt_picker_span = $('<span/>');\n",
"\n",
" var fmt_picker = $('<select/>');\n",
" fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
" fmt_picker_span.append(fmt_picker);\n",
" nav_element.append(fmt_picker_span);\n",
" this.format_dropdown = fmt_picker[0];\n",
"\n",
" for (var ind in mpl.extensions) {\n",
" var fmt = mpl.extensions[ind];\n",
" var option = $(\n",
" '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
" fmt_picker.append(option)\n",
" }\n",
"\n",
" // Add hover states to the ui-buttons\n",
" $( \".ui-button\" ).hover(\n",
" function() { $(this).addClass(\"ui-state-hover\");},\n",
" function() { $(this).removeClass(\"ui-state-hover\");}\n",
" );\n",
"\n",
" var status_bar = $('<span class=\"mpl-message\"/>');\n",
" nav_element.append(status_bar);\n",
" this.message = status_bar[0];\n",
"}\n",
"\n",
"mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
" // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
" // which will in turn request a refresh of the image.\n",
" this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
"}\n",
"\n",
"mpl.figure.prototype.send_message = function(type, properties) {\n",
" properties['type'] = type;\n",
" properties['figure_id'] = this.id;\n",
" this.ws.send(JSON.stringify(properties));\n",
"}\n",
"\n",
"mpl.figure.prototype.send_draw_message = function() {\n",
" if (!this.waiting) {\n",
" this.waiting = true;\n",
" this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
" }\n",
"}\n",
"\n",
"\n",
"mpl.figure.prototype.handle_save = function(fig, msg) {\n",
" var format_dropdown = fig.format_dropdown;\n",
" var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
" fig.ondownload(fig, format);\n",
"}\n",
"\n",
"\n",
"mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
" var size = msg['size'];\n",
" if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
" fig._resize_canvas(size[0], size[1]);\n",
" fig.send_message(\"refresh\", {});\n",
" };\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
" var x0 = msg['x0'] / mpl.ratio;\n",
" var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
" var x1 = msg['x1'] / mpl.ratio;\n",
" var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
" x0 = Math.floor(x0) + 0.5;\n",
" y0 = Math.floor(y0) + 0.5;\n",
" x1 = Math.floor(x1) + 0.5;\n",
" y1 = Math.floor(y1) + 0.5;\n",
" var min_x = Math.min(x0, x1);\n",
" var min_y = Math.min(y0, y1);\n",
" var width = Math.abs(x1 - x0);\n",
" var height = Math.abs(y1 - y0);\n",
"\n",
" fig.rubberband_context.clearRect(\n",
" 0, 0, fig.canvas.width, fig.canvas.height);\n",
"\n",
" fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
" // Updates the figure title.\n",
" fig.header.textContent = msg['label'];\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
" var cursor = msg['cursor'];\n",
" switch(cursor)\n",
" {\n",
" case 0:\n",
" cursor = 'pointer';\n",
" break;\n",
" case 1:\n",
" cursor = 'default';\n",
" break;\n",
" case 2:\n",
" cursor = 'crosshair';\n",
" break;\n",
" case 3:\n",
" cursor = 'move';\n",
" break;\n",
" }\n",
" fig.rubberband_canvas.style.cursor = cursor;\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_message = function(fig, msg) {\n",
" fig.message.textContent = msg['message'];\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
" // Request the server to send over a new figure.\n",
" fig.send_draw_message();\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
" fig.image_mode = msg['mode'];\n",
"}\n",
"\n",
"mpl.figure.prototype.updated_canvas_event = function() {\n",
" // Called whenever the canvas gets updated.\n",
" this.send_message(\"ack\", {});\n",
"}\n",
"\n",
"// A function to construct a web socket function for onmessage handling.\n",
"// Called in the figure constructor.\n",
"mpl.figure.prototype._make_on_message_function = function(fig) {\n",
" return function socket_on_message(evt) {\n",
" if (evt.data instanceof Blob) {\n",
" /* FIXME: We get \"Resource interpreted as Image but\n",
" * transferred with MIME type text/plain:\" errors on\n",
" * Chrome. But how to set the MIME type? It doesn't seem\n",
" * to be part of the websocket stream */\n",
" evt.data.type = \"image/png\";\n",
"\n",
" /* Free the memory for the previous frames */\n",
" if (fig.imageObj.src) {\n",
" (window.URL || window.webkitURL).revokeObjectURL(\n",
" fig.imageObj.src);\n",
" }\n",
"\n",
" fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
" evt.data);\n",
" fig.updated_canvas_event();\n",
" fig.waiting = false;\n",
" return;\n",
" }\n",
" else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
" fig.imageObj.src = evt.data;\n",
" fig.updated_canvas_event();\n",
" fig.waiting = false;\n",
" return;\n",
" }\n",
"\n",
" var msg = JSON.parse(evt.data);\n",
" var msg_type = msg['type'];\n",
"\n",
" // Call the \"handle_{type}\" callback, which takes\n",
" // the figure and JSON message as its only arguments.\n",
" try {\n",
" var callback = fig[\"handle_\" + msg_type];\n",
" } catch (e) {\n",
" console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
" return;\n",
" }\n",
"\n",
" if (callback) {\n",
" try {\n",
" // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
" callback(fig, msg);\n",
" } catch (e) {\n",
" console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
" }\n",
" }\n",
" };\n",
"}\n",
"\n",
"// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
"mpl.findpos = function(e) {\n",
" //this section is from http://www.quirksmode.org/js/events_properties.html\n",
" var targ;\n",
" if (!e)\n",
" e = window.event;\n",
" if (e.target)\n",
" targ = e.target;\n",
" else if (e.srcElement)\n",
" targ = e.srcElement;\n",
" if (targ.nodeType == 3) // defeat Safari bug\n",
" targ = targ.parentNode;\n",
"\n",
" // jQuery normalizes the pageX and pageY\n",
" // pageX,Y are the mouse positions relative to the document\n",
" // offset() returns the position of the element relative to the document\n",
" var x = e.pageX - $(targ).offset().left;\n",
" var y = e.pageY - $(targ).offset().top;\n",
"\n",
" return {\"x\": x, \"y\": y};\n",
"};\n",
"\n",
"/*\n",
" * return a copy of an object with only non-object keys\n",
" * we need this to avoid circular references\n",
" * http://stackoverflow.com/a/24161582/3208463\n",
" */\n",
"function simpleKeys (original) {\n",
" return Object.keys(original).reduce(function (obj, key) {\n",
" if (typeof original[key] !== 'object')\n",
" obj[key] = original[key]\n",
" return obj;\n",
" }, {});\n",
"}\n",
"\n",
"mpl.figure.prototype.mouse_event = function(event, name) {\n",
" var canvas_pos = mpl.findpos(event)\n",
"\n",
" if (name === 'button_press')\n",
" {\n",
" this.canvas.focus();\n",
" this.canvas_div.focus();\n",
" }\n",
"\n",
" var x = canvas_pos.x * mpl.ratio;\n",
" var y = canvas_pos.y * mpl.ratio;\n",
"\n",
" this.send_message(name, {x: x, y: y, button: event.button,\n",
" step: event.step,\n",
" guiEvent: simpleKeys(event)});\n",
"\n",
" /* This prevents the web browser from automatically changing to\n",
" * the text insertion cursor when the button is pressed. We want\n",
" * to control all of the cursor setting manually through the\n",
" * 'cursor' event from matplotlib */\n",
" event.preventDefault();\n",
" return false;\n",
"}\n",
"\n",
"mpl.figure.prototype._key_event_extra = function(event, name) {\n",
" // Handle any extra behaviour associated with a key event\n",
"}\n",
"\n",
"mpl.figure.prototype.key_event = function(event, name) {\n",
"\n",
" // Prevent repeat events\n",
" if (name == 'key_press')\n",
" {\n",
" if (event.which === this._key)\n",
" return;\n",
" else\n",
" this._key = event.which;\n",
" }\n",
" if (name == 'key_release')\n",
" this._key = null;\n",
"\n",
" var value = '';\n",
" if (event.ctrlKey && event.which != 17)\n",
" value += \"ctrl+\";\n",
" if (event.altKey && event.which != 18)\n",
" value += \"alt+\";\n",
" if (event.shiftKey && event.which != 16)\n",
" value += \"shift+\";\n",
"\n",
" value += 'k';\n",
" value += event.which.toString();\n",
"\n",
" this._key_event_extra(event, name);\n",
"\n",
" this.send_message(name, {key: value,\n",
" guiEvent: simpleKeys(event)});\n",
" return false;\n",
"}\n",
"\n",
"mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
" if (name == 'download') {\n",
" this.handle_save(this, null);\n",
" } else {\n",
" this.send_message(\"toolbar_button\", {name: name});\n",
" }\n",
"};\n",
"\n",
"mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
" this.message.textContent = tooltip;\n",
"};\n",
"mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
"\n",
"mpl.extensions = [\"eps\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\"];\n",
"\n",
"mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
" // Create a \"websocket\"-like object which calls the given IPython comm\n",
" // object with the appropriate methods. Currently this is a non binary\n",
" // socket, so there is still some room for performance tuning.\n",
" var ws = {};\n",
"\n",
" ws.close = function() {\n",
" comm.close()\n",
" };\n",
" ws.send = function(m) {\n",
" //console.log('sending', m);\n",
" comm.send(m);\n",
" };\n",
" // Register the callback with on_msg.\n",
" comm.on_msg(function(msg) {\n",
" //console.log('receiving', msg['content']['data'], msg);\n",
" // Pass the mpl event to the overriden (by mpl) onmessage function.\n",
" ws.onmessage(msg['content']['data'])\n",
" });\n",
" return ws;\n",
"}\n",
"\n",
"mpl.mpl_figure_comm = function(comm, msg) {\n",
" // This is the function which gets called when the mpl process\n",
" // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
"\n",
" var id = msg.content.data.id;\n",
" // Get hold of the div created by the display call when the Comm\n",
" // socket was opened in Python.\n",
" var element = $(\"#\" + id);\n",
" var ws_proxy = comm_websocket_adapter(comm)\n",
"\n",
" function ondownload(figure, format) {\n",
" window.open(figure.imageObj.src);\n",
" }\n",
"\n",
" var fig = new mpl.figure(id, ws_proxy,\n",
" ondownload,\n",
" element.get(0));\n",
"\n",
" // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
" // web socket which is closed, not our websocket->open comm proxy.\n",
" ws_proxy.onopen();\n",
"\n",
" fig.parent_element = element.get(0);\n",
" fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
" if (!fig.cell_info) {\n",
" console.error(\"Failed to find cell for figure\", id, fig);\n",
" return;\n",
" }\n",
"\n",
" var output_index = fig.cell_info[2]\n",
" var cell = fig.cell_info[0];\n",
"\n",
"};\n",
"\n",
"mpl.figure.prototype.handle_close = function(fig, msg) {\n",
" var width = fig.canvas.width/mpl.ratio\n",
" fig.root.unbind('remove')\n",
"\n",
" // Update the output cell to use the data from the current canvas.\n",
" fig.push_to_output();\n",
" var dataURL = fig.canvas.toDataURL();\n",
" // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
" // the notebook keyboard shortcuts fail.\n",
" IPython.keyboard_manager.enable()\n",
" $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
" fig.close_ws(fig, msg);\n",
"}\n",
"\n",
"mpl.figure.prototype.close_ws = function(fig, msg){\n",
" fig.send_message('closing', msg);\n",
" // fig.ws.close()\n",
"}\n",
"\n",
"mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
" // Turn the data on the canvas into data in the output cell.\n",
" var width = this.canvas.width/mpl.ratio\n",
" var dataURL = this.canvas.toDataURL();\n",
" this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
"}\n",
"\n",
"mpl.figure.prototype.updated_canvas_event = function() {\n",
" // Tell IPython that the notebook contents must change.\n",
" IPython.notebook.set_dirty(true);\n",
" this.send_message(\"ack\", {});\n",
" var fig = this;\n",
" // Wait a second, then push the new image to the DOM so\n",
" // that it is saved nicely (might be nice to debounce this).\n",
" setTimeout(function () { fig.push_to_output() }, 1000);\n",
"}\n",
"\n",
"mpl.figure.prototype._init_toolbar = function() {\n",
" var fig = this;\n",
"\n",
" var nav_element = $('<div/>')\n",
" nav_element.attr('style', 'width: 100%');\n",
" this.root.append(nav_element);\n",
"\n",
" // Define a callback function for later on.\n",
" function toolbar_event(event) {\n",
" return fig.toolbar_button_onclick(event['data']);\n",
" }\n",
" function toolbar_mouse_event(event) {\n",
" return fig.toolbar_button_onmouseover(event['data']);\n",
" }\n",
"\n",
" for(var toolbar_ind in mpl.toolbar_items){\n",
" var name = mpl.toolbar_items[toolbar_ind][0];\n",
" var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
" var image = mpl.toolbar_items[toolbar_ind][2];\n",
" var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
"\n",
" if (!name) { continue; };\n",
"\n",
" var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
" button.click(method_name, toolbar_event);\n",
" button.mouseover(tooltip, toolbar_mouse_event);\n",
" nav_element.append(button);\n",
" }\n",
"\n",
" // Add the status bar.\n",
" var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
" nav_element.append(status_bar);\n",
" this.message = status_bar[0];\n",
"\n",
" // Add the close button to the window.\n",
" var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
" var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
" button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
" button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
" buttongrp.append(button);\n",
" var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
" titlebar.prepend(buttongrp);\n",
"}\n",
"\n",
"mpl.figure.prototype._root_extra_style = function(el){\n",
" var fig = this\n",
" el.on(\"remove\", function(){\n",
"\tfig.close_ws(fig, {});\n",
" });\n",
"}\n",
"\n",
"mpl.figure.prototype._canvas_extra_style = function(el){\n",
" // this is important to make the div 'focusable\n",
" el.attr('tabindex', 0)\n",
" // reach out to IPython and tell the keyboard manager to turn it's self\n",
" // off when our div gets focus\n",
"\n",
" // location in version 3\n",
" if (IPython.notebook.keyboard_manager) {\n",
" IPython.notebook.keyboard_manager.register_events(el);\n",
" }\n",
" else {\n",
" // location in version 2\n",
" IPython.keyboard_manager.register_events(el);\n",
" }\n",
"\n",
"}\n",
"\n",
"mpl.figure.prototype._key_event_extra = function(event, name) {\n",
" var manager = IPython.notebook.keyboard_manager;\n",
" if (!manager)\n",
" manager = IPython.keyboard_manager;\n",
"\n",
" // Check for shift+enter\n",
" if (event.shiftKey && event.which == 13) {\n",
" this.canvas_div.blur();\n",
" // select the cell after this one\n",
" var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n",
" IPython.notebook.select(index + 1);\n",
" }\n",
"}\n",
"\n",
"mpl.figure.prototype.handle_save = function(fig, msg) {\n",
" fig.ondownload(fig, null);\n",
"}\n",
"\n",
"\n",
"mpl.find_output_cell = function(html_output) {\n",
" // Return the cell and output element which can be found *uniquely* in the notebook.\n",
" // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
" // IPython event is triggered only after the cells have been serialised, which for\n",
" // our purposes (turning an active figure into a static one), is too late.\n",
" var cells = IPython.notebook.get_cells();\n",
" var ncells = cells.length;\n",
" for (var i=0; i<ncells; i++) {\n",
" var cell = cells[i];\n",
" if (cell.cell_type === 'code'){\n",
" for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
" var data = cell.output_area.outputs[j];\n",
" if (data.data) {\n",
" // IPython >= 3 moved mimebundle to data attribute of output\n",
" data = data.data;\n",
" }\n",
" if (data['text/html'] == html_output) {\n",
" return [cell, data, j];\n",
" }\n",
" }\n",
" }\n",
" }\n",
"}\n",
"\n",
"// Register the function which deals with the matplotlib target/channel.\n",
"// The kernel may be null if the page has been refreshed.\n",
"if (IPython.notebook.kernel != null) {\n",
" IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
"}\n"
],
"text/plain": [
"<IPython.core.display.Javascript object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAYAAADVKCZpAAAgAElEQVR4Xu2dX4xcV33Hf+fuLnHipIlVVFIHCqx3zowRtozypyjtgys/ETkFVNom5AWkEIyaSEUgqkp5IEECNYK+JKgkTSKRqC0QFxIFp6RN1ISQqk1MRUO2nXtm7KxQtOKhNYpMvd717r3V0FlYb3Y9c+6dc+d37nz8grJzfuf8zud72E/OnfHECH8gAAEIQAACERIwEfZMyxCAAAQgAAFBYBwCCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizI2moYABCAAAQTGGYAABCAAgSgJILAoY6NpCEAAAhBAYJwBCEAAAhCIkgACizK285teOpfnw25j17W3DzuUcSUJ/Ozl+0rOQLkWAhfPGH5XagljQx+EojAU35YQmC+xasYjsGo4V7EKAquCsv8aCMyfmboKBKYukl80hMB05lKkKwRWhFr4GgQWnnHwFRBYcMSFFkBghbCpLEJgKmMRBKYzF6+uEJgXrsoGI7DKUAdfCIEFR1xoAQRWCJuuIgSmK4/1bhCYzlyKdIXAilALX4PAwjMOvgICC4640AIIrBA2lUUITGUsPELUGYtfVwjMj1dVoxFYVaTDr4PAwjMusgI3sCLUlNUgMGWB9NtBYDpzKdIVAitCLXwNAgvPOPgKCCw44kILILBC2FQWITCVsfAIUWcsfl0hMD9eVY1GYFWRDr8OAgvPuMgK3MCKUFNWg8CUBcIjRJ2BlOgKgZWAF7AUgQWEW9XUCKwq0n7rcAPz46V5NALTmQ4CC5hLs9k8muf5IRF52jl30549e94xNTX15IYlrYh81Dn3+MY2ms3mc3mev01Elns/d84duFCbCCxgiCWmRmAl4CkrRWDKAum3g8AC5tJqtQ6urq5eliTJLT2BbVxq9+7dl1x66aULO3bsePcrr7zyv1sI7Ihzrj1MewhsGErVj0Fg1TMPtSICC0W23LwIrBy/gdU9iWVZ1pPReQJrNpsfEZGb0jTt/e95f/o3MAQ2kK7uAQhMdz4+3SEwH1rVjUVggVlvJzBr7TeNMX+fpum3thHYLhFZE5GHnHNf5RFi4KACTI/AAkAd05QIbEzgByyLwALnspXA3v72t198ySWXLPz85z9/9+Li4pktBLY7TdPFffv27VpeXn4qSZI/b7fbz23XKo8QA4dYcHoEVhCcwjIEpjAUEf4eWOhYthJYo9H4A2PMHzvn/mjQ+tbaz+R5bjqdzpcR2CBaul5HYLryKNMNAitDL1wtN7BwbH8x81YCs9Z+wxhzNE3To5uXP3jw4PTi4uIVzrn/npubuyhJku+KyF865/4BgQUOa8TTI7ARAx3jdAhsjPAvsDQCC5iLtfaYiFwnIjtF5FSWZTeePXu2fckll7x25syZd7/++utL68tbax80xnztoosu+q/l5eXv53k+IyKJiDzmnLvrQm3yCDFgiCWmRmAl4CkrRWDKAum3g8B05uLVFQLzwlXZYARWGergCyGw4IgLLYDACmHTVYTAdOWx3g0C05lLka4QWBFq4WsQWHjGwVfwEZhPM7uuvd1nOGM3EUBgbz4SPmdKEz8EpvP/3ghMZy5eXSEwL1yVDdb0C7iyTQ9YCIFpSaIefSCwGuSIwHSGiMC4gek8mfXpCoHVIEsEpjNEBIbAdJ7M+nSFwGqQJQLTGSICQ2A6T2Z9ukJgNcgSgekMEYEhMJ0nsz5dIbAaZInAdIaIwBCYzpNZn64QWA2yRGA6Q0RgCEznyaxPVwisBlkiMJ0hIjAEpvNk1qcrBFaDLBGYzhARGALTeTLr0xUCq0GWCExniAgMgek8mfXpCoHVIMtQAguJxucbGUL24Ts3UjqfWKw5+uZ+9kdf5XelL7QKxhNKBZBDL4HAQhP+1fwIDIFVd9pYaRABBDaIUASvI7DqQkJgCKy608ZKgwggsEGEIngdgVUXEgJDYNWdNlYaRACBDSIUwesIrLqQEBgCq+60sdIgAghsEKEIXkdg1YWEwBBYdaeNlQYRQGCDCEXwOgKrLiQEhsCqO22sNIgAAhtEqMTrzWbzaJ7nh0TkaefcTb2prLULIvKGiOQisuicu2HzEnNzc3uSJPmmiFxhjHkmTdNP9cdv2Q0CKxGSZykCQ2CeR4bhAQkgsIBwW63WwdXV1cuSJLllo8BWVlZaCwsLZ7dbuic+EXkwTdPvNRqNx0Tk651O57vbjUdgAUPcNDUCQ2DVnTZWGkQAgQ0iVPL1nsSyLDviITBjrf2Jc+4dvaXn5uZuTJLksHPukwisZBgjKEdgCGwEx4gpRkQAgY0I5HbTbCGw14wxp/I8XzPG3JOmae+29cs/1tq3isgzzrkDfYG9L0mSLzjnDiOwwGENMT0CQ2BDHBOGVEQAgQUGvVlgzWZzd5qmi9baq0Tk2bW1tcMnTpzorreBwAIHUnJ6BIbASh4hykdIAIGNEOZWU20W2Kbb1r3GmOc33cK2eoR4o3PuNm5ggcMaYvpJENikfL/hEHH/cgjfhehDq7qxCCww640C279//87l5eUkTdPTs7Ozl09PT7+QZdnN3W53fpPYvm2MeWD9Qxx5nj/S7XafRGCBwxpiegQ2BKQaDkFgOkNFYAFzsdYeE5HrRGSniJxKkuSDWZY93Fsyz/NERO7rdDr39/7ZWvugMeZraZoe37t3byPLsm/keX5F7zGjc+6IiGQILGBYQ06NwIYEVbNhCExnoAhMZy5eXfExei9cpQYjsFL4oi1GYDqjQ2A6c/HqCoF54So1GIGVwhdtMQLTGR0C05mLV1cIzAtXqcEIrBS+aIsRmM7oEJjOXLy6QmBeuEoNRmCl8EVbjMB0RofAdObi1RUC88JVajACK4Uv2mIEpjM6BKYzF6+uEJgXrlKDEVgpfNEWIzCd0SEwnbl4dYXAvHCVGozASuGLthiB6YwOgenMxasrBOaFq9RgBFYKX7TFCExndAhMZy5eXSEwL1ylBiOwUviiLUZgOqNDYDpz8eoKgXnhKjUYgZXCF20xAtMZHQLTmYtXVwjMC1epwQisFL5oixGYzugQmM5cvLpCYF64Sg1GYKXwRVuMwHRGh8B05uLVFQLzwlVqMAIrhS/aYgSmMzoEpjMXr64QmBeuUoMRWCl80RYjMJ3RITCduXh1hcC8cJUajMBK4Yu2GIHpjA6B6czFqysE5oWr1GAEVgpftMUITGd0CExnLl5dITAvXKUGI7BS+KItRmA6o0NgOnPx6gqBeeEqNRiBlcIXbTEC0xkdAtOZi1dXCMwLV6nBCKwUvmiLEZjO6BCYzly8ukJgXrhKDUZgpfBFW4zAdEaHwHTm4tXVjgN/knsVBBrs88t917W3B+oi7LQ+ewzbSbjZY80mHBERBBaSbvG5EVhxdgMrm83m0TzPD4nI0865m1qt1q9nWfYtEfnNPM/XkiS5O03TxzZPZK1dEJE3RKQnpkXn3A0XWgyBDYxiZAMQ2MhQRjURAtMZFwILmEur1Tq4urp6WZIkt/QEtm/fvl3Ly8tzzrmX9+zZ8xtTU1P/fubMmcbrr7++tLGNnsBWVlZaCwsLZ4dpD4ENQ2k0YxDYaDjGNgsC05kYAgucS09iWZYd6Qlsi5vWfxhjPpCm6SICCxzEiKZHYCMCGdk0CExnYAgscC7bCWxubu59SZI86px77xZie80Yc6r3mNEYc0+apkd5hBg4qCGnR2BDgqrZMASmM1AEFjiXrQT2rne964qZmZkXpqamPtlut/9lcwvNZnN371Zmrb1KRJ5dW1s7fOLEie52rfIIMXCIG6ZHYNWx1rQSAtOUxq96QWCBc9kssKuvvnrm9OnT3zPG/E2apg8PWt5ae68x5vkL3cIQ2CCKo3sdgY2OZUwzITCdaSGwwLlsFpi19q/zPD/V6XT+bKul9+/fv3N5eTlJ0/T07Ozs5dPT0y9kWXZzt9ud5wYWOKwhpkdgQ0Cq4RAEpjNUBBYwF2vtMRG5TkR2isgpEblZRL4vIq/0PyIvWZZ9tNvt/qe19iljzK15nu8Qke/02srzPBGR+zqdzv28BxYwKI+pEZgHrBoNRWA6w0RgOnPx6opHiF64Sg1GYKXwRVuMwHRGh8B05uLVFQLzwlVqMAIrhS/aYgSmMzoEpjMXr65CCUzLL+uQX22kZY+9wEPu0+tAMfhNBBCYzkOBwHTm4tUVAvPCdd5gBFac3SRVIjCdaSMwnbl4dYXAvHAhsOK4JrYSgemMHoHpzMWrKwTmhQuBFcc1sZUITGf0CExnLl5dITAvXAisOK6JrURgOqNHYDpz8eoKgXnhQmDFcU1sJQLTGT0C05mLV1cIzAsXAiuOa2IrEZjO6BGYzly8ukJgXrgQWHFcE1uJwHRGj8B05uLVFQLzwoXAiuOa2EoEpjN6BKYzF6+uEJgXLgRWHNfEViIwndEjMJ25eHWFwLxwIbDiuCa2EoHpjB6B6czFqysE5oULgRXHNbGVCExn9AhMZy5eXfkITNNXJ3ltcgIG812IekNGYDqzQWA6c/HqCoF54VI7GIGpjUYQmM5sEJjOXLy6QmBeuNQORmBqo0FgSqNBYEqD8WkLgfnQ0jsWgenNhhuYzmwQmM5cvLpCYF641A5GYGqj4QamNBoEpjQYn7YQmA8tvWMRmN5suIHpzAaB6czFqysE5oVL7WAEpjYabmBKo0FgAYNpNptH8zw/JCJPO+du6i1lrb1ORB4WkR0i8ohz7u7NLczNze1JkuSbInKFMeaZNE0/JSL5dq0isIAhVjg1AqsQtudS3MA8gVU0HIEFBN1qtQ6urq5eliTJLRsE9lKWZR/vdrtta+2LSZJ8ot1u/3hjGz3xiciDaZp+r9FoPCYiX+90Ot9FYAHDUjA1AlMQwjYtIDCd2SCwwLn0JJZl2ZGewJrN5u48z59wzl3bW7bRaNxhjLnUOfelDW0Ya+1PnHPv6P1sbm7uxiRJDjvnPonAAoc15ukR2JgDuMDyCExnNggscC6bBHZNnud3Ouc+1H+c+GEROeScu329DWvtW0XkGefcgb7A3pckyRecc4cRWOCwxjw9AhtzAAhMbwDbdIbAAkdWhcCWzuXbvj8WeHvRT480oo+wkg1wA6sEs/ciCMwbmV/BEI8QL3POfXHAI8QbnXO3bbcyAvPLZONoBFac3SRVIjCdaSOwwLlsFFj/seHLWZZ9rP8hjh8kSXLb5g9xWGu/bYx5YP1DHHmeP9Ltdp9EYKMPC4GNnmkdZ0RgOlNFYAFzsdYeE5Hex+Z3isipLMtuFJGLkiR5qP8x+kedc5/vi+1BY8zX0jQ9vnfv3kaWZd/I8/wKEXnWOXdERDIENvqwENjomdZxRgSmM1UEpjMXr654hOiF67zBCKw4u0mqRGA600ZgOnPx6gqBeeFCYMVxTWwlAtMZPQLTmYtXVwjMCxcCK45rYisRmM7oEZjOXLy6QmBeuBBYcVwTW4nAdEaPwHTm4tUVAvPChcCK45rYSgSmM3oEpjMXr64QmBcuBFYc18RWIjCd0SMwnbl4dYXAvHAhsOK4JrYSgemMHoHpzMWrKwTmhQuBFcc1sZUITGf0CExnLl5dIbDzcfF3u7yOz0QO/tnL93nt++IZw+9KL2LVDCaUajgHXQWBIbCgB6yGkyOweoSKwGqQIwJDYDU4xpVuAYFVijvYYggsGNrqJkZgCKy601aPlRBYPXJEYDXIEYEhsBoc40q3gMAqxR1sMQQWDG11EyMwBFbdaavHSgisHjkisBrkiMAQWA2OcaVbQGCV4g62GAILhra6iREYAqvutNVjJQRWjxwRWA1yRGAIrAbHuNItILBKcQdbDIEFQ1vdxAgMgVV32uqxEgKrR44IrAY5IjAEVoNjXOkWEFiluIMthsCCoa1uYgSGwKo7bfVYCYHVI0cENoYcG43Gbxtj7t+w9HvzPL+m0+n8aP1n1toFEXlDRHIRWXTO3bBdqwgMgY3hGEe9JAKLOr5fNo/Axpzj3r1737m2tva8c+5dG1vpCWxlZaW1sLBwdlCLCAyBDTojvH4+AQRWjxOBwMacY6PR+Kwx5krn3GcR2GjC4NvoR8OxzrMgsHqki8DGnKO19t9E5A7n3EubBPaaMeZUnudrxph70jQ9yiPE4cJCYMNxmuRRCKwe6SOwMebYf3z4z8652c1tNJvN3WmaLlprrxKRZ9fW1g6fOHGiu1W7PELkEeIYj3GUSyOwKGN7U9MIbIw5Wms/Y4x5W5qmn7tQG9bae40xz293C0NgCGyMxzjKpRFYlLEhME2xWWv/1Rhze5qmxzf2tX///p3Ly8tJmqanZ2dnL5+enn4hy7Kbu93uPDewwQnyCHEwo0kfgcDqcQK4gY0px9nZ2d+anp7uPT7cs96CtfYpY8yteZ7vEJHv9H6e53kiIvd1Op2NH7s/r2tuYNzAxnSMo10WgUUb3XmNI7Aa5IjAEFgNjnGlW0BgleIOthgCC4a2uokRGAKr7rTVYyUEVo8cEVgNckRgxUPk/bLi7GKuRGAxp/er3hFYDXJEYMVDRGDF2cVcicBiTg+B1SO9/i4QWPE4EVhxdjFXIrCY00Ng9UgPgZXOEYGVRhjlBAgsytje1DSPEGuQIzew4iEisOLsYq5EYDGnxw2sHulxAyudIwIrjTDKCRBYlLFxA6tHbOfvghtY8VQRWHF2MVcisJjT4wZWj/S4gZXOEYGVRhjlBAgsyti4gdUjNm5go8oRgY2KZFzzILC48tquWz7EUYMceYRYPEQEVpxdzJUILOb0eIRYj/R4hFg6RwRWGmGUEyCwKGPjEWI9YuMR4qhyRGCjIhnXPAgsrrx4hFiPvLbcBY8Qi4eLwIqzi7kSgcWcHo8Q65EejxBL54jASiOMcgIEFmVsPEKsR2w8QhxVjghsVCTjmgeBxZUXjxDrkRePEEecIwIbMdBIpkNgkQQ1oE0+Rl+DHHkPrHiICKw4u5grEVjM6fEeWD3S4z2w0jkisNIIo5wAgUUZG++B1SM23gMbVY4IbFQk45oHgcWVF++BKcvLWntOROb7bR13zt26sUVr7XUi8rCI7BCRR5xzd2+3BR4hFg8XgRVnF3MlAos5PR4hjj09a+1PnXNXbteItfalLMs+3u1229baF5Mk+US73f7xVuMRWPE4EVhxdjFXIrCY00NgY0/vQgJrNpu78zx/wjl3ba/RRqNxhzHmUufclxDYaKNDYKPlGctsCCyWpC7cJ59CHFOO1tplY8yreZ4vJUlyZ7vdfm69lWazeU2e53c65z7U+5m19sMicsg5dzsCG21gCGy0PGOZDYHFkhQCU5lU75aVpuni3Nzce5IkOba6unrg5MmTb/SaRWDVRYbAqmPtu5KvZHzn9xl/8YzhX/Z9gFU0llAqAn2hZay1Txpj7krT9HhfYFs9QrzMOfdFbmCjDQyBjZbnKGdDYKOkWc+5ENgYct23b9+upaWlM91ud7n/fteLMzMzV8/Pz59ab8da+3KWZR/rf4jjB0mS3MaHOEYfFgIbPdNRzYjARkWyvvMgsDFk22q1rs+y7AERWRORTETucs49bq19yhhza//R4vuTJHmo/zH6R51zn9+uVT6FWDxEBFacXehKBBaacPzzI7D4MxQEVjxEBFacXehKBBaacPzzI7D4M0RgJTJEYCXgBS5FYIEB12B6BFaDELmBFQ8RgRVnF7oSgYUmHP/8CCz+DLmBlcgQgZWAF7gUgQUGXIPpEVgNQuQGVjxEBFacXehKBBaacPzzI7D4M+QGViJDBFYCXuBSBBYYcA2mR2A1CJEbWPEQEVhxdqErEVhowvHPj8Diz5AbWIkMEVgJeIFLEVhgwDWYHoHVIERuYMVDRGDF2YWuRGChCcc/PwKLP0NuYCUyRGAl4AUuRWCBAddgegRWgxC5gRUPEYEVZxe6EoGFJhz//Ags/gy5gZXIEIGVgBe4FIEFBlyD6RFYDULkBlY8RARWnF3oSgQWmnD88yOw+DPkBlYiQwRWAl7gUgQWGHANpkdgNQiRG1jxEBFYcXahKxFYaMLxz4/A4s+QG1iJDBFYCXiBSxFYYMA1mB6B1SBEbmDFQ0RgxdmFrkRgoQnHPz8Ciz9DbmAlMkRgJeAFLkVggQHXYHoEVoMQuYEVDxGBFWcXuhKBhSYc//wILP4MuYFtyhAp1eBQiwgCq0eOIXeBwELS3WbuZrPZzPP8YRG5XERWsiz7dLfbfX7jcGvtORGZ7//suHPu1u1a5QZ2PhkENoZDHWBJBBYAas2mRGBjCHTv3r3vzLJsR5qmqbW2JSJPOucamwT2U+fclcO0h8AQ2DDnJLYxCCy2xKrvF4FVz3zzisZauy6rfP3FDT8b2CECQ2ADD0mEAxBYhKFV3DICqxj45uWazeYH8zw/4pz7wKYb2LIx5tU8z5eSJLmz3W4/xyPE4cLiEeJwnLSPQmDaExp/fwhsjBn0HiWurq7+U5ZlN5w4caK7sZVms7k7TdPFubm59yRJcmx1dfXAyZMn39iqXW5g3MDGeIyDLY3AgqGtzcQIbExRzs3N/VqSJN/Psuxz3W73Hy/UhrX2SWPMXWmaHkdggwPjBjaYUQwjEFgMKY23RwQ2Hv5TjUbjWJIkT6Rp+lebW9i3b9+upaWlM91ud7l3E8vz/MWZmZmr5+fnTyGwwYEhsMGMYhiBwGJIabw9IrAx8G80GoeNMd/Z8DF5WVlZOfiWt7zlOefcgVardX2WZQ+IyJqIZCJyl3Pucd4DGy4sBDYcJ+2jEJj2hMbfHwIbfwalO+A9MN4DK32IFE6AwBSGoqwlBKYskCLtIDAEVuTcaK9BYNoTGn9/CGz8GZTuAIEhsNKHSOEECExhKMpaQmDKAinSjhaB8d5TkfSo2Y4AAuNsDCKAwAYRiuB1BBZBSLToTQCBeSObuAIEVoPIEVgNQmQLbyKAwDgUgwggsEGEIngdgUUQEi16E0Bg3sgmrgCB1SByBFaDENkCNzDOgDcBBOaNTF8BAtOXCR2VJ8ANrDzDus+AwGqQMAKrQYhsgRsYZ8CbAALzRqavAIHpy4SOyhPgBlaeYd1nQGA1SBiB1SBEtsANjDPgTQCBeSPTV4DA9GVCR+UJcAMrz7DuMyCwGiSMwGoQIlvgBsYZ8CaAwLyR6SsIJTC+Gkpf1usd+d5OYszSd48h07p4xvC7MiTggnMTSkFwmsoQmKY0qunF95c7AiuXCwIrxy9UNQILRbbCeRFYhbCVLIXAqg0CgVXLe9jVENiwpBSPQ2CKwwnUGgILBHabaRFYtbyHXQ2BDUtK8TgEpjicQK0hsEBgEVi1YEuuhsBKAtRQjsA0pFBtDwisWt7cwKrlPexqCGxYUorHITDF4QRqDYEFAssNrFqwJVdDYCUBFilvNBqHjTFfEZFERP7COffgxnmstdeJyMMiskNEHnHO3X2hdRBYkRTirkFg1ebHDaxa3sOuhsCGJTWicQcPHpxeXFx8dWpq6veSJDl97ty540mS/E673f6f9SWstS9lWfbxbrfbtta+mCTJJ9rt9o+3awGBjSiciKZBYNWGhcCq5T3saghsWFIjGtdqta5fW1v7dKfT+cPelNbar+R5frzT6fxd75+bzebuPM+fcM5d2/vnRqNxhzHmUufclxDYiEKowTQIrNoQEVi1vIddDYENS2pE45rN5keyLPvdTqfzp32BfTrP86lOp/PlvsCuyfP8Tufch/qvf1hEDjnnbkdgIwqhBtMgsGpDRGDV8h52NQQ2LKkRjQshsBG1xjQQgAAEoiKAwCqOa5tHiD/sdDp/e4FHiJc5575YcassBwEIQEA1AQRWcTz9D3HMT01NHex/iOOHSZJcv+lDHC9nWfax/oc4fpAkyW0X+hBHxVtgOQhAAAIqCCCwMcTQbDZ/P8uyLxtjeh+jv8c594C19iljzK1pmi7Ozc29P0mSh/ofo3/UOff5MbTJkhCAAARUE0BgquOhOQhAAAIQ2I4AAuNsQAACEIBAlAQQWJSx/X/Tg77RI+Ktnde6tfaciMz3f3jcOXdr7HtrNptH8zw/JCJPO+du6u3H9xtYtDPYZo8LIvKGiOQisuicu0H7PuhPLwEEpjebC3Y2zDd6RLq1N7Vtrf2pc+7Kuuynt49Wq3VwdXX1siRJbtkgMK9vYNHOY5s9LqysrLQWFhbOau+f/vQTQGD6M9qyw0Hf6BHptrZsu44CW5dYlmVHegIr8g0sMWTck9j6Hvu3TAQWQ3CR9IjAIglqc5uD/kJ0pNvaTmDLxphX8zxfSpLkzna7/Vwd9rfxl3uz2fT+BpYYGGwhsNeMMafyPF8zxtyTpunRGPZBjzoJIDCduQzsapIE1rud9P96wXuSJDm2urp64OTJk733UaL+M4kCW8/SWnuViDy7trZ2+MSJE92og6T5sRFAYGNDX27hQd/oUW52vdXW2ieNMXelaXpcb5fDdbZJYFt9iXP038Cy+Qa2kYy19l5jzPPcwoY7L4x6MwEEFumpGOYbPSLd2nlt79u3b9fS0tKZbre73H+f6MWZmZmr5+fnT8W+vy0er9XuG1g27nH//v07l5eXkzRNT8/Ozl4+PT39QpZlN3e73fVPmMYeKf1XTACBVQx8lE1CVn0AAAFrSURBVMtt9Y0eo5xfw1y9m2aWZQ+IyJqIZCJyl3PucQ29lenBWntMRHr/4dKdInIqy7IbReSiOn0Dy+Y9JknywSzLev+hVsnzvPctNPd1Op37y3CkdrIJILDJzp/dQwACEIiWAAKLNjoahwAEIDDZBBDYZOfP7iEAAQhESwCBRRsdjUMAAhCYbAIIbLLzZ/cQgAAEoiWAwKKNjsYhAAEITDYBBDbZ+bN7CEAAAtESQGDRRkfjEIAABCabAAKb7PzZPQQgAIFoCSCwaKOjcQhAAAKTTQCBTXb+7B4CEIBAtAQQWLTR0TgEIACBySaAwCY7f3YPAQhAIFoCCCza6GgcAhCAwGQTQGCTnT+7hwAEIBAtAQQWbXQ0DgEIQGCyCSCwyc6f3UMAAhCIlgACizY6GocABCAw2QQQ2GTnz+4hAAEIREsAgUUbHY1DAAIQmGwCCGyy82f3EIAABKIlgMCijY7GIQABCEw2gf8DfGT05NKMSmcAAAAASUVORK5CYII=\" width=\"432\">"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# イジングモデルのアニメーション\n",
"%matplotlib nbagg\n",
"from matplotlib import animation\n",
"L = 20\n",
"theta = 0.90\n",
"init = [1 - 2*np.random.randint(0, 2) for i in range(L*L)]\n",
"samples = np.array(mcmc_gibbs_ising(init, 10000, 2000, theta))\n",
"\n",
"# animation \n",
"fig = plt.figure()\n",
"plt.axes().set_aspect('equal')\n",
"ims = []\n",
"for sample in samples[::20]:\n",
" sample = np.reshape(sample, [L, L])\n",
" ims.append((plt.pcolor(np.arange(L), np.arange(L), sample, cmap=plt.cm.Blues),))\n",
"\n",
"ani = animation.ArtistAnimation(fig, ims, interval=5, blit=True)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.3"
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"toc_cell": false,
"toc_position": {},
"toc_section_display": "block",
"toc_window_display": false
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment