EMF_RBM_Energy_Densities_Ominglot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  EMF RBM Energy Densities\n",
    "\n",
    "Look at different ways of computing the\n",
    "\n",
    "- energy density\n",
    "- free energy density\n",
    "- entropy density\n",
    "\n",
    "We can\n",
    "\n",
    "- use the TAP S and MF Energy *U_naive   \n",
    "    - which takes an equlibration step first  \n",
    "    - S(mv, mh), U_naive(mv, mh)\n",
    "- use a simple h = means_hiddens() step\n",
    "    - S(v, h), E_MF(v, h)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T23:33:50.367876",
     "start_time": "2016-11-08T23:33:50.361534"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import h5py\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline \n",
    "\n",
    "from sklearn import linear_model, datasets, metrics, preprocessing \n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "import psutil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T23:33:50.512461",
     "start_time": "2016-11-08T23:33:50.507597"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'ignore'}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.datasets import fetch_mldata\n",
    "from sklearn.utils.validation import assert_all_finite\n",
    "from scipy.sparse import csc_matrix, csr_matrix, lil_matrix\n",
    "from sklearn.preprocessing import Binarizer\n",
    "np.seterr(all='warn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import rbm_datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_x, train_t, train_char, test_x, test_t, test_char = rbm_datasets.load_omniglot_iwae()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "X = binarize(train_x, threshold=0.001, copy=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.io import loadmat\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "def show_image(image): \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    imgplot = ax.imshow(image, cmap=mpl.cm.Greys)\n",
    "    imgplot.set_interpolation('nearest')\n",
    "    ax.xaxis.set_ticks_position('top')\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFdCAYAAADSR9wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEUBJREFUeJzt3X+sZGV9x/H3pxWDruESsdlNikF0JcE/NL1rpRRX12KC\n0gT9CzNqCPUPQtDG3KTRmJhC4Q8jxF5azTamSVGDmpBYijXAag3+IBYwIEYES6SLKD9WkAQSfijC\n0z9mts5edu/uzJ2z35lz36/kJMyZh3u+z32GD+c+c55z0lpDklTjj6oLkKTNzBCWpEKGsCQVMoQl\nqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQnMZwkk+nGRvkmeS3JLkz6trmoUkFyd5Yc12d3Vd\n00iyM8nXkzw46sc5B2lzaZKHkjyd5FtJtlfUOo3D9S/JVQcZy+ur6j1SST6R5LYkTybZl+TaJKcc\npN1Cjt2R9G/exm7uQjjJ+4DPABcDfwb8GNiT5FWlhc3OXcBWYNtoe2ttOVPbAtwJXAS86AYkST4O\nfAS4AHgL8BTDcXzp0SxyA9bt38gNHDiWg6NT2obsBD4LnAa8EzgG+GaSl+1vsOBjd9j+jczP2LXW\n5moDbgH+aex1gF8BH6uubQZ9uxi4o7qODvr1AnDOmn0PAStjr48DngHOra53Rv27Cvj36tpm0LdX\njfr31p6O3cH6N1djN1dnwkmOAXYA396/rw1/a/8FnF5V14y9fvQn7n1Jrk7y6uqCZi3JyQzPLsbH\n8UngVvozjgC7Rn/y/izJ7iSvrC5oCsczPNN/HHo5dgf0b8zcjN1chTDD/2v9MbBvzf59DD8Yi+4W\n4HzgLOBC4GTge0m2VBbVgW0MP/h9HUcY/jl7HvBXwMeAtwPXJ0lpVRMY1XolcHNrbf93E70Zu0P0\nD+Zs7F5ScdDNqrW2Z+zlXUluA34BnMvwTyQtiNbaNWMvf5rkJ8B9wC7gppKiJrcbeANwRnUhHTlo\n/+Zt7ObtTPgx4HmGE+bjtgKPHP1yutVaewK4F1iIb54n8AjDufxNMY4ArbW9DD+/CzGWST4HnA3s\naq09PPZWL8Zunf69SPXYzVUIt9aeA24Hzty/b/QnwpnAD6rq6kqSVzAc+HU/JItm9KF+hAPH8TiG\n31j3bhwBkpwInMACjOUooN4DvKO19sD4e30Yu/X6d4j2pWM3j9MR/wh8IcntwG3ACvBy4AuVRc1C\nkiuA/2Q4BfGnwD8AzwFfraxrGqN57O0Mz5oAXpvkTcDjrbVfMpyL+2SSnwP3A5cxvMrluoJyJ7Ze\n/0bbxcDXGAbWduDTDP+q2fPinzY/kuxmeDnWOcBTSfaf8T7RWnt29M8LO3aH699oXOdr7KovzzjE\nZSUXMRz8Z4D/Bt5cXdOM+vVVhh/mZ4AHgK8AJ1fXNWVf3s7w0p/n12z/NtbmEoaXOz3N8AO+vbru\nWfQPOBa4keF/xM8C/wv8C/An1XUfQb8O1qfngfPWtFvIsTtc/+Zx7DIqTJJUYK7mhCVpszGEJamQ\nISxJhQxhSSpkCEtSIUNYkgqVL9ZIcgLDG9rcz/C6PUladMcCrwH2tNZ+s17DzkI4yYeBv2N456Uf\nA3/bWvvhQZqeBXy5qzokqdAHGC7KOqROQnjs6RgX8Ielx3uSnNJae2xN8/sBrr76ak499dQD3lhZ\nWWF1dbWLEsv1uW/Q7/7Zt8V1tPp3zz338MEPfhBG+baers6EV4DPt9a+BJDkQuCvgQ8Bl69p+yzA\nqaeeyvLy8gFvLC0tvWhfX/S5b9Dv/tm3xVXQv8NOsc78i7lN8nQMSZqJLq6O6PvTMSRpZrxETZIK\ndTEnPNXTMVZWVlhaWjpg30knnTTz4ubFYLAIT0efXp/7Z98W1zz2r5NbWSa5Bbi1tfbR0eswvH/u\nP7fWrljTdhm4/fbbb+/1FwKSNo877riDHTt2AOxord2xXtuuro7o7dMxJGmWOgnh1to1SV4FXMpw\nGuJO4KzW2qNdHE+SFlVnK+Zaa7sZPnJaknQIXh0hSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAk\nFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iS\nChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJ\nhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSCs08hJNcnOSF\nNdvdsz6OJPXBSzr6uXcBZwIZvf59R8eRpIXWVQj/vrX2aEc/W5J6o6s54dcneTDJfUmuTvLqjo4j\nSQutixC+BTgfOAu4EDgZ+F6SLR0cS5IW2synI1pre8Ze3pXkNuAXwLnAVbM+niQtsq7mhP9fa+2J\nJPcC29drt7KywtLS0gH7BoMBg8Ggy/IkqVTnIZzkFQwD+EvrtVtdXWV5ebnrciRprnRxnfAVSd6W\n5KQkfwlcCzwHfHXWx5KkRdfFmfCJwFeAE4BHgZuBv2it/aaDY0nSQuviizkncSXpCHnvCEkqZAhL\nUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgp1fitLaZElOXyjBdVaqy5B\neCYsSaUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSrksmUdNX1eAryIJhkP\nlzh3xzNhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1Ihly3rRVxerLVc\n4twdz4QlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYVctrzAXF7cvUVb\ngutnYvFMfCacZGeSryd5MMkLSc45SJtLkzyU5Okk30qyfTblSlK/TDMdsQW4E7gIeNFpQpKPAx8B\nLgDeAjwF7Eny0g3UKUm9NPF0RGvtRuBGgBz8b5+PApe11r4xanMesA94L3DN9KVKUv/M9Iu5JCcD\n24Bv79/XWnsSuBU4fZbHkqQ+mPXVEdsYTlHsW7N/3+g9SdKYubk6YmVlhaWlpQP2DQYDBoNBUUWS\n1L1Zh/AjQICtHHg2vBX40Xr/4urqKsvLyzMuR5Lm20ynI1prexkG8Zn79yU5DjgN+MEsjyVJfTDx\nmXCSLcB2hme8AK9N8ibg8dbaL4ErgU8m+TlwP3AZ8CvguplULEk9Ms10xJuBmxh+AdeAz4z2fxH4\nUGvt8iQvBz4PHA98H3h3a+13M6hXknplmuuEv8thpjFaa5cAl0xX0ua2aMtOF21Zb99NMh5dfdZ8\nMvNkvIGPJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJanQ3NxPuM/mYSny\noi0PnYffGSze720Si7bEGfo5Hp4JS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhL\nUiFDWJIKuWx5CvOypFbS4vNMWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZ\nwpJUyBCWpELeO2IKkz52u6t7TfTx8d/zZpKxczw0Dc+EJamQISxJhQxhSSpkCEtSIUNYkgoZwpJU\nyBCWpEKGsCQVMoQlqZAhLEmFJg7hJDuTfD3Jg0leSHLOmvevGu0f366fXcnS0dNaO+Jt0SQ54q0r\nk/x+F/F3fCSmORPeAtwJXAQc6rdyA7AV2DbaBlNVJ0k9N/ENfFprNwI3AuTQ/4v8bWvt0Y0UJkmb\nQVdzwruS7EvysyS7k7yyo+NI0kLr4laWNwBfA/YCrwM+BVyf5PTW10kdSZrSzEO4tXbN2MufJvkJ\ncB+wC7hp1seTpEXW+U3dW2t7kzwGbGedEF5ZWWFpaemAfYPBgMHA7/Qk9VfnIZzkROAE4OH12q2u\nrrK8vNx1OZI0VyYO4SRbGJ7V7r8y4rVJ3gQ8PtouZjgn/Mio3aeBe4E9syhYkvpkmjPhNzOcVmij\n7TOj/V9keO3wG4HzgOOBhxiG79+31p7bcLWS1DPTXCf8Xda/tO1d05cjSZuLT1teYH1+ivM81DCp\neXgyc5dLjI/UIo5dJW/gI0mFDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkq\n5LLlo2CSZZzzsOx0HpbfLqJFG2fNB8+EJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQV\nMoQlqZAhLEmFXLY8ZxZt6es81KDpzMMTn1327pmwJJUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKS\nVMgQlqRChrAkFTKEJamQy5YX2KItcVb35mEZ8DzUsEg8E5akQoawJBUyhCWpkCEsSYUMYUkqZAhL\nUiFDWJIKGcKSVMgQlqRChrAkFXLZ8iYxD0tJXTr9Bz7pWPtNdCac5BNJbkvyZJJ9Sa5NcspB2l2a\n5KEkTyf5VpLtsytZkvpj0umIncBngdOAdwLHAN9M8rL9DZJ8HPgIcAHwFuApYE+Sl86kYknqkYmm\nI1prZ4+/TnI+8GtgB3DzaPdHgctaa98YtTkP2Ae8F7hmg/VKUq9s9Iu544EGPA6Q5GRgG/Dt/Q1a\na08CtwKnb/BYktQ7U4dwht8AXAnc3Fq7e7R7G8NQ3rem+b7Re5KkMRu5OmI38AbgjFkUsrKywtLS\n0gH7BoMBg8FgFj9ekubSVCGc5HPA2cDO1trDY289AgTYyoFnw1uBH633M1dXV1leXp6mHElaWBNP\nR4wC+D3AO1prD4y/11rbyzCIzxxrfxzDqyl+sLFSJal/JjoTTrIbGADnAE8l2Tp664nW2rOjf74S\n+GSSnwP3A5cBvwKum0nFktQjk05HXMjwi7fvrNn/N8CXAFprlyd5OfB5hldPfB94d2vtdxsrVZL6\nZ9LrhI9o+qK1dglwyRT1qMdcJts9f8eLxxv4SFIhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUM\nYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKG\nsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFD\nWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEIThXCSTyS5LcmT\nSfYluTbJKWvaXJXkhTXb9bMtW5L6YdIz4Z3AZ4HTgHcCxwDfTPKyNe1uALYC20bbYIN1SlIvvWSS\nxq21s8dfJzkf+DWwA7h57K3fttYe3XB1ktRzG50TPh5owONr9u8aTVf8LMnuJK/c4HEkqZcmOhMe\nlyTAlcDNrbW7x966AfgasBd4HfAp4Pokp7fW2kaKlaS+mTqEgd3AG4Azxne21q4Ze/nTJD8B7gN2\nATcd6oetrKywtLR0wL7BYMBg4HSypP6aKoSTfA44G9jZWnt4vbattb1JHgO2s04Ir66usry8PE05\nkrSwJg7hUQC/B3h7a+2BI2h/InACsG5YS9JmNOl1wruBDwDvB55KsnW0HTt6f0uSy5OcluSkJGcC\n/wHcC+yZdfGStOgmvTriQuA44DvAQ2PbuaP3nwfeCFwH/A/wr8APgbe11p6bQb2S1CuTXie8bmi3\n1p4F3rWhiiRpE/HeEZJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoaw\nJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhSZ62nJHjgW45557\nquuQpJkYy7NjD9c2rbVuqzlcAcn7gS+XFiFJ3fhAa+0r6zWYhxA+ATgLuB94trQYSZqNY4HXAHta\na79Zr2F5CEvSZuYXc5JUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1Kh/wMMMzPnhb48HAAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f769b2dcf90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFdCAYAAADSR9wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEPZJREFUeJzt3W+sZHV9x/H3pwUDruESodlNikF0JcEHkt61UqqrazFB\naYI+wowaQn1ACNKYmzQaElIoPDBC6KWVbEOaFDWoCYmlWAOs1uAfYgEDYvhbAl1E+bOCJJDwRxF+\nfTCz6dzL7t2duXP2OzP3/Uom4Z457HwP5+6bc889ZyatNSRJNf6oegBJ2siMsCQVMsKSVMgIS1Ih\nIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFZrKCCf5XJLdSV5OcnuSP6+eaRKSXJzk9VWPB6rnGkeS\n7Um+k+SJwXacuY91Lk3yZJKXknw/ydaKWcdxoO1Lcu0+9uVNVfMerCQXJrkzyQtJ9iS5IcmJ+1hv\nJvfdwWzftO27qYtwkk8CVwIXA38G/ALYleTY0sEm5z5gM7Bl8PhA7Thj2wTcA5wPvOENSJJ8EbgA\nOBd4H/Ai/f34pkM55DqsuX0DN7NyX/YOzWjrsh34CnAK8BHgcOB7SY7cu8KM77sDbt/A9Oy71tpU\nPYDbgX8a+jrAr4EvVM82gW27GLi7eo4Otut14MxVy54Eloa+Pgp4GTiret4Jbd+1wL9XzzaBbTt2\nsH0fmNN9t6/tm6p9N1VHwkkOB7YBP9i7rPX/q/0XcGrVXBP2rsGPuI8muS7J26oHmrQkJ9A/uhje\njy8AdzA/+xFgx+BH3oeS7Ezy1uqBxnA0/SP952Au992K7RsyNftuqiJM//9afwzsWbV8D/1vjFl3\nO3AOcDpwHnAC8OMkmyqH6sAW+t/487ofof/j7NnAXwFfAD4E3JQkpVONYDDrVcBtrbW9v5uYm323\nn+2DKdt3h1W86EbVWts19OV9Se4EfgmcRf9HJM2I1tr1Q1/en+Re4FFgB3BryVCj2wm8G3h/9SAd\n2ef2Tdu+m7Yj4WeB1+ifMB+2GXj60I/Trdba88DDwEz85nkET9M/l78h9iNAa203/e/fmdiXSa4G\nzgB2tNaeGnpqLvbdGtv3BtX7bqoi3Fp7FbgLOG3vssGPCKcBP62aqytJ3kJ/x6/5TTJrBt/UT7Ny\nPx5F/zfWc7cfAZIcBxzDDOzLQaA+Dny4tfb48HPzsO/W2r79rF+676bxdMQ/Al9NchdwJ7AEvBn4\nauVQk5DkCuA/6Z+C+FPgH4BXgW9VzjWOwXnsrfSPmgDekeRk4LnW2q/on4u7KMkjwGPAZfSvcrmx\nYNyRrbV9g8fFwLfpB2sr8GX6P9XseuOfNj2S7KR/OdaZwItJ9h7xPt9ae2XwzzO77w60fYP9Ol37\nrvryjP1cVnI+/Z3/MvDfwHurZ5rQdn2L/jfzy8DjwDeBE6rnGnNbPkT/0p/XVj3+bWidS+hf7vQS\n/W/wrdVzT2L7gCOAW+j/JX4F+F/gX4A/qZ77ILZrX9v0GnD2qvVmct8daPumcd9lMJgkqcBUnROW\npI3GCEtSISMsSYWMsCQVMsKSVMgIS1Kh8ps1khxD/w1tHqN/3Z4kzbojgLcDu1prv11rxc4inORz\nwN/Rf+elXwB/21r72T5WPR34RldzSFKhT9O/KWu/Oonw0KdjnMv/33q8K8mJrbVnV63+GMB1113H\nSSedtOKJpaUllpeXuxix3DxvG8z39rlts+tQbd+DDz7IZz7zGRj0bS1dHQkvAde01r4OkOQ84K+B\nzwKXr1r3FYCTTjqJxcXFFU8sLCy8Ydm8mOdtg/nePrdtdhVs3wFPsU78F3Mb5NMxJGkiurg6Yt4/\nHUOSJsZL1CSpUBfnhMf6dIylpSUWFhZWLDv++OMnPty06PVm4dPRxzfP2+e2za5p3L5O3soyye3A\nHa21zw++Dv33z/3n1toVq9ZdBO6666675voXApI2jrvvvptt27YBbGut3b3Wul1dHTG3n44hSZPU\nSYRba9cnORa4lP5piHuA01trz3TxepI0qzq7Y661tpP+R05LkvbDqyMkqZARlqRCRliSChlhSSpk\nhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGW\npEJGWJIKGWFJKmSEJalQZ5+2rOmSpHoEjam1Vj2COuSRsCQVMsKSVMgIS1IhIyxJhYywJBUywpJU\nyAhLUiEjLEmFjLAkFTLCklTI25alKdfVLefeDj0dPBKWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZ\nYUkqZIQlqZARlqRCRliSChlhSSo08QgnuTjJ66seD0z6dSStT5KDfqg7Xb2Bz33AacDevfeHjl5H\nkmZaVxH+Q2vtmY7+bEmaG12dE35XkieSPJrkuiRv6+h1JGmmdRHh24FzgNOB84ATgB8n2dTBa0nS\nTJv46YjW2q6hL+9LcifwS+As4NpJv54kzbLOP1mjtfZ8koeBrWutt7S0xMLCwoplvV6PXq/X5XiS\nVKrzCCd5C/0Af32t9ZaXl1lcXOx6HEmaKl1cJ3xFkg8mOT7JXwI3AK8C35r0a0nSrOviSPg44JvA\nMcAzwG3AX7TWftvBa0nSTOviF3OexJWkg+RH3s+wabid1I9N79407Gd1xzfwkaRCRliSChlhSSpk\nhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgp527IOmWm4/XYWb7MeZeau/ht39efO4v6YNI+E\nJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpELetjxlvLVX6zENtziPossZ\nZuX72CNhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZC3LQ+McvvkqLdD\nTsPtobNyC6cOnVm7xXlUXf6dniSPhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQl\nqZARlqRC3rY84C2c3fPW6e5Nw37WaEY+Ek6yPcl3kjyR5PUkZ+5jnUuTPJnkpSTfT7J1MuNK0nwZ\n53TEJuAe4HzgDYc2Sb4IXACcC7wPeBHYleRN65hTkubSyKcjWmu3ALcAZN8/+3weuKy19t3BOmcD\ne4BPANePP6okzZ+J/mIuyQnAFuAHe5e11l4A7gBOneRrSdI8mPTVEVvon6LYs2r5nsFzkqQhU3N1\nxNLSEgsLCyuW9Xo9er1e0USS1L1JR/hpIMBmVh4NbwZ+vta/uLy8zOLi4oTHkaTpNtHTEa213fRD\nfNreZUmOAk4BfjrJ15KkeTDykXCSTcBW+ke8AO9IcjLwXGvtV8BVwEVJHgEeAy4Dfg3cOJGJJWmO\njHM64r3ArfR/AdeAKwfLvwZ8trV2eZI3A9cARwM/AT7WWvv9BOaVpLkyznXCP+IApzFaa5cAl4w3\nkuaVt9RqtWm5lb3yk5l9Ax9JKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWp\n0NS8n/AsmcVbLaVDZVr+foyicmaPhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQl\nqZARlqRC3rY8w2bx9lBJK3kkLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYyw\nJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUaOcJJ\ntif5TpInkrye5MxVz187WD78uGlyI0vS/BjnSHgTcA9wPtD2s87NwGZgy+DRG2s6SZpzh436L7TW\nbgFuAUiS/az2u9baM+sZTJI2gq7OCe9IsifJQ0l2JnlrR68jSTNt5CPhg3Az8G1gN/BO4EvATUlO\nba3t7/SFJG1IE49wa+36oS/vT3Iv8CiwA7h10q8nSbOsiyPhFVpru5M8C2xljQgvLS2xsLCwYlmv\n16PX83d6kuZX5xFOchxwDPDUWustLy+zuLjY9TiSNFVGjnCSTfSPavdeGfGOJCcDzw0eF9M/J/z0\nYL0vAw8DuyYxsCTNk3GOhN9L/7RCGzyuHCz/Gv1rh98DnA0cDTxJP75/31p7dd3TStKcGec64R+x\n9qVtHx1/HEnaWHzvCEkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZak\nQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZ\nYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQl\nqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKjRShJNcmOTOJC8k2ZPkhiQn7mO9S5M8meSl\nJN9PsnVyI0vS/Bj1SHg78BXgFOAjwOHA95IcuXeFJF8ELgDOBd4HvAjsSvKmiUwsSXPksFFWbq2d\nMfx1knOA3wDbgNsGiz8PXNZa++5gnbOBPcAngOvXOa8kzZX1nhM+GmjAcwBJTgC2AD/Yu0Jr7QXg\nDuDUdb6WJM2dsSOcJMBVwG2ttQcGi7fQj/KeVavvGTwnSRoy0umIVXYC7wbeP4lBlpaWWFhYWLGs\n1+vR6/Um8cdL0lQaK8JJrgbOALa31p4aeuppIMBmVh4NbwZ+vtafuby8zOLi4jjjSNLMGvl0xCDA\nHwc+3Fp7fPi51tpu+iE+bWj9o+hfTfHT9Y0qSfNnpCPhJDuBHnAm8GKSzYOnnm+tvTL456uAi5I8\nAjwGXAb8GrhxIhNL0hwZ9XTEefR/8fbDVcv/Bvg6QGvt8iRvBq6hf/XET4CPtdZ+v75RJWn+jHqd\n8EGdvmitXQJcMsY8krSh+N4RklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhL\nUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmF\njLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLC\nklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUaKcJJLkxyZ5IXkuxJckOSE1etc22S11c9\nbprs2JI0H0Y9Et4OfAU4BfgIcDjwvSRHrlrvZmAzsGXw6K1zTkmaS4eNsnJr7Yzhr5OcA/wG2Abc\nNvTU71prz6x7Okmac+s9J3w00IDnVi3fMThd8VCSnUneus7XkaS5NNKR8LAkAa4CbmutPTD01M3A\nt4HdwDuBLwE3JTm1tdbWM6wkzZuxIwzsBN4NvH94YWvt+qEv709yL/AosAO4dX9/2NLSEgsLCyuW\n9Xo9ej1PJ0uaX2NFOMnVwBnA9tbaU2ut21rbneRZYCtrRHh5eZnFxcVxxpGkmTVyhAcB/jjwodba\n4wex/nHAMcCasZakjWjU64R3Ap8GPgW8mGTz4HHE4PlNSS5PckqS45OcBvwH8DCwa9LDS9KsG/Xq\niPOAo4AfAk8OPc4aPP8a8B7gRuB/gH8FfgZ8sLX26gTmlaS5Mup1wmtGu7X2CvDRdU0kSRuI7x0h\nSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQV\nMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFRvq05Y4cAfDggw9WzyFJEzHUsyMOtG5aa91O\nc6ABkk8B3ygdQpK68enW2jfXWmEaInwMcDrwGPBK6TCSNBlHAG8HdrXWfrvWiuURlqSNzF/MSVIh\nIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFfo/H3nOHf02K54AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f76e0528fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFdCAYAAADSR9wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEPhJREFUeJzt3W+sZHV9x/H3pwWDruESsdlNikF0JcEHmt61UqqrazFB\naYI+wowaQn1ACNqYmzQaE1MoPDBi7KXVbGOaFDWoCYmlWAMs1uAfYgEDYkSwRLqI8mcFSSDhjyL8\n+mBm69zL7l1m7hy+M+e+X8kkzJnDnd/vnuHNuWfOmUlrDUlSjT+qHoAkbWVGWJIKGWFJKmSEJamQ\nEZakQkZYkgoZYUkqZIQlqZARlqRCRliSCs1lhJN8OMn+JE8luSnJn1ePaRaSXJjkuXW3O6vHNY0k\nu5N8I8n9o3mcdYh1Lk7yQJInk3wryc6KsU7jSPNLcvkhtuU1VeN9oZJ8IsktSR5PciDJVUlOPsR6\nC7ntXsj85m3bzV2Ek7wP+CxwIfBnwI+BfUleWTqw2bkD2A7sGN3eWjucqW0DbgcuAJ73ASRJPg58\nBDgPeDPwBMPt+JIXc5CbsOH8Rq5l7bYcvDhD25TdwOeAU4F3AkcD1yd56cEVFnzbHXF+I/Oz7Vpr\nc3UDbgL+aex+gF8BH6se2wzmdiFwW/U4OpjXc8BZ65Y9AKyM3T8WeAo4u3q8M5rf5cC/V49tBnN7\n5Wh+b+3ptjvU/OZq283VnnCSo4FdwLcPLmvD39p/AadVjWvGXjf6E/eeJFckeVX1gGYtyUkM9y7G\nt+PjwM30ZzsC7Bn9yfuzJHuTvKJ6QFM4juGe/qPQy223Zn5j5mbbzVWEGf5f64+BA+uWH2D4wlh0\nNwHnAmcA5wMnAd9Lsq1yUB3YwfCF39ftCMM/Z88B/gr4GPB24JokKR3VBEZjvQy4sbV28L2J3my7\nw8wP5mzbHVXxpFtVa23f2N07ktwC/AI4m+GfSFoQrbUrx+7+NMlPgHuAPcANJYOa3F7g9cBbqgfS\nkUPOb9623bztCT8CPMvwgPm47cBDL/5wutVaewy4G1iId54n8BDDY/lbYjsCtNb2M3z9LsS2TPJ5\n4ExgT2vtwbGHerHtNpjf81Rvu7mKcGvtGeBW4PSDy0Z/IpwO/KBqXF1J8nKGG37DF8miGb2oH2Lt\ndjyW4TvWvduOAElOAI5nAbblKFDvAd7RWrtv/LE+bLuN5neY9Uu33TwejvhH4ItJbgVuAVaAlwFf\nrBzULCT5DPCfDA9B/CnwD8AzwNcqxzWN0XHsnQz3mgBek+SNwKOttV8yPBb3ySQ/B+4FLmF4lsvV\nBcOd2EbzG90uBL7OMFg7gU8z/Ktm3/N/2vxIspfh6VhnAU8kObjH+1hr7enRPy/stjvS/Ebbdb62\nXfXpGYc5reQChhv/KeC/gTdVj2lG8/oawxfzU8B9wFeBk6rHNeVc3s7w1J9n193+bWydixie7vQk\nwxf4zupxz2J+wDHAdQz/I34a+F/gX4A/qR73C5jXoeb0LHDOuvUWctsdaX7zuO0yGpgkqcBcHROW\npK3GCEtSISMsSYWMsCQVMsKSVMgIS1Kh8os1khzP8ANt7mV43p4kLbpjgFcD+1prv9loxc4inOTD\nwN8x/OSlHwN/21r74SFWPQP4SlfjkKRCH2B4UdZhdRLhsW/HOI8/XHq8L8nJrbVH1q1+L8AVV1zB\nKaecsuaBlZUVVldXuxhiuT7PDfo9P+e2uF6s+d1111188IMfhFHfNtLVnvAK8IXW2pcBkpwP/DXw\nIeDSdes+DXDKKaewvLy85oGlpaXnLeuLPs8N+j0/57a4CuZ3xEOsM39jbot8O4YkzUQXZ0f0/dsx\nJGlmPEVNkgp1cUx4qm/HWFlZYWlpac2yE088ceaDmxeDwSJ8O/r0+jw/57a45nF+nXyUZZKbgJtb\nax8d3Q/Dz8/959baZ9atuwzceuutt/b6DQFJW8dtt93Grl27AHa11m7baN2uzo7o7bdjSNIsdRLh\n1tqVSV4JXMzwMMTtwBmttYe7eD5JWlSdXTHXWtvL8CunJUmH4dkRklTICEtSISMsSYWMsCQVMsKS\nVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1Ih\nIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYyw\nJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBWaeYST\nXJjkuXW3O2f9PJLUB0d19HPvAE4HMrr/+46eR5IWWlcR/n1r7eGOfrYk9UZXx4Rfl+T+JPckuSLJ\nqzp6HklaaF1E+CbgXOAM4HzgJOB7SbZ18FyStNBmfjiitbZv7O4dSW4BfgGcDVw+6+eTpEXW1THh\n/9daeyzJ3cDOjdZbWVlhaWlpzbLBYMBgMOhyeJJUqvMIJ3k5wwB/eaP1VldXWV5e7no4kjRXujhP\n+DNJ3pbkxCR/CVwFPAN8bdbPJUmLros94ROArwLHAw8DNwJ/0Vr7TQfPJUkLrYs35jyIK0kvUOfH\nhKVpJDnySgustVY9BM0JP8BHkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJ\nKuRly9qUvl9e3JV5+L11den0JHPz8m33hCWplBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkq\nZIQlqZARlqRCXra8RczDZbKaL129JrwUeTLuCUtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUy\nwpJUyAhLUiEjLEmFvGx5znh58ZCXvv7Bor0muhxvH18X7glLUiEjLEmFjLAkFTLCklTICEtSISMs\nSYWMsCQVMsKSVMgIS1IhIyxJhbxs+UWwaJedzsOloX4T8B90NeZFe1321cR7wkl2J/lGkvuTPJfk\nrEOsc3GSB5I8meRbSXbOZriS1C/THI7YBtwOXAA873/RST4OfAQ4D3gz8ASwL8lLNjFOSeqliQ9H\ntNauA64DyKH/nvkocElr7Zujdc4BDgDvBa6cfqiS1D8zfWMuyUnADuDbB5e11h4HbgZOm+VzSVIf\nzPrsiB0MD1EcWLf8wOgxSdKYuTk7YmVlhaWlpTXLBoMBg8GgaESS1L1ZR/ghIMB21u4Nbwd+tNG/\nuLq6yvLy8oyHI0nzbaaHI1pr+xmG+PSDy5IcC5wK/GCWzyVJfTDxnnCSbcBOhnu8AK9J8kbg0dba\nL4HLgE8m+TlwL3AJ8Cvg6pmMWJJ6ZJrDEW8CbmD4BlwDPjta/iXgQ621S5O8DPgCcBzwfeDdrbXf\nzWC8ktQr05wn/F2OcBijtXYRcNF0Q1K1SS5nXcTLgDU/fP34AT6SVMoIS1IhIyxJhYywJBUywpJU\nyAhLUiEjLEmFjLAkFTLCklTICEtSobn5POE+6/O35XqJ83yZh9eE23ky7glLUiEjLEmFjLAkFTLC\nklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhbxseYFNcnnoPFzOuohjmIdLcOfh9zYPv4e+ck9Y\nkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQnx2xRSza50zMC38X\n6pp7wpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1KhiSOcZHeSbyS5P8lz\nSc5a9/jlo+Xjt2tmN2R1rbVWftN03B6LZ5o94W3A7cAFwOG2zrXAdmDH6DaYanSS1HMTf4BPa+06\n4DqAHP7TTX7bWnt4MwOTpK2gq2PCe5IcSPKzJHuTvKKj55GkhdbFR1leC3wd2A+8FvgUcE2S05oH\nlyRpjZlHuLV25djdnyb5CXAPsAe4YdbPJ0mLrPMPdW+t7U/yCLCTDSK8srLC0tLSmmWDwYDBwPf0\nJPVX5xFOcgJwPPDgRuutrq6yvLzc9XAkaa5MHOEk2xju1R48M+I1Sd4IPDq6XcjwmPBDo/U+DdwN\n7JvFgCWpT6bZE34Tw8MKbXT77Gj5lxieO/wG4BzgOOABhvH9+9baM5serST1zDTnCX+XjU9te9f0\nw5GkrcVvW9Zc8mxGbRV+gI8kFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJU\nyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEj\nLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAk\nFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhSaKcJJPJLklyeNJDiS5KsnJh1jv4iQPJHky\nybeS7JzdkCWpPybdE94NfA44FXgncDRwfZKXHlwhyceBjwDnAW8GngD2JXnJTEYsST1y1CQrt9bO\nHL+f5Fzg18Au4MbR4o8Cl7TWvjla5xzgAPBe4MpNjleSemWzx4SPAxrwKECSk4AdwLcPrtBaexy4\nGThtk88lSb0zdYSTBLgMuLG1dudo8Q6GUT6wbvUDo8ckSWMmOhyxzl7g9cBbZjGQlZUVlpaW1iwb\nDAYMBoNZ/HhJmktTRTjJ54Ezgd2ttQfHHnoICLCdtXvD24EfbfQzV1dXWV5enmY4krSwJj4cMQrw\ne4B3tNbuG3+stbafYYhPH1v/WIZnU/xgc0OVpP6ZaE84yV5gAJwFPJFk++ihx1prT4/++TLgk0l+\nDtwLXAL8Crh6JiOWpB6Z9HDE+QzfePvOuuV/A3wZoLV2aZKXAV9gePbE94F3t9Z+t7mhSlL/THqe\n8As6fNFauwi4aIrxSNKW4mdHSFIhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMs\nSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQV\nMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgI\nS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklRooggn+USSW5I8nuRAkquSnLxuncuTPLfu\nds1shy1J/TDpnvBu4HPAqcA7gaOB65O8dN161wLbgR2j22CT45SkXjpqkpVba2eO309yLvBrYBdw\n49hDv22tPbzp0UlSz232mPBxQAMeXbd8z+hwxc+S7E3yik0+jyT10kR7wuOSBLgMuLG1dufYQ9cC\nXwf2A68FPgVck+S01lrbzGAlqW+mjjCwF3g98Jbxha21K8fu/jTJT4B7gD3ADYf7YSsrKywtLa1Z\nNhgMGAw8nCypv6aKcJLPA2cCu1trD260bmttf5JHgJ1sEOHV1VWWl5enGY4kLayJIzwK8HuAt7fW\n7nsB658AHA9sGGtJ2oomPU94L/AB4P3AE0m2j27HjB7fluTSJKcmOTHJ6cB/AHcD+2Y9eEladJOe\nHXE+cCzwHeCBsdvZo8efBd4AXA38D/CvwA+Bt7XWnpnBeCWpVyY9T3jDaLfWngbetakRSdIW4mdH\nSFIhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJ\nhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSoYm+bbkjxwDcdddd1eOQpJkY69kxR1o3rbVu\nR3OkASTvB75SOghJ6sYHWmtf3WiFeYjw8cAZwL3A06WDkaTZOAZ4NbCvtfabjVYsj7AkbWW+MSdJ\nhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTo/wAe+Gxykwh3+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f769912ef50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFdCAYAAADSR9wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEVNJREFUeJzt3X+sZGV9x/H3pxWDruESsdlNikF0JcE/NL1rpVRX12KC\n0gT9i2bUEOofhKCNuUmjMTGFwh9GjL20mm1Mk6IGNSGxFGuAxRr8QSxgQIwIlkgXUX6sIAkk/FCE\np3/MbDv37u7dnblz9jtz9v1KTrJz5mHO97nP8LnnPnOeM2mtIUmq8QfVBUjSscwQlqRChrAkFTKE\nJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEJzGcJJPpxkb5Jnk9ya5E+ra5qFJJckeXHddk91\nXdNIsjPJN5I8NOrHuQdpc1mSh5M8k+RbSbZX1DqNw/UvyVUHGcvrq+o9Ukk+keT2JE8l2Zfk2iSn\nHaTdQo7dkfRv3sZu7kI4yV8BnwUuAf4E+DGwJ8mrSgubnbuBrcC20fa22nKmtgW4C7gYOOAGJEk+\nDnwEuBB4C/A0w3F86dEschM27N/IDawdy8HRKW1TdgKfA84A3gUcB9yU5GX7Gyz42B22fyPzM3at\ntbnagFuBfxx7HOBXwMeqa5tB3y4B7qyuo4N+vQicu27fw8DK2OMTgGeB86rrnVH/rgL+rbq2GfTt\nVaP+va2nY3ew/s3V2M3VmXCS44AdwLf372vDn9p/AmdW1TVjrx/9iXt/kquTvLq6oFlLcirDs4vx\ncXwKuI3+jCPArtGfvD9LsjvJK6sLmsKJDM/0n4Bejt2a/o2Zm7GbqxBm+FvrD4F96/bvY/jGWHS3\nAhcAZwMXAacC30uypbKoDmxj+Mbv6zjC8M/Z84G/AD4GvAO4PklKq5rAqNYrgVtaa/s/m+jN2B2i\nfzBnY/eSioMeq1pre8Ye3p3kduAXwHkM/0TSgmitXTP28KdJfgLcD+wCbi4panK7gTcAb60upCMH\n7d+8jd28nQk/DrzAcMJ83Fbg0aNfTrdaa08C9wEL8cnzBB5lOJd/TIwjQGttL8P370KMZZLPA+cA\nu1prj4w91Yux26B/B6geu7kK4dba88AdwFn7943+RDgL+EFVXV1J8gqGA7/hm2TRjN7Uj7J2HE9g\n+Il178YRIMnJwEkswFiOAuq9wDtbaw+OP9eHsduof4doXzp28zgd8Q/AF5PcAdwOrAAvB75YWdQs\nJPkM8B8MpyD+GPh74Hnga5V1TWM0j72d4VkTwGuTvAl4orX2S4ZzcZ9M8nPgAeByhle5XFdQ7sQ2\n6t9ouwT4OsPA2g58muFfNXsOfLX5kWQ3w8uxzgWeTrL/jPfJ1tpzo38v7Ngdrn+jcZ2vsau+POMQ\nl5VczHDwnwX+C3hzdU0z6tfXGL6ZnwUeBL4KnFpd15R9eQfDS39eWLf961ibSxle7vQMwzf49uq6\nZ9E/4HjgRob/Ez8H/A/wz8AfVdd9BP06WJ9eAM5f124hx+5w/ZvHscuoMElSgbmaE5akY40hLEmF\nDGFJKmQIS1IhQ1iSChnCklSofLFGkpMY3tDmAYbX7UnSojseeA2wp7X2m40adhbCST4M/C3DOy/9\nGPib1toPD9L0bOArXdUhSYU+wHBR1iF1EsJj345xIf+/9HhPktNaa4+va/4AwNVXX83pp5++5omV\nlRVWV1e7KLFcn/sG/e6ffVtcR6t/9957Lx/84AdhlG8b6epMeAX4QmvtywBJLgL+EvgQcMW6ts8B\nnH766SwvL695Ymlp6YB9fdHnvkG/+2ffFldB/w47xTrzD+aOkW/HkKSZ6OLqiL5/O4YkzYyXqElS\noS7mhKf6doyVlRWWlpbW7DvllFNmXty8GAwW4dvRp9fn/tm3xTWP/evkVpZJbgVua619dPQ4DO+f\n+0+ttc+sa7sM3HHHHXf0+gMBSceOO++8kx07dgDsaK3duVHbrq6O6O23Y0jSLHUSwq21a5K8CriM\n4TTEXcDZrbXHujieJC2qzlbMtdZ2M/zKaUnSIXh1hCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSChnC\nklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxh\nSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoaw\nJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSo08xBOckmS\nF9dt98z6OJLUBy/p6HXvBs4CMnr8+46OI0kLrasQ/n1r7bGOXluSeqOrOeHXJ3koyf1Jrk7y6o6O\nI0kLrYsQvhW4ADgbuAg4Ffheki0dHEuSFtrMpyNaa3vGHt6d5HbgF8B5wFWzPp4kLbKu5oT/T2vt\nyST3Ads3areyssLS0tKafYPBgMFg0GV5klSq8xBO8gqGAfzljdqtrq6yvLzcdTmSNFe6uE74M0ne\nnuSUJH8OXAs8D3xt1seSpEXXxZnwycBXgZOAx4BbgD9rrf2mg2NJ0kLr4oM5J3El6Qh57whJKmQI\nS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKdX4rS/VbksM30gFaa9Ul\naE54JixJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKuSyZR3Apcjdm+Rn\n3NUS53moQZ4JS1IpQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKuWx5zrhk\neHJdLqmdh/HoqgaXIs8Hz4QlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEs\nSYVctnwUzMPS1670fenrPHzTcVf8tuX5MPGZcJKdSb6R5KEkLyY59yBtLkvycJJnknwryfbZlCtJ\n/TLNdMQW4C7gYuCAX49JPg58BLgQeAvwNLAnyUs3Uack9dLE0xGttRuBGwFy8L9nPgpc3lr75qjN\n+cA+4H3ANdOXKkn9M9MP5pKcCmwDvr1/X2vtKeA24MxZHkuS+mDWV0dsYzhFsW/d/n2j5yRJY+bm\n6oiVlRWWlpbW7BsMBgwGg6KKJKl7sw7hR4EAW1l7NrwV+NFG/+Hq6irLy8szLkeS5ttMpyNaa3sZ\nBvFZ+/clOQE4A/jBLI8lSX0w8Zlwki3AdoZnvACvTfIm4InW2i+BK4FPJvk58ABwOfAr4LqZVCxJ\nPTLNdMSbgZsZfgDXgM+O9n8J+FBr7YokLwe+AJwIfB94T2vtdzOoV5J6ZZrrhL/LYaYxWmuXApdO\nV9L8m4clp5Ny2el8mWQ85uH95jc+d8cb+EhSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJ\nKmQIS1IhQ1iSCs3N/YSrzcPS0Ekt2pJPl75OZ9GWOE9i0nr7ONaeCUtSIUNYkgoZwpJUyBCWpEKG\nsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSCvV62fKiLeGE7pao9nG5pw7U5yXO0M/3vGfCklTIEJak\nQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFer1vSPmhfeD0GYs4j0edOQ8\nE5akQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFZo4hJPsTPKNJA8leTHJ\nueuev2q0f3y7fnYlL54kR7y11o54WzST9G0Rfw6TjPMk26LpapwX9X1xONOcCW8B7gIuBg7V0xuA\nrcC20TaYqjpJ6rmJb+DTWrsRuBEgh/41/dvW2mObKUySjgVdzQnvSrIvyc+S7E7yyo6OI0kLrYtb\nWd4AfB3YC7wO+BRwfZIz2yJN1EjSUTDzEG6tXTP28KdJfgLcD+wCbp718SRpkXV+U/fW2t4kjwPb\n2SCEV1ZWWFpaWrNvMBgwGPiZnqT+6jyEk5wMnAQ8slG71dVVlpeXuy5HkubKxCGcZAvDs9r9V0a8\nNsmbgCdG2yUM54QfHbX7NHAfsGcWBUtSn0xzJvxmhtMKbbR9drT/SwyvHX4jcD5wIvAww/D9u9ba\n85uuVpJ6ZprrhL/Lxpe2vXv6ciTp2LJw37bc1bcRz8vyUL9teWhexmPR9Pk90VfewEeSChnCklTI\nEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVGjhli13tRS570ucNV9cXqz9PBOWpEKG\nsCQVMoQlqZAhLEmFDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBVauGXLk5iHJc6Tvrbmi8uL1TXP\nhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhXq9bHkSXS5PXbSlr/Ow\nzHrRfmbStDwTlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVctmyDuCS\nYenomehMOMknktye5Kkk+5Jcm+S0g7S7LMnDSZ5J8q0k22dXsiT1x6TTETuBzwFnAO8CjgNuSvKy\n/Q2SfBz4CHAh8BbgaWBPkpfOpGJJ6pGJpiNaa+eMP05yAfBrYAdwy2j3R4HLW2vfHLU5H9gHvA+4\nZpP1SlKvbPaDuROBBjwBkORUYBvw7f0NWmtPAbcBZ27yWJLUO1OHcIY3nb0SuKW1ds9o9zaGobxv\nXfN9o+ckSWM2c3XEbuANwFtnUcjKygpLS0tr9g0GAwaDwSxeXpLm0lQhnOTzwDnAztbaI2NPPQoE\n2Mras+GtwI82es3V1VWWl5enKUeSFtbE0xGjAH4v8M7W2oPjz7XW9jIM4rPG2p/A8GqKH2yuVEnq\nn4nOhJPsBgbAucDTSbaOnnqytfbc6N9XAp9M8nPgAeBy4FfAdTOpWJJ6ZNLpiIsYfvD2nXX7/xr4\nMkBr7YokLwe+wPDqie8D72mt/W5zpUpS/0x6nfARTV+01i4FLp2iHkk6pngDH0kqZAhLUiFDWJIK\nGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmF\nDGFJKmQIS1IhQ1iSChnCklTIEJakQoawJBUyhCWpkCEsSYUMYUkqZAhLUiFDWJIKGcKSVMgQlqRC\nhrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKGsCQVMoQlqZAhLEmFDGFJKmQIS1Ih\nQ1iSChnCklRoohBO8okktyd5Ksm+JNcmOW1dm6uSvLhuu362ZUtSP0x6JrwT+BxwBvAu4DjgpiQv\nW9fuBmArsG20DTZZpyT10ksmadxaO2f8cZILgF8DO4Bbxp76bWvtsU1XJ0k9t9k54ROBBjyxbv+u\n0XTFz5LsTvLKTR5HknppojPhcUkCXAnc0lq7Z+ypG4CvA3uB1wGfAq5PcmZrrW2mWEnqm6lDGNgN\nvAF46/jO1to1Yw9/muQnwP3ALuDmQ73YysoKS0tLa/YNBgMGA6eTJfXXVCGc5PPAOcDO1tojG7Vt\nre1N8jiwnQ1CeHV1leXl5WnKkaSFNXEIjwL4vcA7WmsPHkH7k4GTgA3DWpKORZNeJ7wb+ADwfuDp\nJFtH2/Gj57ckuSLJGUlOSXIW8O/AfcCeWRcvSYtu0qsjLgJOAL4DPDy2nTd6/gXgjcB1wH8D/wL8\nEHh7a+35GdQrSb0y6XXCG4Z2a+054N2bqkiSjiHeO0KSChnCklTIEJakQoawJBUyhCWpkCEsSYUM\nYUkqZAhLUiFDWJIKGcKSVMgQlqRChrAkFTKEJamQISxJhQxhSSpkCEtSIUNYkgoZwpJUyBCWpEKG\nsCQVMoQlqdBE37bckeMB7r333uo6JGkmxvLs+MO1TWut22oOV0DyfuArpUVIUjc+0Fr76kYN5iGE\nTwLOBh4AnistRpJm43jgNcCe1tpvNmpYHsKSdCzzgzlJKmQIS1IhQ1iSChnCklTIEJakQoawJBUy\nhCWp0P8CwZe9q9rpNkYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7699073890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFdCAYAAADSR9wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAERxJREFUeJzt3W+sZHV9x/H3pwWDruESsdlNikF0JcEHmt61UqqrazFB\naYI+ohk1hPqAELQxN2k0JqZQeGDE2Eur2cY0KWpQExJLsQZYrcE/xAIGxIhgiXQR5c8KkkDCH0X4\n9cHMtnMvu3d35s7hO3Pu+5VMwpw53vmd+xvfe+6Zc2bSWkOSVOMPqgcgSVuZEZakQkZYkgoZYUkq\nZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJzGeEkH0qyP8nTSW5O8qfVY5qFJBcneX7d7a7qcU0j\nye4kX0/ywGg7zjnEOpcmeTDJU0m+lWRnxVincaTtS3LlIebyuqrxHq0kH09ya5InkhxIck2SUw+x\n3kLO3dFs37zN3dxFOMlfAZ8BLgb+BPgxsC/JK0sHNjt3AtuBHaPbW2uHM7VtwB3ARcALPoAkyceA\nDwMXAG8GnmQ4jy95MQe5CRtu38j1rJ3LwYsztE3ZDXwWOB14J3As8M0kLz24woLP3RG3b2R+5q61\nNlc34GbgH8fuB/gV8NHqsc1g2y4Gbq8eRwfb9TxwzrplDwIrY/ePB54Gzq0e74y270rg36rHNoNt\ne+Vo+97a07k71PbN1dzN1Z5wkmOBXcC3Dy5rw9/afwJnVI1rxl43+hP33iRXJXlV9YBmLckpDPcu\nxufxCeAW+jOPAHtGf/L+LMneJK+oHtAUTmC4p/8Y9HLu1mzfmLmZu7mKMMN/tf4QOLBu+QGGL4xF\ndzNwPnAWcCFwCvC9JNsqB9WBHQxf+H2dRxj+OXse8BfAR4G3A9clSemoJjAa6xXATa21g+9N9Gbu\nDrN9MGdzd0zFk25VrbV9Y3fvTHIr8AvgXIZ/ImlBtNauHrv70yQ/Ae4F9gA3lgxqcnuB1wNvqR5I\nRw65ffM2d/O2J/wo8BzDA+bjtgMPv/jD6VZr7XHgHmAh3nmewMMMj+VviXkEaK3tZ/j6XYi5TPI5\n4GxgT2vtobGHejF3G2zfC1TP3VxFuLX2LHAbcObBZaM/Ec4EflA1rq4keTnDid/wRbJoRi/qh1k7\nj8czfMe6d/MIkOQk4EQWYC5HgXoP8I7W2v3jj/Vh7jbavsOsXzp383g44h+ALyS5DbgVWAFeBnyh\nclCzkOTTwH8wPATxx8DfA88CX60c1zRGx7F3MtxrAnhNkjcCj7XWfsnwWNwnkvwcuA+4jOFZLtcW\nDHdiG23f6HYx8DWGwdoJfIrhXzX7XvjT5keSvQxPxzoHeDLJwT3ex1trz4z+e2Hn7kjbN5rX+Zq7\n6tMzDnNayUUMJ/9p4L+AN1WPaUbb9VWGL+angfuBrwCnVI9rym15O8NTf55bd/vXsXUuYXi601MM\nX+A7q8c9i+0DjgNuYPh/4meA/wH+Gfij6nEfxXYdapueA85bt95Czt2Rtm8e5y6jgUmSCszVMWFJ\n2mqMsCQVMsKSVMgIS1IhIyxJhYywJBUqv1gjyYkMP9DmPobn7UnSojsOeDWwr7X2m41W7CzCST4E\n/C3DT176MfA3rbUfHmLVs4AvdzUOSSr0foYXZR1WJxEe+3aMC/j/S4/3JTm1tfboutXvA7jqqqs4\n7bTT1jywsrLC6upqF0Ms1+dtg35vn9u2uF6s7bv77rv5wAc+AKO+baSrPeEV4POttS8BJLkQ+Evg\ng8Dl69Z9BuC0005jeXl5zQNLS0svWNYXfd426Pf2uW2Lq2D7jniIdeZvzG2Rb8eQpJno4uyIvn87\nhiTNjKeoSVKhLo4JT/XtGCsrKywtLa1ZdvLJJ898cPNiMFiEb0efXp+3z21bXPO4fZ18lGWSm4Fb\nWmsfGd0Pw8/P/afW2qfXrbsM3Hbbbbf1+g0BSVvH7bffzq5duwB2tdZu32jdrs6O6O23Y0jSLHUS\n4dba1UleCVzK8DDEHcBZrbVHung+SVpUnV0x11rby/ArpyVJh+HZEZJUyAhLUiEjLEmFjLAkFTLC\nklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtS\nISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWM\nsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYWMsCQVmnmE\nk1yc5Pl1t7tm/TyS1AfHdPRz7wTOBDK6//uOnkeSFlpXEf59a+2Rjn62JPVGV8eEX5fkgST3Jrkq\nyas6eh5JWmhdRPhm4HzgLOBC4BTge0m2dfBckrTQZn44orW2b+zunUluBX4BnAtcOevnk6RF1tUx\n4f/TWns8yT3Azo3WW1lZYWlpac2ywWDAYDDocniSVKrzCCd5OcMAf2mj9VZXV1leXu56OJI0V7o4\nT/jTSd6W5OQkfw5cAzwLfHXWzyVJi66LPeGTgK8AJwKPADcBf9Za+00HzyVJC62LN+Y8iCtJR8nP\njpCkQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKtT5R1lq8SQ58kp6\ngdZa9RC0gNwTlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgp52fIW4aXI\n3evqd+zl0P3mnrAkFTLCklTICEtSISMsSYWMsCQVMsKSVMgIS1IhIyxJhYywJBUywpJUyMuWXwST\nXM7qJapaz8uh+809YUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQly2P\ndHlp8STrL9q3Invp63QWbZ4n4WX6k5l4TzjJ7iRfT/JAkueTnHOIdS5N8mCSp5J8K8nO2QxXkvpl\nmsMR24A7gIuAF/wzluRjwIeBC4A3A08C+5K8ZBPjlKRemvhwRGvtBuAGgBz6746PAJe11r4xWuc8\n4ADwXuDq6YcqSf0z0zfmkpwC7AC+fXBZa+0J4BbgjFk+lyT1wazPjtjB8BDFgXXLD4wekySNmZuz\nI1ZWVlhaWlqzbDAYMBgMikYkSd2bdYQfBgJsZ+3e8HbgRxv9D1dXV1leXp7xcCRpvs30cERrbT/D\nEJ95cFmS44HTgR/M8rkkqQ8m3hNOsg3YyXCPF+A1Sd4IPNZa+yVwBfCJJD8H7gMuA34FXDuTEUtS\nj0xzOOJNwI0M34BrwGdGy78IfLC1dnmSlwGfB04Avg+8u7X2uxmMV5J6ZZrzhL/LEQ5jtNYuAS6Z\nbkg1+nxp8aT6fCnppHPX1e9iHl5v8/AtzvMyH5X8AB9JKmSEJamQEZakQkZYkgoZYUkqZIQlqZAR\nlqRCRliSChlhSSpkhCWp0Nx8nrC6NQ+Xe87D5d7z8HuY1Dxc4jwJv215Mu4JS1IhIyxJhYywJBUy\nwpJUyAhLUiEjLEmFjLAkFTLCklTICEtSISMsSYW8bHmBdXU5q5eSLq4+X+IM/XxtuicsSYWMsCQV\nMsKSVMgIS1IhIyxJhYywJBUywpJUyAhLUiEjLEmFjLAkFTLCklTIz45YYIv2eRCL9rkGXer79h2t\neXhdVnNPWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSo0cYST7E7y9SQP\nJHk+yTnrHr9ytHz8dt3shlyvtTbRTZOb9Hfc1Vwk6eS2aOZlPvpomj3hbcAdwEXA4X6L1wPbgR2j\n22Cq0UlSz038AT6ttRuAGwBy+H/Sf9tae2QzA5OkraCrY8J7khxI8rMke5O8oqPnkaSF1sVHWV4P\nfA3YD7wW+CRwXZIzmgeBJGmNmUe4tXb12N2fJvkJcC+wB7hx1s8nSYus8w91b63tT/IosJMNIryy\nssLS0tKaZYPBgMHA9/Qk9VfnEU5yEnAi8NBG662urrK8vNz1cCRprkwc4STbGO7VHjwz4jVJ3gg8\nNrpdzPCY8MOj9T4F3APsm8WAJalPptkTfhPDwwptdPvMaPkXGZ47/AbgPOAE4EGG8f271tqzmx6t\nJPXMNOcJf5eNT2171/TDkaStxW9b3iIW8VJZdcszRueDH+AjSYWMsCQVMsKSVMgIS1IhIyxJhYyw\nJBUywpJUyAhLUiEjLEmFjLAkFfKy5RfBJJeHennx1uAlwzrIPWFJKmSEJamQEZakQkZYkgoZYUkq\nZIQlqZARlqRCRliSChlhSSpkhCWpkJctzxkvZ5W2FveEJamQEZakQkZYkgoZYUkqZIQlqZARlqRC\nRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlh\nSSpkhCWpkBGWpEITRTjJx5PcmuSJJAeSXJPk1EOsd2mSB5M8leRbSXbObsiS1B+T7gnvBj4LnA68\nEzgW+GaSlx5cIcnHgA8DFwBvBp4E9iV5yUxGLEk9cswkK7fWzh6/n+R84NfALuCm0eKPAJe11r4x\nWuc84ADwXuDqTY5Xknpls8eETwAa8BhAklOAHcC3D67QWnsCuAU4Y5PPJUm9M3WEkwS4AriptXbX\naPEOhlE+sG71A6PHJEljJjocsc5e4PXAW2YxkJWVFZaWltYsGwwGDAaDWfx4SZpLU0U4yeeAs4Hd\nrbWHxh56GAiwnbV7w9uBH230M1dXV1leXp5mOJK0sCY+HDEK8HuAd7TW7h9/rLW2n2GIzxxb/3iG\nZ1P8YHNDlaT+mWhPOMleYACcAzyZZPvoocdba8+M/vsK4BNJfg7cB1wG/Aq4diYjlqQemfRwxIUM\n33j7zrrlfw18CaC1dnmSlwGfZ3j2xPeBd7fWfre5oUpS/0x6nvBRHb5orV0CXDLFeCRpS/GzIySp\nkBGWpEJGWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJG\nWJIKGWFJKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJ\nKmSEJamQEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQ\nEZakQkZYkgoZYUkqNFGEk3w8ya1JnkhyIMk1SU5dt86VSZ5fd7tutsOWpH6YdE94N/BZ4HTgncCx\nwDeTvHTdetcD24Edo9tgk+OUpF46ZpKVW2tnj99Pcj7wa2AXcNPYQ79trT2y6dFJUs9t9pjwCUAD\nHlu3fM/ocMXPkuxN8opNPo8k9dJEe8LjkgS4AriptXbX2EPXA18D9gOvBT4JXJfkjNZa28xgJalv\npo4wsBd4PfCW8YWttavH7v40yU+Ae4E9wI2H+2ErKyssLS2tWTYYDBgMPJwsqb+minCSzwFnA7tb\naw9ttG5rbX+SR4GdbBDh1dVVlpeXpxmOJC2siSM8CvB7gLe31u4/ivVPAk4ENoy1JG1Fk54nvBd4\nP/A+4Mkk20e340aPb0tyeZLTk5yc5Ezg34F7gH2zHrwkLbpJz464EDge+A7w4Njt3NHjzwFvAK4F\n/hv4F+CHwNtaa8/OYLyS1CuTnie8YbRba88A79rUiCRpC/GzIySpkBGWpEJGWJIKGWFJKmSEJamQ\nEZakQkZYkgoZYUkqZIQlqZARlqRCRliSChlhSSpkhCWpkBGWpEJGWJIKGWFJKmSEJamQEZakQkZY\nkgoZYUkqZIQlqdBE37bckeMA7r777upxSNJMjPXsuCOtm9Zat6M50gCS9wFfLh2EJHXj/a21r2y0\nwjxE+ETgLOA+4JnSwUjSbBwHvBrY11r7zUYrlkdYkrYy35iTpEJGWJIKGWFJKmSEJamQEZakQkZY\nkgoZYUkq9L/QOo2epN7O6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7698f9db10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    show_image(X[i+100].reshape(28,28))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(24345, 784)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Helper methods for simple case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T22:38:39.742812",
     "start_time": "2016-11-08T22:38:39.737075"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.utils.fixes import expit    \n",
    "from sklearn.utils.extmath import safe_sparse_dot\n",
    "\n",
    "def sig_means(x, b, W):\n",
    "    a = safe_sparse_dot(x, W.T) + b\n",
    "    return expit(a, out=a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:56:14.156132",
     "start_time": "2016-10-20T17:56:14.153682"
    }
   },
   "source": [
    "### Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import emf_rbm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T22:42:32.963306",
     "start_time": "2016-11-08T22:42:32.959982"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sswap(total=0, used=0, free=0, percent=0.0, sin=0, sout=0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psutil.swap_memory()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### mean Entropy, Energy, and Free Energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T22:55:48.509468",
     "start_time": "2016-11-08T22:45:04.953105"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EMF_RBM(batch_size=100, decay=0.01, learning_rate=0.005, momentum=0.5,\n",
      "    monitor=True, n_components=256, n_iter=100, neq_steps=3,\n",
      "    random_state=<mtrand.RandomState object at 0x7f76d81289b0>,\n",
      "    sigma=0.001, thresh=1e-08, verbose=True, weight_decay='L1')\n",
      "[EMF_RBM] Iteration 1, pseudo-likelihood = -261.78, time = 9.99s\n",
      "computing TAP Free Energies\n",
      "monitor:  -322.220443321 425.601949856 103.390377782\n",
      "[EMF_RBM] Iteration 2, pseudo-likelihood = -255.23, time = 23.88s\n",
      "computing TAP Free Energies\n",
      "monitor:  -319.513760446 422.970978189 103.523571456\n",
      "[EMF_RBM] Iteration 3, pseudo-likelihood = -236.98, time = 23.89s\n",
      "computing TAP Free Energies\n",
      "monitor:  -307.996179814 396.015426782 88.7889970446\n",
      "[EMF_RBM] Iteration 4, pseudo-likelihood = -221.07, time = 23.79s\n",
      "computing TAP Free Energies\n",
      "monitor:  -297.705911069 367.702882371 71.2889738607\n",
      "[EMF_RBM] Iteration 5, pseudo-likelihood = -212.67, time = 24.18s\n",
      "computing TAP Free Energies\n",
      "monitor:  -292.385953058 347.00263872 56.5277631836\n",
      "[EMF_RBM] Iteration 6, pseudo-likelihood = -201.19, time = 23.94s\n",
      "computing TAP Free Energies\n",
      "monitor:  -290.412873794 329.523787112 41.4522623482\n",
      "[EMF_RBM] Iteration 7, pseudo-likelihood = -198.25, time = 23.79s\n",
      "computing TAP Free Energies\n",
      "monitor:  -288.176431463 318.111617682 32.6591872351\n",
      "[EMF_RBM] Iteration 8, pseudo-likelihood = -189.13, time = 23.71s\n",
      "computing TAP Free Energies\n",
      "monitor:  -287.131358546 309.130764056 25.1704556729\n",
      "[EMF_RBM] Iteration 9, pseudo-likelihood = -184.82, time = 24.16s\n",
      "computing TAP Free Energies\n",
      "monitor:  -288.061831869 300.866291336 16.4625499447\n",
      "[EMF_RBM] Iteration 10, pseudo-likelihood = -180.39, time = 24.17s\n",
      "computing TAP Free Energies\n",
      "monitor:  -290.140907747 293.241650231 7.22989422779\n",
      "[EMF_RBM] Iteration 11, pseudo-likelihood = -172.95, time = 24.20s\n",
      "computing TAP Free Energies\n",
      "monitor:  -292.577006403 286.325885665 -1.69461087517\n",
      "[EMF_RBM] Iteration 12, pseudo-likelihood = -173.14, time = 23.75s\n",
      "computing TAP Free Energies\n",
      "monitor:  -294.818755643 280.37581557 -9.49969958907\n",
      "[EMF_RBM] Iteration 13, pseudo-likelihood = -167.37, time = 23.95s\n",
      "computing TAP Free Energies\n",
      "monitor:  -296.841672536 275.345648003 -16.1772344188\n",
      "[EMF_RBM] Iteration 14, pseudo-likelihood = -164.32, time = 23.97s\n",
      "computing TAP Free Energies\n",
      "monitor:  -298.967906647 270.884656634 -22.3945474038\n",
      "[EMF_RBM] Iteration 15, pseudo-likelihood = -161.35, time = 23.79s\n",
      "computing TAP Free Energies\n",
      "monitor:  -301.249982612 266.821716971 -28.3791092254\n",
      "[EMF_RBM] Iteration 16, pseudo-likelihood = -160.01, time = 23.65s\n",
      "computing TAP Free Energies\n",
      "monitor:  -303.644178279 263.04705621 -34.2010786311\n",
      "[EMF_RBM] Iteration 17, pseudo-likelihood = -152.33, time = 23.90s\n",
      "computing TAP Free Energies\n",
      "monitor:  -306.103820203 259.478278096 -39.9009683855\n",
      "[EMF_RBM] Iteration 18, pseudo-likelihood = -151.73, time = 24.14s\n",
      "computing TAP Free Energies\n",
      "monitor:  -308.569199597 256.084701724 -45.4516295781\n",
      "[EMF_RBM] Iteration 19, pseudo-likelihood = -148.98, time = 23.85s\n",
      "computing TAP Free Energies\n",
      "monitor:  -310.999946004 252.857278439 -50.8207891841\n",
      "[EMF_RBM] Iteration 20, pseudo-likelihood = -144.36, time = 23.60s\n",
      "computing TAP Free Energies\n",
      "[EMF_RBM] Iteration 22, pseudo-likelihood = -142.05, time = 23.85s\n",
      "computing TAP Free Energies\n",
      "monitor:  -318.063535897 244.151856893 -65.8130475602\n",
      "[EMF_RBM] Iteration 23, pseudo-likelihood = -141.84, time = 23.61s\n",
      "computing TAP Free Energies\n",
      "monitor:  -320.382116766 241.541516022 -70.5076430116\n",
      "[EMF_RBM] Iteration 24, pseudo-likelihood = -136.46, time = 23.84s\n",
      "computing TAP Free Energies\n",
      "monitor:  -322.688867994 239.049391895 -75.0854539272\n",
      "[EMF_RBM] Iteration 25, pseudo-likelihood = -138.74, time = 24.10s\n",
      "computing TAP Free Energies\n",
      "monitor:  -324.969774076 236.666486476 -79.5429870105\n",
      "[EMF_RBM] Iteration 26, pseudo-likelihood = -134.60, time = 24.09s\n",
      "computing TAP Free Energies\n",
      "monitor:  -327.199866937 234.394815146 -83.8537523361\n",
      "[EMF_RBM] Iteration 27, pseudo-likelihood = -132.60, time = 23.94s\n",
      "computing TAP Free Energies\n",
      "monitor:  -329.35034324 232.24333219 -87.9789199122\n",
      "[EMF_RBM] Iteration 28, pseudo-likelihood = -134.20, time = 23.60s\n",
      "computing TAP Free Energies\n",
      "monitor:  -331.399978177 230.220996205 -91.885856582\n",
      "[EMF_RBM] Iteration 29, pseudo-likelihood = -135.85, time = 23.79s\n",
      "computing TAP Free Energies\n",
      "monitor:  -333.34119133 228.331297174 -95.5603360863\n",
      "[EMF_RBM] Iteration 30, pseudo-likelihood = -134.14, time = 23.85s\n",
      "computing TAP Free Energies\n",
      "monitor:  -335.1826792 226.570410148 -99.0118124991\n",
      "[EMF_RBM] Iteration 31, pseudo-likelihood = -131.11, time = 23.72s\n",
      "computing TAP Free Energies\n",
      "monitor:  -336.944416626 224.928540008 -102.2675241\n",
      "[EMF_RBM] Iteration 32, pseudo-likelihood = -132.53, time = 23.76s\n",
      "computing TAP Free Energies\n",
      "monitor:  -338.651975748 223.391992455 -105.364984076\n",
      "[EMF_RBM] Iteration 33, pseudo-likelihood = -127.21, time = 24.17s\n",
      "computing TAP Free Energies\n",
      "monitor:  -340.330011563 221.945445412 -108.343254473\n",
      "[EMF_RBM] Iteration 34, pseudo-likelihood = -126.88, time = 24.05s\n",
      "computing TAP Free Energies\n",
      "monitor:  -341.998984147 220.574653575 -111.236770112\n",
      "[EMF_RBM] Iteration 35, pseudo-likelihood = -126.13, time = 23.83s\n",
      "computing TAP Free Energies\n",
      "monitor:  -343.67320731 219.266823654 -114.072877852\n",
      "[EMF_RBM] Iteration 36, pseudo-likelihood = -124.31, time = 23.75s\n",
      "computing TAP Free Energies\n",
      "monitor:  -345.361835723 218.011421013 -116.871791403\n",
      "[EMF_RBM] Iteration 37, pseudo-likelihood = -122.23, time = 23.87s\n",
      "computing TAP Free Energies\n",
      "monitor:  -347.067529882 216.800897999 -119.644357156\n",
      "[EMF_RBM] Iteration 38, pseudo-likelihood = -120.76, time = 23.83s\n",
      "computing TAP Free Energies\n",
      "monitor:  -348.788561394 215.629354547 -122.395384524\n",
      "[EMF_RBM] Iteration 39, pseudo-likelihood = -122.84, time = 23.68s\n",
      "computing TAP Free Energies\n",
      "monitor:  -350.521521859 214.492827501 -125.126019759\n",
      "[EMF_RBM] Iteration 40, pseudo-likelihood = -123.04, time = 23.96s\n",
      "computing TAP Free Energies\n",
      "monitor:  -352.260556504 213.388279008 -127.833921598\n",
      "[EMF_RBM] Iteration 41, pseudo-likelihood = -121.75, time = 24.41s\n",
      "computing TAP Free Energies\n",
      "monitor:  -354.000634082 212.313492571 -130.516675273\n",
      "[EMF_RBM] Iteration 42, pseudo-likelihood = -121.22, time = 24.14s\n",
      "computing TAP Free Energies\n",
      "monitor:  -355.734282901 211.267341212 -133.168200373\n",
      "[EMF_RBM] Iteration 43, pseudo-likelihood = -119.20, time = 23.67s\n",
      "computing TAP Free Energies\n",
      "monitor:  -357.456156451 210.249073063 -135.784028498\n",
      "[EMF_RBM] Iteration 44, pseudo-likelihood = -118.24, time = 23.88s\n",
      "computing TAP Free Energies\n",
      "monitor:  -359.160682605 209.257950017 -138.359386004\n",
      "[EMF_RBM] Iteration 45, pseudo-likelihood = -115.19, time = 23.95s\n",
      "computing TAP Free Energies\n",
      "monitor:  -360.843711501 208.293620478 -140.890446071\n",
      "[EMF_RBM] Iteration 46, pseudo-likelihood = -117.34, time = 23.99s\n",
      "computing TAP Free Energies\n",
      "monitor:  -362.501128723 207.355651366 -143.3734372\n",
      "[EMF_RBM] Iteration 47, pseudo-likelihood = -120.80, time = 23.73s\n",
      "computing TAP Free Energies\n",
      "monitor:  -364.130952311 206.443821249 -145.806428898\n",
      "[EMF_RBM] Iteration 48, pseudo-likelihood = -116.08, time = 23.99s\n",
      "computing TAP Free Energies\n",
      "monitor:  -365.731455057 205.557950494 -148.187672793\n",
      "[EMF_RBM] Iteration 49, pseudo-likelihood = -114.45, time = 24.16s\n",
      "computing TAP Free Energies\n",
      "monitor:  -367.301401896 204.697511043 -150.516222663\n",
      "[EMF_RBM] Iteration 50, pseudo-likelihood = -117.60, time = 24.11s\n",
      "computing TAP Free Energies\n",
      "monitor:  -368.840980774 203.861958719 -152.792587048\n",
      "[EMF_RBM] Iteration 51, pseudo-likelihood = -111.79, time = 23.87s\n",
      "computing TAP Free Energies\n",
      "monitor:  -370.350524526 203.050639116 -155.017553437\n",
      "[EMF_RBM] Iteration 52, pseudo-likelihood = -113.26, time = 23.77s\n",
      "computing TAP Free Energies\n",
      "monitor:  -371.831761362 202.262941588 -157.19323874\n",
      "[EMF_RBM] Iteration 53, pseudo-likelihood = -113.25, time = 23.91s\n",
      "computing TAP Free Energies\n",
      "monitor:  -373.284616714 201.498138424 -159.320101753\n",
      "[EMF_RBM] Iteration 54, pseudo-likelihood = -110.93, time = 23.77s\n",
      "computing TAP Free Energies\n",
      "monitor:  -374.710522842 200.755496389 -161.400119231\n",
      "[EMF_RBM] Iteration 55, pseudo-likelihood = -113.79, time = 23.71s\n",
      "computing TAP Free Energies\n",
      "monitor:  -376.111330562 200.03436878 -163.435609562\n",
      "[EMF_RBM] Iteration 56, pseudo-likelihood = -115.80, time = 23.88s\n",
      "computing TAP Free Energies\n",
      "monitor:  -377.488383854 199.333769167 -165.428765983\n",
      "[EMF_RBM] Iteration 57, pseudo-likelihood = -109.06, time = 23.96s\n",
      "computing TAP Free Energies\n",
      "monitor:  -378.842678404 198.653008579 -167.381124802\n",
      "[EMF_RBM] Iteration 58, pseudo-likelihood = -110.97, time = 23.66s\n",
      "computing TAP Free Energies\n",
      "monitor:  -380.175555443 197.991406844 -169.294571765\n",
      "[EMF_RBM] Iteration 59, pseudo-likelihood = -109.25, time = 23.67s\n",
      "computing TAP Free Energies\n",
      "monitor:  -381.488122617 197.348215756 -171.170853592\n",
      "[EMF_RBM] Iteration 60, pseudo-likelihood = -108.53, time = 24.30s\n",
      "computing TAP Free Energies\n",
      "monitor:  -382.781401291 196.722927852 -173.011369097\n",
      "[EMF_RBM] Iteration 61, pseudo-likelihood = -113.90, time = 24.48s\n",
      "computing TAP Free Energies\n",
      "monitor:  -384.056490553 196.114768645 -174.817891145\n",
      "[EMF_RBM] Iteration 62, pseudo-likelihood = -109.54, time = 24.23s\n",
      "computing TAP Free Energies\n",
      "monitor:  -385.314876782 195.522878692 -176.592678688\n",
      "[EMF_RBM] Iteration 63, pseudo-likelihood = -108.18, time = 23.71s\n",
      "computing TAP Free Energies\n",
      "monitor:  -386.557276471 194.946625813 -178.336992151\n",
      "[EMF_RBM] Iteration 64, pseudo-likelihood = -107.99, time = 24.00s\n",
      "computing TAP Free Energies\n",
      "monitor:  -387.784380225 194.385620647 -180.051804749\n",
      "[EMF_RBM] Iteration 65, pseudo-likelihood = -104.52, time = 23.88s\n",
      "computing TAP Free Energies\n",
      "monitor:  -388.997297142 193.8392409 -181.73877442\n",
      "[EMF_RBM] Iteration 66, pseudo-likelihood = -107.90, time = 23.82s\n",
      "computing TAP Free Energies\n",
      "monitor:  -390.197506899 193.306700389 -183.400101927\n",
      "[EMF_RBM] Iteration 67, pseudo-likelihood = -104.73, time = 23.71s\n",
      "computing TAP Free Energies\n",
      "monitor:  -391.385633618 192.787498214 -185.036847889\n",
      "[EMF_RBM] Iteration 68, pseudo-likelihood = -105.98, time = 23.82s\n",
      "computing TAP Free Energies\n",
      "monitor:  -392.562641112 192.2810863 -186.650477143\n",
      "[EMF_RBM] Iteration 69, pseudo-likelihood = -106.99, time = 24.12s\n",
      "computing TAP Free Energies\n",
      "monitor:  -393.729058 191.786755432 -188.242198288\n",
      "[EMF_RBM] Iteration 70, pseudo-likelihood = -105.57, time = 23.79s\n",
      "computing TAP Free Energies\n",
      "monitor:  -394.885408832 191.303936069 -189.813070408\n",
      "[EMF_RBM] Iteration 71, pseudo-likelihood = -101.47, time = 23.62s\n",
      "computing TAP Free Energies\n",
      "monitor:  -396.032556024 190.832069552 -191.364490788\n",
      "[EMF_RBM] Iteration 72, pseudo-likelihood = -105.17, time = 23.68s\n",
      "computing TAP Free Energies\n"
     ]
    }
   ],
   "source": [
    "rbm = emf_rbm.EMF_RBM(verbose=True, monitor=True, n_iter=100)\n",
    "print rbm\n",
    "rbm= rbm.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-11-08T22:57:27.868434",
     "start_time": "2016-11-08T22:57:27.697130"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.plot(rbm.entropies-rbm.entropies[0], color='red')\n",
    "plt.plot(rbm.mean_field_energies-rbm.mean_field_energies[0], color='blue')\n",
    "plt.plot(rbm.free_energies-rbm.free_energies[0], color='green')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rbm = emf_rbm.EMF_RBM(verbose=True, monitor=True, n_iter=100)\n",
    "print rbm\n",
    "rbm= rbm.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(rbm.entropies-rbm.entropies[0], color='red')\n",
    "plt.plot(rbm.mean_field_energies-rbm.mean_field_energies[0], color='blue')\n",
    "plt.plot(rbm.free_energies-rbm.free_energies[0], color='green')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}