{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# MKS Introduction\n",
      "\n",
      "The following notebook is an introductory tutorial describing the calculation of influence coefficients for the MKS in both real and Fourier space. In the MKS, the goal is to quickly calculate a response $p_{a,s}$ from a microstructure $m_{a,s}$ where $a$ represents different samples and $s$ represents the spatial position. The first step in the MKS is to bin (or discretize) the state space, $m_{a, s}$, in terms of a new representation, $m_{a, s}^h$, given by,\n",
      "\n",
      "$$ m_{a,s} = \\sum\\limits_{h=0}^{H-1} m^h_{a, s} \\chi^h $$\n",
      "\n",
      "and\n",
      "\n",
      "$$ \\sum\\limits_{h=0}^{H-1} m^h_{a, s} = 1 $$\n",
      "\n",
      "where $\\chi^h$ is the basis representation of the state space (microstructure space). $H$ is the size of the state space discretization.\n",
      "\n",
      "Why do we discretize the microstructure?\n",
      "\n",
      "Can we represent \n",
      "\n",
      "When the state space is discretized, the relationship between the response and microstructure can be written as,\n",
      "\n",
      "$$ p_{a,s} = \\sum\\limits_{h=0}^{H - 1} \\sum\\limits_{t=0}^{S-1} \\alpha_t^h m_{a,s + t}^h $$\n",
      "\n",
      "where the $p_{a,s}$ are the responses, $\\alpha_t^h$ are the influence coefficients and $m_{a,s}^h$ is the microstructure. $S$ is the size of the spatial discretization The $s$ and $t$ indices run over the spatial discretization, the $a$ index runs over the number of samples and the $h$ index runs over the microstructure discretization. The goal of the notebook is to introduce the method for calculating the $\\alpha_t^h$ by\n",
      "\n",
      " * creating a response function using the Cahn-Hilliard equation,\n",
      " \n",
      " * demonstrate Cahn-Hilliard evolution,\n",
      " \n",
      " * create some sample data,\n",
      " \n",
      " * solve MKS regression in real space for demonstration purposes,\n",
      " \n",
      " * solve MKS regression in Fourier space,\n",
      " \n",
      " * test that MKS regression seems to work with a single test sample and\n",
      " \n",
      " * show that the calculated coefficients are not quite right."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Questions for Tony\n",
      "\n",
      "Why do we discretize the microstructure?\n",
      "\n",
      "Which type of non-linearities does the first order MKS handle?\n",
      "\n",
      "Can we represent any PDE in one variable with a high enough order MKS relationship?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Calculating the Response"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the MKS a sample set of microstructures and responses are required. In this example we will use the Cahn-Hilliard equation to provide the example response. [FiPy](http://www.ctcms.nist.gov/fipy/) is used to solve the governing equation, which is given by,\n",
      "\n",
      "$$ \\frac{\\partial \\phi}{\\partial t} = \\nabla \\cdot D \\nabla \\left( \\frac{\\partial f}{\\partial \\phi}   - \\epsilon^2 \\nabla^2 \\phi \\right).$$\n",
      "\n",
      "where the free energy is given by,\n",
      "\n",
      "$$ f = (a^2/2) \\phi^2 (1 - \\phi)^2 $$\n",
      "\n",
      "In this example $D = 1$, $a = \\sqrt{200}$ and $\\epsilon=0.1$. See [the FiPy CH example](http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html) for further details.\n",
      "\n",
      "The `fipy_response` function takes an initial field, $\\phi_0$ (the microstructure), and yields $\\frac{\\partial \\phi}{\\partial t}$ (the response). \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "%load_ext autoreload\n",
      "%autoreload 2\n",
      "\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import fipy as fp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The autoreload extension is already loaded. To reload it, use:\n",
        "  %reload_ext autoreload\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fipy_response(phi0, dt, N):\n",
      "    from fipy.solvers.scipy.linearLUSolver import LinearLUSolver\n",
      "    dx = 0.005\n",
      "    a = np.sqrt(200.)\n",
      "    epsilon = 0.1\n",
      "    nx = ny = N\n",
      "    mesh = fp.PeriodicGrid2D(nx=nx, ny=ny, dx=dx, dy=dx)\n",
      "    phi = fp.CellVariable(name=r\"$\\phi$\", mesh=mesh, value=phi0.copy())\n",
      "    PHI = phi.arithmeticFaceValue\n",
      "    D = 1.\n",
      "    eq = (fp.TransientTerm()\n",
      "      == fp.DiffusionTerm(coeff=D * a**2 * (1 - 6 * PHI * (1 - PHI)))\n",
      "      - fp.DiffusionTerm(coeff=(D, epsilon**2)))\n",
      "    \n",
      "    eq.solve(phi, dt=dt, solver=LinearLUSolver())\n",
      "    \n",
      "    return (np.array(phi) - phi0) / dt\n",
      "    \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Demonstrate Cahn-Hilliard Evolution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following cell iterates the `fipy_response` function to demonstrate the evolution of the microstructure for an initially uniform random field. Using the `fipy_response` function is quite an inefficient method of using [FiPy](http://www.ctcms.nist.gov/fipy/), but useful for demonstration purposes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import clear_output\n",
      "import time\n",
      "\n",
      "N = 21\n",
      "np.random.seed(0)\n",
      "phi0 = np.random.random(N * N)\n",
      "dt = 1e-8\n",
      "\n",
      "fig = plt.figure()\n",
      "\n",
      "for i in range(30):\n",
      "    response = fipy_response(phi0, dt=dt, N=N)\n",
      "    #Euler backward in this case since implicit\n",
      "    phi0 = response * dt + phi0\n",
      "    #print phi0\n",
      "    plt.contourf(phi0.reshape((N,N)))\n",
      "    time.sleep(1)\n",
      "    clear_output()\n",
      "    display(fig)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD7CAYAAABzGc+QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGLRJREFUeJzt3U9vm1d2x/EfnUDdyYzkbTsDOt4aMuxsZhxMKkp5Azai\nZG8KKOBdbFANECDTAHEU2zsDbUTvAxviK4gIjJE0XVSRAW+NUJ0CXYp2BHQjjKMu3IemJIp8/tx7\nzzn3/j5AAEuxqEfSw6+OD//VDg8PD0FEROackT4AIiIqhwEnIjKKASciMooBJyIy6u1pf6HT6QAA\nfvnlF3z99dcAgG63i3q9jn6/j1ar5fcIiYhorIkTeK/Xw9LSElqtFvr9Pnq9Hp4+fQoAaDabADB8\nm4iIwpoY8H6/j62tLQBAo9FAv9/Ho0ePUK/Xh+/L/j8REYU1cYUyuh7Z2dnBysoKfv75Z8zPzw/f\nv7e35+/oiIjoVLluxNzZ2cHly5dx6dIlAAAf+0NEJG/qjZjA6134nTt3AAD1eh2DwQAA8OLFiyPT\neKZWqzk8RCKidBQZkKcGfGNjA7dv3wbwOuQrKyvY3t5Gs9nE7u4ulpeXx37cr3+byX0Qk8y2D5xc\njmVf/AR88Qfpo3hjf93Nz/Y099665eyy1gftE+97tX4Hb7X/ufBlHdycdXFIk333zP/ncOpfAfyT\nm4v65KKbyxlj5sG+08trz607vbzMl7WvCv39iSuUra0trK2t4d1338Xc3BxqtdpwjdLr9VCv17Gw\nsFD+aInoKI8RU83Y1z1uMJAwcQJfWloarktGZTduZnclpHRYn75N+OSiwUlct4Obs86ncA34SEwD\nPvh76SMIw2W8J6n98WqQz1OJmYn0ivQBiNEwIKgPuO+Jz4JUAu7SpCvXmavvBzySCkxE/D3pA8jN\nx20Y0hFXH3DSw+cv01DTtzkmIl6R8a9RMuIMOIlzHW/pqci5Ty6aj5wWQe5JFBADTrlwlaUAI66W\n1NDAgJMoTt8FxRZxga/H1xQuce4x4DQVp29lYot4REJHnAEnMZy+K7C+Fxc+fp+78JDnYa7nQiFy\nzdq9TtTe+HU8gpofAKTsF47PB/esD9reHm4/ykTA99dn+JwoQnysT3zEO6npexKNQVcW7lBCRNxE\nwEkG4x0BqaAbibbvh9j7jriZHThvSKNxGO+Cst2zr8Aa3M37Xo/5PEc5gXv27D5w8VPpoyjOwvTN\neFeUN7TTpnZjwR7H6iReO/Tw8jq1Ws3Z84EfZ2kX/uz+mz9bi7jrgFuOt9obMMk5389YOC3iX9a+\nKvSCDmZWKBkrq5TReI97WzPGm1JlbZ1iLuAWnBbrZ/f1h9zKL0giXyxF3GTANUcmT6C1R9wlTt9k\nUYiIuziXvQXc2gM1XCgSZo3TOFcnJ3H/na4QP/uq57TXCdxnxDVP4RRHvIlCqHJue1+hpDiJW+Ty\nF2Is8eb0TaHOgbLnuMkdOKWD8SZpmiMeJOC+pnCuUdzQOn0z3qSF1nOCEzg5w3hTzDSeG8ECzl34\nSRoenanxXzGMN2ml7RwJOoH7iLjGAKXI1c9W6q6C2q6YpJemc4UrlIS5+uVn+V9Xmq6MZIeW8yZ4\nwC1f2ekky3tvLVdCsknD+RPFBG5xjSK9/9b2PWO8ySLp80gk4JzC42B17y19paO4SJ5PYhM4Iy7H\nxfRtMd68sZJ8kTqvolihAPpWAjQdX5SBYiJxjkUTcEsk999apm/Gm6g6Btwh6RsmrQgVb65MKHai\nAece3J6qP7OQ8SaKHV+VPiFV1ycW4s1wU0q4QgmMaxZ/GG+SFvoc5AROuWievhluShUn8ERI3s2S\n8SbygxM4TaXxxmaGm4gBD0pq/x3b9J16vGce7I99f+rflxSJr1A0TndVxHYjpbafT8qRmnmwf2q8\nR///tL9HfoU8R8UDfuvVPelDiFpM03fq8S7zMQx5eCG/5+IBJ700Td8px7sqhjxeDHggsa1WppF6\nXUs6HUPuX+jvLwPugZZYx/IMjZy+3WLI4yEacO6/9aqyPnE5fTPe/jDibkl8P3MFvN1uj3270+m4\nPyIiCoYRd0Pq+zg14BsbG+h2u0fe1+l0cOHCBZw/f97bgRU12z6QPoQjRtcoEiuVGNYnnL7JAslf\nglMDvrq6ikajceR9nU4Hz58/x+LiorcDIzla1icUBqfw8qS/d6V24IPBAL1eD3fv3i39ibn/pkk4\nfYclHSIqp9RD6VutFgDg+++/R6/XQ7PZdHpQsUjtofOupm/GmyzQ8Euv8ATe6XSGO/H5+Xn0+33n\nB0VyND14h8LSECQrtHyvCk/gjUYDV65cAQDs7e1heXl57N978sUPwz//7oN/wO8/+N3wba5P4sPp\nm1LgOty//fgDDv/9x9IfPzXgm5ub2N7exsOHD3Hjxg00m83hBH7u3DksLCyM/bg/ffF+6YOi8mK4\n9wmddHBz1vvUx1+ek/n4/p+5+j5w9U0rX33zdaGPnxrw69ev4/r160fed+3atUKfZJSP6VvbXQit\nkl6fMCCkkZZ1yTh8KD2REfwFF57meAOBA87dt1/W731C0/mKOH85HGXl+WL4ijykAgNCGliI9qhg\nEzinb7+qTt/S+2/Kz/UvO/7yfE063u259cIfwwmcyKDR6EqHxzrp71+ZcGeCBJzTd7y4/5Y3boKW\njpIFGr5HVeINBAi473jzLoT27/vNf8K7lzfqqX3vNUQ7UzXeAFcoBO6/U5FarEdpCjfgJt6A54Bz\ndeKf9embyCdt4QbcxRuIYALfX5/hGkWIi/13ylMh+aEx2oDbcGe8BZzTt38upm+uTygWKYU7E8VD\n6VNcI6T4NRONo/lRkz7jDUSwQkmRq3hz+iartAY74zvcmWgCnsouXMvkzft/U2jao50JFW8gooBT\nMZy+yQIr0QbChjsTVcBjn8K1TN9EPlmKNiAT7kxUAY+Zy3hz+iZtrEU7IxlvIMKAxziFa5u8uf8m\nF6xGG5APdya6gAPxRFxbuImqsBzsjJZwZ6IMuHU+w831CYUUQ7QBfeHORBvw0QhamMZDTNsu4s31\nCU0SS7AzWsOdiTbgo8bFUUPUuSKhGMQWbUB/uDNJBHyc4/EMEXSpYHNtQj4w3PKSDfhxeaZ0ixOz\ny3hzfUIAw60JAz6BxWBnOHWTawy3Pgx4hHzEm9N3uhhuvRjwiPiauhnvNDHcYYy+dsKXBT+WAY8E\nVybkCsMdhosXvWHAjfMdbk7f6WC4w3D5amUMuGGcusmV2OIde7gzDLhBocLN6TtusUUb0Bdu368N\nXDs8PDx0fqG1Gn79m9274GkVcuIOGW++Mn1YDLdfVaJ99u0DFEkyJ3DlJNYknLzjw2j75XvSPg0n\ncCU07LMlw80p3A+G2x8f0eYErpyGUI/DqTseMUYb0BFuqUn7NAx4RVqDXATjHQeG2x9t4c4kvUKJ\nIb5VaAs31yjlMNx+SESbKxQwzHloizcVE2u0Adlwa520T2N+Amesi9Eebk7hkzHcfmgJd/QTOINd\nnvZ40+kYbj+0hLss1RM4Y+2GtXBzCn+D4XZPc7SLTuAqA85wV2Mt2OOkHPGYow0w3JOYDjjDXU4M\nwT5N7CGPPdajJMJtIdqjot+Bpy7mWI8z82A/moinFOsMp22/1EzgnL7Ti3NRVkKeYqiP47RdjskV\nSozxZozDCh13Rno8hrsacwG3Em8GOQ6joWeE3WG43fAS8Ha7jfX1Nz+gbreLer2Ofr+PVqt18kJz\nBlxjvBlqonwYbfec34i5sbGBbrc7DPjOzg4AoNlsot/v4+nTp7h06VLhA9USbwabqBiGW4+pAV9d\nXcXm5ubw7cePH+PDDz8EADQaDWxtbRUOuFS8GWui8hhufQrfjfDly5eYm5sbvr23t1fo42N9WTCi\nWDHcJ822D079f/vr4Z6JtdT9wMve7skX4yWyIeX7b0+Kc9GP9x3zwgGv1+sYDAYAgBcvXmB+fj7X\nx4WIN8NNVE1q4a4a6yKX7yPmhQO+srKC7e1tNJtN7O7uYnl5eezfu/Mvr4Z//p9//Bi//6D0MU7F\ncBNVk9KaxHe0J33e4xH/4clv+PFJ+XtyT70b4ebmJlZXV/HNN9/gxo0bAIBOp4NGo5H7boQ+p2/G\nm6iclKZtqWiPM2kSV/lAHl8BZ7yJiksl3JqiPcplwM0+mRXjTZRfKi+aoDXao8atUsryPoH7mL41\nx7voc3Lw4dzkSyovCmwh2uOMi3j0E7h0vF0/aVLey2PoKQ/paAMMd14uJnFTAQ8db01PX8onYaLT\naIg2wHCXUTXiXlcoLtcnoeKtKdp5MOZp0hJtgOF2IYt4lCuUEPG2Fu4MJ/M0aAp2JkS4Y472qLKT\nuPqA+4631XCPc/xrYdBt0hjrUQy3Ht5WKJ8fflb5cnzFO6Zo58GQ66Y92BmG27/a/WLPNaV2AvcR\n79TCnTm4OcuIK2Al1Mcx3HqpDbhrqcab/LEa5CJ8x5vhrkZlwF1P34w3ZVKIrgucum1QGXCXGO90\nMM7VMdy2qAu4y+mb8X4tlv03A+1XivF+dv/o2xc/lTmOstQF3BXG2z4GO4yUwn082Kf9fyshjzbg\nZBOjHU4K4Z4W7Gkfpz3kqgLuan3C6dsWRjusWB/6XjbW0y5Tc8RVBdwFxvsorftvRju8GMPtI9qn\nfQ6NIVcTcOmniSW/GGw5Mb6gQohwn/Y5NYVcTcBd4PStE+MtI8ZwAzLxHvf5NYQ8moAz3idJr08Y\nbhkxv/akdLxHadiPqwg41ydxYbjlSLzieyia4g3IxxsAzkgfAPkhNX0z3nIk4h3zztsCFRM4uSUR\nb4ZbFuOdpigCzv33G6HjzXDLY7zTJR5w7r/dYbzTEvO+O0X76zPA/WK/GMUDTm6EjDfDLU8y3py+\n9WDAIxAq3gy3LA0TN+Oti2jAXaxPuP8Og/EOT0OwRzHe+nACNy7E9G053nkjeO+tW56PJB9t0Q4t\n1Xjvr8+U+jgG3DDf8dYYbl+Bq3q5VX4BWIi29NPC0nimA57y+iSleFsI3KRjPC3uFr6ukFKdvqsQ\nCzjvPlheCvGOKW7Z1zIa8pi+Plcufmon4i4fRl92fQIYn8BTFHO8Y4+a5a9vf30myBolC6PWkGt4\n/pNRDLghMcbbctTIH43TuI94V5m+AaGAc31SXEzxZrRtCjWFZzRFXGO8AU7gJviMN8NNmmmIuLa1\nySizTyebyj1QYoj3rVf3GO9IuJgai5IMqK/P7er7GHwC5/okP+vxZrTJFYkbN7XHGzA8gceO8Sat\nJKbwTKhpXPPaZBR34ApZjjfDnYbQN2iO8rkX9x1u17/8gk7gXJ9Mx3gTTecjtFam7lGcwBWxGm+G\nO02SUzhQbBLXEGcfqycGXAnGmyzSEHELfN1uECzgXJ+cjvFOm/VnMszixGcsHM/njb6cwIVZjLeG\naFjm8rnHNT1BFkN+ku977AQJOKfv8Rjv9Ph84QgtMR+NVsoxD3F3S07gQhjv9IR81Z/jn0vqZ5fq\nVB7qvvLeA87p+yTGOz3SL9kmPZ2nFPKQD3QqdT/wdvt1lDudjtODycvy86Aw3umRjvdx9966Nfwv\ntP31GdFHcvoW+murHR4eHhb9oLm5OczPz+Pbb7/F4uLiyQut1fD54Wfepm+rAWe806Mt3qeR+jnH\nMJG7jPbZtw9QJMmlViidTgfXrl0r86HJYrzTYyXewJtjDf0zt3yDp4Z/SZQK+GAwQK/Xw87ODm7f\nvu36mKLj+8UYfGC8q7EU71GSu3IrMdcQ7kypFUpmbW0Ny8vLaDabRy+0VsPM3q+VD24SK2sUi6+k\nw3iXZzXck0ifDxpiHira3lconU4Hc3NzuHbtGubn59Hv908EHABerd8Z/rn2x6s4c/X9op/KPMY7\nLTHGG9BzD5ZMyKD7DvcPT37Dj08O8VPtD9l7Cn184Qm81+vhypUrOHv2LNbW1vDxxx9jYWHh6IUG\nmMAB3VM4452WWOM9iZbzxXXQQ69IRs+dL2tf+Z3Am80mut0uAODcuXMn4p06i/tuqibFeAPyk3lG\n0066qKrnTqUd+KkXGmgCB3RN4aHizelbh1TDPQ3PpXzGnT9FJ3DzAQfkIx5y6ma8dWC88+P5ddSk\nc8f7CoWOYrzTw3gXo2XVooHrcyeKgM882A8+hXPXnSZN8c4e6RziRapdSTnmPs6dKFYomVARl4g3\np295GuLt4+kpNPwCiP1czHvuJLkDz/gOuNTUzXjLk4x36Gf0ZNDdKnLuJL0D97lKiSnepJv0UzAf\n//wS52AMq5YQv/SjmsAzLiMuvevm9C0vxBVROtpFSA8VFs7fsudM0iuU47KQS0e4LMZbns94W4r2\nJJJB13g+VzlnGPBI8Cli5THexUlP54DMOe7qXGHAI8B4y/MV71jDPY6GmI+qev6HWKUlfSNmDBjv\neKUUb+Do16sh5qcFePS6oeGuokUw4IpoOMnJz5U4tXgfpy3mo6xFexQDroTPk5rTd36ur8yph3sc\nzTG3hgFXgPHWIfZ4n3b3Wsl7aTHm1TDgwnjS6mAt3i4f6zDusiSizpgXx4AL8n2Sapu+Z9sHpp98\nPy+f8Q71fD/SUWfM82HAhaQYb61cTd8xhLvoMYSIusVnXQzFW8Dbc+vqdoBa8ETUQ3O8NUR7mtFj\n9B1zTuUneZ3As28yQ/5aqJNO2/StldZ4Wwj3OIx5eN4eifn54Wcn3p9yyBlvfVwE3OU5bTXck1h+\nnVgJah5KPy7go1KJecgTi/HOj/EOiyHPx0zAM7GGPPSJxHjnpynesYf7OIZ8MnMBH2U95lInDeOd\nn5Z4pxbucXzH3GLETT+ZlcUbPaVPEsY7LMbbHd/P15/C3Q9VBTxz/BuuLehaTgjGuxgNT1rEeJ90\ncHPW6zS+Pmiruc66pmqFUlTIsGs7ARjvYjSsTkTj/d2zfH/vk4t+j2OK1NcqpnfgrpS5omn/wY5i\nvItJKt55Q12EQNR9hlzzdZ0BjxzjXVzVgKuNt49Y5xEo6ClO46ZvxKTJGO/ipPfeXuItFe7TPr+n\noPNGzunOSB8A5cN4F6dhdeLUd8/k4z2O5+M6uDnrdQWl6mdcEAOu3K1X9xhvIWpXJ1p5/uXCiJ/E\nFYpSjHY10a1ONE7e43z3TPyeLCnhBK5INm0z3tVEuTqxxPM6xRdVP/OcOIErwGDHxWlkrMU7kx23\nh2nc9wN/LOEELoTTth/S0zfjfYynr8HXJG5tCucEHhiDTbnEEO8M9+LeMOABMNphRDN9xxTvjIeI\nc5XCFYo3XJHYo+KfzzHGO+Pha/OxSlFxHuTECbwiBlqHKO42GHO8M5zEnfIW8Kphk75CHsdQx018\ndZJCvDMGduJWnoJW7QTOYFJe2n7ZF5ZSvDOO72aY6hSuNuCkw2z7YPjn/fUZwSPxR3z6JiqJN2LS\nEbPtgyP/aWd++k6dw399uP5lauHGTE7giSsS6dn2QbRTOJFFDHhiLEzVeUVxv+8U99/kDAMeOdfB\n5hROmqV2Y2apgHe7XdTrdfT7fbRaLdfHRCXFNF1PIz19k0MG7lYYwq1X9/BlwY8pfCPmzs4OAKDZ\nbAIAnj59WvQiqKAfnvw29v3WbnDU4r/+8tfKl8H1SeY/pQ8gaYUD/vjxY7zzzjsAgEajga2tLecH\nRUf9+OT1i5xqCbb1XxZ//ct/Sx9CRLalD8CrUP9KK/u4l8IrlJcvX2Jubm749t7eXqlPTJONRvLv\nfgJm//eV4NHowvVJhByuUVLag5fagRd52XuazOo0m/KNmXzwDrlU5VHnhQNer9cxGAwAAC9evMD8\n/PzYv3f2bZth0urP/yF9BGPcl/oZf1X9Iv58tvpl0P/7NzcX852biwGAA4eXVfSGxZCXXzjgKysr\n2N7eRrPZxO7uLpaXl0/8HU7oRET+Fb4R89KlSwCAXq+Her2OhYUF5wdF5FK7fXTf3e120ev10Ol0\nhI7ItuPfz+xtfj/DK/VcKK1WC81m88R9wHnFcItXjOo2NjbQ7XaHb/NusNUc/34Cr8/PCxcu4Pz5\n80JHZVOn00Gn08Ha2trwfUUb6uzJrHjFcI9XjOpWV1fRaDSGb/NusNUc/34Cr8/T58+fY3FxUeio\n7On1elhaWkKr1UK/30ev1xs2s0hDnQWcVwz3eMVwj3eDdW8wGKDX6+Hu3bvSh2JGv98fNrLRaKDf\n7+PRo0eo1+vD9+VpqLOA84rhHq8YfvBGdreylere3h56vZ704ZjQarWGK+idnR1cuXIFL1++PHKv\nvjwNdfp84LxiuMUrhnt57wZL+XQ6neFOfH5+Hv1+X/iIbNnZ2cHly5eHdw4p2lBnAecVwy1eMfxY\nWVkZfi9Puxss5ddoNLC0tATg9cT43nvvCR+RLb1eD3fu3AFQrqHOAs4rhlu8YrixubmJ7e1tPHz4\nEADvBlvV8e9ns9nE1tYWut0uzp07x+9nARsbG7h9+zaA1+djmYbWDh3uPTqdznAhz6eZrS6bwHd3\nd3HrFl86jCgWW1tb+OijjzA3N4fBYIDNzU0sLi4WbqjTgBMRUTh8UWMiIqMYcCIioxhwIiKjGHAi\nIqMYcCIioxhwIiKj/g9yxy13VDYzBAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x522b5d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD7CAYAAABzGc+QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGLRJREFUeJzt3U9vm1d2x/EfnUDdyYzkbTsDOt4aMuxsZhxMKkp5Azai\nZG8KKOBdbFANECDTAHEU2zsDbUTvAxviK4gIjJE0XVSRAW+NUJ0CXYp2BHQjjKMu3IemJIp8/tx7\nzzn3/j5AAEuxqEfSw6+OD//VDg8PD0FEROackT4AIiIqhwEnIjKKASciMooBJyIy6u1pf6HT6QAA\nfvnlF3z99dcAgG63i3q9jn6/j1ar5fcIiYhorIkTeK/Xw9LSElqtFvr9Pnq9Hp4+fQoAaDabADB8\nm4iIwpoY8H6/j62tLQBAo9FAv9/Ho0ePUK/Xh+/L/j8REYU1cYUyuh7Z2dnBysoKfv75Z8zPzw/f\nv7e35+/oiIjoVLluxNzZ2cHly5dx6dIlAAAf+0NEJG/qjZjA6134nTt3AAD1eh2DwQAA8OLFiyPT\neKZWqzk8RCKidBQZkKcGfGNjA7dv3wbwOuQrKyvY3t5Gs9nE7u4ulpeXx37cr3+byX0Qk8y2D5xc\njmVf/AR88Qfpo3hjf93Nz/Y099665eyy1gftE+97tX4Hb7X/ufBlHdycdXFIk333zP/ncOpfAfyT\nm4v65KKbyxlj5sG+08trz607vbzMl7WvCv39iSuUra0trK2t4d1338Xc3BxqtdpwjdLr9VCv17Gw\nsFD+aInoKI8RU83Y1z1uMJAwcQJfWloarktGZTduZnclpHRYn75N+OSiwUlct4Obs86ncA34SEwD\nPvh76SMIw2W8J6n98WqQz1OJmYn0ivQBiNEwIKgPuO+Jz4JUAu7SpCvXmavvBzySCkxE/D3pA8jN\nx20Y0hFXH3DSw+cv01DTtzkmIl6R8a9RMuIMOIlzHW/pqci5Ty6aj5wWQe5JFBADTrlwlaUAI66W\n1NDAgJMoTt8FxRZxga/H1xQuce4x4DQVp29lYot4REJHnAEnMZy+K7C+Fxc+fp+78JDnYa7nQiFy\nzdq9TtTe+HU8gpofAKTsF47PB/esD9reHm4/ykTA99dn+JwoQnysT3zEO6npexKNQVcW7lBCRNxE\nwEkG4x0BqaAbibbvh9j7jriZHThvSKNxGO+Cst2zr8Aa3M37Xo/5PEc5gXv27D5w8VPpoyjOwvTN\neFeUN7TTpnZjwR7H6iReO/Tw8jq1Ws3Z84EfZ2kX/uz+mz9bi7jrgFuOt9obMMk5389YOC3iX9a+\nKvSCDmZWKBkrq5TReI97WzPGm1JlbZ1iLuAWnBbrZ/f1h9zKL0giXyxF3GTANUcmT6C1R9wlTt9k\nUYiIuziXvQXc2gM1XCgSZo3TOFcnJ3H/na4QP/uq57TXCdxnxDVP4RRHvIlCqHJue1+hpDiJW+Ty\nF2Is8eb0TaHOgbLnuMkdOKWD8SZpmiMeJOC+pnCuUdzQOn0z3qSF1nOCEzg5w3hTzDSeG8ECzl34\nSRoenanxXzGMN2ml7RwJOoH7iLjGAKXI1c9W6q6C2q6YpJemc4UrlIS5+uVn+V9Xmq6MZIeW8yZ4\nwC1f2ekky3tvLVdCsknD+RPFBG5xjSK9/9b2PWO8ySLp80gk4JzC42B17y19paO4SJ5PYhM4Iy7H\nxfRtMd68sZJ8kTqvolihAPpWAjQdX5SBYiJxjkUTcEsk999apm/Gm6g6Btwh6RsmrQgVb65MKHai\nAece3J6qP7OQ8SaKHV+VPiFV1ycW4s1wU0q4QgmMaxZ/GG+SFvoc5AROuWievhluShUn8ERI3s2S\n8SbygxM4TaXxxmaGm4gBD0pq/x3b9J16vGce7I99f+rflxSJr1A0TndVxHYjpbafT8qRmnmwf2q8\nR///tL9HfoU8R8UDfuvVPelDiFpM03fq8S7zMQx5eCG/5+IBJ700Td8px7sqhjxeDHggsa1WppF6\nXUs6HUPuX+jvLwPugZZYx/IMjZy+3WLI4yEacO6/9aqyPnE5fTPe/jDibkl8P3MFvN1uj3270+m4\nPyIiCoYRd0Pq+zg14BsbG+h2u0fe1+l0cOHCBZw/f97bgRU12z6QPoQjRtcoEiuVGNYnnL7JAslf\nglMDvrq6ikajceR9nU4Hz58/x+LiorcDIzla1icUBqfw8qS/d6V24IPBAL1eD3fv3i39ibn/pkk4\nfYclHSIqp9RD6VutFgDg+++/R6/XQ7PZdHpQsUjtofOupm/GmyzQ8Euv8ATe6XSGO/H5+Xn0+33n\nB0VyND14h8LSECQrtHyvCk/gjUYDV65cAQDs7e1heXl57N978sUPwz//7oN/wO8/+N3wba5P4sPp\nm1LgOty//fgDDv/9x9IfPzXgm5ub2N7exsOHD3Hjxg00m83hBH7u3DksLCyM/bg/ffF+6YOi8mK4\n9wmddHBz1vvUx1+ek/n4/p+5+j5w9U0rX33zdaGPnxrw69ev4/r160fed+3atUKfZJSP6VvbXQit\nkl6fMCCkkZZ1yTh8KD2REfwFF57meAOBA87dt1/W731C0/mKOH85HGXl+WL4ijykAgNCGliI9qhg\nEzinb7+qTt/S+2/Kz/UvO/7yfE063u259cIfwwmcyKDR6EqHxzrp71+ZcGeCBJzTd7y4/5Y3boKW\njpIFGr5HVeINBAi473jzLoT27/vNf8K7lzfqqX3vNUQ7UzXeAFcoBO6/U5FarEdpCjfgJt6A54Bz\ndeKf9embyCdt4QbcxRuIYALfX5/hGkWIi/13ylMh+aEx2oDbcGe8BZzTt38upm+uTygWKYU7E8VD\n6VNcI6T4NRONo/lRkz7jDUSwQkmRq3hz+iartAY74zvcmWgCnsouXMvkzft/U2jao50JFW8gooBT\nMZy+yQIr0QbChjsTVcBjn8K1TN9EPlmKNiAT7kxUAY+Zy3hz+iZtrEU7IxlvIMKAxziFa5u8uf8m\nF6xGG5APdya6gAPxRFxbuImqsBzsjJZwZ6IMuHU+w831CYUUQ7QBfeHORBvw0QhamMZDTNsu4s31\nCU0SS7AzWsOdiTbgo8bFUUPUuSKhGMQWbUB/uDNJBHyc4/EMEXSpYHNtQj4w3PKSDfhxeaZ0ixOz\ny3hzfUIAw60JAz6BxWBnOHWTawy3Pgx4hHzEm9N3uhhuvRjwiPiauhnvNDHcYYy+dsKXBT+WAY8E\nVybkCsMdhosXvWHAjfMdbk7f6WC4w3D5amUMuGGcusmV2OIde7gzDLhBocLN6TtusUUb0Bdu368N\nXDs8PDx0fqG1Gn79m9274GkVcuIOGW++Mn1YDLdfVaJ99u0DFEkyJ3DlJNYknLzjw2j75XvSPg0n\ncCU07LMlw80p3A+G2x8f0eYErpyGUI/DqTseMUYb0BFuqUn7NAx4RVqDXATjHQeG2x9t4c4kvUKJ\nIb5VaAs31yjlMNx+SESbKxQwzHloizcVE2u0Adlwa520T2N+Amesi9Eebk7hkzHcfmgJd/QTOINd\nnvZ40+kYbj+0hLss1RM4Y+2GtXBzCn+D4XZPc7SLTuAqA85wV2Mt2OOkHPGYow0w3JOYDjjDXU4M\nwT5N7CGPPdajJMJtIdqjot+Bpy7mWI8z82A/moinFOsMp22/1EzgnL7Ti3NRVkKeYqiP47RdjskV\nSozxZozDCh13Rno8hrsacwG3Em8GOQ6joWeE3WG43fAS8Ha7jfX1Nz+gbreLer2Ofr+PVqt18kJz\nBlxjvBlqonwYbfec34i5sbGBbrc7DPjOzg4AoNlsot/v4+nTp7h06VLhA9USbwabqBiGW4+pAV9d\nXcXm5ubw7cePH+PDDz8EADQaDWxtbRUOuFS8GWui8hhufQrfjfDly5eYm5sbvr23t1fo42N9WTCi\nWDHcJ822D079f/vr4Z6JtdT9wMve7skX4yWyIeX7b0+Kc9GP9x3zwgGv1+sYDAYAgBcvXmB+fj7X\nx4WIN8NNVE1q4a4a6yKX7yPmhQO+srKC7e1tNJtN7O7uYnl5eezfu/Mvr4Z//p9//Bi//6D0MU7F\ncBNVk9KaxHe0J33e4xH/4clv+PFJ+XtyT70b4ebmJlZXV/HNN9/gxo0bAIBOp4NGo5H7boQ+p2/G\nm6iclKZtqWiPM2kSV/lAHl8BZ7yJiksl3JqiPcplwM0+mRXjTZRfKi+aoDXao8atUsryPoH7mL41\nx7voc3Lw4dzkSyovCmwh2uOMi3j0E7h0vF0/aVLey2PoKQ/paAMMd14uJnFTAQ8db01PX8onYaLT\naIg2wHCXUTXiXlcoLtcnoeKtKdp5MOZp0hJtgOF2IYt4lCuUEPG2Fu4MJ/M0aAp2JkS4Y472qLKT\nuPqA+4631XCPc/xrYdBt0hjrUQy3Ht5WKJ8fflb5cnzFO6Zo58GQ66Y92BmG27/a/WLPNaV2AvcR\n79TCnTm4OcuIK2Al1Mcx3HqpDbhrqcab/LEa5CJ8x5vhrkZlwF1P34w3ZVKIrgucum1QGXCXGO90\nMM7VMdy2qAu4y+mb8X4tlv03A+1XivF+dv/o2xc/lTmOstQF3BXG2z4GO4yUwn082Kf9fyshjzbg\nZBOjHU4K4Z4W7Gkfpz3kqgLuan3C6dsWRjusWB/6XjbW0y5Tc8RVBdwFxvsorftvRju8GMPtI9qn\nfQ6NIVcTcOmniSW/GGw5Mb6gQohwn/Y5NYVcTcBd4PStE+MtI8ZwAzLxHvf5NYQ8moAz3idJr08Y\nbhkxv/akdLxHadiPqwg41ydxYbjlSLzieyia4g3IxxsAzkgfAPkhNX0z3nIk4h3zztsCFRM4uSUR\nb4ZbFuOdpigCzv33G6HjzXDLY7zTJR5w7r/dYbzTEvO+O0X76zPA/WK/GMUDTm6EjDfDLU8y3py+\n9WDAIxAq3gy3LA0TN+Oti2jAXaxPuP8Og/EOT0OwRzHe+nACNy7E9G053nkjeO+tW56PJB9t0Q4t\n1Xjvr8+U+jgG3DDf8dYYbl+Bq3q5VX4BWIi29NPC0nimA57y+iSleFsI3KRjPC3uFr6ukFKdvqsQ\nCzjvPlheCvGOKW7Z1zIa8pi+Plcufmon4i4fRl92fQIYn8BTFHO8Y4+a5a9vf30myBolC6PWkGt4\n/pNRDLghMcbbctTIH43TuI94V5m+AaGAc31SXEzxZrRtCjWFZzRFXGO8AU7gJviMN8NNmmmIuLa1\nySizTyebyj1QYoj3rVf3GO9IuJgai5IMqK/P7er7GHwC5/okP+vxZrTJFYkbN7XHGzA8gceO8Sat\nJKbwTKhpXPPaZBR34ApZjjfDnYbQN2iO8rkX9x1u17/8gk7gXJ9Mx3gTTecjtFam7lGcwBWxGm+G\nO02SUzhQbBLXEGcfqycGXAnGmyzSEHELfN1uECzgXJ+cjvFOm/VnMszixGcsHM/njb6cwIVZjLeG\naFjm8rnHNT1BFkN+ku977AQJOKfv8Rjv9Ph84QgtMR+NVsoxD3F3S07gQhjv9IR81Z/jn0vqZ5fq\nVB7qvvLeA87p+yTGOz3SL9kmPZ2nFPKQD3QqdT/wdvt1lDudjtODycvy86Aw3umRjvdx9966Nfwv\ntP31GdFHcvoW+murHR4eHhb9oLm5OczPz+Pbb7/F4uLiyQut1fD54Wfepm+rAWe806Mt3qeR+jnH\nMJG7jPbZtw9QJMmlViidTgfXrl0r86HJYrzTYyXewJtjDf0zt3yDp4Z/SZQK+GAwQK/Xw87ODm7f\nvu36mKLj+8UYfGC8q7EU71GSu3IrMdcQ7kypFUpmbW0Ny8vLaDabRy+0VsPM3q+VD24SK2sUi6+k\nw3iXZzXck0ifDxpiHira3lconU4Hc3NzuHbtGubn59Hv908EHABerd8Z/rn2x6s4c/X9op/KPMY7\nLTHGG9BzD5ZMyKD7DvcPT37Dj08O8VPtD9l7Cn184Qm81+vhypUrOHv2LNbW1vDxxx9jYWHh6IUG\nmMAB3VM4452WWOM9iZbzxXXQQ69IRs+dL2tf+Z3Am80mut0uAODcuXMn4p06i/tuqibFeAPyk3lG\n0066qKrnTqUd+KkXGmgCB3RN4aHizelbh1TDPQ3PpXzGnT9FJ3DzAQfkIx5y6ma8dWC88+P5ddSk\nc8f7CoWOYrzTw3gXo2XVooHrcyeKgM882A8+hXPXnSZN8c4e6RziRapdSTnmPs6dKFYomVARl4g3\np295GuLt4+kpNPwCiP1czHvuJLkDz/gOuNTUzXjLk4x36Gf0ZNDdKnLuJL0D97lKiSnepJv0UzAf\n//wS52AMq5YQv/SjmsAzLiMuvevm9C0vxBVROtpFSA8VFs7fsudM0iuU47KQS0e4LMZbns94W4r2\nJJJB13g+VzlnGPBI8Cli5THexUlP54DMOe7qXGHAI8B4y/MV71jDPY6GmI+qev6HWKUlfSNmDBjv\neKUUb+Do16sh5qcFePS6oeGuokUw4IpoOMnJz5U4tXgfpy3mo6xFexQDroTPk5rTd36ur8yph3sc\nzTG3hgFXgPHWIfZ4n3b3Wsl7aTHm1TDgwnjS6mAt3i4f6zDusiSizpgXx4AL8n2Sapu+Z9sHpp98\nPy+f8Q71fD/SUWfM82HAhaQYb61cTd8xhLvoMYSIusVnXQzFW8Dbc+vqdoBa8ETUQ3O8NUR7mtFj\n9B1zTuUneZ3As28yQ/5aqJNO2/StldZ4Wwj3OIx5eN4eifn54Wcn3p9yyBlvfVwE3OU5bTXck1h+\nnVgJah5KPy7go1KJecgTi/HOj/EOiyHPx0zAM7GGPPSJxHjnpynesYf7OIZ8MnMBH2U95lInDeOd\nn5Z4pxbucXzH3GLETT+ZlcUbPaVPEsY7LMbbHd/P15/C3Q9VBTxz/BuuLehaTgjGuxgNT1rEeJ90\ncHPW6zS+Pmiruc66pmqFUlTIsGs7ARjvYjSsTkTj/d2zfH/vk4t+j2OK1NcqpnfgrpS5omn/wY5i\nvItJKt55Q12EQNR9hlzzdZ0BjxzjXVzVgKuNt49Y5xEo6ClO46ZvxKTJGO/ipPfeXuItFe7TPr+n\noPNGzunOSB8A5cN4F6dhdeLUd8/k4z2O5+M6uDnrdQWl6mdcEAOu3K1X9xhvIWpXJ1p5/uXCiJ/E\nFYpSjHY10a1ONE7e43z3TPyeLCnhBK5INm0z3tVEuTqxxPM6xRdVP/OcOIErwGDHxWlkrMU7kx23\nh2nc9wN/LOEELoTTth/S0zfjfYynr8HXJG5tCucEHhiDTbnEEO8M9+LeMOABMNphRDN9xxTvjIeI\nc5XCFYo3XJHYo+KfzzHGO+Pha/OxSlFxHuTECbwiBlqHKO42GHO8M5zEnfIW8Kphk75CHsdQx018\ndZJCvDMGduJWnoJW7QTOYFJe2n7ZF5ZSvDOO72aY6hSuNuCkw2z7YPjn/fUZwSPxR3z6JiqJN2LS\nEbPtgyP/aWd++k6dw399uP5lauHGTE7giSsS6dn2QbRTOJFFDHhiLEzVeUVxv+8U99/kDAMeOdfB\n5hROmqV2Y2apgHe7XdTrdfT7fbRaLdfHRCXFNF1PIz19k0MG7lYYwq1X9/BlwY8pfCPmzs4OAKDZ\nbAIAnj59WvQiqKAfnvw29v3WbnDU4r/+8tfKl8H1SeY/pQ8gaYUD/vjxY7zzzjsAgEajga2tLecH\nRUf9+OT1i5xqCbb1XxZ//ct/Sx9CRLalD8CrUP9KK/u4l8IrlJcvX2Jubm749t7eXqlPTJONRvLv\nfgJm//eV4NHowvVJhByuUVLag5fagRd52XuazOo0m/KNmXzwDrlU5VHnhQNer9cxGAwAAC9evMD8\n/PzYv3f2bZth0urP/yF9BGPcl/oZf1X9Iv58tvpl0P/7NzcX852biwGAA4eXVfSGxZCXXzjgKysr\n2N7eRrPZxO7uLpaXl0/8HU7oRET+Fb4R89KlSwCAXq+Her2OhYUF5wdF5FK7fXTf3e120ev10Ol0\nhI7ItuPfz+xtfj/DK/VcKK1WC81m88R9wHnFcItXjOo2NjbQ7XaHb/NusNUc/34Cr8/PCxcu4Pz5\n80JHZVOn00Gn08Ha2trwfUUb6uzJrHjFcI9XjOpWV1fRaDSGb/NusNUc/34Cr8/T58+fY3FxUeio\n7On1elhaWkKr1UK/30ev1xs2s0hDnQWcVwz3eMVwj3eDdW8wGKDX6+Hu3bvSh2JGv98fNrLRaKDf\n7+PRo0eo1+vD9+VpqLOA84rhHq8YfvBGdreylere3h56vZ704ZjQarWGK+idnR1cuXIFL1++PHKv\nvjwNdfp84LxiuMUrhnt57wZL+XQ6neFOfH5+Hv1+X/iIbNnZ2cHly5eHdw4p2lBnAecVwy1eMfxY\nWVkZfi9Puxss5ddoNLC0tATg9cT43nvvCR+RLb1eD3fu3AFQrqHOAs4rhlu8YrixubmJ7e1tPHz4\nEADvBlvV8e9ns9nE1tYWut0uzp07x+9nARsbG7h9+zaA1+djmYbWDh3uPTqdznAhz6eZrS6bwHd3\nd3HrFl86jCgWW1tb+OijjzA3N4fBYIDNzU0sLi4WbqjTgBMRUTh8UWMiIqMYcCIioxhwIiKjGHAi\nIqMYcCIioxhwIiKj/g9yxy13VDYzBAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x522b5d0>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Create Samples\n",
      "\n",
      "Using the `fipy_response` function, we can now create a sample set of microstructures and responses. We create `Nsample` microstrucures over a 2D space of $N \\times N$. We choose a very small system to first demonstrate the linear regression in real space."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 10\n",
      "Nbin = 6\n",
      "Nsample = 5\n",
      "\n",
      "np.random.seed(1)\n",
      "microstructures = np.random.random((Nsample, N**2))\n",
      "responses = np.array([fipy_response(M, dt=dt, N=N) for M in microstructures])\n",
      "print microstructures.shape\n",
      "print responses.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(5, 100)\n",
        "(5, 100)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Bin the Microstructure"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The function `bin`, discretizes the original microstructure, $m_{a,s}$, into the binned microstructure, $m_{a,s}^h$, given by\n",
      "\n",
      "  $$ m_{a, s} = \\sum\\limits_{s=0}^{S-1} m_{a, s}^h \\chi^h $$ \n",
      "  \n",
      "The `bin` function takes $m_{a,s}$ and returns $m_{a,s}^h$. Don't look at it's code yet, this is the next exercise."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymks import bin\n",
      "?bin"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can examine three points to see a graphical representation of $M_{a, s}$ in terms of $m_{a, s}^h$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymks import draw_microstructure_discretization\n",
      "for s in range(3):\n",
      "    plt.figure()\n",
      "    draw_microstructure_discretization(microstructures, a=0, s=s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAACmCAYAAADpuSUWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE8xJREFUeJzt3T9s22afB/Cv3BuKAm9sy8DdFV1KyUvb4ZJY7lDciwNq\nK92uRV25fVF0aywtt8Wh8g5vkxvqP+lQ4JZIfKemSxJpOSBDLTLA4b2ptqS7od1MdUp7B0Smckun\nlzc8pSRKssQ/ovTQ+X4AISL1kP6RIX969JB8noRt2zaIiEgaC/MOgIiI3JiYiYgkw8RMRCQZJmYi\nIskwMRMRSYaJmYhIMn8z7wDoArt9e/R7IhorwfuYKTKJRO89DzMiz1hjpuh88cW8I5iNu3eBVApo\nt8X09evely0UgHv3htcHAMfHwPo6sLs7nTgpNlhjJgpDVYH33gPefVdMF4simW5teVvWMICTk968\nYhHY3+9NZzLAxx8zOb9gePGPKAxN6yVlAMhmgVJp8nKm6W7qAYBOB1hZcc/L54G9vfBxUqwwMRMF\n1WgMz1teBnR98rKGIZJ4v2fPRC36p5/c67OsUGFS/LCNmeanXBbJ7eCgl8xqNVFLBESbrWkCp6fu\nn/eyaLeBZNI9b2lJ/Pv8OXDp0ujlDAPY3hZtyP1SKbE/Xn+9N69WG07gdOExMdN8GIZoO63XgVwO\nODrqfZbLiUR95QqwsQEsLESXmAsFb+UODoDFRfc8y+pd8HM4ibrdPj8xW9bwuhyXL7vLPXo0umZO\nFxoTM0Vn3H3MqZRITvU68Oc/9+afnooaolNrNE1Rtp+miXmmCWxuAooy/Le9lAGG74jww6kd93MS\n9WBN2lGterswCIha9ZMn7ho0vRDYxkzRuXOn9xrkJErTdNcSB3+667p72mna2NgQt6Wp6vC6vZSZ\nhmRyuP3XmR5VW261RifzUYpF8erfN/TCYI2Z5kfXxe1g/ep1UcN1lErivt5OR9SwKxUgne59Pupn\nvpcyjjBNGVevDifadvv8NuFGQ3xpOPEcH4tE/tVXohbtfFlVq8C1a727PZpN0axDLwwmZpqfWk0k\nIEejIZof+mubrZZIUJomar/ttjvpAsMX2ryUcYRpygCAnR1384Su9y5eAiIRN5vi88EmjHJZfH7j\nRm+erov4Nzd7bdgPHjAxv2DYlEHz02oBH33Um67XgU8+cZfZ3haJb329N2/wgtsoXspMw/6+SK7V\nqqjZr64CH37Y+9wwRAIepGmiZt9q9X4RWJb4osrnxW1yyaRYX6s1m20hafDJP4pOFH1l3L0rmg+c\nx56TyeEk7KUMkcRYY6bofPFF7zUtH30kataAqGH2t+c6F97GlSGKAdaYKX6cW+EaDXHPs3M7WTLZ\ne0DjvDJEMcDETBdHpyOaLM67Z5koJtiUQReHaTIp04XAGjMRkWRYY6b4MU1xG9nz5/OOhCgSfMCE\nohPVmH+p1PCDKEQXCJsyKDpRjflnmqIfiYcPp7fOaeCj0zTIMESfLT6xxkzxo+vi1jjDEI91y9BX\nc6MhntDrT8xhxgKc5bLOAzmW5X4wp1+lIobA8ruvg8QyrrwTKyDi7R9yq1IRX9q5nHhyUtPEPe1e\nLghHub+DJGebKCqinixe05TL2bZhiPf5vG2b5nTX79fZmYij382bvRht27ZV1bYrFW/rm+WyOzu2\nXa32pnM529b13rSu2/bhoW1ns7ZdKHiLIWgsk8ofHrrLNxrueaWSbScS4rW87N6uacbpd1lVtW3L\n8ra+3zAxU3SiSsxra6Pfz8vBgW03m+55y8vuaV0Xyc2LWS17diaSWL9KZXR5VR3+8plmLF7Kj/q/\nzuV678tl2+50bLvVijZOv8uapth/PvCuDIoXy+p1nG+aotvQZnO+Mem6u9/kMGMBznLZ/tG5HYoy\ner5ffmPxUj6ZFJ1adTpiulod7vTq0iV/T3nOYn8riu9RaNjGTNGZ1EdGkDH/DKN3MiYSor1xcLTp\nWTLN0X0yBxkLcNbLnjfKyjQGf/Ubi5fypZLo90RRgFu3xBd0f09+gGhXdtZjmu426GnEGXTZVEpc\ng/D4ABQTM0Vn3C1yQcf86+/TWFHCX/gL01E+MHroq6BjAc562atXxb/OIARAr7Y8KSlN4jcWL+UV\nRXxp12piVJqbN93Hw+AQYoVCrx/vacUZdNl0WlRCmJhJamHG/CuXxYE+7zH/AJHUVlbc84KMBTiv\nZUslsT+dmuW4obH88BuLl/KqCvzhDyJWwxBf4KbZu21y8P84mxXLjEvMs9rfS0siVo/YxkzzEWbM\nP9OUY8w/YHQ/z37HApznstevi5pztSpe6fTw6C9B+I1lUvlGQzRZOcfKxoZoGnCawCxL/LLqfxp0\ncXFyMpzX/9UErDHT/MR9zD9AnJynp+55fscCnPey/ffYqup0vsj8xjKp/NnZ8C+TxcXesZJIiKaN\n/oRompO/ZGa1v/svWnvAGjPNj58x/5yfq6MuuAz2meGljOPePW+vUUkZEPGOuljmjAXoGDUWYKUy\nep2zXHZ1tXdXi2WJJoLPPx9e73lPboaNpf/zceU3NsTx0q8/2S0uDifuSkV8oU6KNcr97Xj2zNcv\nET6STdGZ1FfG9jZweNhrT9Y0UTPuH5y0UBA1kHRa/IwtFsV7p91wdVUk9P5k7qXMNF275r546XCe\nCDNNcRtVf8I7PBRx/vWvo9cZ1bLOWIPffdebdp74Oz0F/vhH935qNkWyKZVErbVYFLVU5wnHacYy\nqXyrJeJYWenVVPvbjzsd0V7ubMvbb7vv2hgXa1T723HeMXIOJmaKzosy5l+hcH5TxzgB+1EIvey0\nyRTLJEFjDbONliUSu48LzWzKoHiRccw/VQX29vwvF+Z+4WncazwtMsUySdBYw2yjpnm/lvEbXvyj\neFEUYG1N1GCch1Mczhh/48pEFVM67esBArRavi4GTW3ZaZMplkmCxhpmG50HkPrvPPKATRkUnai6\n/TzPvMf8m/QwA714Ah4TTMwUnVknZvaHTBcEmzIoOpP6ypg2JmW6IFhjJiKSDO/KICKSDBMzEZFk\nmJiJiGTja7yTAUdHR7aiKDYAvvjiiy++fnspimIfHR0Fzq2hLv6lUim0Wq2gi9MF90Xf+ztzi4Jo\nPhRFgemjD+Z+oRJzYp5D+pD0+g8sHin0IgqaXtnGTEQkGT5gMmV/+/t/nncIM/e/f/l31/RHt/5B\nvNn77+F5F1Slb1sB4L/+9G9zimR+Lv/rv7im/+P7/5tTJPPzT2//birrYY2ZiEgyTMxERJJhUwZF\n5uE//t28QyCKJSZmisyj3//9vEMgiiU2ZRARSYaJmYhIMkzMRESSYWImIpIML/5RZHJ/+aX7nhcC\nibxjYqbIbP/n/3TfMzETecemDCIiyTAxExFJhomZiEgyTMxERJLhxT+KDPvKIAqGiZkiwzsxiIJh\nUwYRkWSYmImIJMPETEQkGSZmIiLJ8OIfRYZ9ZRAFw8RMkWFfGUTBsCmDiEgyTMxERJJhYiYikgwT\nMxGRZHjxjyLDvjKIgmFipsjwTgyiYNiUQUQkGSZmIiLJMDETEUmGiZmISDK8+EeRYV8ZRMEwMVNk\n2FcGUTBsyiAikgwTMxGRZJiYiYgkw8RMRCQZXvyjyLCvDKJgmJgpMrwTgygYNmUQEUmGiZmISDJM\nzEREkmFiJiKSDC/+UWTYVwZRMEzMFBn2lUEUDJsyiIgkw8RMRCQZJmYiIskwMRMRSYYX/ygy7CuD\nKBgmZooM78QgCoZNGUREkmFiJiKSDBMzEZFkmJiJiCTDi38UGfaVQRQMEzNFhn1lEAXDpgwiIskw\nMRMRSSZh27YdeOFEYpqx0AXTf2DxSKEXUdD0yhozEZFkQl38UxQFrVZrWrHQBXN73gEQzZGiKMEX\ntkM4OjqyFUWxIX618sUXX3zxBdiKothHR0eBc2uoNmYiIpo+tjETEUmGiZmISDJMzEREkmFiJiKS\nDPvKIE/u3r2LVCqFdrsNALh+/brnZQuFAu7duxdVaEQzValUcHJygv39/Yllg543nhKz35WHOYll\nFWQfAMDx8THW19exu7sbeYxRUVUV7733Ht59910AQLFYRLVaxdbWlqdlT05Oog5xJvweA5ZlYX9/\nH+vr62i328hkMrhy5cosQo1MkPNgaWkJgNgfcT4PDMNAo9FArVZDOp2eWD7MeTPxPuabN2/ahmF0\np1VVtSuVytTKx4HfbVJV1TW9trZmHx4eRhZf1JaXl13Tuq7b2Wx24nKnp6e2qqr22tpaVKHNjN9j\n4OzszLXdh4eHdi6XizTGqPndB4PHfKPRiPV54FBV1c7n8xPLBT1vbNu2JyZmvysPE4ys/GyTZVlD\nB1+5XB5aR1zU6/Wh2Ov1up1IJCYuWy6XbV3XL0Ri9ntc7+zs2JqmueZZlhVJbLPidx+M+n+P+5eT\nbXtLzGHOG9u27bEX/xqNxtC85eVl6Lo+lfJx4Hebnj17BlVV8dNPP7nKW5YVVYiRarfbSCaTrnnO\nT9Pnz5+fu5xhGNje3g7ciYtMghzXmqZhc3PTNW9xcXHqsc1KkH2QTCaxvb2NTqcDAKhWq/jkk08i\ni1EmQc8bx9jE7HflYYORkd9tSqVSaDQaeP3117vzarUastlspHFGxbKsbnuiw9kfg/MHl4tzIurn\n9xgwTRMAcHp6imq1Ck3Tutcc4irIuV0qldBoNKAoSnf7P/zww2gDlUTQ88YxNjH7XXnYYGQUZJsu\nX77sWv7Ro0colUrRBRkh5+Tr52z34Inq8HyBIyb8HgNOYk4kEtja2upeICsWixFHGp0g54GiKMjn\n88hkMlBVFcfHx5HHKYsg502/sYnZ78rDBiOjsNu0vb2NJ0+euGrQcZJMJoeaYZzpS5cuDZVvtVoj\n91mc+T0GnHmZTKY7b2NjA4eHhxFFGL0g54Gqqshmszg6OkKtVkO5XMb29nakccrC73kzaOztcn5X\nHjYYGYXZpmKxiGKx6KpBx83Vq1eHTsp2u31u00yj0YBpmt02yePjY1iWha+++gpbW1vhukKcE7/H\ngLO/+j/r/9kfx3PB7z5oNBpIJBLdY39jYwOtViuW//9B+D1vBo2tMftdedhgZBR0m6rVKq5du9a9\nh7HZbEYWY9R2dnZQrVa707quI5/Pd6dN0+x+vrW1hd3d3e5rc3MTS0tLuHHjRmxPSr/HQCqVwtLS\nkquv8rhXUPzug7OzM6ysrLjmLS4uDl0QjaPzLmj3nwfA5PNmnJdu3759e1yBs7MzdDodvPnmmwCA\ncrmM999/H2+88UY3mCdPnnQ/n1Q+jvzuA13X8fPPPyObzeLXX3/F06dP8c0338T2oNzc3MTjx4/x\nyy+/4PHjx3jttdfw6aefdj+vVCq4f/8+PvvsM9dymqahUqngxx9/xEsvvYS33noLL7/88qzDnwq/\nx8DCwgKePn2KtbU1AMDXX3+NDz74AO+88858NmAK/OyDVCqFO3fuuI4Jy7Lwww8/xPY8aDab+Pbb\nb3H//n3U63UsLCzglVdewauvvgpg+DyYdN6M46k/ZudpH9M0sby8jM8//7z7mXPyfffdd57Kx5XX\nfWBZ1sg2t1wuhwcPHswyZJqyIOeBI5FI4MaNGzONNwp+9kGr1UKpVMLKykq3tn0RngKeBXaUT0Qk\nGfYuR0QkGSZmIiLJMDETEUmGiZmISDJMzDRVpmlidXU1tn2jEMmAiZmmKpVKIZVKxfZBCiIZMDHT\nVJmmOfW+MuL61KRhGPMOgWKKiZmmStd1JJNJGIYxld7UnL43BhmGgWw2i4WFBRQKhW4SbDabyOVy\nQ/OjYllWd8ggTdOGvkSYnCmQAB34E50rl8t1hx/K5/O2aZqB13V2djZ2pIiDgwM7k8kMzS+VSvbq\n6mrgv+uVl+GjVFWN/cglNHusMdNUmabZ7bjp5OQkVMdF5XIZhULh3M9rtdrIfhfOmz9tqqq64tvd\n3YWmaa4y+Xwee3t7kcdCFwsTM02NZVlIpVIARILOZDKh2od1XR/bZaphGFhfXx+a32w2Z9KjoZfh\noxRFGTksE9E4Y/tjJnJomobT01NYloV79+65PqtWq9jc3IRhGN0x3RKJBJaWlpBIJAL9vUkXEZ1k\nN5gYLcuCaZqR15j7h4+q1+tot9uwLAu7u7tDZVOp1AvVFzGFx8RMnji9giWTSeTzeVy5cgWAqLWm\nUiksLi66hpNSFAX7+/uB/55pmt3a9yi6rmNpaWloyK7vv/8e6XTa0+1645pJ+h0cHAzVhAeHjwJE\nz2vFYnFou9PpdHfsOyIvmJjJl52dHezt7eHhw4fd5OQk6WnqdDpDHa33q9Vq+Pjjj4dqqIVCwXNt\nebDm78d5w0dlMpmhxLy0tDTyzhKi87CNmXzZ2dlBpVJBq9VCs9nExsZGJH9n0uC99Xp9ZDvyrEYk\nnzR8FFEYrDGTL86TfaVSKVRTxSTJZBKnp6cjP2s0GrAsa2T7cqvV8lxjDtOU0T98lNNEcd7wUf0X\nRYm8YGIm31Kp1NDAnFH8jVqtNvIzXddHPvZ93vzzhGnKAIBbt25B1/Vu+/uDBw9GjoT97NmzWI97\nSbPHEUzIF03TkMlksLGxMbG5Iaxr167h6OioO93pdPDll19C0zSsrKxAVdXu0EZ3795FqVTC2dkZ\nbt26NbNhnLwMHzW4HUSTMDGTZ9VqFdlsFpcuXUIymcSjR48ia2MGRFPDqGaEOHEe2Q5bO6cXCy/+\nkSfOvcpOM8HOzs7QrWrTpqpq7J+a0zTNc1s2kYOJmSZqNpvde5Ud+/v7ME0TmqZF1lGPoihIp9No\ntVqRrD9qzkMy455eJBqFTRkkPU3TYjnsfVzjpvljYiYikgybMoiIJMPETEQkGSZmIiLJMDETEUmG\niZmISDJMzEREkvl/+MoukKEIPn4AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5536d10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAACmCAYAAAC7ttoYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFFZJREFUeJzt3U9sG1d+wPEflR4WAWpTFLZtkIs58mV3D7UsqocARQBT\ndNAeugvLkg0EudmiLr1ZSxkLdO0eYlnKIUAvJien9V4ci5cCabHm0MCit0gk20Ny4zCnJFvA1NC9\n7I09PAzJ4f8hOaQe+f0AhDT/OG84nPnx/Zn3QvV6vS4AAGhmadYJAABgFAQwAICWCGAAAC0RwAAA\nWiKAAQC0RAADAGjpL2adAABT8OhR9/8BjYV4DgxYAKFQ838uecwJcmDAIvjtb2edguk4PhYxDJFq\nVU3fv9973WRS5OBAJBrt/34iIqenIhsbIvv7k0srxkYODMB8SKVEPvpI5MYNNX1woILO1lb39VdX\nRSqVzvmZjMi9e2r7w8Pm/FhM5M4dgtgFQiMOAPPBNJvBS0QkkRBJp3uvn0iIFIsitq1e5bIKgvfu\niTiOyMqKd/1kUuTJk2DSjpEQwADor1jsnLe8LGJZ3dev1VSwunZN5MoV9bIskYcP1fJqVS3/7jvv\n+znOhBOOcVAHBkAVmxWLIk+fNm/6uZzKdYioG7qbS2ktVrsoqlWRSMQ7LxxWf9++Fbl0ybvs8mX1\nchWLqu7MXc8w1LwrV5rr5HIq14YLgwAGLLp8XtXtFAoi29sir141l21vq4C2tiYSj4ssLQUXwPb2\nhlvv6VNv8BFROSO34YbLDWjVamcAa5fJiDx75p137Zr3/V++7J7Tw8wQwIBF0O85MMNQAaFQEPni\ni+b8clnlONxciG2rdVuZpppn2yKbm71b9GUyqkHF2lrvNLYHED/c3FYrN6C158zaWZZq0NHPzo7I\n69feHBlmjjowYBE8ftx8tXODjm17cx3tRWaW5Z12ixTjcdVcPZXqvu9MRr3Oz8c/jl4ikc76KXd6\nUO4rne4fwA4O1Kv1s8GFQA4MgApOsZh3XqGgclWudFo9F1WrqRzbyYn3xt+reG13Vy0b9MTOOEWI\n16935sKq1eHqrLJZkd/8pveymzebrRtLpf65SEwVAQyAym3dvNmcbm/UIKKembpxQxUb3r+vAkR7\nzqVbg4lhjVOEKKICZTbbfO7LspqNUERUjrFU8j4X5ubSuhVBWpY6xs3NZh3bixcEsAuEIkQAKjjd\nvt2cLhRE7t71rrOzowLExkZzXnvDiVk6PFRBKptVOcWrV0Vu3Wouz+dVUWa71dXOejLHUQE9mVTN\n5yMR9X7dHnzGzNATB7AIgugL8fhY5Vzc7poikd4BbW9PBcDWB42BMVGECCyCIPpCvH1b1Ufdv69y\nLK31TY7TWSzHb2VMGAEMWARBDKESjYqsr6uiOfchaFfrg8AnJyJnZyoXGIlQh4SJoQgRwOTVaqo4\nsV9P78CYaMQBYPJsm+CFwJEDAwBoiRwYgNHZtmpe/vbtrFOCBUQjDmAR9OsLcRyG0fnAMzAlFCEC\niyCI58BEVA7s4EDkyy8n956TQJdPesrnVd+aQyIHBmB0lqWaxufzqjuqizBWWLGoesxoDWDHxyqn\n6D5o7T58PQwdt93b89c1l999OY461xsbaptYbPgfDIP25SOIUQcGYHSWpXrYiMfVTW3WXS05juou\nqrW/w1RKPa+2taVuluWy6m5qGLpue3Y23Lqj7MtxVP+Qh4dqG8cRefJkMvuKx9UPoVptuPerA5h/\nquBQvSZpfb37/7Py9Gm9Xip55y0ve6ctq15PJIZ7P922LZfr9VTK37nwu6/d3XrdNL3zHGdy+7Jt\ndQxDIAcGYDSO0xzg0rZVMVKpNNs0WZZ33K5uQ7wsL6v1BtFx23x+uCFkxtmXaXqH2RHpHN5mnH1F\no0OPfE0dGLAIBvWFmMk0u4Nybyi5XHM4kmq1OYClW8+Vzzd7rA+FVN+HrY1Fps22u48J1t7TvLvO\noKFfdNs2n1fFuaenvd973H3ZtvpbLqsRC6pV9UNmf3+y+zIMVRw94GF4AhiwCPo1nc/nRe7cUTek\n7W2RV6+ay7a3VUBbW1P1E0tLzQDWWs8UjY7fgGOcAS1F1M3VzRG63HG8Wrk30Wq1fyDRbVvHGS4n\nNM6+3AAWCjXP//Gxaok66Pz72dfqqvpBRQAD0JdhqBtfoSDyxRfN+eWyKo66ckVNdwsQmYy62di2\nKlbqdcPJZFSLtX4t1cYd0LJWE1lZ8c7rNlClexNtzw2002nb1oE8/fC7L3de6+jd8biaHhTA/Owr\nHG4Gyz6oAwMWnRt0bNtbf5TLeetTLMs7bdvqFY+rFmWpVPf3z2TU6/x88mlv1W0sskikOeqyy50e\n9PC1LttWKt2DwzD87svdT+uy1mLASe5rCOTAAKjg1PqrWkTlyFor69NpVVxUq6kc28mJyn25elW8\n7+6qZYMeoB63CDESUbnGVtevd68XG6ahgy7bFovqh4T7+Z+eqsDw2WcqV9avGM7vvgxDrd9aPzVs\nEPKzr9YGQn0QwACo3NbNm83pYrGzi6hKRY2obJoqx1WtegOYyOAGCv2MW4RoGOo42u3ueovYLKvZ\nOEWkefO/fVuvbUsltby96DCTUcsfPJj8vkREHj5U67gPIL94IXJ0NJl9ud68GSrYU4QILIJHj5qv\nbioV7w2nUGi2MHTt7Kibz8ZGc163YrtZWVvrXm9yeKjmZ7MqB3n1qsitW83lJyfq2Lq5qNvm8ypQ\ntTNN9b6VisqBuQ8ET3Jf+/sqh3R8rF4//ak3WI6zL1ex6C3O7oG+EIFFEERfiMfHqkjI/SUeifQO\naHt76qZ248Zk9t3L3l7vIsZ+fPbBdyG2ndd9OY5q1ThEjpwcGIDR3L6tcmoi6qbTWuTTXlkvMtlO\nhHtJpYbv1qhVt/Re9G3ndV+mOXR9KHVgAEYTjap+7fL55kPQLsNQ865cUUVKZ2cqFxiJBNtLfDSq\n6uWGeAi2oVIZqsHAhdp2XvflPow+RPGhCEWIwGIIajiVXmo1VZw4bBCZNLehCfTi87wRwIBFMO0A\nxnhcmAKKEIFFMKgvxEkjeGEKyIEBALREK0QAgJYIYAAALRHAAABaGqsRRy6Xk2QyKZVKZVLpAQDt\nRaNRSafTkvAzOjJ8G6sRh2EYBC9AA61tEB/PLBWLJRqNij3EmFYY3VgBLDTL4cMBDK31IueqnR4a\neQeLOjAAgJZ4kHnC/urv/2nWSZi6//2vf/dM3374tzNKyeycPPkfz/R//8u/zSglPfzrPzf+DSpt\n11r2ISLyx6//L5D9XGQf/t1fzjoJC4UcGABASwQwAICWKEIEFsCPH/7DrJMATBwBDFgAf/rwH2ed\nBGDiKEIEAGiJAAYA0BIBDACgJQIYAEBLNOIAFsBf//E/Gv/ToAPzggAGLIC/+eN/Nv4ngGFeUIQI\nANASAQwAoCUCGABASwQwAICWaMQBLAD6QsQ8IoABC4CWh5hHFCECALREAAMAaIkABgDQEgEMAKAl\nGnEAC4C+EDGPCGDAAqAvRMwjihABAFoigAEAtEQAAwBoiQAGANASjTiABUBfiJhHBDBgAdDyEPOI\nIkQAgJYIYAAALRHAAABaIoABALREIw5gAdAXIuYRAQxYAPSFiHlEESIAQEsEMACAlghgAAAtEcAA\nAFqiEQewAOgLEfOIAAYsAFoeYh5RhAgA0BIBDACgJQIYAEBLBDAAgJZoxAEsAPpCxDwigAELgL4Q\nMY8oQgQAaIkABgDQEgEMAKAlAhgAQEs04gAWAH0hYh4RwIAFQMtDzCOKEAEAWiKAAQC0RAADAGiJ\nAAYA0BKNOIAFQF+ImEcEMGAB0Bci5hFFiAAALRHAAABaCtXr9frIG4dCk0wLgIC0XuRctdMzxu0V\nQyAHBgDQ0liNOKLRqFQqlUmlBUBAHs06AQsoGo3OOglzb6wcWDqd5iQBGnjc8kLwotGopNPpWSdj\n7o1VBwYAwKxQBwYA0BIBDACgJQIYAEBLBDAAgJboCxFDOT4+FsMwpFqtiojI/fv3h952b29Pnj17\nFlTSgKk6OTmRs7MzOTw8HLjuONcNBhsqgPk9CfN40kb5DERETk9PZWNjQ/b39wNPY1BSqZR89NFH\ncuPGDREROTg4kGw2K1tbW0Nte3Z2FnQSp8Lvd8BxHDk8PJSNjQ2pVqsSi8VkbW1tGkkNzCjXQTgc\nFhH1eeh8HeTzeSkWi5LL5WR1dXXg+uNcNxhSfYBf//rX9Xw+35hOpVL1k5OTia2vA7/HlEqlPNPr\n6+v1o6OjwNIXtOXlZc+0ZVn1RCIxcLtyuVxPpVL19fX1oJI2NX6/A+fn557jPjo6qm9vbweaxqD5\n/Qzav/PFYlHr68CVSqXqyWRy4HqjXjcY3sAA5vckzONJ83NMjuN0XKSZTKbjPXRRKBQ60l4oFOqh\nUGjgtplMpm5Z1lwEML/f693d3bppmp55juMEkrZp8fsZdDvvugfxen24ADbOdYPh9W3EUSwWO+Yt\nLy+LZVkTWV8Hfo/pzZs3kkql5LvvvvOs7zhOUEkMVLValUgk4pnnFgm9ffu253b5fF52dnbmojPT\nUb7XpmnK5uamZ97ly5cnnrZpGeUziEQisrOzI7VaTUREstms3L17N7A0XiSjXjfwp28A83sS5vGk\n+T0mwzCkWCzKlStXGvNyuZwkEolA0xkUx3Ea9R0u9/Non9++nc437FZ+vwO2bYuISLlclmw2K6Zp\nNupEdTXKtZ1Op6VYLEo0Gm0c/61bt4JN6AUx6nUDf/oGML8nYR5P2ijHdO3aNc/2L1++1LZfNPcm\n1co97vYbmmveKqr9fgfcABYKhWRra6vR0OHg4CDglAZnlOsgGo1KMpmUWCwmqVRKTk9PA0/nRTHK\ndQP/+gYwvydhHk/auMe0s7Mjr1+/9uTIdBKJRDqKP93pS5cudaxfqVS6fmY68/sdcOfFYrHGvHg8\nLkdHRwGlMHijXAepVEoSiYS8evVKcrmcZDIZ2dnZCTSdF4Xf6waj6duM3u9JmMeTNs4xHRwcyMHB\ngSdHppvr16933Lyq1WrPItFisSi2bTfqTE5PT8VxHPnss89ka2tLy9EL/H4H3M+rdVlrcZuO14Lf\nz6BYLEooFGp89+PxuFQqFS3P/yj8XjcYTd8cmN+TMI8nbdRjymazcvPmzcYzIKVSKbA0Bm13d1ey\n2Wxj2rIsSSaTjWnbthvLt7a2ZH9/v/Ha3NyUcDgsDx480Pbm5fc7YBiGhMNhz1h5uv+Q8/sZnJ+f\ny8rKimfe5cuXOxq26KhXw6TW60Bk8HWD8b3z6NGjR/1WOD8/l1qtJj//+c9FRCSTycgvf/lL+dnP\nfiYi6qS9fv26sXzQ+jry+xlYliU//PCDJBIJ+fOf/yzff/+9/O53v9P24t3c3JSvvvpKfvzxR/nq\nq6/k/fffl48//rix/OTkRJ4/fy6ffPKJZzvTNOXk5ES+/fZbeeedd+QXv/iF/OQnP5l28ifC73dg\naWlJvv/+e1lfXxcRkc8//1x+9atfyQcffDCbA5gAP5+BYRjy+PFjz3fCcRz55ptvtL0OSqWS/P73\nv5fnz59LoVCQpaUleffdd+W9994Tkc7rYNB1g/ENNR6Y+/S9bduyvLws9+7dayxzb1J/+MMfhlpf\nV8N+Bo7jdK0T2N7elhcvXkwzyZiwUa4DVygUkgcPHkw1vUHw8xlUKhVJp9OysrLSyL3NQ688uDgY\n0BIAoCV6owcAaIkABgDQEgEMAKAlAhgAQEsEMEyUbdty9epVbfu+BKAPAhgmyjAMMQxD2wd2AeiD\nAIaJsm174n0h6tqLST6fn3USgLlGAMNEWZYlkUhE8vn8RHpfd/tWbJfP5yWRSMjS0pLs7e01gkWp\nVJLt7e2O+UFxHKcxVLxpmh3BliAGBGiWo2li/mxvbzeGnU8mk3Xbtkd+r/Pz874j3z59+rQei8U6\n5qfT6frVq1dH3u+wzs/PPaMOHx0ddYw4nEqltB+JGbioyIFhomzbbnRgfHZ2NlYHvplMRvb29nou\nz+VyXfvV6zV/0lKplCd9+/v7YpqmZ51kMilPnjwJPC3AIiKAYWIcxxHDMEREBbJYLDZW/ZVlWX2H\nosnn87KxsdExv1QqTWUEBNM0OwJl+yjU0Wi0MbQMgMnqOx4Y4DJNU8rlsjiOI8+ePfMsy2azsrm5\nKfl8Xu7evSsiqvPacDgsoVBopP0NagziBoX2AOI4jti2HXgOzK2XK5fLUigUpFqtiuM4sr+/37Gu\nYRgLNRYWMC0EMAzF7UU8EolIMpmUtbU1EVG5IMMw5PLly7K1tdVYPxqNyuHh4cj7s227kZvrxrIs\nCYfDkk6nPfO//vprWV1dHaoZf7/iyVZPnz7tyFm5ASwUCjWO+/j4WA4ODjqOe3V1VYrFIgEMmDAC\nGHzZ3d2VJ0+eyJdfftm4ibvBbJJqtVrHgIitcrmc3LlzpyPHs7e3N3Tuqz0n6Yc7ZE4sFmvMi8fj\nEovFOgJYOBzu2pISwHioA4Mvu7u7cnJyIpVKRUqlksTj8UD2U61W+y4vFApd67lyudxU6r/c4s3W\nnJ47j15IgOkgBwZf3J420un0WEWEg0QiESmXy12XFYtFcRyna/1XpVIZOgc2ThGiYRgSDoc9dVuO\n44iIdBRftjZuATA5BDD4ZhhG42Yd5D5yuVzXZZZlde2uqtf8XsYpQhQRefjwoViW1agffPHihRwd\nHXWs9+bNm6nkCoFFw4jM8MU0TYnFYhKPxwcW843r5s2b8urVq8Z0rVaTTz/9VEzTlJWVFUmlUo0h\n7Y+PjyWdTsv5+bk8fPhQHjx4EGjaXMfHx43/Q6FQ1/22HweAySCAYWjZbFYSiYRcunRJIpGIvHz5\nMrA6MBFVxNet+E4nbldT4+b2AHSiEQeG4j7r5RbP7e7udjRhn7RUKqV9LxamaQ5d1wbAHwIYBiqV\nSo1nvVyHh4di27aYphlYh7XRaFRWV1elUqkE8v5Bcx/G7tebCIDRUYSIC880zUZDCZ3omm5AFwQw\nAICWKEIEAGiJAAYA0BIBDACgJQIYAEBLBDAAgJYIYAAALf0/aXS/5Gc/sxIAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5c05210>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAACmCAYAAADpuSUWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEtJJREFUeJzt3T9s22afB/Cv3BuKAmfLMnBA0SWks7QdLqnlDgUOB8RW\nsl2LOHICFL0ptvwOt8Whghe4Jjc0ttOhwA1vJL5T0yWJtByQIRYV4HBbY0l3Q7uZ6pQWB0Smckun\nlzc8L/WX+sN/1kPn+wGImORD6iFD/vTo4cPnSdi2bYOIiKQxN+sMEBFRPwZmIiLJMDATEUmGgZmI\nSDIMzEREkmFgJiKSzN8E3sPdu+5/ExGRL4nA7ZgTie7fbBJNRBRY8BLz11+HkI0YePAAUFWg1RLz\nW1vB0gddDwClEnB0BOztuedh0noiklLwEvPbQNOAK1eAS5fEfD4PrK4CGxv+0gddX60C9TpQqQDL\ny8Cf/tT/+ZPWE5HUGJinkUp1S66ACHz7+8Dhob/0Qdc78nnAsoCHD93zMWk9EUmJrTImqdeHly0u\nAobhL33Q9UR05gWvY56kWBTBZn+/G1wqFSCXE3+3WoBpAsfHctaFtlqiBNsrmRT/vnkDzM97Sx90\n/eDnEdGZE21grlaB69eBWg3IZvt/imezIlBfvAisrQFzc9EF5p2d6dLt7wMLC/3LLKu/WgHoBs5W\nazhQTkofdD0DM9GZF207ZlUVga5WA/785+7y42MgkwHOnRPzpinS9tJ1scw0gfV1QFGGP1vXxb+1\nmnhg5pYGCFbH6pRWezmBc7BkO036oOuJ6MwLHpjv3ev+PRiYnUBpmsCFC93llQrwhz905w1DBGqH\nU7WxtSVK05ubwJMn/ftuNIB0WpS4VVWUwI+OAh/OkFRKlGJ7OfNupddJ6YOuJ6IzL/o6ZsMQAbRX\nrSZKwY5CQbTbbbdFCbtUEs28HG4PxExTBPiHD4GVFTE/SpCqjE8+GS7Ftlr9XyRe0gddT0RnXvSB\nuVIBLl/uztfrooTbW/prNkWbXV0XpeRWqz8wA8MPvjY2usF9sMQ9KGhzse1toFzutiM2jO7DS0B8\nKTQa3fWT0gdd75jU0pEtIYliKfrA3Gz2V1vUasCNG/1pNjdFIFpd7S4bfADmxindPnnSX4cdtr09\nUaIvl0UQPn8euHq1u75aFaV8J5BOSh90faMhgnW5DJyciC+x9XVRrTPNeiKSmpx9ZTx4IH7OO68h\nD75wMZg2l2P9KxGdGXL2lXHtmqjv3doSD756qyksq1sHWyqJn/3z86KE2FtvTUQUU/K+ku00l6vX\nRYsLp2ldKiWWtVoiEDul6dVV4PnzmWaZiCgM8gbmUdptEYhHtVkmIoq5+PWVYZoMykR0psWvxExE\ndMbJW2J2mom9eTPrnBARnSp5x/xT1eEXUYiI3gJytmMGRIk5nx/uI2PWGg2+qPG2q1ZFHy5EEYn+\nzT+/DEM0hatWxWvdMvTVXK+LNxkHA7PfsfVkGkewVBJfhtms6Jhf10V7cudBay4nvij9PHj1epxx\n2JbBmaJkByXKyWIKUzZr29Wq+DuXs23TDHf/Xp2ciHz0MgzbPjiw7UzGtnd2vO3v9u3u8dm2bWua\nbZdK/tMHXV8o2HYiIabFRdsul/s/X1W763snXQ/3OOOyrabZtmVNtz8ij+QNzCsr7n/Pyv6+bTca\n7us0bThoT7K42D9vGCLA+00fdH2xaNvttm03m+6fn8uJ4282xWSatp3Pj87vtJ8b121NU/y/E0VA\nzlYZltXtON80RbehjcZs82QY/X1KByHrOILz8903LHu122IgggsXxPpz58S2d+6453fafMd5W0Vx\nT0sUguj7yvAz5l+12u2BLpEQfWP0PmQ8babpPrKIX7KOI6jr3XSmCezuir8XFvr7qXbrujWM44zb\ntqoqnjnwhScKWbjN5Qb5HfPP6T4TEBd90Ad/QTrKB9yHvgpCxnEEB4fv2tnp9o89qFicro/rIOMX\nxmHb5WXxJcXATCGLtlVGkDH/ikVx4c96zD9A/JRfWgq2j14yjiM4eO4yGXFOBwOzYQwPYjBKkPEL\n47BtMjl+5Bwin6KtYx435l9vV55uY/6ZpihJb22JADHIGfNva0uUvrPZaI4BmK7Tfi9kG0fQssQv\nlt63LBcW3INOoTB9YA4yfmEctyUKCcf86zWqKiOVEqX8sMg2jmAiAdy+3R94TNM9AJfLwB//OPrY\nvOQ77tv2PqQmClH0rTK8jPnnvOXn9gBmsM+MjQ0RSIHpxvybZnILyoDI72ApqteoNx5NU3zJuHHG\n9XO4jSPYu35S+iDrFxaGq2pKpe75dTjnwO3n/qhjneY4g5yjWWzreP16+l8PRB4EfyV7Ul8Zm5vA\nwUG3PlnXRcn41q1ump0dEViXl0WVRz4v/nbqN8+fFwF91E/JzU1Rhx3lT83Ll/sfXgLdsfUKBTG2\nXj7fP7bewYFY9pe/uO/TebvMNEWTrJs3u+t0XQSO3s7/x6UPur7dFvX6yaT4dfDpp/3jDAIiMK+u\nil88g+d63LGO+9wg52hW2zrcrgmiMARuCR3FCyYHB+KFB8dgg//BtO12eJ89Si7n700vwwg/L7Ly\ne6xBztGstnV7E5QoJHK+YHLtmiiVAe5j/jkGx/yLkqYB9+97325cFchZ4/dYg5yjWW2r69M/uyDy\nSM5OjBRFPNCrVrsvpziccQBbLRGUe8f8i3IwVkUR1SteXihoNt+eh0N+jzXIOZrVts4LR2G9CUo0\nQN5uP0eZ9Zh/o166oLcHrwGKWPwCM/tDJqIzLvq+MsLGoExEZxwHYyUikoycrTKIiN5iDMxERJJh\nYCYikk2Qt1MODw9tRVFsAJw4ceLE6a+Toij24eGh79ga6OGfqqr452azM3/P746IiM4YRVFg+uyv\nO1BgTiQS6N14hoM/ERFJx294ZR0zEZFk5OwrI8b+7h/+adZZOHX/+1//0Td/7c7fzygns1O6/z99\n8//9r/8+o5zMzoV/+5e++f/88f9mlJPZ+cdP/zaU/bDETEQkGQZmIiLJBK7KuBtCJoiIqCtwYGYT\nOSKicLEqg4hIMgzMRESSYWAmIpIMAzMRkWQCP/z7uudvPggkIgou1OZyDMxERMGxKoOISDIMzERE\nkmFgJiKSDAMzEZFk2FcGEZFk2FcGEZFkWJVBRCQZBmYiIskwMBMRSYaBmYhIMuwrg4hIMuwrg4hI\nMqzKICKSDAMzEZFkGJiJiCTDwExEJBn2lUFEJBn2lUFEJBlWZRARSYaBmYhIMgzMRESSYWAmIpIM\n+8ogIpIM+8ogIpIMqzKIiCTDwExEJBkGZiIiyTAwExFJhn1lEBFJhn1lEBFJhlUZRESSYWAmIpIM\nAzMRkWQYmImIJMO+MoiIJMO+MoiIJMOqDCIiyTAwExFJhoGZiEgyDMxERJJhXxlERJJhXxlERJJh\nVQYRkWQYmImIJMPATEQkGQZmIiLJsK8MIiLJsK8MIiLJsCqDiEgyDMxERJJJ2LZt+944kUDvxokQ\nMkREdFb4Da8sMRMRSSbQwz9FUXC32QwrL0REZ4aiKP43tgM4PDy0FUWxAXDixIkTp79OiqLYh4eH\nvmNroDpmIiIKH+uYiYgkw8BMRCQZBmYiIskwMBMRSSZwXxn0dnjw4AFUVUWr1QIAbG1tTb3tzs4O\nHj58GFXWiE5VqVTC0dER9vb2Jqb1e99MFZi97jzITSwrP+cAAF6+fInV1VXs7u5GnseoaJqGK1eu\n4NKlSwCAfD6PcrmMjY2NqbY9OjqKOounwus1YFkW9vb2sLq6ilarhXQ6jYsXL55GViPj5z5IJpMA\nxPmI831QrVZRr9dRqVSwvLw8MX2Q+2ZiO+bbt2/b1Wq1M69pml0qlUJLHwdej0nTtL75lZUV++Dg\nILL8RW1xcbFv3jAMO5PJTNzu+PjY1jTNXllZiSprp8brNXByctJ33AcHB3Y2m400j1Hzeg4Gr/l6\nvR7r+8ChaZqdy+UmpvN739i2bU8MzF53HiQzsvJyTJZlDV18xWJxaB9xUavVhvJeq9XsRCIxcdti\nsWgbhnEmArPX63p7e9vWdb1vmWVZkeTttHg9B27/73H/crLt6QJzkPvGtm177MO/er0+tGxxcRGG\nYYSSPg68HtPr16+haRp++eWXvvSWZUWVxUi1Wi2kUqm+Zc5P0zdv3ozcrlqtYnNz03cnLjLxc13r\nuo719fW+ZQsLC6Hn7bT4OQepVAqbm5tot9sAgHK5jBs3bkSWR5n4vW8cYwOz150HzYyMvB6Tqqqo\n1+s4d+5cZ1mlUkEmk4k0n1GxLKtTn+hwzsfg8sHt4hyIenm9BkzTBAAcHx+jXC5D1/XOM4e48nNv\nFwoF1Ot1KIrSOf6rV69Gm1FJ+L1vHGMDs9edB82MjPwc04ULF/q2f/r0KQqFQnSZjJBz8/Vyjnvw\nRnVM/YAjJrxeA05gTiQS2NjY6Dwgy+fzEec0On7uA0VRkMvlkE6noWkaXr58GXk+ZeHnvuk1NjB7\n3XnQzMgo6DFtbm7ixYsXfSXoOEmlUkPVMM78/Pz8UPpms+l6zuLM6zXgLEun051la2trODg4iCiH\n0fNzH2iahkwmg8PDQ1QqFRSLRWxubkaaT1l4vW8GjW0u53XnQTMjoyDHlM/nkc/n+0rQcfPJJ58M\n3ZStVmtk1Uy9Xodpmp06yZcvX8KyLHz77bfY2NgI1hXijHi9Bpzz1buu92d/HO8Fr+egXq8jkUh0\nrv21tTU0m81Y/v/74fW+GTS2xOx150EzIyO/x1Qul3H58uVOG8ZGoxFZHqO2vb2NcrncmTcMA7lc\nrjNvmmZn/cbGBnZ3dzvT+vo6kskkbt26Fdub0us1oKoqkskkmj19lce9gOL1HJycnGBpaalv2cLC\nwtAD0Tga9UC79z4AJt8347xz9+7du+MSnJycoN1u46OPPgIAFItFfP755/jwww87mXnx4kVn/aT0\nceT1HBiGgV9//RWZTAa///47Xr16he+//z62F+X6+jqePXuG3377Dc+ePcMHH3yAL7/8srO+VCrh\n0aNH+Oqrr/q203UdpVIJP//8M9555x18/PHHePfdd087+6Hweg3Mzc3h1atXWFlZAQB89913+OKL\nL/DZZ5/N5gBC4OUcqKqKe/fu9V0TlmXhp59+iu190Gg08MMPP+DRo0eo1WqYm5vDe++9h/fffx/A\n8H0w6b4ZZ6r+mJ23fUzTxOLiIm7evNlZ59x8z58/nyp9XE17DizLcq1zy2azePz48WlmmULm5z5w\nJBIJ3Lp161TzGwUv56DZbKJQKGBpaalT2j4LbwGfBnaUT0QkGfYuR0QkGQZmIiLJMDATEUmGgZmI\nSDIMzBQq0zRx/vz52PaNQiQDBmYKlaqqUFU1ti9SEMmAgZlCZZpm6H1lxPWtyWq1OussUEwxMFOo\nDMNAKpVCtVoNpTc1p++NQdVqFZlMBnNzc9jZ2ekEwUajgWw2O7Q8KpZldYYM0nV96EuEwZl88dGB\nP9FI2Wy2M/xQLpezTdP0va+Tk5OxI0Xs7+/b6XR6aHmhULDPnz/v+3OnNc3wUZqmxX7kEjp9LDFT\nqEzT7HTcdHR0FKjjomKxiJ2dnZHrK5WKa78Lo5aHTdO0vvzt7u5C1/W+NLlcDvfv3488L3S2MDBT\naCzLgqqqAESATqfTgeqHDcMY22VqtVrF6urq0PJGo3EqPRpOM3yUoiiuwzIRjTO2P2Yih67rOD4+\nhmVZePjwYd+6crmM9fV1VKvVzphuiUQCyWQSiUTC1+dNeojoBLvBwGhZFkzTjLzE3Dt8VK1WQ6vV\ngmVZ2N3dHUqrqupb1RcxBcfATFNxegVLpVLI5XK4ePEiAFFqVVUVCwsLfcNJKYqCvb09359nmman\n9O3GMAwkk8mhIbt+/PFHLC8vT9Vcb1w1Sa/9/f2hkvDg8FGA6Hktn88PHffy8nJn7DuiaTAwkyfb\n29u4f/8+njx50glOTpAOU7vdHupovVelUsH169eHSqg7OztTl5YHS/5ejBo+Kp1ODwXmZDLp2rKE\naBTWMZMn29vbKJVKaDabaDQaWFtbi+RzJg3eW6vVXOuRT2tE8knDRxEFwRIzeeK82VcoFAJVVUyS\nSqVwfHzsuq5er8OyLNf65WazOXWJOUhVRu/wUU4Vxajho3ofihJNg4GZPFNVdWhgzig+o1KpuK4z\nDMP1te9Ry0cJUpUBAHfu3IFhGJ3698ePH7uOhP369etYj3tJp48jmJAnuq4jnU5jbW1tYnVDUJcv\nX8bh4WFnvt1u45tvvoGu61haWoKmaZ2hjR48eIBCoYCTkxPcuXPn1IZxmmb4qMHjIJqEgZmmVi6X\nkclkMD8/j1QqhadPn0ZWxwyIqga3aoQ4cV7ZDlo6p7cLH/7RVJy2yk41wfb29lBTtbBpmhb7t+Z0\nXZ+6LpvIwcBMEzUajU5bZcfe3h5M04Su65F11KMoCpaXl9FsNiPZf9Scl2TGvb1I5IZVGSQ9Xddj\nOex9XPNNs8fATEQkGVZlEBFJhoGZiEgyDMxERJJhYCYikgwDMxGRZBiYiYgk8/+BWbWceiUwNQAA\nAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5fafb10>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Run `bin` for each microstructure and rebuild the array (maybe this operation could be vectorized)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "binnedMicrostructures = np.array([bin(M, Nbin) for M in microstructures])\n",
      "print binnedMicrostructures.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(5, 100, 6)\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The new `binnedMicrostructures`, $m_{a,s}^h$, has a shape of `(Nsample, N*N, Nbin)`. To double check that the binning worked we can evaluate $\\sum\\limits_{h=0}^{H-1} m_{a, s}^h \\chi^h$ and check against the original $m_{a,s}$. The summation is over the last axis (the binning axis)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.linspace(0, 1, Nbin)\n",
      "reconstructedMicrostructure = np.sum(binnedMicrostructures * X[np.newaxis, np.newaxis, :], axis=-1)\n",
      "print np.allclose(reconstructedMicrostructure, microstructures)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print X"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.   0.2  0.4  0.6  0.8  1. ]\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Exercise 02-0\n",
      "\n",
      "Try and create a binning function. Just use loops to do this and then we will discuss how to vectorize.\n",
      "\n",
      "Test it with `np.array((0.2, 0.5, 0.7))` and `Nbin = 4`.\n",
      "\n",
      "If this is too difficult, look at the code, `??bin`, and then reproduce this with loops."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "??bin"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Real space\n",
      "\n",
      "In order to understand how to compute the influence coefficients it is\n",
      "useful to see a demonstration in real space. Although we have a tensor\n",
      "representations of $m_{a, s}^h$, we need to create an intermediate\n",
      "matrix to calculate the convolution, $\\sum\\limits_{h=0}^{H - 1} \\sum\\limits_{t=0}^{S-1} \\alpha_t^h m_{a,s + t}^h$. This matrix representation of $m_{a, s}^h$ is\n",
      "given by the `microstructureMatrix` and has shape `(N*N*Nsample, N*N*Nbin)`. The `microstructureMatrix` is essentially a\n",
      "[circulant matrix](http://en.wikipedia.org/wiki/Circulant_matrix) that\n",
      "isn't square.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rollMatrix(m, N, Nbin):\n",
      "    matrix = np.zeros((N**2, N**2 * Nbin))\n",
      "    for i in range(N**2):\n",
      "        matrix[i] = np.roll(m, -i, axis=0).swapaxes(0,1).flatten()\n",
      "    return matrix\n",
      "\n",
      "microstructureMatrix = np.concatenate([rollMatrix(m, N, Nbin) for m in binnedMicrostructures])\n",
      "\n",
      "print microstructureMatrix.shape\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(500, 600)\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To calculate the influence coefficients, $\\alpha_s^h$, we use `numpy`'s `lstsq` function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "responses = responses.flatten()\n",
      "coefficients = np.linalg.lstsq(microstructureMatrix, responses)[0]\n",
      "print coefficients.shape\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(600,)\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If `Nbin > Nsample` then we can check that the influence coeffiencts exacly reproduce the `responses`. The result below should be `True` for an over-determined system."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.allclose(np.dot(microstructureMatrix, coefficients), responses)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Fourier Space"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Having calculated the coefficients for a small system in real space,\n",
      "we will now calculate a much larger system in Fourier space and\n",
      "then confirm that a sensible response is reproduced. Calculating the regression\n",
      "in Frequency space drastically reduces the\n",
      "computational cost. The size of the regression is reduced\n",
      "from $ \\left(N^2 N_{\\text{sample}} \\times N^2 H \\right)$ to\n",
      "$\\left(N_{\\text{sample}} \\times H \\right)$ for each point in\n",
      "reciprocal space $N^2$. The convolution,\n",
      "\n",
      "  $$ \\sum\\limits_{h=0}^{H-1} \\sum\\limits_{t=0}^{S-1} \\alpha_t^h m_{a,s + t}^h $$\n",
      "  \n",
      "can be deconvolved in Fourier space with,\n",
      "\n",
      "$$ \\begin{split}\n",
      "\\mathcal{F}_k \\left( \\sum\\limits_{h=0}^{H-1} \\sum\\limits_{t=0}^{S-1} \\alpha_t^h m_{a,s + t}^h \\right)\n",
      "=&  \\sum\\limits_{s=0}^{S-1} \\left[ \\sum\\limits_{h=0}^{H-1} \\sum\\limits_{t=0}^{S-1} \\alpha_t^h m_{a,s + t}^h \\right] e^{-2 \\pi i k s / S} \n",
      " \\;\\;\\; & \\text{(definition)}\n",
      " \\\\\n",
      "=&  \\sum\\limits_{h=0}^{H-1} \\left[ \\sum\\limits_{t=0}^{S-1} \\alpha_t^h \\left[ \\sum\\limits_{s=0}^{S-1} m_{a,s + t}^h e^{-2 \\pi i k s / S} \\right] \\right]\n",
      " \\;\\;\\; & \\text{(rearrange)}\n",
      " \\\\\n",
      "=&  \\sum\\limits_{h=0}^{H-1} \\left[ \\sum\\limits_{t=0}^{S-1} \\alpha_t^h \\left[ \\sum\\limits_{v=t}^{S-1+t} m_{a,v}^h e^{-2 \\pi i k \\left(v - t\\right) / S} \\right] \\right]\n",
      " \\;\\;\\; & \\text{($v=s+t$)} \n",
      " \\\\\n",
      " =&  \\sum\\limits_{h=0}^{H-1} \\left[ \\sum\\limits_{t=0}^{S-1} \\alpha_t^h e^{2 \\pi i k t / S} \\left[ \\sum\\limits_{v=t}^{S-1+t} m_{a,v}^h e^{-2 \\pi i k v / S} \\right] \\right]\n",
      " \\;\\;\\; & \\text{(rearrange)} \n",
      "  \\\\\n",
      " =&  \\sum\\limits_{h=0}^{H-1} \\left[ \\sum\\limits_{t=0}^{S-1} \\alpha_t^h e^{2 \\pi i k t / S} \\left[ \\sum\\limits_{v=0}^{S-1} m_{a,v}^h e^{-2 \\pi i k v / S} \\right] \\right]\n",
      " \\;\\;\\; & \\text{(periodicity of $m_v^h$)} \n",
      " \\\\\n",
      " =&  \\sum\\limits_{h=0}^{H-1} \\left[ \\left[ \\mathcal{F}_k \\left( \\alpha_t^h \\right) \\right]^* \\mathcal{F}_k \\left( m_{a,t}^h \\right) \\right]\n",
      " \\;\\;\\; & \\text{} \n",
      "\\end{split}$$\n",
      "  \n",
      " The above is a proof of the [circular convolution theorem](http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Circular_convolution_theorem_and_cross-correlation_theorem). If we write $P_{a,k} =\n",
      " \\mathcal{F}_k \\left( p_{a, s} \\right)$, $M_{a, k}^h = \\mathcal{F}_k\n",
      " \\left( m_{a,s}^h \\right)$ and $\\beta_k^h = \\mathcal{F}_k \\left(\n",
      " \\alpha_s^h \\right)$, then we just need to solve\n",
      " \n",
      " \n",
      "   $$ P_{a,k} = \\sum\\limits_{h=0}^{h-1} \\beta_k^h M_{a, k}^h $$\n",
      "   \n",
      "with a linear regression at each discretization location in $k$ to calculate the $\\beta_k^h$.\n",
      "\n",
      "In the example below we will create a microstructure with 11 bins and\n",
      "160 samples. This is enough to give reasonable influence coefficients."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 20\n",
      "Nbin = 11\n",
      "Nsample = 160\n",
      "\n",
      "np.random.seed(101)\n",
      "microstructures = np.array([np.random.random(N**2) for i in range(Nsample)])\n",
      "responses = np.array([fipy_response(m, dt=dt, N=N) for m in microstructures])\n",
      "binnedMicrostructures = np.array([bin(m, Nbin) for m in microstructures])\n",
      "print microstructures.shape\n",
      "print responses.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(160, 400)\n",
        "(160, 400)\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We use `numpy`'s `fft2` to calculate the $P_{a,k}$ and $M^h_{k,a}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "microstructuresRS = binnedMicrostructures.reshape((Nsample, N, N, Nbin))\n",
      "responsesRS = responses.reshape((Nsample, N, N))\n",
      "Fmicrostructures = np.fft.fft2(microstructuresRS, axes=(1, 2))\n",
      "Fresponses = np.fft.fft2(responsesRS, axes=(1, 2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Fcoeff = np.zeros((N, N, Nbin), dtype=np.complex)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We calculate $\\beta_k^h$ at every point in $k$ space."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for ki in range(N):\n",
      "    for kj in range(N):\n",
      "        Fcoeff[ki,kj,:] = np.linalg.lstsq(Fmicrostructures[:,ki,kj,:], Fresponses[:,ki,kj] )[0]\n",
      "        "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reconstructedResponses = np.sum(Fmicrostructures * Fcoeff[np.newaxis], axis=-1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For a loose check let's see how close the reconstructed responses are to the sample responses."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn import metrics\n",
      "mse = metrics.mean_squared_error\n",
      "responses = np.fft.ifft2(reconstructedResponses, axes=(1, 2)).real\n",
      "MSE = mse(responsesRS, np.fft.ifft2(reconstructedResponses, axes=(1, 2)).real)\n",
      "print 'Mean square error: {0:1.3e}'.format(MSE)\n",
      "print 'Typical values'\n",
      "print responsesRS[0,0,:10]\n",
      "print responses[0,0,:10]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Mean square error: 1.358e+10\n",
        "Typical values\n",
        "[  3134920.24329635 -14859922.85111164  24350855.96847351  16518931.7101288\n",
        " -13399504.96081172 -20608868.32626785  21532496.27938246\n",
        " -24030785.28692817 -11643521.47041003  27758619.45859079]\n",
        "[  3013635.84       -14692726.82641511  24504455.19776576\n",
        "  16802677.09155771 -13414803.84246958 -20575614.08        21450135.01626968\n",
        " -23919938.87934423 -11691192.67329024  27626894.75236416]\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Exercise 02-1"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Test the influence coefficients using one sample.\n",
      " \n",
      " - first create a test microstructure\n",
      " - use `fipy_response` to create a test response\n",
      " - bin the test microstructure\n",
      " - take the FFT of the binned microstructure\n",
      " - take the inner prodcuct of `Fcoeff` with the binned microstructure along the `bin` axis\n",
      " - send the result back to real space to get the calculated response\n",
      " - compare the calculated response with the test response."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.seed(103)\n",
      "test_microstructure = np.random.random(N**2)\n",
      "test_response = fipy_response(test_microstructure, dt=dt, N=N)\n",
      "\n",
      "binned_test_microstructure = bin(test_microstructure, Nbin).reshape((N, N, Nbin))\n",
      "Fm = np.fft.fft2(binned_test_microstructure, axes=(0, 1))\n",
      "Fr = np.sum(Fm * Fcoeff, axis=-1)\n",
      "            \n",
      "calc_response = np.fft.ifft2(Fr, axes=(0, 1)).real.flatten()\n",
      "\n",
      "\n",
      "print test_response[:20]\n",
      "print calc_response[:20]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 10321080.57374739  18170943.91955532  21953730.26038766\n",
        " -27668437.71011165  -1222554.20968084   4850834.021033   -26983622.32046547\n",
        " -28409259.90543952   5996828.80698766  10823111.24422134\n",
        "   6559197.50876078 -30826019.2029601   -8894200.91717253\n",
        "  -6046499.57272069 -12751789.16653129  33920217.9478319    7588871.20941563\n",
        "  10018364.24439061   2116472.25868405 -25814016.94703518]\n",
        "[ 10136893.8         18057017.95942609  21878377.34953976\n",
        " -27708096.82038447  -1299137.4447629    4855990.40000001\n",
        " -26852174.36448368 -28324788.67386933   6094005.82172765\n",
        "  10782186.63208983   6469447.32       -30907790.79029018\n",
        "  -9013656.85033033  -6133981.44410236 -12830583.26266264\n",
        "  33932517.76000001   7771430.67336559  10094396.05258644\n",
        "   2064878.55760657 -25968923.07545603]\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Influence Coefficients in Real Space\n",
      "\n",
      "Let's view the influence coefficients in real space."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coeff = np.fft.ifft2(Fcoeff, axes=(0, 1)).real\n",
      "\n",
      "plt.figure()\n",
      "plt.title(r'$\\beta_s^0$', fontsize=18)\n",
      "plt.contourf(coeff[:,:,0])\n",
      "plt.colorbar()\n",
      "plt.figure()\n",
      "plt.title(r'$\\beta_s^1$', fontsize=18)\n",
      "plt.contourf(coeff[:,:,1])\n",
      "c = plt.colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEQCAYAAAD8jMw7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnc9vHNeV77+MA22eQbaaAmY147CtLJ9GbctePMSI4Sb5\nB7jlKFkGz03nAV49i2h68AZwJoAtSvYuQGKW9oPQ3f4DxG7AHvtlo3Z3YrydrepkFgFmALaaygDB\nEHb0FpxqVRerbt0f5/6oqvsBCKl/1Y9b937r1LnnnrPy+PHjx/B4PB6Pdr5j+wA8Ho+nKnjB9Xg8\nHkN4wfV4PB5DeMH1eDweQ3jB9Xg8HkN81/YBeKrJYDDAycnJ4nW73bZ4NB6PGbyF67HCwcEB2u02\n2u02jo6OlsTX4ykrXnA9xhkMBqjX64vXzz77LA4PDy0ekcdjBi+4HuOEYYharbZ4XavV8ODBA4tH\n5PGYwftwPdqYTqcYDAYAgPl8jkajgXa7jZOTE6yvry99dz6f2zhEj8coXnA9WphOpwjDEJ1OZ/He\n5cuXsbm5mWrRxi1ej6eseJeCRwvj8RitVmvpvTAMsba2hkajsWTRzudzvPjii6YP0eMxjhdcDzmT\nyQTPPffc0nt7e3s4ODgAALRaLcxms8VnDx48wObmptFj9HhssOKzhXmo6ff7aLfbeP7557GysoJ6\nvY7XXnttyb0wHA4XVu7KygpeffVVW4fr8RjD+3A92vjiiy8AACcnJ2i1WkuCm3Q3eDxVwLsUPKSc\nnJyg0Wgsvbe2toZ6vY7hcGjpqDweN/CC6yFlNBqh2Wyee382m/nVZJ7K4wXXY4TxeOwnxjyVxwuu\nh5QwDM+91+v1sLW1hdXVVQtH5PHwEYbhwu11cnKCyWSC4XCI6XRKtg8vuB4ywjA8t6AhDEMcHBzg\no48+Wno/CAIMh0P0+32Th+ipKN1ud+l1v9/HcDhEEARL70WRM6PRaJHzg3IVpBdcDxnT6RRvv/02\ngiBAv99f/Hvv3r0l63YymQA4i1S4f/++rcP1VISDg4OlG/t4PAbwJFIm6o9xl9e1a9fw4MED7Ozs\nnJsEVsGHhXlIWVtbWwr/SqNWq6Hb7eLo6Aj7+/uGjsxTVXZ2dtDr9RavDw8Psb29DQBoNBoYDAZo\nNpuIL0k4PDzEr3/9a5ycnODg4AC7u7skx+ItXA8ZvI9eGxsb+MMf/oCtrS288cYbmo/KU3Ymk8mS\nn3U4HDIjYubz+VJ60OPj48XvRqPRIrRxOBwiDEPSyV5v4XpImE6nXI9ek8kEv/nNb3Dr1i10Op3F\n453HI0uz2UQQBNjc3FxM2q6trTF/k7bANm7F6lqY4y1cDwnj8Tg1/jZJo9HA1tbWYsLCuxQ8FHQ6\nHezv7yMMw1yxrNVqi1weDx8+xNdff40gCLC3t8f83Z07d5SP0wuuhwTemmRra2totVqLpb4+VMxD\nwXA4xGuvvQYAuWFcN27cWFjCR0dH+OlPf4pOp7MUFpZkMBjg6OhI+Ti94Ho8nkITj3rpdDrnCpT2\nej2MRiPcvXsXABZPYsPhEH/5y1/wpz/9CcDZ01daHDlwlmCJAp8tzOPxeABsb2/j9u3buHr16tL7\nk8kEzWYT29vbuHfvntI+vIXr8Xgqz3g8xvPPP39ObAEs5W5WxQuux+OpPMPhEO+999659yeTCWnE\ngpawMCp/h8fjqQaqns21lRU84vzu008/jT//+c+L1/GFDcPhEK1WC/P5HLVaDWEYIgxDHB8fYzab\nLdwLsmiLw/09gCtvpX/2aP/Cuffef+qm1H5ufvu+1O8oWe2eMj//8oMn//8VgP+V8p2stmKR1o5J\nZNs1j/1Zl/l5t/4k3Ct+jZJtFbVNWruw2oTn3AGx83ehL8WJ2uqd3wL/+zO+8+XBVJtk7SfZd07X\n2TGzPDzCmebw8Pf/8R+L/w8GA+zt7eH27duYzWaLFWmbm5sYjUaL6JsgCHBycqJsTBp3KfAOFF7e\nf+rm4s8GeWILZAvHlbee/OnCNRHhvf4UYiuKrT6URXSe/7n1FOl2efuEa31HB5ubm5jNZvj6668x\nm83wyiuvADhLXhOn0+ngq6++SvXxiuDESjOqjh7fjmudJRKQv/ktcOV/qG9Pl+jwkGfdJnn/qZtW\nrodrAirDo/0LwD99S77d6HroaiNe67ZqaBNcnVYbD66K78t/a/sI3OLKW2duhWuJ97KweaOxxQ9+\nqG9OJEt4XRozZcKoS4HSdytC3O1AuT8ed0ISCsEVFR3XBk/a8b/wX/9S3ajLYN1GvPRD/cM03kdU\n+4u3brNxwqVgmqhDuCZEZWJ/1l2aOCsL8Zts2axtivFQphudDowJri3rVhcy1q2HTZ5165rAZfUB\n147TBbx1e0alFz6YFnyKgSi7DSprXnbgsNqaOlJD5bqyfst7k13tnqb+lZ0iG1CmMGLh8lq3ZXwM\nTZ57/HUVBmEWj/YvCJ0/tdVoo69F51s2CzhPaFk36dM3q5UtzjkfbnRxXBbePKHgHVDJ71Ftt8rw\nWFm2H29Xu6eluJbeohVHu+CKWLdpr10W3giqwaPb+r357fuVHyS2xTaiyKIr0oe8dbuMcxZuHN1i\nK+vXNDFQkuJLtc/4ORdFfCmX8Xbr+86IbtGQ6S+s9r7wy0eVE10t+XBXVlbwODYJImPluiq2ZYZ3\nQPEKVvIaprU5lRtFxupi9bGs/kHx5FE0y5bqxpzVb07X15ST16ysrPDnUoB6shxZnBNcEy4EW2Lr\n6uq3LFgDTSR5TZzkebMETIfYilKE66QDnW2alrzGC67KRhOCC7gVh2t6EPGcp6sDW1Z0eQUXOC+6\nIhagiT7k6rVhoeKGMtGm8b5TJcGtXByuycEjspTYdtazLEyLjWtia3I/1PuVcX+YOtdufb8QE+LU\nWBVc04PZtNiq/NZVAY6TNWBkB5KLYmt6f8nrLrPfIsV3V010jbkUAD1uBZHwsaxE2JSTGGX1J4q4\nFvKuha7kKCbQcQ2oXU6yLhpb7fqLlXcr41KwHhamEhsqEt7D6rAUYVemsp4BdoSXdZ1MhlrZtvhZ\n+xe5LrrOQ9a6td2uVcGo4FLGk6qIbVqnlD02Gx1Vp/CyLH/bCydcFwWdx2cribuHFicmzUQ7UprY\n8gowywIQtQ5sC4Du/L5p7ZF1rSI3gg6fnOu+bNfJ69dlc9G4jFEfLpDtT6IIuk8b7CJxn3HyrF0X\nBUC288osPlA9f95jdbGdbZLnGsuC1Z9NtnHa8a9997QyPlxtgitatRdQyzoUwVrdlFYxVqaci44O\nKnNuWahMsGRBLbo8x+jF9jzUgutCHHOVBFebS0FUbItE1R6DbF2zqrUzDyz3StZ1sim20X78zfMM\nZ1wKgH4LFxD309p8FJPJLSEjUrL5DGTPX+YY/YBNh5XzwbbQsogfd5UsXKNRCkW0blkdV/esvegE\nVFEWdsigu6y3q+TFmWdFrLi0lD6N5WN519pxRHS7Xezvp7fxeDzGdDrFbDZDp9NR2o8TUQoReYJh\na1UKywJ05bFX5ThM3ghV28uV9jZB/Alnf9Zd/KXBemz3j/RsDg4O0O/3Mz+/desW2u025vM5JpOJ\n0r6sFpG0gWhplwiblm4eOkXIlesWx4S1m9emJt1JWZ+lGSBxi9eLLB87Ozvo9Xqpn/V6PbzwwgsA\ngN3dXeV9WatpZgrKgHEXRZfq3GRvRBTILrO2ae3qFH3emHKWj9+LLQ2j0QgAMJlMMBgMlEXXKZdC\nEchzL5gUARfy+sbJEwqe4y1ahVvqayC7RJrlbvCocenSJTSbTQBguh540C64Lj6Sqh5TniCYEEId\n+2BVGOZFdNBntWWRhJfqRkshmEUTXtePdX19HRsbGwCAWq2G+/fvK21Pq0tBZsAWxfeUl3tB53lU\naeLIpdLiecei4mbILD+TUvPrwi8fCW3P1RSIpsT2PoCR4G/m8zlqtRquX7++8O/O53O8+OKLSsdi\nPVuYKC4VAbQhurrFVpcvV6aeWdp3bQhvVo4JKuEVEdus91ki7FoFbF3jN2ux1RUA/zP2+tcfLH/e\n6/UwGo1w9+5dvP766wCAzc1NjEYjbGxsoFarod/vYzabKftwtS18OPlGfmCoLIBQKe0ii6m8C65Y\ntnk+XJ6cFoBa++sWXpFj07GYh6KabZYI2xJe3UUksxZbnfvuByVc2lsldPt0TU/GsXDF3aPTvyu6\nbeo4barS4advrqZuKy+ml5qi+ZV14qRLwZQfN8syobJuAHcESjfJIP20MCWKm4YJl0J8Hypl3GWt\n2wu/fKTVwgXMWrkuuQFt46RLATBTolsnZXMjRIiEhKm6FlyYKBOFIs1ohIjo8kykAWyhTV4bamPB\nuxQctXCLDGUndU1sReG1dG0Ia/w6uVpjLcvS5RXXOCJCG3+f8ty8pVtQCxeQmzgD7BUBlME1waV8\n6nBt0QaLtGNVueYmRYdqLFD28eT5ewvXw0VZK/TqIGntmi6IqVq2vmjICm1WyBultVtlS5crSqHb\n7aa+DoKA/ogKgO7sS2UT24i0QaZbzFzLlKVTaLr1/cVfGqxol+TKvqx6dmXtm6bItXCj1GXxXJFB\nEODjjz/Ghx9+qPXgdCA7W+7SoC0yabG6lP7UtG1S4eLqLZ7jkC3Lw7J2AbU2rqqVmyu4aanLgiBA\nu93WdlAuYSvRdtlJm1AD0ttbt69RZuDLzCFQCQxF9Q/RVX5pE5uqboYqii7XpNn29jbu3bu3eB0E\nARqNBsbjcepSN4pJM0B+kkb1rl/W3LYqqEyYsZC1FFVDmIo20EXaSVZov0wseU0uldVVbuoXK+/6\nSTMWUZmJo6MjDIdDtFqtc9/RPSmi4+7I02l0igvLvxaniPGpWWRZunnIDnDZ60cVniUClcgCYkIb\nfz8uuqwcEkVJOmUbYcENggD1eh3tdhvr6+sIwzBVcCNMz0bzEvfl6hba+DZEF2WkDRTZhN2uYsI3\nKnr9eBYdxL9DIb4y565DaJPfSVq6ulwMVUBYcBuNBq5duwYAOD4+xtbWVur3Pn3ns8X/n3n57/D+\ny2aEV6TTmhJaDz+2M1vJLqk9fXNVSHQpzk9lLPGIbfy7vKLLwx8++SP++Mm/Sv226OQKbjJ1WavV\nWmQ9v3TpEq5evZr6ux++89K596jFNulWkHkE83dkt9AltPHtsm6icdHkFV8eodVxXir5Ka68xS+6\naWkPVZ6svvfyM/jey88sXv/Lzz9jfLtcaFtp9o+P/2HpPZ2hWKJ+QNEJF92zyyqzyaZcCrJldXhx\nLWVgkqT42hLZLKgmytIQFVxRI8ZPmhHjkv827VhM+J5cidt0EZttk9w3K4OX7DZNwJorYSWVz7N0\ns5J6e+TQLrgqYssrhLo7uEpEhM5jK8OEmWs3Il7XA+t3NslyM8hU8sgS2zL0O1tUKpcCS/zzxF1G\ndFXjgT30iLiT8sTXFZFNkmXtZomuiD/Xo4ZWwXVJTHiOhcq1QDUQXa9aq+q/1SVYIv0u/l1e8S0K\nWakweURXxZXgw8Oy0Sa4RRNbHvKs3CIOSmp4w6NU2kpX3+IVX5XtykAtXiqWrncnqFF6l4JoZ5e9\nO+sIWi8avGFUpvMVyyAjvjpvBLKiK+rT9ZNkeim14OoYACqxv2WCZemLLgIwJbIsF01eufsIqtps\nougQXVG8dauO81V7ZTuKanQEi7y8oyZwrfPLWrcmxCuZ65X1Hd0VmFVQ2XeaWLvWh6qA84IrA8Wg\n0D2wXHl0NonpGxSPgLJ+5+KkJW+/4U32rkt0q9i/eSiE4IpePIpJBl/RQT++HfQQia1t0S0K/X4f\nw+Ews4JN3uciFMaHK5r7QCVLmStFIdMmNmxXuJXFhvslaitRS7UsAsRbpl4XRQgPG4/HAIBWq4Uw\nDDGZTNBsNhefTyYTNBqNxXvJz0UphIUbR8ba5a1rpaP+VVQHSrajxwe/y0KgOzesCo/2L3C1He/3\nXEf05sZ7Q3LRxaLK4eEhLl68COAsE+JgMDj3naiGYxiGSmILFMjCjSOb6SvL6tUhspQUTQR4xNeG\nOyHN4i1a2wLFsByLwnw+R71eX7w+Pj5e+rzZbGJjYwP1et1tl0La3ZC6c8t2PFdcBkWgqO4EFkUU\nWVF0tbloHtwi3BxYmcPm8zkuX76MIAjQ6XTw3HPPYWNjQ3pfRi1cnkcS0cFgO69tGUVWFpfdCp7y\nk6Udn336V3z+aVxUv138r1arYTabAQAePnyI9fX1pd8GQYA33ngDq6urqNVq6PV6qXUceXHOpcCq\nm+QKRRDZ5A1Id/J3EYrQfmXEVrvbtnJf+uF38NIPn7y+9U9PBPfGjRsYjUZotVqYTqeLCjbz+Ry1\nWg0AsLp6FmMeTayp4JzgRrj66CLbaU36DdPawVZtOepHWxOuKg8blfI6rtFsNjEajTAcDlGr1RYV\nbDY3NzEajbC7u4s7d+6g0WhgNpstCujK4qzgAuW4sKwikDrOLe+mIyK8tn1vIrPnLveTeDtS3fBs\nW42iuHy8kYjGi+GORqPF/1VcCEmcFlxA72BSLc3DgkcsKM+NKmKDEpG2NRlypNvdkrWf6D2XXCqy\n7e76Tc5VtNU0+33K+6qZiHgusEwlXqpE4SqiodJ5KSyHtPNTrfXGG3hPJbYU/YNCDHmuh6n95O3T\ndJ9NO2aqmmYn3/Adz9p3T8td0ywimWtTVIAp7qppIpFXmjtvcFAIhsy5UT6iyVi8rIkzG6ucWG0o\nukIR0Luk3JZP3WMX510KtjEhtvFt2Z4oFN1esgyN7eWklPCKoso1sO1ikKl1Fv1OFFd9uCaxJrgy\n7gWK0sxZVpnMbLoOv2ORfWM2Ui/yoDLQTZyDjOiqLPjhLbuT/I4sXmifUBgLl0JsI3iTiLMGgSvr\nyl2c/aVI0iNDUW9UgJjoql7vrFpnAF0EjWt90hWsJK9xoYxHJLIuiq3o9l2xJAH5YymyWPLAs0jE\npEhlJWpSTeCjIwFUmSiEhUtp3cZxUWzj+7HtzxXBJdGngvqceMIQ8yxdXT572XP14iqGccEVtW51\nia0MrrgRsrAlulTCROVacI24dasiujqvrWmhryqFsHDTKFrGr3hIHO9NpwgTaPH2ovD/lU10s8IQ\nRSdpTQieF1X9GPXhUlm3NsRWVgS+/OB8/LEILvtz08Q2+r9KTTAXbjIU7cjy2+b5dKsgfvuzrnQC\npKJizMKlmihL64h5Cxd40CW2We+LtIdpfy6PBZYltklkrd6yWbpppePz+m0ZRLdqgpqHsy6FtAHK\nEtvo/zKiSy22KhYtFbKimzdAWELLcpsUvcqCKGntmCa60XddS9KehRdQNYzkUqBwJeSJbRKqpDSi\ngisitpQTiFmIiG5ejokssc0756zz5Dkfkfbn7TcAW+Ty+kNyP8l9JNvx9M3VpddZidptia5tET1d\nX6tMLoXCFZGMyOskvJ3IltjKoPMRW8SyFSWrXXRUAEmSddymQgKTYpv1HqBf+CKfafLPYw6tgnvl\nLX3WLQ+qnUl0sEfnG/9jfU/38Yi0m6upFAG+YHzW57w3i6zvJScAKc4/S3R1UhS3RZnR5sN1YTUZ\nD7qThxSlHYDlJc/xwenKooa4qIr6hFm171xY6BKh25+bTDakY7ui/ILsKNzH2UkzFiIdpUgTEmnI\nPlLLPhXoaCsdNx3ZdklOJrKsWp24UHAzutYywlvkMWWTQgquSWyEJ7mUOMQV65YSaouWFRFy4ZeP\nhN0Hpo2EPOH14kqH84JbhlhEHmyV2hGlTLGxScp8bjx4YdWPFcGljsPMim/kxVYS6CrEoxYBGaFl\nLUZRKSHvKTfaBNclMVF9RKN0K+hsF2rrtozuhCSuWrVFm3tQ6St+0sxhsiwHVStXJyZuPlVxvVBB\nLbS202Paogo3ZUoKJ7hJTMUzyli5pqx8HQOdqsKri2Fxpq1amYkzlymbyPb7fdRqNYRhiE6nI/y5\nCIVdaSaKSZ+aatZ8Tz5RZQGRm41KBrOs7enENT/wzW/fL53YjsdjAECr1QIATCYToc9FcVpwRS03\nFSuCwkq0IbSmrFsdyOaGSJ4zTxuYtmqL5H/No4xCG3F4eIiLFy8CABqNBgaDgdDnomhzKaiW7oh+\nKyIoLB8uxQBIcyvYtGSLLLai5J2r7XLjrpHVFrx9piptOZ/PUa/XF6+Pj4+FPhdFuw9XdSAkRTcZ\nchP5xyjElucmQZWxS7VD2xDbZGXXiMhPy1vVgqcNZcuAmxQK1nmw3AG6J3dZbVAVIRUhL3MYZWax\nQk6apYlu1vdk4Fn6KbIN6n3YtmzzhDfvd1lQnFeW6OpeMajjmoj236qLadY1+MMnf8QfP/nX2Duf\nLf5Xq9Uwm80AAA8fPsT6+vrSb/M+F6UQgpvmWsgLLqfyofG6RlQGXPK3RSnmlyW8ad/JQsf5mBBd\n3a4kL7Z0fO/lZ/C9l59ZvP6Xnz8R3Bs3bmA0GqHVamE6nWJrawvAmSuhVqtlfi6LEcGleNTjFV1d\nkxVpwqtL/LKsX1cLZ8r4tlVLAEWwytNkiS5AO4lGNc8gg0zyHR9B84Rms4nRaIThcIharYarV68C\nADY3NzEajTI/l0VbxYd/fPwPS+9R3YFZ2fXLNDNsCpOWEa8wiYZDsa677PnxChZvJZIogoZHcHn7\nMVWsdBomRZmiAkOa5mTxi5V3rVV84LJwu90u9vefdAKZQGCqCY2sSTQvtu6SJ7Sq8aas6y/b72xZ\ngTz92ERKSYqS957z5MbhHhwcoN/vL15TBwJT4MVWDp3Wbd7CBOoSL6ztmHD9mMJ0/l7qxSJVJ9fC\n3dnZQa/XW7w+PDzE9vY2gCeBwM1mk2tnuqxc05TBotYltrqt2bxt55Uctzm5lOdOkHGNpIkhT209\n0SXX3uKlQXjSjDoQWBbTopsUirRSNFWGdS14RJYq+VDezZDqps+y3GUwKbY838sSZFZaSk8+UlEK\nthzOSXSLLs/gKeKkHaWVpyI8yaXYWUuzRYXYlOhSYVpseYi2lSa83tqVR1hwVQOBXevsSVQeeYtg\n9eqKFomQEVqZ7+eJMI/osqCOhZax3nmjECiFNm3bLGsX8MIrgrDg8gYCf/rOk+DiZ17+u6XAY9cQ\nEdn44M8aREW0enlRdR2oYirnsWoRToq2SDNO0mKer7ylV3Sp+ezTv+LzT914SjZNbhxur9fDzs4O\nbt++jddffx0AEAQBGo1GZljYysoKLhyfANATI5m2MosqkD4PHsGN46Lgqli5OgRXxoJlkdfmuq38\nOCJtknXctixd1aXavFQpDlfbwodIcAF10WV1bIqVWLyDQkYYXBRcQF+OCNs5XHXc4PMoi/DyRC7o\ncB9USXCNLO2VCdeREU/b4WIee5iyavO2ndX/RApLZs0FyLgYqCpueD8tDUYs3DiUFp9qzs8I1Yme\nqli5FDczHRawDas2D8qwsbTz07msN44JofUWrkbKsGigLNiIGDF57XlCqnQJSrRv0Sx3aaSNGV5r\nVxZv0erBuIUbR2XwUadLZA0CnjAm3dUmZOC9ufGKbpHcNaIWoOnKypQLJFw8Vx6i9qCwOIti4Vqt\naaZzgsW05eZaVdaobanauChim1V/K8/yM5EvIH5csjfhtOuZdm2i+nppfy5QlP5EjfUikpQJTJJQ\niG6qkP7zl0LbMD1zn7UMOQubnT8SyOSf7LbS4BVTE4laRM/t9M3Vc30wbcyIVjC2SVGOUwfWBTdC\nRJREOq02SzdFdF2zcuO4JLo8wioiwiyrNitsKvpLw1R2rDwrN96fsoQ3iep11N0Pqiy2gGMldlyb\nUKOwbk1jOw42C6oscbzwxqdmhU9RLVtNm6CLhy9mTaJl3byTSX7SQsgoRddE3o0qYXXSLAsTq9OS\ncHf6pOD+5Mq5r6RNoJm4kfCIrUy8quxA4blWOsrOU6csFD0mmWoR8Wu31O+i/sbRz3T3MV0rFP2k\nmcPoHPy5OG7d5g04ngGZbF/TVomux/m8BQCsz0WPycTElKmcEnFU/MQuJ6wyiVOC263vO+VSKCJZ\n7SfSrnnVGkS2Y4MswcsSVarVWMljSB4HqxZa/PosielPrpyzbm09QaniRdegD5eyQ8gG7LOW/soE\npJtAJldE8lyKMBjjUCS5zloEkMysxSu2ssdEbe3qtmx1z6NUffm9Ngs3slZds1qLdJdNHqtIyFTU\n5i61fRLd0QB5YidbZkYVHsFJE9YssaW6xqYMjiKNQWqccimIoHKXLMJy1rwwKF37LRtponvlLT1u\nBFWS1ysusCbFVkZ4RWOYqyq6hRXcMsPTGcveYSmtX8rHetXjoqwW4aJlKyq6RejH/X4fw+EQQRBI\nfR6n0IJbJCuXFycWdRjAl97mw5bYqohwma7teDwGALRaLQDAZDJZ+nwymaDRaKDVaqHRaJz7PIl1\nwVVd2kkpurYfwWUnAqmFN1o6ansSsaxWLguePuiaZZvWHiaWSZvg8PAQFy9eBAA0Gg0MBoNz3+l2\nz9oxDEM0m03m9qwJruxyzTTKMOupKppUopuWh8EF8c1C5NrbFl2X+qloRWpZii668/kc9Xp98fr4\n+Hjp82aziY2NDdTr9aXvZWFccEUtMl6rV0fyDpkQHKrqrDKobIdHVF0WXl5sZcsS6ZssC1b3U5iO\nfCBFF13WqrT5fI7Lly8jCAJ0Oh1Mp1PmtozE4VILCqvziq4D1x0XaLr0i8z5yCTEjpARAJEBSBGT\nm4QqUTfvsVH1L92uhEhs4/ka8uJyRTKx6b7ZZZ3XXz//DI//7+eZv0ub7KrX62i326jVapjNZgCA\nhw8fYn19/dxv33jjDayurqJWq6HX62F3dzdzX9oFV8fEDo/wRp+L7p+5AOInV8iW95qosUVd3DDv\n97Z94CLEB79OC0yl3LqOts0T2zxU+q0J0U3jOz94CfjBS4vX396+tfR5WuXxiBs3bmA0GqHVamE6\nnWJrawvAmWVbq9UAAKurZ23XarUQhiHzWLQlr4knktApLmkdWmZ/ye1kJhMBzkQ3Z7mlSOVV3ZgQ\n3gjq+lt5A1RHmSARAWYdn0s+2yQs4Y335fj1VK2jltVWVDXNeBNmna6vCe0vCAI0Gg2EYbgQ52vX\nrmE0GgFJnFr2AAAUrElEQVQA7ty5g0ajgdlsxhRvwJDgAvpFl2r7rDIoWVYAj9jKHh9lNi1TJc6p\nRNeF6gRZx+i60PK4lvKud/I6qiR4Z7WX64JLiTHBBYoTNyoiujrEVrUD52FCeFVE1wWhzSLvsdgV\nsY0je73zrNsI1fppXnBVN5qTm7IIwpsnuraENg0ZgdItuiIuFRNVdE3gotjGEbnmvGILyD0FxPGC\nq7pRjmTARRNdgM8SiFCtTiCDa8IrKrpFFVvXhTYJ63iTUQm821VxfXnBVd0oZ/b1Mogur58rie1M\nWXFcEd0iUjSxjZN37LJGg+iN0wuu6kYFyl0UYeClia6sVWs6CJxKeL3onqfIYhuHKtJHFi+4qhsV\nEFygGAOPNbhcFNokFMIrK7o6atTZpCxCGyd5Tl5w9eCE4EYUYfDJdEzbYhuHV3i96J7HBaEF+Puc\nrE/f9PXwgqu6UUnBBfQOPqoQpGjgyc7e8pJWaZYqcbaK8FZNdF0RWkD+Bu/yhKQXXNWNKgguQLdI\ngBcdnVH2WHjKeUeoiq8Na5dniaoL4uuSyEZQPE25KLxecFU3qii4gLkY1iSqHVLlOETENsKE6FJN\nplGXGtJBWYU2jmui6wVXdaMEghthI54VEO+U1H7aPPGlcC1QZboSXSIqi27xLaLYyvQ7L7hecHOx\nGd+a10F1TorJlPTOQ9WVkEaa6OrKHqZLeF0TXF1zBF5wveByYzsqwFRaP13oENs4eflTk8Svp8w+\nKcXXJcGVzV3Agxdce4JrJAE5JbwRAroooshG6BZbgN+qTbt+MsnTbYQx2aLIfc9zRuEs3CQqKeOq\nhAmxzUNnjTpdK66S6AiJS6Kad5YHl6xcb+EWiCwLh6qMStEpktCq4Lqly+tqSZ6DaB9OTrZS+f09\nNFgvk05Blli4dBc3yaP9C4s/HnSJrWwJd5VkLDpvHKYLaIqI7ZcfyIUVpuGSL7tslEJwAS+6AIRE\nNoJ6cPFWWebZjiwy52RCZPIEW+ac84RWVoS96OqhNIILmBFdSkuCAlFrNo4OsXUF6nPTbd2KuhJ0\n9cN4u+l+YqgipRJcQM+dOerc8Q6e9p4pVEQWyB9I+7OukMBQWLRZ282DdZxFFQuW2Mr0ubTvi7gr\nitqOLuJciR2qi6s602tjmS0LKitddJlu3kSPatnsOFnnyJPLQTUxjq5cwGnkJa1P66eqN/a0vpnW\n3pRJyXnxUQoEyF4YkVhM1oBLm7UuUuSCDt+zjKiILmRQISuloEx8bhyXIhhkK4RceUv/0xRPG7vU\nlkWk0GFh3fq+Nt9a0iLI6uzUVq1poQXULDibAzC6/jw3BNZxmnpklhXbCN4+mfc7Gcu26nS7Xezv\np/ez8XiM6XSK2WyGTqfD3I6TPlwX76BX3nryF39NgapPloWqr5ZXzKjIegLJ6hOmrG9VRMSWtx/E\n+2RWX/Riq87BwQH6/X7m57du3UK73cZ8PsdkMmFuq9AWri1MJ5GRwaRP0hQ6XAsmxEbGspVxf+X1\nSy+2cuzs7KDX66V+1uv18MILLwAAdnd3c7flpIULuGnlUqHLmo2wkb/WhJULqPcLkwLTre8ruxGo\n8GKrh9FohOPjY0wmE9y5cyf3+84KbpmIuwx0Cy1rRl+3ZSsygFUmL4twM5at6hyHqq+IiG1Rn35s\ncunSJTSbTQBguh4ASZdC5EAOgiDXSayC6iOkDWytbCubCyGvCKJK34hcCzr6FnVJ+LQ2UKnuwNNP\nZCJTomMq4srOIAjOvVev19Fut3N/u76+jo2NDQBArVbD/fv3mb+TEtwgCPDxxx/jww8/lPl5KXCl\nY+kSWtnJqKKEDYmILU80DLXQslBZ8JIF1U1ZtlowBadvrqZ/8G+fAP/+SebvZIzG+XyOWq2G69ev\nL/y78/kcL774IvN3Ui6FIAjw1Vdf4ZVXXpH5uRCuDV7dbgFeeFaL2YLKcsyz5FzpGybFVhbRviLS\nf5LXyblY9795Gfjv7zz5E6TX62E0GuHu3buL9zY3NwEAGxsbqNVq6Pf7mM1mePXVV5nbkrJwZ7MZ\nhsMhxuMx18xcWXBBaAG9cbVFIxI1G66nIggtYMfdZNPSpeb69eu4fv360nuj0Wjx/8hC5nFBSAlu\ntIOjoyMMh0O0Wq1z36Gst2Tbl+tKxzEltBSxrVSuBd6Bq1qqRxQRsY2PBdN9SUVs03y5okVdXRk7\nriAsuEEQLBzK6+vrCMMwVXDLgEudpUhiG2HLn6s71lZWbOOvdfUtmWXxKuQZVl50lxEW3EajgWvX\nrgEAjo+PsbW1lfq9d3775P8v/+3Zn07S7sayA97FDpL36KyyzFnHaq28djeZ18LGYocIyvOkPO7o\nmrP6jEySnTSSovvZp3/F55/aSR5jG6lsYVGs2XQ6xc2b5zvBysoKTr6hEy2RmEHe+EfeTFWuovKo\nqHs5LFVJe93XRLYMvKjvVqWvmbo5JPsMldjGyTpvqmxh+AnnNv55pXxl0qkEV1ZYXFnhoxMRN4Oo\nyJpqr6yBa+oGSJGZLk6e6LootnF4/LaU8zOAF1z1jRIJrursKsVqnyKgM4ewbtIGr0sTS3GiPicr\nuK6LbRpUYhuR1gZecFU3SiC4lKEsVbB2AXvlxVWxOYsfIVIiXWdSdt5jMQG12EYkr7EXXNWNplR8\nEOmEOuIGq2LtAnay9qvgyrJQKheN6QKYOtAlthHxa10lwTWWnpHVkUzEUGZFMcjgmmAlSYtooJrI\n4kFUOE0LbZbA2473LrLYRsnQdZaZKgNWaprJ4lKsqUnR1W0B6grPsm2xJuHxF/M+XelwK9gWXNk6\ngMnKE7yiG7W9t3AdhHIJosk6XbKkdXTqUDbdcbCuuApY55mczLJl5doUW95wNmqquChCm4UbVdCk\nEjZda75lj0+HhavSwXk7rs3EIrZcBzzwWrq8E2cAfx+xIbZ5x8bbdhS1/lY+gLdwqVCJBU3bBjUU\nx6cKhQiyrF/VwcNL3iAzZfHKtKctS9e02PLcBCj6y5cfeH9uGkZdCi6IGwsR64UCndZmJCCmxDa+\nDR7hLYI/WrfomhRbk5Omnmy0uxTiqAiZ6ZSDusJ+stAV4yiybd0WboROC5faLaMzttmU4OoWWz9p\nxo9WwdVhKZoSXtOCm4Sn84sKl6zwUj4amvLjiooH67hEhJEq3pwSSjeCKn6lWYHCwiJ0i67ulUS2\nsfHYaGM2WmXSLImOxSQmBNe22PJc9yoJbmHCwuLwpJbTRdHFFngyCGxaNab2S/WkwPLnutonbIht\n1cK8RLEuuCq5SlVywJadqB1Zg85ETlrbA1DnOboqtID5Y7N9nYuCMcHl7QCiNapsWrsukdVeecKr\ny9p1aQCyRFfkOONWbhnEVudEbRrZY/pdkuMoAs5mCwPEfVy6l/66NMh05YGwkUtBFlHxcyEVZBId\nflyTrgSKyI5frLzrfbguYDuZiGtQtEVe2SHbApRHWhvwllJKWvOun6sMpsRWJAeF5wnfsX0Aedz8\n9n1uC4YiDM21BRnvP3Vz8SdKlsUvuz1b8LSByPkUSWhFntpMPYF5sZXHaQs3Dq+1q2sizbQ7QbYT\nJ8+dtXrOZX+kzPmLFA0tkuhSQl2xgaaatNs+3CAIAAAPHjzArVu3Mr93584d7O7uMrdlTHApYiLL\n7mKgEtms72RZ77bKmacdB8U2XDgX0+h2JeQJbVknrYfDITY3N7GxsYEf/ehHGA6HaLVa5743GAxw\ndHRkX3ApEonE4RFd2cgFW5NllO4Cnt+4Zu1S30SrKrq6qKrYAkAYhgjDEJ1OB41GA2EYpgruysoK\n1/aMRClQV/kE6CMYiiK2JpKwmxQrXU8sRRBc6nPXFYHCE4mg0i9P19cKEaWwvb2N27dv4+rVq0vv\nTyYTNJtNbG9v4969e8xtOD1pxuogogOqW9/PFBlRsbWxNHZ/1i21JVE1bLjGKP3WybHh2mQzNePx\nGM8///w5sQWA2WzGvR1jcbg6k2urPJKLFJekrrgAiFWLpca2latTdFy2cG2ft8xY5I23lemr2i3c\nf/sE+PdPnrz+fz9f2l80KRanXq+j3W4vXmdNiEXWLQAuC9fowgfXRDcNUatWRXR5j1mqE7+5CgC4\n8MtHmd+xKbi2RccWLiSsoRJcgMZgIBNc/J7z238vtL+DgwPs7OwAwGLSbD6fo1arod/vAwCOj49x\ncHCAIAgWApyG0y6FOHmdRCRel7UN0X3rdi+Iiu3pm6sLsVXZdpmjQVRQud6m2jRvPzpdC2lE7rzk\nXxEYDAbY29vD5cuXUa/XF5Njm5ubAIB2u412u31mZJ6c5E6eaU3PKFsFlIVoZ+Ht5Kr+WplOTBPD\neAZLZGWtXECPtahbeHRZuCor1FyrW2YrX3AaZEt7NVm4lGi1cNMuhOlKs5Hly+p8omKbVhlB9Lio\nxFbEopXZj46QrSISv75FKUPDamvTVq7nDO0uBR0DTLbDp4mvrNhSiK4KIkKrIshFw8XyR0W9yYjg\nRZcPbYIbt5ySHY7i7qo6CFhWL69lKyu6KtatqkUruj+ATjBY23E15K3oYktl5fK0gxfdfLRauK6L\nbtr2RNwIWe+rHFeW8KgKrapIqwoHj9i6KrppUNxYTWHKtQDwufCqjHaXAkt0KaASXdZ28qrZUoku\nS2x1wyN2Oq5fVrIdFyhCRAIvtibtvPAuYyQsLEt0q5qxyUV4wnRkB4+J5csuDW7XxDYi67h4xqEf\nqzRYSc8YTy4iU3NKx8VXqX1FWUbcRVSEzITfVpfQmqj5ZpqsxD5ppZa8yNJjbOEDa3DlXdhH+xeW\n/nSRte2yC6quFWdFFts8iixGrGTuJsZZlTG60ozXn2tKYEXIEt2yirHqI3pZxJZqoY2LyYeKVvmj\nDGgT3KyJHpY/1wWBpdq/7fNIkrXaTCR5Dy9FDP9i4dq1pMYLrzms51Io0oVOWrOy1q0rkzsefuKi\nK5rIJX6TcfmGU+XcxKawIrjJTuea6LIsmkhky+BK8NatGFSWrsvt4K1dvVizcF3udACf6FKjK4MS\nK3lNnCKJbdEny1zv/1549WDVpZD055q8wDr3VQSfH7W4l9my5YHHnVBEvPDSYt2Hm+Ze0H2Ro22b\nzBvqEtSuBD8g5XBRjLOOyQsvDVYWPvASv8AUj5AudRhTJd953Ak6xVaHqFD0BVZlaJOwytebPo7k\n/12r7lwGtFm4vH5DXlTusKzfFsHKpWxLnsHNu7rKpRuYCNH5sZIViaDqTrBl6UaxwVn7Z30mOh69\nQJ+hreLDheOTzM9V7+giF0+14kOE6coPJgtHig4GCqHNOj9b5epVkRVdG9atTN/SWf9u7bunlan4\nYMSl4MIjUx5Za8wjWOvqVa3gNPdC1GaUwpt1HfLOPfldnceSxHWhVaEI48JDi9aaZjrQbY2JWLo6\n3A26rF3VbGCm3QdFFFoRK9em2FJauBTXyVu4FSZvUkC3T1eHtcs7uNMsXS+0/GRNhHbr+0vXTsa6\nL6q/3LNMoQTX5GAUecymJtpvlvAmYQmxqCUVP2+Tg7woQqsS3SByLZLtQRXV4mIoWpXQ5lI4+eYC\n+YCVGZSqx2BbCGSPnxXaowvbbaUbXneSzj6num3qpyTVa77aPcXKB6iMS0EqLKzf72M4HCIIAub3\nKAcgq+CjbGgPT8o8249ysmkSu/V9L7YMbF/XLPLasWjtzKJsyd15EBbc8XgMAGi1WgCAyWTC/D5F\n6ROe34tePJEMTpSrbD779K9Sv3O5OB/FMcm2iyjxaylyTZP9i9XfZNqDdW2TbSPb3jzWrYn6eVTx\nz6bo9XoYDof42c9+lvp5EAQIggB7e3u52xIW3MPDQ1y8eBEA0Gg0MBgMUr+XbFAdYkF50UxZu59/\nqv4o45L4qhxDNPBWu6e4/3++0ToQs26aLqQkzPtuWp8xfe0pfL9FE1oAGA6HGA6HaLVaCMMQv/vd\n7859vrm5iU6ngzAMMRwOmdsTFtz5fI56vb54fXx8nPndtAYWFQqWK4HnvTRYncektUuBDvHVMbES\nF9e8gUc5MHmul+z1pDhGlevmwg2XhyIKbUSr1cKvfvUrAMBsNsPVq1eXPg/DcGF0NhoNhGHI3J5U\nlIKowzlq7PhEg6lcAjLwrG9PHruuzi8S96saRpR0s+S1AeucKQaYasyzSBuwolJkz4XVx3n7y2r3\nFPhvT0ntQxRqdwJPu335AekutXBycoKDgwO8/fbb5z7rdDqL/4/HY/z4xz9mb+yxIN1u9/FgMHj8\n+PHjxx999NHj27dvn/sOAP/n//yf/+P+U0VkX08//bTUPra2th6HYZj62RdffPF4b28vdxvCFu6N\nGzcwGo3QarUwnU6xtbV17juPLYVceIpFt9vF/v4TK7rf76NWqyEMwyXLoWok2yV6HQRBpduFhYrm\npEVb1et1tNttjMdjrKysoNls4rnnnkOv18Pu7u657w+HQ7z33nu5+xL24TabzcUOarXakk+DN1ys\ninS7Z4/rvm3OODg4QL/fX7wWjX4pK8l2Ac76zPe//308++yzlo6q3HQ6nXN/7XYbwJnOzWYzAGfz\nV9E1mM/ni98fHBwsRJh80iw6wFardc5/AfgBk4UfNMvs7Oyg0WgsXvNGv5SdZLsAZ33nq6++wiuv\nvGLpqKrLzs4OwjBEEAS4ePEiXn31VQDA5uYmAGAwGGBvbw+XL19GvV7/rwUY2ZAt7T08PMT29jaA\nJwMmsoY9Z4Mmumt6ziMS/VI1ZrMZhsMhxuNx6uOsRx9ra2upbpzRaATgTHgjC5gHsgTkfsCwiQbN\nnTt3bB+Ks3jffzrRE+Xx8XHuI6vHbUgrPvgBk40fNGxqtdrCUnj48CHW19ctH5EbBEGw8Omur6/n\nxnl63IZMcP2AycYPmnxu3LixaJes6Jcq0mg0Fv7C4+NjvPDCC5aPyKMCmeD6AZONHzTn6fV6GI1G\nuHv3LgB29EuVSLZLq9XCYDBAv9/HpUuXKtsuZYE0PWMQBIvlbT5ecJnIwp1Op7h5080Vdh6PRy9a\n8uF6PB6P5zzayqR7PB6PZxkvuB6Px2MIL7gej8djCC+4Ho/HYwgvuB6Px2MIL7gej8djCC+4Ho/H\nY4j/D7x/ISsqtmFAAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5cbbf10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAEQCAYAAAD8jMw7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnc9vHNeV77+MA22eQbaaAmY147CtLJ9GbctePMSI4Sb5\nB5hKlCyD56bzAK+eRTQ9eAM4E8AWJXsXIDFL+0Hobv8BUhdgj/2yUac7Md7OdnUyiwAzAFtNZYBg\nCDt6C061qotVt+6Pc39U1f0AhNS/6sete7916txzz1l5/PjxY3g8Ho9HO9+yfQAej8dTF7zgejwe\njyG84Ho8Ho8hvOB6PB6PIbzgejwejyG84HqsEUUR9vf3bR+Gx2OMb9s+AE89GQwGuH//vu3D8HiM\n4i1cjxV2dnbwgx/8wPZheDxG8YLrsYZfc+OpG96l4NHGdDrFcDgEAMznc7RaLezs7Cw+X1lZsXVo\nHo8VvOB6tDCdThFFEbrd7uK9y5cvY2trC6urqxaPzOOxh3cpeLQwHo/R6XSW3ouiyIutp9Z4wfWQ\nM5lM8Nxzzy29t7+/j8PDw6X3vA/XUze84HrIiaIIGxsbeP7553Ht2jVsb2/j2WefxWuvvbb4ThiG\nODw8RBiGuHv3rsWj9XjMseLTM3qoGQwGS5NjJycn6HQ6GI1GFo/K47GPt3A9pJycnKDVai29t7a2\nhmaziTAMLR2Vx+MGXnA9pIxGI7Tb7XPvz2YznJycWDgij8cdvOB6jDAej7G5uWn7MDweq3jB9ZAS\nRdG59/r9vo+/9ThPFEULt9fJyQkmkwnCMMR0OiXbhxdcDxlRFOGrr746997h4SE+/PDDpfeDIEAY\nhhgMBiYP0VNTer3e0uvBYIAwDBEEwdJ78/kcwJlrbDgcotlsLt6jwAuuh4zpdIq33noLQRBgMBgs\n/r13796SdTuZTAAAnU4HDx48sHW4nppweHi4dGMfj8cAsFiYE/fHpMvr2rVr+Oqrr7C7u3tuElgF\nv7TXQ8ra2trSct4sGo0Ger0e7t+/j4ODA0NH5qkru7u76Pf7i9dHR0fY3t4GALRaLQyHQ7Tb7aWF\nOEdHR/jVr36Fk5MTHB4eYm9vj+RYvIXrIYP30WtjYwN/+MMfsLW1hddff13zUXmqzmQyWfKzhmHI\njIiZz+doNpuL18fHx4vfjUajRWhjGIaIooh0stdbuB4SptMp16PXZDLBr3/9a9y6dQvdbnfxeOfx\nyNJutxEEATY3NxeTtmtra8zfZK33Slqx6TwgVHgL10PCeDzOjL9N02q1sLW1tZiw8C4FDwXdbhcH\nBweIoqhQLBuNBmazGQDg4cOH+PLLLxEEQWG5pzt37igfpxdcDwnJpbws1tbW0Ol00Ol00O12faiY\nh4QwDBcVRIrCuG7cuLGwhO/fv4+f/OQn6Ha7S2FhaYbDIUlJKC+4Ho+n1CSjXrrdLobD4ZIPt9/v\nYzQaLZIkxU9iYRjiL3/5C/70pz8BOHv6yoojB+iS5fvkNR6PxwNge3sbt2/fxtWrV5fen0wmaLfb\n2N7exr1795T24S1cj8dTe8bjMZ5//vlzYgtg4e+lwAuux+OpPWEY4t133z33/mQyIY1Y0BIW5osD\nejweEVQ9m2srK3jE+d2nn34af/7znxevkwsbwjBEp9PBfD5Ho9FAFEWIogjHx8eYzWYL94Is2uJw\nfw/gypvZnz06uHDuvfeeuim1n5vfvCf1O0pWe6fMzz9//8n/fwngf2V8J6+tWGS1YxrZdi3iYNZj\nft5rPgn3Sl6jdFvFbZPVLqw24Tl3QOz8XehLSeK2evs3wP/+lO98eTDVJnn7Sfed03V2zCwPj3Cm\nOTz8/X/8x+L/w+EQ+/v7uH37Nmaz2WJF2ubmJkaj0SL6JggCnJycKBuTxl0KvAOFl/eeurn4s0GR\n2AL5wnHlzSd/unBNRHivP4XYimKrD+URn+d/bj1Ful3ePuFa39HB5uYmZrMZvvzyS8xmM7zyyisA\ncK46SbfbxRdffJHp4xXBiZVmVB09uR3XOkssIH/zG+DK/1Dfni7R4aHIuk3z3lM3rVwP1wRUhkcH\nF4B/+oZ8u/H10NVGvNZt3dAmuDqtNh5cFd+X/9b2EbjFlTfP3ArXUu/lYfNGY4vvfV/fnEie8Lo0\nZqqEUZcCpe9WhKTbgXJ/PO6ENBSCKyo6rg2erON/4b/+pbpRV8G6jXnp+/qHabKPqPYXb93m44RL\nwTRxh3BNiKrEway3NHFWFZI32apZ2xTjoUo3Oh0YE1xb1q0uZKxbD5si69Y1gcvrA64dpwt46/aM\nWi98MC34FANRdhtU1rzswGG1NXWkhsp1Zf2W9ya72jvN/Ks6ZTagTGHEwuW1bqv4GJo+9+TrOgzC\nPB4dXBA6f2qr0UZfi8+3ahZwkdCybtKnb9QrW5xzPtz44rgsvEVCwTug0t+j2m6d4bGybD/ervZO\nK3EtvUUrjnbBFbFus167LLwxVINHt/V785v3aj9IbIttTJlFV6QPeet2Gecs3CS6xVbWr2lioKTF\nl2qfyXMui/hSLuPtNQ+cEd2yIdNfWO194RePaie6WvLhrqys4HFiEkTGynVVbKsM74DiFaz0Ncxq\ncyo3iozVxepjef2D4smjbJYt1Y05r9+crq8pJ69ZWVnhz6UA9WQ5sjgnuCZcCLbE1tXVb3mwBppI\n8pok6fNmCZgOsRWlDNdJBzrbNCt5jRdclY2mBBdwKw7X9CDiOU9XB7as6PIKLnBedEUsQBN9yNVr\nw0LFDWWiTZN9p06CW7s4XJODR2Qpse2sZ3mYFhvXxNbkfqj3K+P+MHWuveZBKSbEqbEquKYHs2mx\nVfmtqwKcJG/AyA4kF8XW9P7S111mv2WK766b6BpzKQB63Aoi4WN5ibApJzGq6k8UcS0UXQtdyVFM\noOMaULucZF00ttr15yvv1MalYD0sTCU2VCS8h9VhKcKuTGU9A+wIL+s6mQy1sm3xs/Yvcl10nYes\ndWu7XeuCUcGljCdVEdusTil7bDY6qk7hZVn+thdOuC4KOo/PVhJ3Dy1OTJqJdqQsseUVYJYFIGod\n2BYA3fl9s9oj71rFbgQdPjnXfdmuU9Svq+aicRmjPlwg359EEXSfNdhF4j6TFFm7LgqAbOeVWXyg\nev68x+piO9ukyDWWB6s/m2zjrONf+/ZpbXy42gRXtGovoJZ1KIa1uimrYqxMORcdHVTm3PJQmWDJ\ng1p0eY7Ri+15qAXXhTjmOgmuNpeCqNiWibo9Btm6ZnVrZx5Y7pW862RTbOP9+JvnGc64FAD9Fi4g\n7qe1+Sgmk1tCRqRk8xnInr/MMfoBmw0r54NtoWWRPO46WbhGoxTKaN2yOq7uWXvRCaiyLOyQQXdZ\nb1cpijPPi1hxaSl9FsvH8o6144jp9Xo4OMhu4/F4jOl0itlshm63q7QfJ6IUYooEw9aqFJYF6Mpj\nr8pxmLwRqraXK+1tguQTzsGst/jLgvXY7h/p2RweHmIwGOR+fuvWLezs7GA+n2MymSjty2oRSRuI\nlnaJsWnpFqFThFy5bklMWLtFbWrSnZT3WZYBkrR4vcjysbu7i36/n/lZv9/HCy+8AADY29tT3pe1\nmmamoAwYd1F0qc5N9kZEgewya5vWrk7R540pZ/n4vdjSMBqNAACTyQTD4VBZdJ1yKZSBIveCSRFw\nIa9vkiKh4DneslW4pb4GskukWe4GjxqXLl1Cu90GAKbrgQftguviI6nqMRUJggkh1LEPVoVhXkQH\nfV5blkl4qW60FIJZNuF1/VjX19exsbEBAGg0Gnjw4IHS9rS6FGQGbFl8T0W5F3SeR50mjlwqLV50\nLCpuhtzyMxk1vy784pHQ9lxNgWhKbB8AGAn+Zj6fo9Fo4Pr16wv/7nw+x4svvqh0LNazhYniUhFA\nG6KrW2x1+XJl6pllfdeG8OblmKASXhGxzXufJcKuVcDWNX7zFltdAfA/E69/9f7y5/1+H6PRCHfv\n3sVrr70GANjc3MRoNMLGxgYajQYGgwFms5myD1fbwoeTr+UHhsoCCJXSLrKYyrvgimVb5MPlyWkB\nqLW/buEVOTYdi3koqtnmibAt4dVdRDJvsdW5775fwaW9dUK3T9f0ZBwLV9w9Ov27otumjtOmKh1+\n+sZq5raKYnqpKZtfWSdOuhRM+XHzLBMq6wZwR6B0kw7SzwpTorhpmHApJPehUsZd1rq98ItHWi1c\nwKyV65Ib0DZOuhQAMyW6dVI1N0KMSEiYqmvBhYkyUSjSjMaIiC7PRBrAFtr0taE2FrxLwVELt8xQ\ndlLXxFYUXkvXhrAmr5OrNdbyLF1ecU0iIrTJ9ynPzVu6JbVwAbmJM8BeEUAZXBNcyqcO1xZtsMg6\nVpVrblJ0qMYCZR9Pn7+3cD1cVLVCrw7S1q7pgpiqZevLhqzQ5oW8UVq7dbZ0uaIUer1e5usgCOiP\nqATozr5UNbGNyRpkusXMtUxZOoWm1zxY/GXBinZJr+zLq2dX1b5pikILN05dlswVGQQBPvroI3zw\nwQdaD04HsrPlLg3aMpMVq0vpT83aJhUurt7iOQ7ZsjwsaxdQa+O6WrmFgpuVuiwIAuzs7Gg7KJew\nlWi76mRNqAHZ7a3b1ygz8GXmEKgEhqL6h+gqv6yJTVU3Qx1Fl2vSbHt7G/fu3Vu8DoIArVYL4/E4\nc6kbxaQZID9Jo3rXr2puWxVUJsxYyFqKqiFMZRvoIu0kK7Sfp5a8ppfK6io39fOVd/ykGYu4zMT9\n+/cRhiE6nc657+ieFNFxd+TpNDrFheVfS1LG+NQ88izdImQHuOz1owrPEoFKZAExoU2+nxRdVg6J\nsiSdso2w4AZBgGaziZ2dHayvryOKokzBjTE9G81L0perW2iT2xBdlJE1UGQTdruKCd+o6PXjWXSQ\n/A6F+Mqcuw6hTX8nbenqcjHUAWHBbbVauHbtGgDg+PgYW1tbmd/75O1PF/9/5uW/w3svmxFekU5r\nSmg9/NjObCW7pPb0jVUh0aU4P5WxxCO2ye/yii4Pf/j4j/jjx/8q9duyUyi46dRlnU5nkfX80qVL\nuHr1aubvvv/2S+feoxbbtFtB5hHM35HdQpfQJrfLuokmRZNXfHmEVsd5qeSnuPImv+hmpT1UebL6\nzsvP4DsvP7N4/S8/+5Tx7WqhbaXZPz7+h6X3dIZiifoBRSdcdM8uq8wmm3IpyJbV4cW1lIFp0uJr\nS2TzoJooy0JUcEWNGD9pRoxL/tusYzHhe3IlbtNFbLZNet+sDF6y2zQBa66ElVS+yNLNS+rtkUO7\n4KqILa8Q6u7gKhEROo+tChNmrt2IeF0PrN/ZJM/NIFPJI09sq9DvbFGrXAos8S8SdxnRVY0H9tAj\n4k4qEl9XRDZNnrWbJ7oi/lyPGloF1yUx4TkWKtcC1UB0vWqtqv9Wl2CJ9Lvkd3nFtyzkpcLkEV0V\nV4IPD8tHm+CWTWx5KLJyyzgoqeENj1JpK119i1d8VbYrA7V4qVi63p2gRuVdCqKdXfburCNovWzw\nhlGZzlcsg4z46rwRyIquqE/XT5LppdKCq2MAqMT+VgmWpS+6CMCUyLJcNEXl7mOoarOJokN0RfHW\nrTrOV+2V7Siq0REsivKOmsC1zi9r3ZoQr3SuV9Z3dFdgVkFl31li7VofqgPOC64MFINC98By5dHZ\nJKZvUDwCyvqdi5OWvP2GN9m7LtGtY//moRSCK3rxKCYZfEUH/fh20EMstrZFtywMBgOEYZhbwabo\ncxFK48MVzX2gkqXMlaKQWRMbtivcymLD/RK3lailWhUB4i1Tr4syhIeNx2MAQKfTQRRFmEwmaLfb\ni88nkwlardbivfTnopTCwk0iY+3y1rXSUf8qrgMl29GTg99lIdCdG1aFRwcXuNqO93uuI3pz470h\nuehiUeXo6AgXL14EcJYJcTgcnvtOXMMxiiIlsQVKZOEmkc30lWf16hBZSsomAjzia8OdkGXxlq1t\ngXJYjmVhPp+j2WwuXh8fHy993m63sbGxgWaz6bZLIetuSN25ZTueKy6DMlBWdwKLMoqsKLraXDQP\nbhluDqzMYfP5HJcvX0YQBOh2u3juueewsbEhvS+jFi7PI4noYLCd17aKIiuLy24FT/XJ045PP/kr\nPvskKarfLP7XaDQwm80AAA8fPsT6+vrSb4MgwOuvv47V1VU0Gg30+/3MOo68OOdSYNVNcoUyiGz6\nBqQ7+bsIZWi/KmKr3W1buS99/1t46ftPXt/6pyeCe+PGDYxGI3Q6HUyn00UFm/l8jkajAQBYXT2L\nMY8n1lRwTnBjXH10ke20Jv2GWe1gq7Yc9aOtCVeVh41KeR3XaLfbGI1GCMMQjUZjUcFmc3MTo9EI\ne3t7uHPnDlqtFmaz2aKArizOCi5QjQvLKgKp49yKbjoiwmvb9yYye+5yP0m2I9UNz7bVKIrLxxuL\naLIY7mg0WvxfxYWQxmnBBfQOJtXSPCx4xILy3KgiNigRaVuTIUe63S15+4nfc8mlItvurt/kXEVb\nTbPfZ7yvmomI5wLLVOKlShSuIhoqnZfCcsg6P9Vab7yB91RiS9E/KMSQ53qY2k/RPk332axjpqpp\ndvI13/Gsffu02jXNYtK5NkUFmOKumiUSRaW5iwYHhWDInBvlI5qMxcuaOLOxyonVhqIrFAG9S8pt\n+dQ9dnHepWAbE2Kb3JbtiULR7aXL0NheTkoJryiqXAPbLgaZWmfx70Rx1YdrEmuCK+NeoCjNnGeV\nycym6/A7ltk3ZiP1Ig8qA93EOciIrsqCH96yO+nvyOKF9gmlsXApxDaGN4k4axC4sq7cxdlfiiQ9\nMpT1RgWIia7q9c6rdQbQRdC41iddwUryGhfKeMQi66LYim7fFUsSkD+WMoslDzyLREyKVF6iJtUE\nPjoSQFWJUli4lNZtEhfFNrkf2/5cEVwSfSqoz4knDLHI0tXls5c9Vy+uYhgXXFHrVpfYyuCKGyEP\nW6JLJUxUrgXXSFq3KqKr89qaFvq6UgoLN4uyZfxKhsTx3nTKMIGWbC8K/1/VRDcvDFF0ktaE4HlR\n1Y9RHy6VdWtDbGVF4PP3z8cfi+CyPzdLbOP/q9QEc+EmQ9GOLL9tkU+3DuJ3MOtJJ0AqK8YsXKqJ\nsqyOWLRwgQddYpv3vkh7mPbn8lhgeWKbRtbqrZqlm1U6vqjfVkF06yaoRTjrUsgaoCyxjf8vI7rU\nYqti0VIhK7pFA4QltCy3SdmrLIiS1Y5Zoht/17Uk7Xl4AVXDSC4FCldCkdimoUpKIyq4ImJLOYGY\nh4joFuWYyBPbonPOO0+e8xFpf95+A7BFrqg/pPeT3ke6HU/fWF16nZeo3Zbo2hbR0/W12uRSKF0R\nyZiiTsLbiWyJrQw6H7FFLFtR8tpFRwWQNHnHbSokMC22ee8B+oUv9pmm/zzm0Cq4V97UZ93yoNqZ\nRAd7fL7JP9b3dB+PSLu5mkoR4AvGZ33Oe7PI+156ApDi/PNEVydlcVtUGW0+XBdWk/GgO3lIWdoB\nWF7ynBycrixqSIqqqE+YVfvOhYUuMbr9uelkQzq2K8rPyY7CfZydNGMh0lHKNCGRhewjtexTgY62\n0nHTkW2X9GQiy6rViQsFN+NrLSO8ZR5TNiml4JrERniSS4lDXLFuKaG2aFkRIRd+8UjYfWDaSCgS\nXi+udDgvuFWIReTBVqkdUaoUG5umyufGgxdW/VgRXOo4zLz4Rl5sJYGuQzxqGZARWtZiFJUS8p5q\no01wXRIT1Uc0SreCznahtm6r6E5I46pVW7a5B5W+4ifNHCbPclC1cnVi4uZTF9cLFdRCazs9pi3q\ncFOmpHSCm8ZUPKOMlWvKytcx0KkqvLoYFmfaqpWZOHOZqonsYDBAo9FAFEXodrvCn4tQ2pVmopj0\nqalmzfcUE1cWELnZqGQwy9ueTlzzA9/85r3Kie14PAYAdDodAMBkMhH6XBSnBVfUclOxIiisRBtC\na8q61YFsboj0OfO0gWmrtkz+1yKqKLQxR0dHuHjxIgCg1WphOBwKfS6KNpeCaumO+LcigsLy4VIM\ngCy3gk1LtsxiK0rRudouN+4aeW3B22fq0pbz+RzNZnPx+vj4WOhzUbT7cFUHQlp00yE3sX+MQmx5\nbhJUGbtUO7QNsU1Xdo2J/bS8VS142lC2DLhJoWCdB8sdoHtyl9UGdRFSEYoyh1FmFivlpFmW6OZ9\nTwaepZ8i26Deh23Ltkh4i36XB8V55Ymu7hWDOq6JaP+tu5jmXYM/fPxH/PHjf0288+nif41GA7PZ\nDADw8OFDrK+vL/226HNRSiG4Wa6FouByKh8ar2tEZcClf1uWYn55wpv1nTx0nI8J0dXtSvJiS8d3\nXn4G33n5mcXrf/nZE8G9ceMGRqMROp0OptMptra2AJy5EhqNRu7nshgRXIpHPV7R1TVZkSW8usQv\nz/p1tXCmjG9btQRQDKs8TZ7oArSTaFTzDDLIJN/xETRPaLfbGI1GCMMQjUYDV69eBQBsbm5iNBrl\nfi6LtooP//j4H5beo7oDs7LrV2lm2BQmLSNeYRINh2Jdd9nz4xUs3kokcQQNj+Dy9mOqWOksTIoy\nRQWGLM3J4+cr71ir+MBl4fZ6PRwcPOkEMoHAVBMaeZNoXmzdpUhoVeNNWddftt/ZsgJ5+rGJlJIU\nJe895ymMwz08PMRgMFi8pg4EpsCLrRw6rduihQnUJV5Y2zHh+jGF6fy91ItF6k6hhbu7u4t+v794\nfXR0hO3tbQBPAoHb7TbXznRZuaapgkWtS2x1W7NF2y4qOW5zcqnInSDjGskSQ57aeqJLrr3FS4Pw\npBl1ILAspkU3LRRZpWjqDOta8IgsVfKhopsh1U2fZbnLYFJseb6XJ8istJSeYqSiFGw5nNPoFl2e\nwVPGSTtKK09FeNJLsfOWZosKsSnRpcK02PIQbytLeL21K4+w4KoGArvW2dOoPPKWwerVFS0SIyO0\nMt8vEmEe0WVBHQstY73zRiFQCm3WtlnWLuCFVwRhweUNBP7k7SfBxc+8/HdLgceuISKyycGfN4jK\naPXyouo6UMVUzmPVIpwUbZFlnGTFPF95U6/oUvPpJ3/FZ5+48ZRsmsI43H6/j93dXdy+fRuvvfYa\nACAIArRardywsJWVFVw4PgGgJ0Yya2UWVSB9ETyCm8RFwVWxcnUIrowFy6KozXVb+UlE2iTvuG1Z\nuqpLtXmpUxyutoUPseAC6qLL6tgUK7F4B4WMMLgouIC+HBG2c7jquMEXURXh5Ylc0OE+qJPgGlna\nKxOuIyOetsPFPPYwZdUWbTuv/4kUlsybC5BxMVBV3PB+WhqMWLhJKC0+1ZyfMaoTPXWxciluZjos\nYBtWbRGUYWNZ56dzWW8SE0LrLVyNVGHRQFWwETFi8trzhFTpEpR436JZ7rLIGjO81q4s3qLVg3EL\nN4nK4KNOl8gaBDxhTLqrTcjAe3PjFd0yuWtELUDTlZUpF0i4eK48xO1BYXGWxcK1WtNM5wSLacvN\ntaqscdtStXFZxDav/laR5WciX0DyuGRvwlnXM+vaxPX1sv5coCz9iRrrRSQpE5ikoRDdTCH958+F\ntmF65j5vGXIeNjt/LJDpP9ltZcErpiYStYie2+kbq+f6YNaYEa1gbJOyHKcOrAtujIgoiXRabZZu\nhui6ZuUmcUl0eYRVRIRZVm1e2FT8l4Wp7FhFVm6yP+UJbxrV66i7H9RZbAHHSuy4NqFGYd2axnYc\nbB5UWeJ44Y1PzQufolq2mjVBlwxfzJtEy7t5p5P8ZIWQUYquibwbdcLqpFkeJlanpeHu9GnB/fGV\nc1/JmkAzcSPhEVuZeFXZgcJzrXSUnadOWSh6TDLVIpLXbqnfxf2No5/p7mO6Vij6STOH0Tn4C3Hc\nui0acDwDMt2+pq0SXY/zRQsAWJ+LHpOJiSlTOSWSqPiJXU5YZRKnBLfXPHDKpVBG8tpPpF2LqjWI\nbMcGeYKXJ6pUq7HSx5A+DlYttOT1WRLTH185Z93aeoJSxYuuQR8uZYeQDdhnLf2VCUg3gUyuiPS5\nlGEwJqFIcp23CCCdWYtXbGWPidra1W3Z6p5Hqfvye20Wbmytuma1lukumz5WkZCpuM1davs0uqMB\nisROtsyMKjyCkyWseWJLdY1NGRxlGoPUOOVSEEHlLlmG5axFYVC69ls1skT3ypt63AiqpK9XUmBN\niq2M8IrGMNdVdEsruFWGpzNWvcNSWr+Uj/Wqx0VZLcJFy1ZUdMvQjweDAcIwRBAEUp8nKbXglsnK\n5cWJRR0G8KW3+bAltioiXKVrOx6PAQCdTgcAMJlMlj6fTCZotVrodDpotVrnPk9jXXBVl3ZSiq7t\nR3DZiUBq4Y2XjtqeRKyqlcuCpw+6ZtlmtYeJZdImODo6wsWLFwEArVYLw+Hw3Hd6vbN2jKII7Xab\nuT1rgiu7XDOLKsx6qoomlehm5WFwQXzzELn2tkXXpX4qWpFalrKL7nw+R7PZXLw+Pj5e+rzdbmNj\nYwPNZnPpe3kYF1xRi4zX6tWRvEMmBIeqOqsMKtvhEVWXhZcXW9myRPomy4LV/RSmIx9I2UWXtSpt\nPp/j8uXLCIIA3W4X0+mUuS0jcbjUgsLqvKLrwHXHBZou/SJzPjIJsWNkBEBkAFLE5KahStTNe2xU\n/Uu3KyEW22S+hqK4XJFMbLpvdnnn9dfPPsXj//tZ7u+yJruazSZ2dnbQaDQwm80AAA8fPsT6+vq5\n377++utYXV1Fo9FAv9/H3t5e7r60C66OiR0e4Y0/F90/cwHEj6+QLe81UWOLurhh0e9t+8BFSA5+\nnRaYSrl1HW1bJLZFqPRbE6Kbxbe+9xLwvZcWr7+5fWvp86zK4zE3btzAaDRCp9PBdDrF1tYWgDPL\nttFoAABWV8/artPpIIoi5rFoS16TTCShU1yyOrTM/tLbyU0mApyJbsFyS5HKq7oxIbwx1PW3igao\njjJBIgLMOj6XfLZpWMKb7MvJ66laRy2vrahqmvEmzDpdXxPaXxAEaLVaiKJoIc7Xrl3DaDQCANy5\ncwetVguz2Ywp3oAhwQX0iy7V9lllUPKsAB6xlT0+ymxapkqcU4muC9UJ8o7RdaHlcS0VXe/0dVRJ\n8M5qL9eXzC70AAAUa0lEQVQFlxJjgguUJ25URHR1iK1qBy7ChPCqiK4LQptH0WOxK2KbRPZ6F1m3\nMar107zgqm60IDdlGYS3SHRtCW0WMgKlW3RFXComquiawEWxTSJyzXnFFpB7CkjiBVd1oxzJgMsm\nugCfJRCjWp1ABteEV1R0yyq2rgttGtbxpqMSeLer4vrygqu6Uc7s61UQXV4/VxrbmbKSuCK6ZaRs\nYpuk6NhljQbRG6cXXNWNCpS7KMPAyxJdWavWdBA4lfB60T1PmcU2CVWkjyxecFU3KiC4QDkGHmtw\nuSi0aSiEV1Z0ddSos0lVhDZJ+py84OrBCcGNKcPgk+mYtsU2Ca/wetE9jwtCC/D3OVmfvunr4QVX\ndaOSggvoHXxUIUjxwJOdveUlq9IsVeJsFeGtm+i6IrSA/A3e5QlJL7iqG1UQXIBukQAvOjqj7LHw\nlPOOURVfG9YuzxJVF8TXJZGNoXiaclF4veCqblRRcAFzMaxpVDukynGIiG2MCdGlmkyjLjWkg6oK\nbRLXRNcLrupGCQQ3xkY8KyDeKan9tEXiS+FaoMp0JbpEVBbd4ltGsZXpd15wveAWYjO+taiD6pwU\nkynpXYSqKyGLLNHVlT1Ml/C6Jri65gi84HrB5cZ2VICptH660CG2SYryp6ZJXk+ZfVKKr0uCK5u7\ngAcvuPYE10gCckp4IwR0UUaRjdEttgC/VZt1/WSSp9sIY7JFmfue54zSWbhpVFLG1QkTYluEzhp1\nulZcpdEREpdGNe8sDy5Zud7CLRF5Fg5VGZWyUyahVcF1S5fX1ZI+B9E+nJ5spfL7e2iwXiadgjyx\ncOkubpJHBxcWfzzoElvZEu4qyVh03jhMF9AUEdvP35cLK8zCJV921aiE4AJedAEIiWwM9eDirbLM\nsx1ZZM7JhMgUCbbMORcJrawIe9HVQ2UEFzAjupSWBAWi1mwSHWLrCtTnptu6FXUl6OqHyXbT/cRQ\nRyoluICeO3PcuZMdPOs9U6iILFA8kA5mPSGBobBo87ZbBOs4yyoWLLGV6XNZ3xdxV5S1HV3EuRI7\nVBdXdabXxjJbFlRWuugy3aKJHtWy2UnyzpEnl4NqYhxduYCzKEpan9VPVW/sWX0zq70pk5Lz4qMU\nCJC9MCKxmKwBlzVrXabIBR2+ZxlREV3IoEJeSkGZ+NwkLkUwyFYIufKm/qcpnjZ2qS3LSKnDwnrN\nA22+tbRFkNfZqa1a00ILqFlwNgdgfP15bgis4zT1yCwrtjG8fbLodzKWbd3p9Xo4OMjuZ+PxGNPp\nFLPZDN1ul7kdJ324Lt5Br7z55C/5mgJVnywLVV8tr5hRkfcEktcnTFnfqoiILW8/SPbJvL7oxVad\nw8NDDAaD3M9v3bqFnZ0dzOdzTCYT5rZKbeHawnQSGRlM+iRNocO1YEJsZCxbGfdXUb/0YivH7u4u\n+v1+5mf9fh8vvPACAGBvb69wW05auICbVi4VuqzZGBv5a01YuYB6vzApML3mgbIbgQovtnoYjUY4\nPj7GZDLBnTt3Cr/vrOBWiaTLQLfQsmb0dVu2IgNYZfKyDDdj2arOSaj6iojYlvXpxyaXLl1Cu90G\nAKbrAZB0KcQO5CAICp3EKqg+QtrA1sq2qrkQioogqvSN2LWgo29Rl4TPagOV6g48/UQmMiU+pjKu\n7AyC4Nx7zWYTOzs7hb9dX1/HxsYGAKDRaODBgwfM30kJbhAE+Oijj/DBBx/I/LwSuNKxdAmt7GRU\nWcKGRMSWJxqGWmhZqCx4yYPqpixbLZiC0zdWsz/4t4+Bf/8493cyRuN8Pkej0cD169cX/t35fI4X\nX3yR+Tspl0IQBPjiiy/wyiuvyPxcCNcGr263AC88q8VsQWU5FllyrvQNk2Iri2hfEek/6evkXKz7\n37wM/Pe3n/wJ0u/3MRqNcPfu3cV7m5ubAICNjQ00Gg0MBgPMZjO8+uqrzG1JWbiz2QxhGGI8HnPN\nzFUFF4QW0BtXWzZiUbPheiqD0AJ23E02LV1qrl+/juvXry+9NxqNFv+PLWQeF4SU4MY7uH//PsIw\nRKfTOfcdynpLtn25rnQcU0JLEdtK5VrgHbiqpXpEERHb5Fgw3ZdUxDbLlyta1NWVseMKwoIbBMHC\noby+vo4oijIFtwq41FnKJLYxtvy5umNtZcU2+VpX35JZFq9CkWHlRXcZYcFttVq4du0aAOD4+Bhb\nW1uZ33v7N0/+//Lfnv3pJOtuLDvgXewgRY/OKsucdazWKmp3k3ktbCx2iKE8T8rjjq85q8/IJNnJ\nIi26n37yV3z2iZ3kMbaRyhYWx5pNp1PcvHm+E6ysrODkazrREokZ5I1/5M1U5Soqj4q6l8NSlbTX\nfU1ky8CL+m5V+pqpm0O6z1CJbZK886bKFoYfc27jn1eqVyadSnBlhcWVFT46EXEziIqsqfbKG7im\nboAUmemSFImui2KbhMdvSzk/A3jBVd8okeCqzq5SrPYpAzpzCOsma/C6NLGUJO5zsoLruthmQSW2\nMVlt4AVXdaMEgksZylIHaxewV15cFZuz+DEiJdJ1JmXnPRYTUIttTPoae8FV3WhGxQeRTqgjbrAu\n1i5gJ2u/Cq4sC6Vy0ZgugKkDXWIbk7zWdRJcY+kZWR3JRAxlXhSDDK4JVpqsiAaqiSweRIXTtNDm\nCbzteO8yi22cDF1nmakqYKWmmSwuxZqaFF3dFqCu8CzbFmsaHn8x79OVDreCbcGVrQOYrjzBK7px\n23sL10EolyCarNMlS1ZHpw5l0x0H64qrgHWe6cksW1auTbHlDWejpo6LIrRZuHEFTSph07XmW/b4\ndFi4Kh2ct+PaTCxiy3XAA6+lyztxBvD3ERtiW3RsvG1HUetv5X14C5cKlVjQrG1QQ3F8qlCIIMv6\nVR08vBQNMlMWr0x72rJ0TYstz02Aor98/r7352Zh1KXggrixELFeKNBpbcYCYkpsk9vgEd4y+KN1\ni65JsTU5aerJR7tLIYmKkJlOOagr7CcPXTGOItvWbeHG6LRwqd0yOmObTQmubrH1k2b8aBVcHZai\nKeE1LbhpeDq/qHDJCi/lo6EpP66oeLCOS0QYqeLNKaF0I6jiV5qVKCwsRrfo6l5JZBsbj402ZqNV\nJs3S6FhMYkJwbYstz3Wvk+CWJiwsCU9qOV2UXWyBJ4PAplVjar9UTwosf66rfcKG2NYtzEsU64Kr\nkqtUJQds1YnbkTXoTOSktT0AdZ6jq0ILmD8229e5LBgTXN4OIFqjyqa16xJ57VUkvLqsXZcGIEt0\nRY4zaeVWQWx1TtRmkT+m3yE5jjLgbLYwQNzHpXvpr0uDTFceCBu5FGQRFT8XUkGm0eHHNelKoIjs\n+PnKO96H6wK2k4m4BkVbFJUdsi1ARWS1AW8ppbQ17/q5ymBKbEVyUHie8C3bB1DEzW/e47ZgKMLQ\nXFuQ8d5TNxd/ouRZ/LLbswVPG4icT5mEVuSpzdQTmBdbeZy2cJPwWru6JtJMuxNkO3H63Fmr51z2\nR8qcv0jR0DKJLiXUFRtoqkm77cMNggAA8NVXX+HWrVu537tz5w729vaY2zImuBQxkVV3MVCJbN53\n8qx3W+XMs46DYhsunItpdLsSioS2qpPWYRhic3MTGxsb+OEPf4gwDNHpdM59bzgc4v79+/YFlyKR\nSBIe0ZWNXLA1WUbpLuD5jWvWLvVNtK6iq4u6ii0ARFGEKIrQ7XbRarUQRVGm4K6srHBtz0iUAnWV\nT4A+gqEsYmsiCbtJsdL1xFIGwaU+d10RKDyRCCr98nR9rRRRCtvb27h9+zauXr269P5kMkG73cb2\n9jbu3bvH3IbTk2asDiI6oHrNg1yRERVbG0tjD2a9SlsSdcOGa4zSb50eG65NNlMzHo/x/PPPnxNb\nAJjNZtzbMRaHqzO5tsojuUhxSeqKC4BYtVhqbFu5OkXHZQvX9nnLjEXeeFuZvqrdwv23j4F///jJ\n6//3s6X9xZNiSZrNJnZ2dhav8ybEYusWAJeFa3Thg2uim4WoVasiurzHLNWJ31gFAFz4xaPc79gU\nXNuiYwsXEtZQCS5AYzCQCS5+z/ntvxfa3+HhIXZ3dwFgMWk2n8/RaDQwGAwAAMfHxzg8PEQQBAsB\nzsJpl0KSok4iEq/L2obovnW7F0TF9vSN1YXYqmy7ytEgKqhcb1NtWrQfna6FLGJ3XvqvDAyHQ+zv\n7+Py5ctoNpuLybHNzU0AwM7ODnZ2ds6MzJOTwskzrekZZauAshDtLLydXNVfK9OJaWIYz2CJrKyV\nC+ixFnULjy4LV2WFmmt1y2zlC86CbGmvJguXEq0WbtaFMF1pNrZ8WZ1PVGyzKiOIHheV2IpYtDL7\n0RGyVUaS17csZWhYbW3ayvWcod2loGOAyXb4LPGVFVsK0VVBRGhVBLlsuFj+qKw3GRG86PKhTXCT\nllO6w1HcXVUHAcvq5bVsZUVXxbpVtWhF9wfQCQZrO66GvJVdbKmsXJ528KJbjFYL13XRzdqeiBsh\n732V48oTHlWhVRVpVeHgEVtXRTcLihurKUy5FgA+F16d0e5SYIkuBVSiy9pOUTVbKtFlia1ueMRO\nx/XLS7bjAmWISODF1qSdF95ljISF5YluXTM2uQhPmI7s4DGxfNmlwe2a2MbkHRfPOPRjlQYr6RmT\nyUVkak7puPgqta8oy4i7iIqQmfDb6hJaEzXfTJOX2Cer1JIXWXqMLXxgDa6iC/vo4MLSny7ytl11\nQdW14qzMYltEmcWIlczdxDirM0ZXmvH6c00JrAh5oltVMVZ9RK+K2FIttHEx+VDZKn9UAW2CmzfR\nw/LnuiCwVPu3fR5p8labiSTv4aWM4V8sXLuW1HjhNYf1XAplutBpa1bWunVlcsfDT1J0RRO5JG8y\nLt9w6pyb2BRWBDfd6VwTXZZFE4tsFVwJ3roVg8rSdbkdvLWrF2sWrsudDuATXWp0ZVBiJa9JUiax\nLftkmev93wuvHqy6FNL+XJMXWOe+yuDzoxb3Klu2PPC4E8qIF15arPtws9wLui9yvG2TeUNdgtqV\n4AekHC6Kcd4xeeGlwcrCB16SF5jiEdKlDmOq5DuPO0Gn2OoQFYq+wKoMbRJW+XrTx5H+v2vVnauA\nNguX12/Ii8odlvXbMli5lG3JM7h5V1e5dAMTIT4/VrIiEVTdCbYs3Tg2OG//rM9Ex6MX6DO0VXy4\ncHyS+7nqHV3k4qlWfIgxXfnBZOFI0cFAIbR552erXL0qsqJrw7qV6Vs669+tffu0NhUfjLgUXHhk\nKiJvjXkMa129qhWc5V6I24xSePOuQ9G5p7+r81jSuC60KpRhXHho0VrTTAe6rTERS1eHu0GXtaua\nDcy0+6CMQiti5doUW0oLl+I6eQu3xhRNCuj26eqwdnkHd5al64WWn7yJ0F7zYOnayVj3ZfWXe5Yp\nleCaHIwij9nUxPvNE940LCEWtaSS521ykJdFaFWiG0SuRbo9qKJaXAxFqxPaXAonX18gH7Ayg1L1\nGGwLgezxs0J7dGG7rXTD607S2edUt039lKR6zVd7p1h5H7VxKUiFhQ0GA4RhiCAImN+jHICsgo+y\noT08KfNsP8rJpknsNQ+82DKwfV3zKGrHsrUzi6old+dBWHDH4zEAoNPpAAAmkwnz+xSlT3h+L3rx\nRDI4Ua6y+fSTv0r9zuXifBTHJNsuoiSvpcg1TfcvVn+TaQ/WtU23jWx781i3JurnUcU/m6Lf7yMM\nQ/z0pz/N/DwIAgRBgP39/cJtCQvu0dERLl68CABotVoYDoeZ30s3qA6xoLxopqzdzz5Rf5RxSXxV\njiEeeKu9Uzz4P19rHYh5N00XUhIWfTerz5i+9hS+37IJLQCEYYgwDNHpdBBFEX73u9+d+3xzcxPd\nbhdRFCEMQ+b2hAV3Pp+j2WwuXh8fH+d+N6uBRYWC5UrgeS8LVucxae1SoEN8dUysJMW1aOBRDkye\n6yV7PSmOUeW6uXDD5aGMQhvT6XTwy1/+EgAwm81w9erVpc+jKFoYna1WC1EUMbcnFaUg6nCOGzs5\n0WAql4AMPOvb08euq/OLxP2qhhGl3SxFbcA6Z4oBphrzLNIGrKgU2XNh9XHe/rLaOwX+21NS+xCF\n2p3A026fv0+6Sy2cnJzg8PAQb7311rnPut3u4v/j8Rg/+tGP2Bt7LEiv13s8HA4fP378+PGHH374\n+Pbt2+e+A8D/+T//5/+4/1QR2dfTTz8ttY+tra3HURRlfvbb3/728f7+fuE2hC3cGzduYDQaodPp\nYDqdYmtr69x3HlsKufCUi16vh4ODJ1b0YDBAo9FAFEVLlkPdSLdL/DoIglq3CwsVzcmKtmo2m9jZ\n2cF4PMbKygra7Taee+459Pt97O3tnft+GIZ49913C/cl7MNtt9uLHTQajSWfBm+4WB3p9c4e133b\nnHF4eIjBYLB4LRr9UlXS7QKc9Znvfve7ePbZZy0dVbXpdrvn/nZ2dgCc6dxsNgNwNn8VX4P5fL74\n/eHh4UKEySfN4gPsdDrn/BeAHzB5+EGzzO7uLlqt1uI1b/RL1Um3C3DWd7744gu88sorlo6qvuzu\n7iKKIgRBgIsXL+LVV18FAGxubgIAhsMh9vf3cfnyZTSbzf9agJEP2dLeo6MjbG9vA3gyYGJr2HM2\naOK7puc8ItEvdWM2myEMQ4zH48zHWY8+1tbWMt04o9EIwJnwxhYwD2QJyP2AYRMPmjt37tg+FGfx\nvv9s4ifK4+PjwkdWj9uQVnzwAyYfP2jYNBqNhaXw8OFDrK+vWz4iNwiCYOHTXV9fL4zz9LgNmeD6\nAZOPHzTF3LhxY9EuedEvdaTVai38hcfHx3jhhRcsH5FHBTLB9QMmHz9oztPv9zEajXD37l0A7OiX\nOpFul06ng+FwiMFggEuXLtW2XaoCaXrGIAgWy9t8vOAysYU7nU5x86abK+w8Ho9etOTD9Xg8Hs95\ntJVJ93g8Hs8yXnA9Ho/HEF5wPR6PxxBecD0ej8cQXnA9Ho/HEF5wPR6PxxBecD0ej8cQ/x9hJgWZ\nJOgGDwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x6454d90>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Doesn't look good. The influence coefficients should die off away from the origin.\n"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}