{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Appendix A: A fast version\n",
      "\n",
      "This takes a (very long) while to run unless the tools in osmosis.parallel are used to parallelize the computation of all the model params and the rRMSE values.\n",
      "\n",
      "For completeness, the notebook AppendixA will do that for you, but god only knows how long that might take. This here, on the other hand, works FAST!\n",
      "\n",
      "We import the needed modules: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as ni\n",
      "import os\n",
      "import osmosis.model.analysis as oza\n",
      "import osmosis.model.sparse_deconvolution as ssd\n",
      "import osmosis.model.dti as dti\n",
      "import osmosis.viz.mpl as mpl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following points the data-path variable to where the data is stored, relative to the install location"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import osmosis as oz\n",
      "import osmosis.io as oio\n",
      "oio.data_path = os.path.join(oz.__path__[0], 'data')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Set ht subject to be analyzed:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subject = 'SUB1'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We sweep through particular regularization parameter combinations in each b value: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alphas = [0.0001,0.0005, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05]\n",
      "l1_ratios = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the following loop, we collect the median value of the rRMSE from pre-computed rRMSE files for each combination of regularization parameters: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rmse_matrix = {}\n",
      "rel_matrix = {}\n",
      "for b in [1000, 2000, 4000]: # \n",
      "    ad_rd = oio.get_ad_rd(subject, b)\n",
      "    rmse_matrix[b] = np.zeros((len(l1_ratios), len(alphas)))\n",
      "    rel_matrix[b] = np.zeros((len(l1_ratios), len(alphas)))\n",
      "    for l1_idx, l1_ratio in enumerate(l1_ratios):\n",
      "        for alpha_idx, this_alpha in enumerate(alphas):\n",
      "            \n",
      "            rrmse_fname = \"%s/%s/RMSE_b%s_l1ratio%s_alpha%s.nii.gz\"%(oio.data_path,\n",
      "                                                             subject,\n",
      "                                                             b,\n",
      "                                                             l1_ratio,\n",
      "                                                             this_alpha)\n",
      "            \n",
      "            rrmse_fname = \"%s/%s/PDD_angle_b%s_l1ratio%s_alpha%s.nii.gz\"%(oio.data_path,\n",
      "                                                             subject,\n",
      "                                                             b,\n",
      "                                                             l1_ratio,\n",
      "                                                             this_alpha)\n",
      "            \n",
      "            rrmse = ni.load(oio.data_path + '/%s/'%subject + 'RMSE_b%s_l1ratio%s_alpha%s.nii.gz'%(b, l1_ratio, this_alpha)).get_data() \n",
      "            median_rmse = np.median(rrmse[np.isfinite(rrmse)])\n",
      "            rmse_matrix[b][l1_idx, alpha_idx] = median_rmse\n",
      "            rel = ni.load(oio.data_path + '/%s/'%subject + 'PDD_angle_b%s_l1ratio%s_alpha%s.nii.gz'%(b, l1_ratio, this_alpha)).get_data() \n",
      "            median_rel = np.median(rel[np.isfinite(rel)])\n",
      "            rel_matrix[b][l1_idx, alpha_idx] = median_rel\n",
      "            \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The results are visualized: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vmax = 1.5\n",
      "vmin = 0.5\n",
      "for b in [1000, 2000, 4000]:\n",
      "    fig, ax = plt.subplots(1)\n",
      "    cax = ax.imshow(rmse_matrix[b], interpolation='nearest', cmap=cm.RdYlGn_r, vmax=vmax, vmin=vmin)\n",
      "    cbar = fig.colorbar(cax, ticks=[0.5, 1., 1.5])\n",
      "    cax.axes.set_xticks([0,1,2,3,4,5,6,7,8])\n",
      "    cax.axes.set_xticklabels([str('%1.4f'%this) for this in alphas])\n",
      "    cax.axes.set_xlabel(r'$\\lambda$')\n",
      "    cax.axes.set_yticks([0,1,2,3,4,5])\n",
      "    cax.axes.set_yticklabels([str('%1.1f'%(1-this)) for this in l1_ratios])\n",
      "\n",
      "    cax.axes.set_ylabel(r'$\\alpha$')\n",
      "    im = ax.matshow(rmse_matrix[b], cmap=matplotlib.cm.RdYlGn_r, vmax=vmax, vmin=vmin)\n",
      "    fig.set_size_inches([8,6])\n",
      "    fig.savefig('figures/AppendixA_acc_b%s.svg'%b)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAFlCAYAAACJCrOeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHchJREFUeJzt3X9wVdW99/HPDgkXAxRJaVNzTvpEkoMJKCE0lNqOFe1Q\nCJXoWPtM8F6lmNIYpbYzjONM/3hG6890Wmd6zVOMM1Xr2CKaokENpw4+QouURBRxrqQYlQwnQTMi\noFxpTThZzx+9nDFCkg0n+6y9kveLOTPZydp7f8+anXz5rrXO3p4xxggAAIwoy3YAAAC4gqQJAIBP\nJE0AAHwiaQIA4BNJEwAAn0iaAAD4RNIEAIxJN9xwg/Lz83XRRRed9udbt27VtGnTVFFRoYqKCt11\n110jHjN7tIMEACAMVq1apZ/85Ce6/vrrh2xz6aWXatOmTb6PSaUJABiTLrnkEk2fPn3YNmd6fx+S\nJgDAqry8PHmel9Zr6tSpZ3xez/O0Y8cOlZeXa9myZdq7d++I+zA8CwCw6siRI/qDZqV1jH//77fO\neJ/58+crkUgoNzdXmzdv1lVXXaW33hr+OFSaAIBxaerUqcrNzZUkVVVVqb+/X4cPHx52H5ImAMC6\nrKz0Xmejt7c3NafZ3t4uY4zy8vKG3YfhWQCAdWeb+FIGTv3WihUrtG3bNh06dEiFhYW644471N/f\nL0mqq6tTc3Oz1q1bp+zsbOXm5uqJJ54Y8TQejwYDANjkeZ6empjenOYP+t4645WwZ4PhWQAAfGJ4\nFgBgXdrDsxlC0gQAWJfl2Y7AH5ImAMA6VypNR8IEAMA+Kk0AgHWuVJokTQCAdSRNAAB8ciVpOhIm\nAAD2UWkCAKxzpdIkaQIArCNpAgDgE0kTAACfPM+NWwI5ktsBALCPShMAYB3DswAA+ETSBADAJ1eS\npiNhAgBgH5UmAMA6VypNkiYAwDqSJgAAPrmSNB0JEwAA+6g0AQDWuVJpkjQBANaRNAEA8CnLjVvP\nMqcJAIBfVJoAAOsYngUAwCeSJgAAPrmSNB0JEwAA+6g0AQDWZTmyfNapSjMej6u0tFSxWEwNDQ2n\nbXPLLbcoFoupvLxcu3fvHnHfw4cPa/HixZo1a5a++93v6ujRo5Kkrq4unXPOOaqoqFBFRYVuuumm\nYN9cBgTRf0899ZTmzJmjCRMm6LXXXht0rHvvvVexWEylpaV64YUXgnlTGRJE3916660qKytTeXm5\nrr76an300UeSxt61F0Tf3X777YpGo6k+2rx5c+pnY+m6k4Lpv5qamlTfnX/++aqoqJBk99rzJnhp\nvTLGOOLEiROmuLjY7N+/3/T19Zny8nKzd+/eQW2ef/55U1VVZYwxZufOnWbhwoUj7nvrrbeahoYG\nY4wx9913n7ntttuMMcbs37/fXHjhhZl6e4ELqv86OjrMvn37zKJFi8yrr76aOtabb75pysvLTV9f\nn9m/f78pLi42yWQyQ+92dAXVdy+88EKqT2677bYxee0F1Xe33367+fWvf33K+cbSdWdMcP33WWvX\nrjV33nmnMcbetSfJvPX12Wm9MpXOnKk029vbVVJSoqKiIuXk5KimpkYtLS2D2mzatEkrV66UJC1c\nuFBHjx7V+++/P+y+n91n5cqVeuaZZzL7xjIkqP4rLS3VrFmzTjlfS0uLVqxYoZycHBUVFamkpETt\n7e3Bv9EABNV3ixcvVtb/rH5YuHChuru7M/vGMiCovpOkf/2tHWwsXXdSsP0n/asPn3zySa1YsSJj\n72koXpaX1itTnEmaPT09KiwsTG1Ho1H19PT4anPw4MEh9+3t7VV+fr4kKT8/X729val2+/fvV0VF\nhRYtWqTt27cH8r4yJaj+G8rBgwcVjUbPaJ+wykTfPfzww1q2bFlqe6xce0H23QMPPKDy8nLV1tam\nplXG0nUnBX/t/fWvf1V+fr6Ki4tT3xsr115QnEmanufvfxKn+9/n6dqc7nie56W+X1BQoEQiod27\nd+v+++/Xtddeq2PHjp1Z0CEymv0XdAxhE3Tf3X333Zo4caKuvfZaSWPr2guq7+rr67V//369/vrr\nOu+887R27dq0YwijoK+99evXp647ye6158qcpjOrZyORiBKJRGo7kUgM+h/l6dp0d3crGo2qv7//\nlO9HIhFJ/6ou33//fX3lK1/Re++9py9/+cuSpIkTJ2rixImSpPnz56u4uFidnZ2aP39+YO8xSKPZ\nf6fbd6TzfbbPXRNk3z366KNqbW3Viy++mPreWLr2guq7k7+nkvSjH/1Iy5cvH/JYrl53UrDX3okT\nJ/T0008PWsBn89rL5BBrWjIyczoK+vv7zcyZM83+/fvNp59+OuKE+N/+9rfUhPhw+956663mvvvu\nM8YYc++996YWY3zwwQfmxIkTxhhj3nnnHROJRMyRI0cy8l6DEFT/nbRo0SKza9eu1PbJBRmffvqp\neffdd83MmTPNwMBAwO8yGEH13ebNm83s2bPNBx98MOhYY+naC6rvDh48mNr//vvvNytWrDDGjK3r\nzphgf283b95sFi1aNOhYtq49SebdRXPTemUqnTmTNI0xprW11cyaNcsUFxebe+65xxhjzIMPPmge\nfPDBVJubb77ZFBcXm7lz5w5azXm6fY0x5sMPPzTf+c53TCwWM4sXL05dIH/605/MnDlzzLx588z8\n+fPNc889l6F3GZwg+m/jxo0mGo2aSZMmmfz8fLN06dLUz+6++25TXFxsLrjgAhOPxzPwDoMTRN+V\nlJSYr371q2bevHlm3rx5pr6+3hhjTHNz85i69oLou+uuu85cdNFFZu7cuebKK68077//fupnY+m6\nMyaY/jPGmB/+8Iemqalp0Pds/d1zKWl6/xMwAABWeJ6nrsvL0zpG0f/bE+iajJOcmdMEAIxdrsxp\nkjQBANZl9K4+aQh10nR5qTgAjEXjfUYv1ElTkh76rxtH7VjP/naXlt9UOWrHk6SBgVE9XGg9t26X\nrqgf3b5LjqPfvefX7dL3Rrn/xosg+i5nwqgebtTlZo9OwbCx8RVdvWbBqBxLkv6jdN2oHevzPEee\nDRb6pAkAGPsYngUAwCdXFgK5UQ+PklkLCmyH4KxZlfRdOmL031mj785e2dfpu9E2rirNC0iaZ43/\ncKSH/jt79N3ZK/u6O7cQdKXSHFdJEwAQTsxpAgDgUxaVJgAA/rgyPDuuFgIBAJAOKk0AgHXMaQIA\n4JMrw7MkTQCAda5UmsxpAgDgE5UmAMA6btgOAIBPzGkCAOBT1nif04zH4yotLVUsFlNDQ8MpPz90\n6JCWLl2qefPm6cILL9Sjjz4aVCgAAIyKQJJmMpnUmjVrFI/HtXfvXq1fv14dHR2D2jQ2NqqiokKv\nv/66tm7dqrVr1+rEiRNBhAMACDkvy0vrlSmBJM329naVlJSoqKhIOTk5qqmpUUtLy6A25513nj7+\n+GNJ0scff6wvfvGLys5mtBgAxiNvgpfWK1MCyVI9PT0qLCxMbUejUbW1tQ1qs3r1al1++eUqKCjQ\nsWPH9OSTT572WM/+dlfq61kLCni8FwBkSEd7jzraD2bkXON6IZDnjfzm77nnHs2bN09bt27VO++8\no8WLF2vPnj2aOnXqoHbLb6oMIkQAwAjKvh4Z9EzOjY27hmk9PgQyPBuJRJRIJFLbiURC0Wh0UJsd\nO3boBz/4gSSpuLhY559/vvbt2xdEOACAsJvgpffKkECSZmVlpTo7O9XV1aW+vj5t2LBB1dXVg9qU\nlpZqy5YtkqTe3l7t27dPM2fODCIcAEDIubIQKJDh2ezsbDU2NmrJkiVKJpOqra1VWVmZmpqaJEl1\ndXX6+c9/rlWrVqm8vFwDAwP65S9/qby8vCDCAQCE3QQ37gjkGWOM7SCG4nmeHvqvG22HMayBAdsR\nuCsZ2isPY13OBNsRDC83O5yLYv6jdJ2CSBme5+mT25amdYzJDfFAYvs8N1I7AAAhwAcjAQDWufJo\nMJImAMC+8fw5TQAAzogjlSZzmgAA+ESlCQCwblzfRg8AgDPiyOc0SZoAAOtcqTTdSO0AAIQAlSYA\nwD5HVs+SNAEA9pE0AQDwhzlNAADGGCpNAIB9DM+OjrA/eqt/gOdbna2ckA/HTPu3cA/E5E2aYjuE\nYc2YNNV2CEOadyzcf/r2X/+47RAyLyvcv28nhfvKAQCMCzzlBAAAv0I+8nSSG/UwAAAhQKUJALCP\n4VkAAPxx5XOaJE0AgH2OPOXEjSgBAAgBKk0AgH3MaQIA4A9zmgAA+OVIpcmcJgAAPlFpAgDsY3gW\nAAB/uPcsAAB+OfKUEzeiBAAgBKg0AQD2OTKnGVilGY/HVVpaqlgspoaGhtO22bp1qyoqKnThhRdq\n0aJFQYUCAAi7rKz0XhkSSKWZTCa1Zs0abdmyRZFIRAsWLFB1dbXKyspSbY4ePaqbb75Zf/7znxWN\nRnXo0KEgQgEAuGA8z2m2t7erpKRERUVFysnJUU1NjVpaWga1+eMf/6jvf//7ikajkqQZM2YEEQoA\nwAVZXnqvDAmk0uzp6VFhYWFqOxqNqq2tbVCbzs5O9ff367LLLtOxY8f005/+VNddd90px3pu3a7U\n17MqCzRrQUEQIQMAPqfto0/U9vEntsMIlUCSpueNnPX7+/v12muv6cUXX9Tx48d18cUX6xvf+IZi\nsdigdlfUVwYRIgBgBAunTdbCaZNT2/+Z+CC4kzkyPBtI0oxEIkokEqntRCKRGoY9qbCwUDNmzNA5\n55yjc845R9/+9re1Z8+eU5ImAGAccCRpBhJlZWWlOjs71dXVpb6+Pm3YsEHV1dWD2lx55ZXavn27\nksmkjh8/rra2Ns2ePTuIcAAAYTee5zSzs7PV2NioJUuWKJlMqra2VmVlZWpqapIk1dXVqbS0VEuX\nLtXcuXOVlZWl1atXkzQBAKHmGWOM7SCG4nmeHnzjRtthDKt/ILTdF3o5If8w87R/C/dwUd6kKbZD\nGNaMSVNthzCkecfCfV+X/dc/bjuE0yrZ8aaCSBme5yn50k/TOsaEy34TSGyfF+4rBwAwPjgyp0nS\nBADYF/KRp5PcSO0AAIQAlSYAwD6GZwEA8ImkCQCAP37uJBcGbqR2AABCgEoTAGAfw7MAAPhE0gQA\nwCdHPqdJ0gQA2OdIpelGlAAAhACVJgDAPoZnAQDwyZHh2dAnzSRP3kpLVog/MDxlYnhjk6SZ075s\nO4Rhnf+FYtshDGvGhDzbIQzp4A3/x3YIw2rb0W87hMxzJGm6ESUAACEQ+koTADAOMKcJAIBPjgzP\nkjQBAPY5kjTdiBIAgBCg0gQA2MecJgAAPjkyPEvSBADY57mRNN2IEgCAEKDSBADY50ilSdIEANjH\nnCYAAD5RaQIA4JMjSdONKAEACAEqTQCAfeO90ozH4yotLVUsFlNDQ8OQ7V555RVlZ2dr48aNQYUC\nAAi7rKz0XpkKM4iDJpNJrVmzRvF4XHv37tX69evV0dFx2na33Xabli5dKmN42jQAjFteVnqvDAnk\nTO3t7SopKVFRUZFycnJUU1OjlpaWU9o98MADuuaaa/SlL30piDAAABhVgSTNnp4eFRYWpraj0ah6\nenpOadPS0qL6+npJkue5cbNeAEAAHKk0A1kI5CcB/uxnP9N9990nz/NkjBlyePb5dbtSX8cqCzRr\nQcGoxQkAGNpeHVeH/pGZkzmyECiQpBmJRJRIJFLbiURC0Wh0UJtXX31VNTU1kqRDhw5p8+bNysnJ\nUXV19aB236uvDCJEAMAIZitXs5Wb2t6oD4M72XhOmpWVlers7FRXV5cKCgq0YcMGrV+/flCbd999\nN/X1qlWrtHz58lMSJgAAYRJI0szOzlZjY6OWLFmiZDKp2tpalZWVqampSZJUV1cXxGkBAK4a7/ee\nraqqUlVV1aDvDZUsH3nkkaDCAAA4wPMm2A7BF+4IBACwz5E5TTeiBAAgBKg0AQD2OVJpkjQBAPaR\nNAEA8Gm8r54FAMA3RypNN6IEACAEqDQBAPY5UmmSNAEA9pE0AQDwyZGFQG5ECQBACJA0AQD2BfAQ\n6htuuEH5+fm66KKLhjztLbfcolgspvLycu3evXvEMEmaAAD7Akiaq1atUjweH/KUra2tevvtt9XZ\n2amHHnpI9fX1I4Y5YtI8cOCAWlpatH37dklSc3PziAcFAOCMBJA0L7nkEk2fPn3IU27atEkrV66U\nJC1cuFBHjx5Vb2/vsGGOuBDok08+UXNzs1555RVNnjxZ3/zmN3XNNdeMtBsAAIHZunWPtm7bk9Yx\nenp6VFhYmNqORqPq7u5Wfn7+kPuMmDSffvpp/eY3v1FeXp6OHz+ul156Ka0gkVkTPNsRDG3axHNs\nhzCsr079X7ZDGNaM/nAvfjd//5vtEIa0PX7Mdgj4vDP8yMmiyyq06LKK1PYdv3jsrE5rjBkchjf8\nH80Rf+sKCwuVl5cnScrNzT3lBAAApMtY+A9+JBJRIpFIbXd3dysSiQy7z4ipfcaMGaqpqdGzzz6r\nPXv2qKOjI/1IAQD4DGMG0nqdjerqaj322L8q1J07d+rcc88ddmhW8lFpVlVVaebMmfr973+vbdu2\nafXq1WcVHAAAmbRixQpt27ZNhw4dUmFhoe644w719/dLkurq6rRs2TK1traqpKREkydP1iOPPDLi\nMX1NilxwwQW655570oseAIAhDJxltTic9evXj9imsbHxjI4Z7pUEAIBxwWj0k2YQSJoAAOuCqDSD\nQNIEAFjnSqXJbfQAAPCJShMAYB3DswAA+HS2n7XMNJImAMC6AeY0AQAYW6g0AQDWMTwLAIBPLAQC\nAMCncf85zXg8rtLSUsViMTU0NJzy8z/84Q8qLy/X3Llz9a1vfUtvvPFGUKEAADAqAqk0k8mk1qxZ\noy1btigSiWjBggWqrq5WWVlZqs3MmTP1l7/8RdOmTVM8HtePf/xj7dy5M4hwAAAh58rwbCCVZnt7\nu0pKSlRUVKScnBzV1NSopaVlUJuLL75Y06ZNkyQtXLhQ3d3dQYQCAHCAUTKtV6YEUmn29PSosLAw\ntR2NRtXW1jZk+9/97ndatmzZaX/2/Lpdqa9jlQWataBg9AIFAAxpr46rQ//IyLlcqTQDSZqe5/lu\n+9JLL+nhhx/Wyy+/fNqff6++crTCAgCcgdnK1WzlprY36kOL0YRDIEkzEokokUikthOJhKLR6Cnt\n3njjDa1evVrxeFzTp08PIhQAgANc+ZxmIHOalZWV6uzsVFdXl/r6+rRhwwZVV1cPanPgwAFdffXV\nevzxx1VSUhJEGAAARwyk+S9TAqk0s7Oz1djYqCVLliiZTKq2tlZlZWVqamqSJNXV1ekXv/iFjhw5\novr6eklSTk6O2tvbgwgHABByrlSanjHG2A5iKJ7n6f/uudF2GMMaCG/3SZJysvzPL2daZEruyI0s\nqsy/yHYIw8o/Mcl2CMMy77xpO4QhPbXgSdshDCuZucWgZ+Tf9ZaCSBme5ylxrCmtYxROrQskts/j\nhu0AAPjEbfQAANa5chs9kiYAwLpx/TlNAADOhCuVJnOaAAD4RKUJALCO4VkAAHxy5XOaJE0AgHWZ\nvKtPOpjTBADAJypNAIB1DM8CAOBT2G9JehJJEwBgnSurZ5nTBADAJypNAIB1DM+OkrB3ZJYX3kdv\nSVLOhPDGlzdpiu0QhpU/scB2CMMy7+ywHcKwku2dtkMYUlgfvTWeDSjcf+tPCn3SBACMfWEvkE4i\naQIArGMhEAAAYwyVJgDAOoZnAQDwiaQJAIBPzGkCADDGUGkCAKxjeBYAAJ+4uQEAAD4xpwkAwBhD\npQkAsI45TQAAfCJpAgDgkytJkzlNAAB8otIEAFg37lfPxuNxlZaWKhaLqaGh4bRtbrnlFsViMZWX\nl2v37t1BhQIACLkBY9J6ZUogSTOZTGrNmjWKx+Pau3ev1q9fr46OjkFtWltb9fbbb6uzs1MPPfSQ\n6uvrgwgFAOCAcZ0029vbVVJSoqKiIuXk5KimpkYtLS2D2mzatEkrV66UJC1cuFBHjx5Vb29vEOEA\nAEJuIM1/mRLInGZPT48KCwtT29FoVG1tbSO26e7uVn5+/qB2rQ/uSn0dqyxQrLIgiJABAJ+zV8fV\noX/YDiNUAkmanuf5amc+V1Kfbr9lN1aOSkwAgDMzW7mardzU9kZ9GNi5XPnISSBJMxKJKJFIpLYT\niYSi0eiwbbq7uxWJRIIIBwAQcq4kzUDmNCsrK9XZ2amuri719fVpw4YNqq6uHtSmurpajz32mCRp\n586dOvfcc08ZmgUAjA+uLAQKpNLMzs5WY2OjlixZomQyqdraWpWVlampqUmSVFdXp2XLlqm1tVUl\nJSWaPHmyHnnkkSBCAQBg1AR2c4OqqipVVVUN+l5dXd2g7cbGxqBODwBwiCs3N+COQAAA61yZ0yRp\nAgCscyVpcsN2AAB8otIEAFjnSqVJ0gQAWJd0I2eSNAEA9rlSaTKnCQCAT1SaAADrGJ4FAMCnAZIm\nAAD+JB2Z0yRpAgCsc6XSZCEQAAA+UWkCAKxjIRAAAD65MjxL0hzjckI8AJ83abLtEIZ3pNt2BMN7\n+4DtCIb1320HbYcAh7iyECjEf1IBAAgXKk0AgHUDbjyDmqQJALCPhUAAAPjkykIg5jQBAPCJShMA\nYJ0rq2dJmgAA61wZniVpAgCsc2UhEHOaAAD4RKUJALBugDlNAAD8cWV4lqQJALCOhUAAAPjkSqXJ\nQiAAAHyi0gQAWDfgyPgsSRMAYJ0rw7MkTQCAdY4UmsxpAgDgF5UmAMA6hmdHSeuDu1JfxyoLFKss\nsBgNAIwfe3VcHfpHRs7FHYFGybIbK22HAADj0mzlarZyU9sb9WFg53Kl0mROEwAAn0JfaQIAxj4e\nQg0AgE/JAdsR+EPSBABY50qlyZwmAAA+UWkCAKxLOnJLIJImAMA6V4ZnSZoAAOtYCAQAgE+uVJos\nBAIAwCcqTQCAdSwEAgDAJ1eGZ0maAADrXFkINK7mNDt3HbQdgrM62ntsh+C0rS//3XYIztr+3lHb\nIThrr47bDmHMIWnCl7+303fp2LZjn+0QnPXy+x/ZDsFZmXoW5mhIGpPWK1MYngUAWMecJgAAPrmy\netYzJrzp3fM82yEAAD4jiJTheZ4ufXJFWsfY9r/XBxLb54W60gxxPgcAjKKkI3/uQ500AQDjgyvD\nsyRNAIB1riwEGlcfOQEAIB1UmgAA6xieBQDAJxYCAQDgkyuVJnOaQICOHDmia6+9VocPH7YdChBq\nrtxGj6QJBGj69Om6/PLL1dzcbDsUAKOApAkEbPny5WppabEdBhBqrlSazGkCAcvPz9cnn3yijz/+\nWF/4whdshwOEEs/TBCBJ+uc//6kpU6bo+eeftx0KEFquVJokTSBAyWRSt99+u+68804988wztsMB\nkCaGZ4EArV27VqtXr9acOXN04MAB9fX1aeLEibbDAkKHj5wA41xzc7O+9rWvac6cOZKkK664Qq2t\nrZajAsLJleHZUD9PEwAw9nmep/PWXZXWMd6rfyYjj5Ok0gQAwCfmNAEA1hlH5jRJmgAA60iaAAD4\n5ErSZE4TAACfqDQBANa58kEOkiYAwDpXhmdJmgAA60iaAAD45ErSZCEQAAA+UWkCAKxzpdIkaQIA\nrCNpAgDgkytJkzlNAAB8otIEAFjnSqVJ0gQAWEfSBADAJ1duo8ecJgAAPlFpAgCsY3gWAACfSJoA\nAPjkStJkThMAAJ+oNAEA1rlSaZI0AQDWkTQBAPCJpAkAgE+uJE0WAgEA4BOVJgDAOlduo0fSBABY\nx/AsAAA+mQGT1mso8XhcpaWlisViamhoOOXnW7du1bRp01RRUaGKigrdddddw8ZJpQkAGJOSyaTW\nrFmjLVu2KBKJaMGCBaqurlZZWdmgdpdeeqk2bdrk65hUmgAA64KoNNvb21VSUqKioiLl5OSopqZG\nLS0tp577DOZTqTQBANaZ3+5Ma/8pU6ac8r2enh4VFhamtqPRqNra2ga18TxPO3bsUHl5uSKRiH71\nq19p9uzZQ56HpAkAsCqolbOe543YZv78+UokEsrNzdXmzZt11VVX6a233hqyPcOzAIAxKRKJKJFI\npLYTiYSi0eigNlOnTlVubq4kqaqqSv39/Tp8+PCQxyRpAgDGpMrKSnV2dqqrq0t9fX3asGGDqqur\nB7Xp7e1NVbrt7e0yxigvL2/IYzI8CwAYk7Kzs9XY2KglS5YomUyqtrZWZWVlampqkiTV1dWpublZ\n69atU3Z2tnJzc/XEE08Me0zPuHIbBgAALGN4FgAAn0iaAAD4RNIEAMAnkiYAAD6RNAEA8ImkCQCA\nT/8fPW8C/NVS8YgAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4107d90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAFlCAYAAACJCrOeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHe5JREFUeJzt3X9wVdXd7/HPDjlQgoik2FTOSRtJDiaChEC41HasaAdD\neGp0rL0TvGMpppgnldrOMI4z/Uvrz/S2zu01TzXOVKz9gdhUDdZw9MErtNRCBBFnCo9ESS4nQfOI\nQPWKknCy7h99PGOEJBtO9ll7Je8Xc2ayk3X2/p41m3zzXWvtvT1jjBEAABhRju0AAABwBUkTAACf\nSJoAAPhE0gQAwCeSJgAAPpE0AQDwiaQJABiTbrrpJhUUFOiSSy457c+3bNmiadOmqaKiQhUVFbr7\n7rtH3GfuaAcJAEAYrFq1Sj/4wQ/0ne98Z8g2l19+uTZu3Oh7n1SaAIAx6bLLLtP06dOHbXOm9/ch\naQIArMrPz5fneRm9pk6desbH9TxPL7/8ssrLy7V8+XLt3bt3xPcwPAsAsOro0aP6nWZntI//8f/2\nn/F7FixYoGQyqby8PG3atEnXXnut9u8ffj9UmgCAcWnq1KnKy8uTJFVXV6u/v19HjhwZ9j0kTQCA\ndTk5mb3ORm9vb3pOs729XcYY5efnD/sehmcBANadbeJLGzj1WytWrNDWrVt1+PBhFRYW6s4771R/\nf78kqb6+Xi0tLXrooYeUm5urvLw8PfHEEyMexuPRYAAAmzzP0x8mZjan+e2+/We8EvZsMDwLAIBP\nDM8CAKzLeHg2S0iaAADrcjzbEfhD0gQAWOdKpelImAAA2EelCQCwzpVKk6QJALCOpAkAgE+uJE1H\nwgQAwD4qTQCAda5UmiRNAIB1JE0AAHwiaQIA4JPnuXFLIEdyOwAA9lFpAgCsY3gWAACfSJoAAPjk\nStJ0JEwAAOyj0gQAWOdKpUnSBABYR9IEAMAnV5KmI2ECAGAflSYAwDpXKk2SJgDAOpImAAA+5bhx\n61nmNAEA8ItKEwBgHcOzAAD4RNIEAMAnV5KmI2ECAGAflSYAwLocR5bPOlVpJhIJlZaWKh6Pq7Gx\n8bRtbr31VsXjcZWXl2v37t0jvvfIkSNaunSpZs+erauuukrHjh2TJHV1dWny5MmqqKhQRUWFvv/9\n7wf74bIgiP77wx/+oDlz5mjChAl69dVXB+3rvvvuUzweV2lpqV544YVgPlSWBNF3t912m8rKylRe\nXq7rrrtO//jHPySNvXMviL674447FIvF0n20adOm9M/G0nknBdN/tbW16b678MILVVFRIcnuuedN\n8DJ6ZY1xxMmTJ01xcbHp7Ow0fX19pry83Ozdu3dQm+eee85UV1cbY4zZvn27Wbx48Yjvve2220xj\nY6Mxxpj777/f3H777cYYYzo7O83cuXOz9fECF1T/7du3z7zxxhtmyZIlZteuXel9/f3vfzfl5eWm\nr6/PdHZ2muLiYpNKpbL0aUdXUH33wgsvpPvk9ttvH5PnXlB9d8cdd5if//znpxxvLJ13xgTXf5+2\ndu1ac9dddxlj7J17ksz+/3ZxRq9spTNnKs329naVlJSoqKhIkUhEtbW1am1tHdRm48aNWrlypSRp\n8eLFOnbsmN55551h3/vp96xcuVLPPPNMdj9YlgTVf6WlpZo9e/Ypx2ttbdWKFSsUiURUVFSkkpIS\ntbe3B/9BAxBU3y1dulQ5/7X6YfHixeru7s7uB8uCoPpOkv75u3awsXTeScH2n/TPPnzyySe1YsWK\nrH2moXg5XkavbHEmafb09KiwsDC9HYvF1NPT46vNoUOHhnxvb2+vCgoKJEkFBQXq7e1Nt+vs7FRF\nRYWWLFmibdu2BfK5siWo/hvKoUOHFIvFzug9YZWNvnv00Ue1fPny9PZYOfeC7LsHH3xQ5eXlqqur\nS0+rjKXzTgr+3PvLX/6igoICFRcXp783Vs69oDiTND3P318Sp/vr83RtTrc/z/PS3585c6aSyaR2\n796tBx54QDfccIM++OCDMws6REaz/4KOIWyC7rt77rlHEydO1A033CBpbJ17QfVdQ0ODOjs79dpr\nr+mCCy7Q2rVrM44hjII+99avX58+7yS7554rc5rOrJ6NRqNKJpPp7WQyOegvytO16e7uViwWU39/\n/ynfj0ajkv5ZXb7zzjv64he/qLfffltf+MIXJEkTJ07UxIkTJUkLFixQcXGxOjo6tGDBgsA+Y5BG\ns/9O996RjvfpPndNkH332GOPqa2tTS+++GL6e2Pp3Auq7z75fypJ3/ve93T11VcPuS9Xzzsp2HPv\n5MmTevrppwct4LN57mVziDUjWZk5HQX9/f1m1qxZprOz05w4cWLECfG//e1v6Qnx4d572223mfvv\nv98YY8x9992XXozx7rvvmpMnTxpjjHnrrbdMNBo1R48ezcpnDUJQ/feJJUuWmJ07d6a3P1mQceLE\nCXPgwAEza9YsMzAwEPCnDEZQfbdp0yZz8cUXm3fffXfQvsbSuRdU3x06dCj9/gceeMCsWLHCGDO2\nzjtjgv1/u2nTJrNkyZJB+7J17kkyB5bMy+iVrXTmTNI0xpi2tjYze/ZsU1xcbO69915jjDEPP/yw\nefjhh9NtbrnlFlNcXGzmzZs3aDXn6d5rjDHvvfee+cY3vmHi8bhZunRp+gT54x//aObMmWPmz59v\nFixYYP70pz9l6VMGJ4j+e+qpp0wsFjOf+9znTEFBgVm2bFn6Z/fcc48pLi42F110kUkkEln4hMEJ\nou9KSkrMl770JTN//nwzf/5809DQYIwxpqWlZUyde0H03Y033mguueQSM2/ePHPNNdeYd955J/2z\nsXTeGRNM/xljzHe/+13T3Nw86Hu2fu+5lDS9/woYAAArPM9T15XlGe2j6P/sCXRNxiecmdMEAIxd\nrsxpkjQBANZl9a4+GQh10nR5qTgAjEXjfUYv1ElTkv5tz7+O2r6ee2in/qWhctT2J0knUuPjBHq+\neaeq6ke378ZJ10mSXmjeqatGuf/CbrQKhyDOvUjIr1CflDs6nffsL3fq6u+PXt/dPPfhUdvXZ3mO\nPBss9EkTADD2MTwLAIBPriwEcqMeHiXxypm2Q3BW8UL6LhP039mj787e7EX03WgbV5UmJ9DZK+EP\njowU039njXPv7F3k0O88VyrNcZU0AQDhxJwmAAA+5VBpAgDgjyvDs+NqIRAAAJmg0gQAWMecJgAA\nPrkyPEvSBABY50qlyZwmAAA+UWkCAKzjhu0AAPjEnCYAAD7ljPc5zUQiodLSUsXjcTU2Np7y88OH\nD2vZsmWaP3++5s6dq8ceeyyoUAAAGBWBJM1UKqU1a9YokUho7969Wr9+vfbt2zeoTVNTkyoqKvTa\na69py5YtWrt2rU6ePBlEOACAkPNyvIxe2RJI0mxvb1dJSYmKiooUiURUW1ur1tbWQW0uuOACvf/+\n+5Kk999/X5///OeVm8toMQCMR94EL6NXtgSSpXp6elRYWJjejsVi2rFjx6A2q1ev1pVXXqmZM2fq\ngw8+0JNPPnnafT330M701/HKmTzeCwCy5I1XDmn/K4eycqxxvRDI80b+8Pfee6/mz5+vLVu26K23\n3tLSpUu1Z88eTZ06dVC7f2moDCJEAMAILlo0c9AzOZ/95c5hWo8PgQzPRqNRJZPJ9HYymVQsFhvU\n5uWXX9a3v/1tSVJxcbEuvPBCvfHGG0GEAwAIuwleZq8sCSRpVlZWqqOjQ11dXerr69OGDRtUU1Mz\nqE1paak2b94sSert7dUbb7yhWbNmBREOACDkXFkIFMjwbG5urpqamlRVVaVUKqW6ujqVlZWpublZ\nklRfX68f//jHWrVqlcrLyzUwMKCf/vSnys/PDyIcAEDYTXDjjkCeMcbYDmIonufp3/b8q+0whnUi\nFdruCz26bmwL87XqkZD/fp6UG87Ou3nuwwoiZXiepw9vX5bRPqY0JgKJ7bNCfuoAABAeXBgJALDO\nlUeDkTQBAPaN5+s0AQA4I45UmsxpAgDgE5UmAMC6cX0bPQAAzogj12mSNAEA1rlSabqR2gEACAEq\nTQCAfY6sniVpAgDsI2kCAOAPc5oAAIwxVJoAAPsYnh0dYX/0Vv9AuOMLs0jIh2POmRju+KZPitgO\nYVj5k6bYDmFIV5xXYTuEYX103zrbIZzWzUHuPMeNgc/QJ00AwNjHU04AAPAr5CNPn3CjHgYAIASo\nNAEA9jE8CwCAP65cp0nSBADY58hTTtyIEgCAEKDSBADYx5wmAAD+MKcJAIBfjlSazGkCAOATlSYA\nwD6GZwEA8Id7zwIA4JcjTzlxI0oAAEKAShMAYJ8jc5qBVZqJREKlpaWKx+NqbGw8bZstW7aooqJC\nc+fO1ZIlS4IKBQAQdjk5mb2yJJBKM5VKac2aNdq8ebOi0agWLVqkmpoalZWVpdscO3ZMt9xyi55/\n/nnFYjEdPnw4iFAAAC4Yz3Oa7e3tKikpUVFRkSKRiGpra9Xa2jqoze9//3t961vfUiwWkyTNmDEj\niFAAAC7I8TJ7ZUkglWZPT48KCwvT27FYTDt27BjUpqOjQ/39/briiiv0wQcf6Ic//KFuvPHGU/b1\nfPPO9NfFC2eqpHJmECEDAD7jz//3iP5y8IjtMEIlkKTpeSNn/f7+fr366qt68cUXdfz4cV166aX6\nyle+ong8PqhdVX1lECECAEbw9S/n6+tfzk9v37vtzeAO5sjwbCBJMxqNKplMpreTyWR6GPYThYWF\nmjFjhiZPnqzJkyfr61//uvbs2XNK0gQAjAOOJM1AoqysrFRHR4e6urrU19enDRs2qKamZlCba665\nRtu2bVMqldLx48e1Y8cOXXzxxUGEAwAIu/E8p5mbm6umpiZVVVUplUqprq5OZWVlam5uliTV19er\ntLRUy5Yt07x585STk6PVq1eTNAEAoeYZY4ztIIbieZ4eeLXedhjD6h8IbfeFXiTkFzOfMzHc8U2f\nFLEdwrDyJ02xHcKQrjivwnYIw/rovnW2QzitKfclFETK8DxPqZd+mNE+Jlzxi0Bi+yzuCAQAsM+R\nOU2SJgDAvpCPPH3CjdQOAEAIUGkCAOxjeBYAAJ9ImgAA+OPnTnJh4EZqBwAgBKg0AQD2MTwLAIBP\nJE0AAHxy5DpNkiYAwD5HKk03ogQAIASoNAEA9jE8CwCAT44Mz4Y+aaZ48lZGckJ8wXBebnhjk6TC\nc6baDmFYs6Z9wXYIwyrIK7IdwpD6/vfvbYcwrF2/OWA7hOxzJGm6ESUAACEQ+koTADAOMKcJAIBP\njgzPkjQBAPY5kjTdiBIAgBCg0gQA2MecJgAAPjkyPEvSBADY57mRNN2IEgCAEKDSBADY50ilSdIE\nANjHnCYAAD5RaQIA4JMjSdONKAEACAEqTQCAfeO90kwkEiotLVU8HldjY+OQ7V555RXl5ubqqaee\nCioUAEDY5eRk9spWmEHsNJVKac2aNUokEtq7d6/Wr1+vffv2nbbd7bffrmXLlskYnjYNAOOWl5PZ\nK0sCOVJ7e7tKSkpUVFSkSCSi2tpatba2ntLuwQcf1PXXX6/zzz8/iDAAABhVgSTNnp4eFRYWprdj\nsZh6enpOadPa2qqGhgZJkue5cbNeAEAAHKk0A1kI5CcB/uhHP9L9998vz/NkjBlyePaF5p3pr4sX\nzlRx5cxRixMAMLTdJ47rtRMfZedgjiwECiRpRqNRJZPJ9HYymVQsFhvUZteuXaqtrZUkHT58WJs2\nbVIkElFNTc2gdlfVVwYRIgBgBBWT8lQxKS+9ve7994I72HhOmpWVlero6FBXV5dmzpypDRs2aP36\n9YPaHDhwIP31qlWrdPXVV5+SMAEACJNAkmZubq6amppUVVWlVCqluro6lZWVqbm5WZJUX18fxGEB\nAK4a7/eera6uVnV19aDvDZUs161bF1QYAAAHeN4E2yH4wh2BAAD2OTKn6UaUAACEAJUmAMA+RypN\nkiYAwD6SJgAAPo331bMAAPjmSKXpRpQAAIQAlSYAwD5HKk2SJgDAPpImAAA+ObIQyI0oAQAIAZIm\nAMC+AB5CfdNNN6mgoECXXHLJkIe99dZbFY/HVV5ert27d48YJkkTAGBfAElz1apVSiQSQx6yra1N\nb775pjo6OvTII4+ooaFhxDBHTJoHDx5Ua2urtm3bJklqaWkZcacAAJyRAJLmZZddpunTpw95yI0b\nN2rlypWSpMWLF+vYsWPq7e0dNswRFwJ9+OGHamlp0SuvvKIpU6boq1/9qq6//vqR3gYAQGC2bNmj\nLVv3ZLSPnp4eFRYWprdjsZi6u7tVUFAw5HtGTJpPP/20fvGLXyg/P1/Hjx/XSy+9lFGQyK4Jnu0I\nhjZtUrhnB4rOnWE7hGEVTZ1rO4RhRY4csh3CkDo3vmk7hGF9/LHtCCw4w0tOllxRoSVXVKS37/zJ\n42d1WGPM4DC84X9pjpg0CwsLlZ+fL0nKy8s75QAAAGTKWPgDPxqNKplMpre7u7sVjUaHfc+IqX3G\njBmqra3Vs88+qz179mjfvn2ZRwoAwKcYM5DR62zU1NTo8cf/WaFu375d55133rBDs5KPSrO6ulqz\nZs3Sr3/9a23dulWrV68+q+AAAMimFStWaOvWrTp8+LAKCwt15513qr+/X5JUX1+v5cuXq62tTSUl\nJZoyZYrWrVs34j593RHooosu0r333ptZ9AAADGHgLKvF4axfv37ENk1NTWe0T26jBwCwzmj0k2YQ\nSJoAAOuCqDSDQNIEAFjnSqUZ7gvlAAAIESpNAIB1DM8CAODT2V5rmW0kTQCAdQPMaQIAMLZQaQIA\nrGN4FgAAn1gIBACAT+P+Os1EIqHS0lLF43E1Njae8vPf/e53Ki8v17x58/S1r31Nr7/+elChAAAw\nKgKpNFOplNasWaPNmzcrGo1q0aJFqqmpUVlZWbrNrFmz9Oc//1nTpk1TIpHQzTffrO3btwcRDgAg\n5FwZng2k0mxvb1dJSYmKiooUiURUW1ur1tbWQW0uvfRSTZs2TZK0ePFidXd3BxEKAMABRqmMXtkS\nSKXZ09OjwsLC9HYsFtOOHTuGbP+rX/1Ky5cvP+3PXmjemf66eOFMFVfOHL1AAQBD2tN/XK/3f5SV\nY7lSaQaSND3P8932pZde0qOPPqq//vWvp/35VfWVoxUWAOAMlEfyVB7JS2//5vh7FqMJh0CSZjQa\nVTKZTG8nk0nFYrFT2r3++utavXq1EomEpk+fHkQoAAAHuHKdZiBzmpWVlero6FBXV5f6+vq0YcMG\n1dTUDGpz8OBBXXfddfrtb3+rkpKSIMIAADhiIMN/2RJIpZmbm6umpiZVVVUplUqprq5OZWVlam5u\nliTV19frJz/5iY4ePaqGhgZJUiQSUXt7exDhAABCzpVK0zPGGNtBDMXzPP3PXfW2wxjWQHi7T5IU\nyfE/v5xtF0yZYDuEYc0//8u2QxjWrHPn2Q5hWJEjh2yHMKTOb/0v2yEM6619J2yHcFpXHd6vIFKG\n53lKftCc0T4Kp9YHEttnccN2AAB84jZ6AADrXLmNHkkTAGDduL5OEwCAM+FKpcmcJgAAPlFpAgCs\nY3gWAACfXLlOk6QJALAum3f1yQRzmgAA+ESlCQCwjuFZAAB8CvstST9B0gQAWOfK6lnmNAEA8IlK\nEwBgHcOzo8SVjgyrED8ZTOdOnGw7hGGdPzlmO4RhRf7xn7ZDGJY5sN92CEPqTYbz0Vvj2YDc+F0f\n+qQJABj7XCmQSJoAAOtYCAQAwBhDpQkAsI7hWQAAfCJpAgDgE3OaAACMMVSaAADrGJ4FAMAnbm4A\nAIBPzGkCADDGUGkCAKxjThMAAJ9ImgAA+ORK0mROEwAAn6g0AQDWjfvVs4lEQqWlpYrH42psbDxt\nm1tvvVXxeFzl5eXavXt3UKEAAEJuwJiMXtkSSNJMpVJas2aNEomE9u7dq/Xr12vfvn2D2rS1tenN\nN99UR0eHHnnkETU0NAQRCgDAAeM6aba3t6ukpERFRUWKRCKqra1Va2vroDYbN27UypUrJUmLFy/W\nsWPH1NvbG0Q4AICQG8jwX7YEMqfZ09OjwsLC9HYsFtOOHTtGbNPd3a2CgoJB7f79kV3pr2ctvEDF\nC2cGETIA4DP29B/X6/0f2Q4jVAJJmp7n+WpnPlNSn+59S29eOCoxAQDOTHkkT+WRvPT2b46/F9ix\nXLnkJJCkGY1GlUwm09vJZFKxWGzYNt3d3YpGo0GEAwAIOVeSZiBzmpWVlero6FBXV5f6+vq0YcMG\n1dTUDGpTU1Ojxx9/XJK0fft2nXfeeacMzQIAxgdXFgIFUmnm5uaqqalJVVVVSqVSqqurU1lZmZqb\nmyVJ9fX1Wr58udra2lRSUqIpU6Zo3bp1QYQCAMCoCezmBtXV1aqurh70vfr6+kHbTU1NQR0eAOAQ\nV25uwB2BAADWuTKnSdIEAFjnStLkhu0AAPhEpQkAsM6VSpOkCQCwLuVGziRpAgDsc6XSZE4TAACf\nqDQBANYxPAsAgE8DJE0AAPxJOTKnSdIEAFjnSqXJQiAAAHyi0gQAWMdCIAAAfHJleJakmaEcz7Md\nwrAiOeGN79yJk22HMKz8SV+0HcKwTOdfbYcwrIF9SdshDOn4cdsR4LNcWQjEnCYAAD5RaQIArBtw\n4xnUJE0AgH0sBAIAwCdXFgIxpwkAgE9UmgAA61xZPUvSBABY58rwLEkTAGCdKwuBmNMEAMAnKk0A\ngHUDzGkCAOCPK8OzJE0AgHUsBAIAwCdXKk0WAgEA4BOVJgDAugFHxmdJmgAA61wZniVpAgCsc6TQ\nZE4TAAC/qDQBANYxPDtK/v2RXemvZy28QMULZ1qMBgDGjz39x/V6/0dZORZ3BBolS29eaDsEABiX\nyiN5Ko/kpbd/c/y9wI7lSqXJnCYAAD6FvtIEAIx9PIQaAACfUgO2I/CHpAkAsM6VSpM5TQAAfKLS\nBABYl3LklkAkTQCAda4Mz5I0AQDWsRAIAACfXKk0WQgEAIBPVJoAAOtYCAQAgE+uDM+SNAEA1rmy\nEGhczWm+teuQ7RCctf8V+i4TW7a8ZjsEZ239j17bIThrT/9x2yGMOeMqaR7Y9bbtEJzVsZOkmYmt\nW/fYDsFZW//jP22H4KxsPQtzNKSMyeiVLQzPAgCsY04TAACfXFk96xkT3vTueZ7tEAAAnxJEyvA8\nT5c/uSKjfWz97+sDie2zQl1phjifAwBGUcqRX/ehTpoAgPHBleFZkiYAwDpXFgKNq0tOAADIBJUm\nAMA6hmcBAPCJhUAAAPjkSqXJnCYQoKNHj+qGG27QkSNHbIcChJort9EjaQIBmj59uq688kq1tLTY\nDgXAKCBpAgG7+uqr1draajsMINRcqTSZ0wQCVlBQoA8//FDvv/++zj33XNvhAKHE8zQBSJI+/vhj\nnXPOOXruuedshwKEliuVJkkTCFAqldIdd9yhu+66S88884ztcABkiOFZIEBr167V6tWrNWfOHB08\neFB9fX2aOHGi7bCA0OGSE2Cca2lp0cKFCzVnzhxJ0je/+U21tbVZjgoIJ1eGZ0P9PE0AwNjneZ4u\neOjajPbxdsMzWXmcJJUmAAA+MacJALDOODKnSdIEAFhH0gQAwCdXkiZzmgAA+ESlCQCwzpULOUia\nAADrXBmeJWkCAKwjaQIA4JMrSZOFQAAA+ESlCQCwzpVKk6QJALCOpAkAgE+uJE3mNAEA8IlKEwBg\nnSuVJkkTAGAdSRMAAJ9cuY0ec5oAAPhEpQkAsI7hWQAAfCJpAgDgkytJkzlNAAB8otIEAFjnSqVJ\n0gQAWEfSBADAJ5ImAAA+uZI0WQgEAIBPVJoAAOtcuY0eSRMAYB3DswAA+GQGTEavoSQSCZWWlioe\nj6uxsfGUn2/ZskXTpk1TRUWFKioqdPfddw8bJ5UmAGBMSqVSWrNmjTZv3qxoNKpFixappqZGZWVl\ng9pdfvnl2rhxo699UmkCAKwLotJsb29XSUmJioqKFIlEVFtbq9bW1lOPfQbzqVSaAADrzC+3Z/T+\nc84555Tv9fT0qLCwML0di8W0Y8eOQW08z9PLL7+s8vJyRaNR/exnP9PFF1885HFImgAAq4JaOet5\n3ohtFixYoGQyqby8PG3atEnXXnut9u/fP2R7hmcBAGNSNBpVMplMbyeTScVisUFtpk6dqry8PElS\ndXW1+vv7deTIkSH3SdIEAIxJlZWV6ujoUFdXl/r6+rRhwwbV1NQMatPb25uudNvb22WMUX5+/pD7\nZHgWADAm5ebmqqmpSVVVVUqlUqqrq1NZWZmam5slSfX19WppadFDDz2k3Nxc5eXl6Yknnhh2n55x\n5TYMAABYxvAsAAA+kTQBAPCJpAkAgE8kTQAAfCJpAgDgE0kTAACf/j8ZpROSxWvmCwAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4127090>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAFlCAYAAACJCrOeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHlhJREFUeJzt3X9wVdXd7/HPDjmIQURS2gjnxAkksQkoIRCktmONOsiP\nqdGx9k7wjqWYQh4qtZ1hHGf6xx1/a3pb5+mY55EwU7VOW8SmamgNpxavQa1CDCLOLTwQlZSToBkp\nUHn01oSTdf+wnjFATjac7LP2St4v58xkJ+vs/T3LnXz5rrX23p4xxggAAAwrx3YAAAC4gqQJAIBP\nJE0AAHwiaQIA4BNJEwAAn0iaAAD4RNIEAIxKt956qwoKCnTppZee9udtbW2aPHmyKisrVVlZqfvu\nu2/YfeaOdJAAAITBypUr9cMf/lDf/e53h2xz5ZVXavPmzb73SaUJABiVrrjiCk2ZMiVtmzO9vw9J\nEwBgVX5+vjzPy+g1adKkMz6u53l67bXXVFFRoWXLlmnPnj3DvofhWQCAVUePHtVvdHFG+/if/73/\njN8zb948JRIJ5eXlacuWLbrhhhu0f3/6/VBpAgDGpEmTJikvL0+StHTpUvX39+vIkSNp30PSBABY\nl5OT2ets9Pb2puY029vbZYxRfn5+2vcwPAsAsO5sE1/KwKnfWr58ubZt26bDhw+rsLBQd999t/r7\n+yVJ9fX1am5u1qOPPqrc3Fzl5eXpqaeeGvYwHo8GAwDY5Hmefjc+sznN7/TtP+OVsGeD4VkAAHxi\neBYAYF3Gw7NZQtIEAFiX49mOwB+SJgDAOlcqTUfCBADAPipNAIB1rlSaJE0AgHUkTQAAfHIlaToS\nJgAA9lFpAgCsc6XSJGkCAKwjaQIA4BNJEwAAnzzPjVsCOZLbAQCwj0oTAGAdw7MAAPhE0gQAwCdX\nkqYjYQIAYB+VJgDAOlcqTZImAMA6kiYAAD65kjQdCRMAAPuoNAEA1rlSaZI0AQDWkTQBAPApx41b\nzzKnCQCAX1SaAADrGJ4FAMAnkiYAAD65kjQdCRMAAPuoNAEA1uU4snzWqUozHo+rrKxMpaWlamho\nOG2b22+/XaWlpaqoqNCuXbuGfe+RI0e0aNEiXXzxxbr22mt17NgxSVJXV5fOPfdcVVZWqrKyUj/4\nwQ+C/XBZEET//e53v9Ps2bM1btw4vfnmm4P29eCDD6q0tFRlZWV64YUXgvlQWRJE391xxx0qLy9X\nRUWFbrzxRv3jH/+QNPrOvSD67q677lIsFkv10ZYtW1I/G03nnRRM/9XW1qb6bsaMGaqsrJRk99zz\nxnkZvbLGOOLEiROmuLjYHDhwwPT19ZmKigqzZ8+eQW2ef/55s3TpUmOMMdu3bzcLFy4c9r133HGH\naWhoMMYY89BDD5k777zTGGPMgQMHzCWXXJKtjxe4oPpv7969Zt++faa6utrs3Lkzta+//vWvpqKi\nwvT19ZkDBw6Y4uJik0wms/RpR1ZQfffCCy+k+uTOO+8cledeUH131113mZ///OenHG80nXfGBNd/\nX7Ru3Tpz7733GmPsnXuSzP7LZmX0ylY6c6bSbG9vV0lJiYqKihSJRFRbW6uWlpZBbTZv3qwVK1ZI\nkhYuXKhjx47pgw8+SPveL75nxYoVeu6557L7wbIkqP4rKyvTxRdffMrxWlpatHz5ckUiERUVFamk\npETt7e3Bf9AABNV3ixYtUs6/Vj8sXLhQ3d3d2f1gWRBU30nSZ39rBxtN550UbP9Jn/Xh008/reXL\nl2ftMw3Fy/EyemWLM0mzp6dHhYWFqe1YLKaenh5fbQ4dOjTke3t7e1VQUCBJKigoUG9vb6rdgQMH\nVFlZqerqar366quBfK5sCar/hnLo0CHFYrEzek9YZaPvHnvsMS1btiy1PVrOvSD77pFHHlFFRYXq\n6upS0yqj6byTgj/3XnnlFRUUFKi4uDj1vdFy7gXFmaTpef7+JXG6f32ers3p9ud5Xur706dPVyKR\n0K5du/Twww/r5ptv1vHjx88s6BAZyf4LOoawCbrv7r//fo0fP14333yzpNF17gXVd2vWrNGBAwf0\n1ltvadq0aVq3bl3GMYRR0Ofexo0bU+edZPfcc2VO05nVs9FoVIlEIrWdSCQG/YvydG26u7sVi8XU\n399/yvej0aikz6rLDz74QBdeeKHef/99feUrX5EkjR8/XuPHj5ckzZs3T8XFxers7NS8efMC+4xB\nGsn+O917hzveF/vcNUH23RNPPKHW1la9+OKLqe+NpnMvqL77/PdUkr7//e/ruuuuG3Jfrp53UrDn\n3okTJ/Tss88OWsBn89zL5hBrRrIyczoC+vv7zcyZM82BAwfMp59+OuyE+Ouvv56aEE/33jvuuMM8\n9NBDxhhjHnzwwdRijA8//NCcOHHCGGPMu+++a6LRqDl69GhWPmsQguq/z1VXV5uOjo7U9ucLMj79\n9FPz3nvvmZkzZ5qBgYGAP2Uwguq7LVu2mFmzZpkPP/xw0L5G07kXVN8dOnQo9f6HH37YLF++3Bgz\nus47Y4L9vd2yZYuprq4etC9b554k8171nIxe2UpnziRNY4xpbW01F198sSkuLjYPPPCAMcaY9evX\nm/Xr16fa3Hbbbaa4uNjMmTNn0GrO073XGGP+/ve/m2uuucaUlpaaRYsWpU6Q3//+92b27Nlm7ty5\nZt68eeaPf/xjlj5lcILov2eeecbEYjEzYcIEU1BQYJYsWZL62f3332+Ki4vNV7/6VROPx7PwCYMT\nRN+VlJSYiy66yMydO9fMnTvXrFmzxhhjTHNz86g694Lou1tuucVceumlZs6cOeb66683H3zwQepn\no+m8MyaY/jPGmO9973umqalp0Pds/d1zKWl6/woYAAArPM9T19UVGe2j6P/sDnRNxuecmdMEAIxe\nrsxpkjQBANZl9a4+GQh10nR5qTgAjEZjfUYv1ElTkv59V/2I7Su+vkNL/q1qxPYnSf0DI7q7ETcw\nQif4nzfs1KLV80dkX2NREP2XM0b+UflCU4eurR/Z39sJ40Z0dyMuLzIyl9C3/Mcbuv62BSOyL0mq\nm/3oiO3rZJ4jzwYLfdIEAIx+DM8CAOCTKwuB3KiHR0hJ1XTbIThr5vxptkNwGv139orn83t7tr66\ngL4baWOq0iRpnj3+cGWG/jt7xfzenrWyy9y5haArleaYSpoAgHBiThMAAJ9yqDQBAPDHleHZMbUQ\nCACATFBpAgCsY04TAACfXBmeJWkCAKxzpdJkThMAAJ+oNAEA1nHDdgAAfGJOEwAAn3LG+pxmPB5X\nWVmZSktL1dDQcMrPDx8+rCVLlmju3Lm65JJL9MQTTwQVCgAAIyKQpJlMJrV27VrF43Ht2bNHGzdu\n1N69ewe1aWxsVGVlpd566y21tbVp3bp1OnHiRBDhAABCzsvxMnplSyBJs729XSUlJSoqKlIkElFt\nba1aWloGtZk2bZo++ugjSdJHH32kL33pS8rNZbQYAMYib5yX0StbAslSPT09KiwsTG3HYjHt2LFj\nUJtVq1bp6quv1vTp03X8+HE9/fTTp91XfH1H6uuSquk83gsAsuS/2nu0741DWTnWmF4I5HnDf/gH\nHnhAc+fOVVtbm959910tWrRIu3fv1qRJkwa1W/JvVUGECAAYRtll0UHP5Nz8nx1pWo8NgQzPRqNR\nJRKJ1HYikVAsFhvU5rXXXtN3vvMdSVJxcbFmzJihffv2BREOACDsxnmZvbIkkKRZVVWlzs5OdXV1\nqa+vT5s2bVJNTc2gNmVlZdq6daskqbe3V/v27dPMmTODCAcAEHKuLAQKZHg2NzdXjY2NWrx4sZLJ\npOrq6lReXq6mpiZJUn19vX7yk59o5cqVqqio0MDAgH76058qPz8/iHAAAGE3zo07AnnGGGM7iKF4\nnqd/31VvO4y0+gdsR5DeQHj/9yJDOT7WDuD0JoyzHUF6eZFwJpC62Y8qiJTheZ4+vnNJRvuY2BAP\nJLaThfP/DAAAIcSFkQAA61x5NBhJEwBg31i+ThMAgDPiSKXJnCYAAD5RaQIArBvTt9EDAOCMOHKd\nJkkTAGCdK5WmG6kdAIAQoNIEANjnyOpZkiYAwD6SJgAA/jCnCQDAKEOlCQCwj+HZkRH2R2/1D4T7\n0VthPg8jIR+OmTIh3AMx+edMsB1CWheck2c7hCFdMa3adghpmX2v2w7htOqC3HlOuH/fPhf6pAkA\nGP14ygkAAH6FfOTpc27UwwAAhACVJgDAPoZnAQDwx5XrNEmaAAD7HHnKiRtRAgAQAlSaAAD7mNME\nAMAf5jQBAPDLkUqTOU0AAHyi0gQA2MfwLAAA/nDvWQAA/HLkKSduRAkAQAhQaQIA7HNkTjOwSjMe\nj6usrEylpaVqaGg4bZu2tjZVVlbqkksuUXV1dVChAADCLicns1eWBFJpJpNJrV27Vlu3blU0GtWC\nBQtUU1Oj8vLyVJtjx47ptttu05/+9CfFYjEdPnw4iFAAAC4Yy3Oa7e3tKikpUVFRkSKRiGpra9XS\n0jKozW9/+1t9+9vfViwWkyRNnTo1iFAAAC7I8TJ7ZUkglWZPT48KCwtT27FYTDt27BjUprOzU/39\n/brqqqt0/Phx/ehHP9Itt9xyyr5eaOpIfV08f7qKq6YHETIA4CRt7Qe07Y0u22GESiBJ0/OGz/r9\n/f1688039eKLL+qTTz7R5Zdfrq997WsqLS0d1O7a+qogQgQADKP6shmqvmxGavue/3gpuIM5Mjwb\nSNKMRqNKJBKp7UQikRqG/VxhYaGmTp2qc889V+eee66++c1vavfu3ackTQDAGOBI0gwkyqqqKnV2\ndqqrq0t9fX3atGmTampqBrW5/vrr9eqrryqZTOqTTz7Rjh07NGvWrCDCAQCE3Vie08zNzVVjY6MW\nL16sZDKpuro6lZeXq6mpSZJUX1+vsrIyLVmyRHPmzFFOTo5WrVpF0gQAhJpnjDG2gxiK53n63zvr\nbYeRVv9AaLtPUrifthMJ+cXMUyaEe7go/5wJtkNI64Jz8myHMKQrplXbDiEts+912yGc1rhZ/0tB\npAzP85R86UcZ7WPcVb8IJLaTcUcgAIB9jsxpkjQBAPaFfOTpc26kdgAAQoBKEwBgH8OzAAD4RNIE\nAMAfP3eSCwM3UjsAACFApQkAsI/hWQAAfCJpAgDgkyPXaZI0AQD2OVJpuhElAAAhQKUJALCP4VkA\nAHxyZHg29ElzILxPLnNCTogvGJ58Trh/SWLnnW87hLQumjTVdghpTTnnQtshDMl07rAdQnp9/bYj\nyD5HkqYbUQIAEAKhrzQBAGMAc5oAAPjkyPAsSRMAYJ8jSdONKAEACAEqTQCAfcxpAgDgkyPDsyRN\nAIB9nhtJ040oAQAIASpNAIB9jlSaJE0AgH3MaQIA4BOVJgAAPjmSNN2IEgCAEKDSBADYN9YrzXg8\nrrKyMpWWlqqhoWHIdm+88YZyc3P1zDPPBBUKACDscnIye2UrzCB2mkwmtXbtWsXjce3Zs0cbN27U\n3r17T9vuzjvv1JIlS2R42DQAjF1eTmavLAnkSO3t7SopKVFRUZEikYhqa2vV0tJySrtHHnlEN910\nk7785S8HEQYAACMqkKTZ09OjwsLC1HYsFlNPT88pbVpaWrRmzRpJkue5cbNeAEAAHKk0A1kI5CcB\n/vjHP9ZDDz0kz/NkjBlyePbPG3amvp45f5qK508fsTgBAENr6+jSto6/ZedgjiwECiRpRqNRJRKJ\n1HYikVAsFhvUZufOnaqtrZUkHT58WFu2bFEkElFNTc2gdotWzw8iRADAMKqrilRdVZTavmf9y8Ed\nbCwnzaqqKnV2dqqrq0vTp0/Xpk2btHHjxkFt3nvvvdTXK1eu1HXXXXdKwgQAIEwCSZq5ublqbGzU\n4sWLlUwmVVdXp/LycjU1NUmS6uvrgzgsAMBVjtx71jMhvtbD8zw1dKy2HUZaydD23mciIX4aev6E\ncP+SXDTpfNshpHXRpKm2Q0hryjkX2g5hSF86mKV5urP1z09tR3Ba4+beG8jlgZ7nyRz7bWb7uODm\nrFy6yB2BAAD2OTKn6UaUAACEAJUmAMA+RypNkiYAwD6SJgAAPjmyepakCQCwz5FK040oAQAIASpN\nAIB9jlSaJE0AgH0kTQAAfHJkIZAbUQIAEAIkTQCAfQE8hPrWW29VQUGBLr300iEPe/vtt6u0tFQV\nFRXatWvXsGGSNAEA9gWQNFeuXKl4PD7kIVtbW/XOO++os7NTGzZs0Jo1a4YNc9ikefDgQbW0tOjV\nV1+VJDU3Nw+7UwAAzkgASfOKK67QlClThjzk5s2btWLFCknSwoULdezYMfX29qYNc9iFQB9//LGa\nm5v1xhtvaOLEifr617+um266abi3AQAQmLa23WrbtjujffT09KiwsDC1HYvF1N3drYKCgiHfM2zS\nfPbZZ/WLX/xC+fn5+uSTT/TSSy9lFORoMy68j6uUFO74zh8f7sXb0ycO/S/UMCjIK7IdQlrnmwm2\nQxiS6XvHdgg42RleclJ9VaWqr6pMbd99z5NnddiTn8Hpeen/aA77V6uwsFD5+fmSpLy8vKw85BMA\nMLYYC//Aj0ajSiQSqe3u7m5Fo9G07xk2tU+dOlW1tbX6wx/+oN27d2vv3r2ZRwoAwBcYM5DR62zU\n1NToySc/q1C3b9+uCy64IO3QrOSj0ly6dKlmzpypX/3qV9q2bZtWrVp1VsEBAJBNy5cv17Zt23T4\n8GEVFhbq7rvvVn9/vySpvr5ey5YtU2trq0pKSjRx4kQ9/vjjw+7TMyEeb/U8Tw0dq22H4bRITngn\nNQsnjbcdQlqz8tMP09gWPa/YdghphXpOc3+H7RDSGzi7yilo4+beG8gUned56k++kNE+IuOuzcr0\nYbhXYgAAxgSjcP5D4WQkTQCAdQNnOS+ZbSRNAIB1rlSa3EYPAACfqDQBANYxPAsAgE9ne61ltpE0\nAQDWDTCnCQDA6EKlCQCwjuFZAAB8YiEQAAA+jfnrNOPxuMrKylRaWqqGhoZTfv6b3/xGFRUVmjNn\njr7xjW/o7bffDioUAABGRCCVZjKZ1Nq1a7V161ZFo1EtWLBANTU1Ki8vT7WZOXOmXn75ZU2ePFnx\neFyrV6/W9u3bgwgHABByrgzPBlJptre3q6SkREVFRYpEIqqtrVVLS8ugNpdffrkmT54sSVq4cKG6\nu7uDCAUA4ACjZEavbAmk0uzp6VFhYWFqOxaLaceOHUO2/+Uvf6lly5ad9md/3rAz9fXM+dNUPH/6\nyAUKABhSW0eXtnX8LSvHcqXSDCRpep7/Zzi+9NJLeuyxx/SXv/zltD9ftHr+SIUFADgD1VVFqq4q\nSm3fs/5le8GERCBJMxqNKpFIpLYTiYRisdgp7d5++22tWrVK8XhcU6ZMCSIUAIADXLlOM5A5zaqq\nKnV2dqqrq0t9fX3atGmTampqBrU5ePCgbrzxRv36179WSUlJEGEAABwxkOF/2RJIpZmbm6vGxkYt\nXrxYyWRSdXV1Ki8vV1NTkySpvr5e99xzj44ePao1a9ZIkiKRiNrb24MIBwAQcq5Ump4xxtgOYiie\n56mhY7XtMJwWyfE/v5xthZPG2w4hrVn5UdshpBU9r9h2CGmdbybYDmFIZn+H7RDSGwhnAhk3914F\nkTI8z1PieFNG+yicVB9IbCfjhu0AAPjEbfQAANa5chs9kiYAwLoxfZ0mAABnwpVKkzlNAAB8otIE\nAFjH8CwAAD65cp0mSRMAYF027+qTCeY0AQDwiUoTAGAdw7MAAPg0EN47ug5C0gQAWOfK6lnmNAEA\n8IlKEwBgHcOzIyQZ8n4M86O3JCkS4rGE88efazuEtKac82XbIaR1vvJsh5De0YO2IxhaSB+9NZYN\nKOR/7P8l9EkTADD6UWkCAOATC4EAABhlqDQBANYxPAsAgE8kTQAAfGJOEwCAUYZKEwBgHcOzAAD4\nxM0NAADwiTlNAABGGSpNAIB1zGkCAOATSRMAAJ9cSZrMaQIA4BOVJgDAujG/ejYej6usrEylpaVq\naGg4bZvbb79dpaWlqqio0K5du4IKBQAQcgPGZPTKlkCSZjKZ1Nq1axWPx7Vnzx5t3LhRe/fuHdSm\ntbVV77zzjjo7O7VhwwatWbMmiFAAAA4Y00mzvb1dJSUlKioqUiQSUW1trVpaWga12bx5s1asWCFJ\nWrhwoY4dO6be3t4gwgEAhNxAhv9lSyBzmj09PSosLExtx2Ix7dixY9g23d3dKigoGNTuxQ07U1/P\nmD9NM+dPDyJkAMBJ2jq6tK3jb7bDCJVAkqbneb7amZNK6tO975rV80ckJgDAmamuKlJ1VVFq+571\nLwd2LFcuOQkkaUajUSUSidR2IpFQLBZL26a7u1vRaDSIcAAAIedK0gxkTrOqqkqdnZ3q6upSX1+f\nNm3apJqamkFtampq9OSTT0qStm/frgsuuOCUoVkAwNjgykKgQCrN3NxcNTY2avHixUomk6qrq1N5\nebmampokSfX19Vq2bJlaW1tVUlKiiRMn6vHHHw8iFAAARkxgNzdYunSpli5dOuh79fX1g7YbGxuD\nOjwAwCGu3NyAOwIBAKxzZU6TpAkAsM6VpMkN2wEA8IlKEwBgnSuVJkkTAGBd0o2cSdIEANjnSqXJ\nnCYAAD5RaQIArGN4FgAAnwZImgAA+JN0ZE6TpAkAsM6VSpOFQAAA+ESlCQCwjoVAAAD45MrwLElz\nlIvkeLZDGNJ5kQm2Q0hr0vh82yGkd/yw7QjSMv/vI9shwCGuLARiThMAAJ+oNAEA1g248QxqkiYA\nwD4WAgEA4JMrC4GY0wQAwCcqTQCAda6sniVpAgCsc2V4lqQJALDOlYVAzGkCAOATlSYAwLoB5jQB\nAPDHleFZkiYAwDoWAgEA4JMrlSYLgQAA8IlKEwBg3YAj47MkTQCAda4Mz5I0AQDWOVJoMqcJAIBf\nVJoAAOsYnh0hL27Ymfp6xvxpmjl/usVoAGDsaOvo0raOv2XlWNwRaIRcs3q+7RAAYEyqripSdVVR\navue9S8HdixXKk3mNAEA8Cn0lSYAYPTjIdQAAPiUHLAdgT8kTQCAda5UmsxpAgDgE5UmAMC6pCO3\nBCJpAgCsc2V4lqQJALCOhUAAAPjkSqXJQiAAAHyi0gQAWMdCIAAAfHJleJakCQCwzpWFQGNqTvO9\nnYdsh+Cs/W/Qd5l4Zdv/tR2Cs9q2v2s7BGe1dXTZDmHUGVNJ88DO922H4Kz9HSTNTLzyMknzbG3b\n8Z7tEJyVrWdhjoSkMRm9soXhWQCAdcxpAgDgkyurZz1jwpvePc+zHQIA4AuCSBme5+nKp5dntI9t\n/2NjILGdLNSVZojzOQBgBCUd+XMf6qQJABgbXBmeJWkCAKxzZSHQmLrkBACATFBpAgCsY3gWAACf\nWAgEAIBPrlSazGkCATp69KhuvvlmHTlyxHYoQKi5chs9kiYQoClTpujqq69Wc3Oz7VAAjACSJhCw\n6667Ti0tLbbDAELNlUqTOU0gYAUFBfr444/10Ucf6fzzz7cdDhBKPE8TgCTpn//8p8477zw9//zz\ntkMBQsuVSpOkCQQomUzqrrvu0r333qvnnnvOdjgAMsTwLBCgdevWadWqVZo9e7YOHjyovr4+jR8/\n3nZYQOhwyQkwxjU3N2v+/PmaPXu2JOlb3/qWWltbLUcFhJMrw7Ohfp4mAGD08zxP0x69IaN9vL/m\nuaw8TpJKEwAAn5jTBABYZxyZ0yRpAgCsI2kCAOCTK0mTOU0AAHyi0gQAWOfKhRwkTQCAda4Mz5I0\nAQDWkTQBAPDJlaTJQiAAAHyi0gQAWOdKpUnSBABYR9IEAMAnV5Imc5oAAPhEpQkAsM6VSpOkCQCw\njqQJAIBPrtxGjzlNAAB8otIEAFjH8CwAAD6RNAEA8MmVpMmcJgAAPlFpAgCsc6XSJGkCAKwjaQIA\n4BNJEwAAn1xJmiwEAgDAJypNAIB1rtxGj6QJALCO4VkAAHwyAyaj11Di8bjKyspUWlqqhoaGU37e\n1tamyZMnq7KyUpWVlbrvvvvSxkmlCQAYlZLJpNauXautW7cqGo1qwYIFqqmpUXl5+aB2V155pTZv\n3uxrn1SaAADrgqg029vbVVJSoqKiIkUiEdXW1qqlpeXUY5/BfCqVJgDAOvOf2zN6/3nnnXfK93p6\nelRYWJjajsVi2rFjx6A2nufptddeU0VFhaLRqH72s59p1qxZQx6HpAkAsCqolbOe5w3bZt68eUok\nEsrLy9OWLVt0ww03aP/+/UO2Z3gWADAqRaNRJRKJ1HYikVAsFhvUZtKkScrLy5MkLV26VP39/Tpy\n5MiQ+yRpAgBGpaqqKnV2dqqrq0t9fX3atGmTampqBrXp7e1NVbrt7e0yxig/P3/IfTI8CwAYlXJz\nc9XY2KjFixcrmUyqrq5O5eXlampqkiTV19erublZjz76qHJzc5WXl6ennnoq7T4948ptGAAAsIzh\nWQAAfCJpAgDgE0kTAACfSJoAAPhE0gQAwCeSJgAAPv1/jOggXtRhE6wAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x410bb10>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also visualize the reliability of the PDD predicted from the SFM model parameters in each combination of regularization parameters:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for b in [1000, 2000, 4000]:\n",
      "    fig, ax = plt.subplots(1)\n",
      "    cax = ax.imshow(rel_matrix[b], interpolation='nearest', cmap=cm.OrRd, vmax=90, vmin=0)\n",
      "    cbar = fig.colorbar(cax)\n",
      "    cax.axes.set_xticks([0,1,2,3,4,5,6,7,8])\n",
      "    cax.axes.set_xticklabels([str('%1.4f'%this) for this in alphas])\n",
      "    cax.axes.set_xlabel(r'$\\lambda$')\n",
      "    cax.axes.set_yticks([0,1,2,3,4,5])\n",
      "    cax.axes.set_yticklabels([str('%1.1f'%(1-this)) for this in l1_ratios])\n",
      "\n",
      "    cax.axes.set_ylabel(r'$\\alpha$')\n",
      "    #ax.matshow(rmse_matrix[b], cmap=matplotlib.cm.hot)\n",
      "    fig.set_size_inches([8,6])\n",
      "    fig.savefig('figures/AppendixA_rel_b%s.svg'%b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}