Estimators/Force.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()\n",
    "import numpy as np\n",
    "from collections import OrderedDict\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance of the bare force diverges in the Monte Carlo integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r = Symbol('r',positive=True)\n",
    "V = 1/r\n",
    "F = - diff(V, r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABQAAAAqBAMAAABB12bjAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzRAiu5mrdu/dZkSJ\nVDL+6OvUAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAhklEQVQYGWNgAAEhMAkiTD7BmCphcCYDOz2Z\nrBk/s2COoAX9Hw7wmn7UVQEqz2PAJgBlchrwfYYyOTYw/4YyGRh4EQHEWQAXPQpncQUAmdzHtI8X\nMOxh2MfAwMMiy9LAJ3vXiYHhLWcC9wOO//8/MDAo6G+A67gPZzEEIpg9cCaS5dwTYKIAWw4uW6ji\nMHsAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{1}{r^{2}}$$"
      ],
      "text/plain": [
       "1 \n",
       "──\n",
       " 2\n",
       "r "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "epos = np.array([[0.0, 1.0, 0.0],\n",
    "                 [0.2, 0.3, 0.0]])\n",
    "npos = np.array([[0.0, 0.0, .0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def compute_bare_force(npos, epos, F, r):\n",
    "    forces = np.zeros_like(npos)\n",
    "    for ion_idx,ion_pos in enumerate(npos):   \n",
    "        for elec in epos:\n",
    "            dr = elec - ion_pos\n",
    "            dr2 = np.dot(dr,dr)\n",
    "            dr_norm = np.sqrt(dr2)\n",
    "            #print dr_norm\n",
    "            dr_hat = dr/dr_norm\n",
    "            force_component = float(F.subs(r, dr_norm))\n",
    "            #print dr_norm, 1.0/dr_norm, force_component\n",
    "            forces[ion_idx] += force_component*dr_hat\n",
    "            #print ''\n",
    "    return forces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4.26692459  7.40038688  0.        ]]\n"
     ]
    }
   ],
   "source": [
    "forces = compute_bare_force(npos, epos, F, r)\n",
    "print forces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "One solution is to smooth the contribution inside some cutoff radius $R_c$.\n",
    "See Chiesa, Ceperley, and Zhang (PRL 94, 036404 (2005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M = Symbol('M', integer=True)\n",
    "k = Symbol('k', integer=True,positive=True)\n",
    "a = IndexedBase('a',(M,))\n",
    "Rc = Symbol('R_c', positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to fit a constant function with a polynomial that is missing the constant term (only powers of $r$ with $k\\gt 0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEkAAAA+BAMAAACVeiuWAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARLt2iatUMs1mIhDv\n3ZnG2WhPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACVElEQVRIDdVWTWgTQRh922yyG0ZrDoUeRBAE\nqSCkVMVLDwuihx5qIEcRV4+iNCgelNYEerFeGgrtSSE92KMJgl4MWKwnFdoinhQSvHiSRGyr2Npx\n/+anQydz7hzyve97b+bbmWXeBgiGlR8G7M1ciHuMGwXgXqWHIKIuBYp5z6Qa6yLT9E0qbxuv3jQM\nKiu3juotgwipximrdtGksjGexYxZNVfFb5PqPsq+0zWo7I/+C/zYNe3RsMqBp2epGL/4blLvOIzA\nBC1F8fXi/OqOoNYFDFGarrFCtuUzSLYZSmJrkxf6KwxafGpSKdMq4/CdIbt9t8BwFLP0L89TDLkL\nR5We+X+ME/EIX5XVmrTNII+TIz7HMbDollIBRpzlklK8vNNQKjjhrHlKzaUVpYLHKKprkdUNSTVw\nZeBMkj67MCV1maA5IfOG/NNxZr1HHhj1Eq6ffuAqx/+MwThrVvENOMxXaImWhPxkMzoNfAKEM9S/\nMAbIdBl+CxJc0mMsxWiNQ/BXTjaCGQTjzrmYtOR7fchLZpAt9FVuYuZ60rPsJ0QY6vxpp/G04GG3\nFJNkRRJJz/Hy7INizlk5GbNuQVYpOHV80olKQwqxJ+0bnrLDgi2OIbNHECVuLX0+BE9qURr+LHCk\nAukYrGWV5HndZ5B02gyqkfBbvHQtucCqJMjdP4+i8TUwgn2uSTKjI2yCCp/YZ7kDXJJd745+H8L1\nHl7VqmTXu61Vya6nV8mup1eFrpcO32TwFdWrZNfTq2TX06tC14s6TvfqKLue/s+C5HrPW0XNgf0H\nqTaxJzFyqH8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\sum_{k=1}^{M} r^{k} a_{k}$$"
      ],
      "text/plain": [
       "  M          \n",
       " ___         \n",
       " ╲           \n",
       "  ╲    k     \n",
       "  ╱   r ⋅a[k]\n",
       " ╱           \n",
       " ‾‾‾         \n",
       "k = 1        "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Equation 4 in the paper\n",
    "fbar = Sum(a[k]*r**k, (k,1,M))\n",
    "fbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADAAAAAxBAMAAABucDcRAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZqu7IjJ23e9UmRCJ\nzUQBtxRlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABnUlEQVQ4Ea2TPUvDUBSG3zZNtflog4OIi7WD\nuAhxEd2C/gCd/JgMiuKmU0fppJPgL7DBQd0s9Ae0sw7qpCjRbgoutRqlCtZ7bz7IRzPZA7n3Pe+T\nc85NQgA3tuquAq96koi+BlmQJRd/vEOlG7UKVRxdMgHwTK0AGOgUCpPABb8WAngH5lXcHNFm/lb8\nJ5C08L0QbpW2yJE++L0ZYMK8M03FHc6VgLKRVpf50IxcHpkiOH1J8gC5lcTq2PobIGvJDZLQ4anT\nlk42bFcwpFFBg53KlpgGcoajIbiC7C3yGHlf7kqeDCgbbubb6WNc6aLPcaTcBEb0lQgQnloKhCkd\n0qYegczY5UtdgaCK7C1GYP9hxLKNnBYD5AaMrih1MljpCnpgdmKiB63/06K6H61+YNZ9BMz9MKvt\nAlGzVfWYgayNiZfRbYA0s6RrJw0DuZ44sz9qqKL2evnSYEUhUF6MaXX+GAOKs2y6ZJq3B6Zp0LvY\nqURLbKfsmuCMRFP8pf8YiSBIqrC6Ak7BsM4KvAph/GvUdpzVaxVwSZJVfM4facN8cl6rjAQAAAAA\nSUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{R_{c}^{k + 1}}{k + 1}$$"
      ],
      "text/plain": [
       "   k + 1\n",
       "R_c     \n",
       "────────\n",
       " k + 1  "
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Integrate individual terms\n",
    "integrate(r**k,(r,0,Rc)).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Still need to get to equations 6-8 from equations 2-4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = Symbol('m')\n",
    "c = IndexedBase('c',(M,))\n",
    "j = Symbol('j',integer=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix and vector that define the coefficients, equation 8 in the paper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC8AAAAzBAMAAAD1HUx/AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZqu7IjJ23e9UmRCJ\nzUQBtxRlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABh0lEQVQ4EaWSsU/CQBTGvxYQSzHtn1A10REm\nY1jozAKJi3EBY8KAA8bJDQcXY4xM6khkNmHAzcE46GSiq0kRE3cNGhJhwLtrD9rrEQde0rzvfb97\n7+5yBSbR4XKDCy/XSe6Sb2vVM/wpTYsLv+NpCdArdJQfKKne0rKJhPYsABhpRL/RTloiKJaAVAER\nU3WcU+dlsnm1ANRsGLbYsQ5oA2CT+O7ml1SR6JHmErBDNT1VrrZCJWK9ctYkeY0WFPCYG0I5gvpx\nQo0Gd0mOHAIDqPc20YEwLETJNuEomuxQYUCukeyjFQZvgNKP1fFQDrLr0SLQLCPeZTcMQlo1cRM2\nqXMst4GfaaCDhBzltnflYHbXGEnja/bJM07Q2G/gG7LgvZpi+Uwio1nhOdsu1zMpAZzzxuo/4Jb9\ntGS12NGYN90ZAlBsoy4FOj5tKQAeqf/kOAeO88qW8FO5VXhz/dedFDqVNpwCEnsCuDrbZxeINwTg\nlrFWvuCBdy+7KX+XCdTjQq2weeOaiz/pxHUuUfWxCAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{R_{c}^{j + 1}}{j + 1}$$"
      ],
      "text/plain": [
       "   j + 1\n",
       "R_c     \n",
       "────────\n",
       " j + 1  "
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Skj = Rc**(m+k+j+1)/(m+k+j+1)\n",
    "hj = Rc**(j+1)/(j+1)\n",
    "hj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grfits = OrderedDict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  2.04800000e-03   6.82666667e-04   2.34057143e-04   8.19200000e-05]\n",
      " [  6.82666667e-04   2.34057143e-04   8.19200000e-05   2.91271111e-05]\n",
      " [  2.34057143e-04   8.19200000e-05   2.91271111e-05   1.04857600e-05]\n",
      " [  8.19200000e-05   2.91271111e-05   1.04857600e-05   3.81300364e-06]]\n",
      "[ 0.08        0.02133333  0.0064      0.002048  ]\n"
     ]
    }
   ],
   "source": [
    "Mval = 4\n",
    "mval = 2\n",
    "S = np.zeros((Mval,Mval))\n",
    "h = np.zeros(Mval)\n",
    "Rcval = 0.4\n",
    "for kval in range(Mval):\n",
    "    for jval in range(Mval):\n",
    "        S[kval,jval] = Skj.subs({M:Mval,m:mval,k:kval+1,j:jval+1,Rc:Rcval})\n",
    "for jval in range(Mval):\n",
    "    h[jval] = hj.subs({j:jval+1, Rc:Rcval})\n",
    "print S\n",
    "print h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   4375.00000011,  -44296.87500127,  147656.25000466,\n",
       "       -161132.81250545])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ck = np.linalg.solve(S,h)\n",
    "ck"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "// for m_exp=2, N_basis=4, Rcut=0.400000\n",
      "coeff[4] = { 4375,-44296.9,147656,-161133 };\n"
     ]
    }
   ],
   "source": [
    "# Dump in C++ format for inclusion in the unit test\n",
    "print '// for m_exp=%d, N_basis=%d, Rcut=%f'%(mval, Mval, Rcval)\n",
    "print 'coeff[%d] = {'%len(ck),\n",
    "print ','.join(['%g'%ck1 for ck1 in ck]),\n",
    "#for c in ck:\n",
    "#    print '%g,'%c,\n",
    "print '};'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   4375.00000011,  -44296.87500127,  147656.25000465,\n",
       "       -161132.81250544])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.linalg.inv(S),h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGMAAAA+BAMAAADNBbiPAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARLt2iatUMs1mIhDv\n3ZnG2WhPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACuElEQVRIDe1Vz2sTQRR+STa7yW4Tc6iIqHRA\n8AcILg3ixcOCKChYCooICq6CNw8B8SBUzEVIBGlU1JOa/6BBwYMEDf4AUaHpwZ5zsNRLaQpWLbUd\nZyYzs7NxU3fOOrC87/veezNvZ2fnAZCR2O8CGMsFiuOOi+MAV8txo1ncURJ+39NKOd4Fs+lrpXir\n8PxVQyclUZiB2mWdDEg19ibqR7RSDBjLQlUz5W4NfmilXIMbvtXVSTE++U9hYV1rx3Tm/x+bwcpw\nxX6k3goUYU3c++zOi63ncEsGzDBUkDwEzq9J+ui7gM4qQ0jwsM3hcSl8FSjRZggJHrYWXpFCvsSh\n0blCJ0Lkefyw2Cnu5jo3T3BdCi5HmQfbaGmIPHdW4F3jANe5SeJeFaq6aZ4xBGA/a8MSfFGdAM70\nelggbGLUh3Sl8r5ShbRrtWC2L2IE+30KjFot+lqIPEOlFLK7lCnDwB8VxuBOq+0RgMgzRdbJIsrU\nMflTZRTfgpNilc2QqdunGn0RTewFyvDZ4aJgSABqjcMHA2ripYB4e/x9gnUEINaehQoc8oSy+Esu\nbPmfYYvQVTvkwWnIFYSUCypzHGVF4ad2kc4a3ML2mlwFzK4aGODXFO6QPPVBQsiXA6wg5xslY5Y4\nOiOlwElqjhzk8JpQvcBLs+R/QoKnCpEZMAfOPKyLuZuuEhWUq4gEGsfugfVmFxdvhp2DWQpNWMyb\nLA8OCnuS7nWDKSca0uFLFAky9TQ7NNlgh7N/SRHzKDu8XWgbWzvY4Zd//p6RuU2Py/YlfgtGhqni\nZIWN29Pkom2rjoE4r1zJWP2kAzP+dYfarywUbzd6/YrGWnPxvgvvV2x6M14K71c6Kbxf6aSwfnWG\nHrMaxCyM9yudVWi/EiPmKqxfscLcuIXxftUrrCXW29DyfkVj7IXlOD/Mb5SdwMtMcvQdAAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$\\sum_{k=1}^{M} r^{k + m} c_{k}$$"
      ],
      "text/plain": [
       "  M              \n",
       " ___             \n",
       " ╲               \n",
       "  ╲    k + m     \n",
       "  ╱   r     ⋅c[k]\n",
       " ╱               \n",
       " ‾‾‾             \n",
       "k = 1            "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gr = Sum(c[k]*r**(k+m),(k,1,M))\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gr2 =  r**6*c[4] + r**5*c[3] + r**4*c[2] + r**3*c[1]\n",
      "0 c[1] 4375.00000011\n",
      "0 r**6*c[4] + r**5*c[3] + r**4*c[2] + 4375.00000010899*r**3\n",
      "1 c[2] -44296.8750013\n",
      "1 r**6*c[4] + r**5*c[3] - 44296.8750012704*r**4 + 4375.00000010899*r**3\n",
      "2 c[3] 147656.250005\n",
      "2 r**6*c[4] + 147656.250004658*r**5 - 44296.8750012704*r**4 + 4375.00000010899*r**3\n",
      "3 c[4] -161132.812505\n",
      "3 -161132.812505447*r**6 + 147656.250004658*r**5 - 44296.8750012704*r**4 + 4375.00000010899*r**3\n"
     ]
    }
   ],
   "source": [
    "gr2 = gr.subs({M:Mval, m:mval}).doit()\n",
    "\n",
    "cc = c.subs(M,Mval)\n",
    "print 'gr2 = ',gr2\n",
    "for kval in range(Mval):\n",
    "    print kval, c[kval+1],ck[kval]\n",
    "    gr2 = gr2.subs(cc[kval+1],ck[kval])\n",
    "    print kval,gr2\n",
    "gr2\n",
    "grfits[Mval] = gr2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def compute_smoothed_force(npos, epos, F, r):\n",
    "    forces = np.zeros_like(npos)\n",
    "    for ion_idx,ion_pos in enumerate(npos):   \n",
    "        for elec in epos:\n",
    "            dr = elec - ion_pos\n",
    "            dr2 = np.dot(dr,dr)\n",
    "            dr_norm = np.sqrt(dr2)\n",
    "            #print dr_norm\n",
    "            dr_hat = dr/dr_norm\n",
    "            #print 'dr_norm',dr_norm\n",
    "            if dr_norm < Rcval:\n",
    "                force_component = float(gr2.subs(r, dr_norm)/dr2)\n",
    "            else:\n",
    "                force_component = float(F.subs(r, dr_norm))\n",
    "            #print dr_norm, 1.0/dr_norm, force_component\n",
    "            forces[ion_idx] += force_component*dr_hat\n",
    "            #print ''\n",
    "    return forces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bare = [[ 4.26692459  7.40038688  0.        ]]\n",
      "smoothed = [[  9.24061106  14.86091659   0.        ]]\n"
     ]
    }
   ],
   "source": [
    "print 'bare =',compute_bare_force(npos, epos, F, r)\n",
    "forces = compute_smoothed_force(npos, epos, F, r)\n",
    "print 'smoothed =',forces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xss = OrderedDict()\n",
    "yss = OrderedDict()\n",
    "xs = []\n",
    "ys = []\n",
    "step = Rcval/50\n",
    "for i in range(50):\n",
    "    rval = i*step + step\n",
    "    xs.append(rval)\n",
    "    ys.append(exp(-rval**2))\n",
    "#xss[0] = xs\n",
    "#yss[0] = ys\n",
    "for Mval,gr2 in grfits.iteritems():\n",
    "    xs = []\n",
    "    ys = []\n",
    "    for i in range(50):\n",
    "        rval = i*step + step\n",
    "        #print rval, gr2.subs(r,rval)\n",
    "        xs.append(rval)\n",
    "        #ys.append(exp(-rval**2)*gr2.subs(r,rval))\n",
    "        ys.append(gr2.subs(r,rval))\n",
    "    xss[Mval] = xs\n",
    "    yss[Mval] = ys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHIFJREFUeJzt3XmYVOWZ/vHvA9gooCAqoqC44J6gUUSNRFvFBMzELWrU\njM4YtyRuSdTEnVaj0agTDDImOHHX4TIuGMA1/iwFFEVAxS2DRhFRIoLN1ggNvL8/nm5pll7rVJ1T\n9d6f6+qLbqg69XCovnn7XS2EgIiIlLd2aRcgIiKFp7AXEYmAwl5EJAIKexGRCCjsRUQioLAXEYlA\nh3wvYGYdgReBirqPx0MIl+V7XRERSY4lMc/ezDqFEGrMrD0wEbgwhDAx7wuLiEgiEunGCSHU1H3a\nse6aXyZxXRERSUYiYW9m7cxsGjAHyIUQ3kniuiIikoykWvarQgjfAnoDB5nZwUlcV0REkpH3AG1D\nIYSFZjYO6A+80PDPzEyb8IiItEEIwfK9Rt4tezPb3My61n2+EXA48Pr6HhtCyNTH0KFDU6+hFGrK\nal2qSTXFUFdSkmjZbwXcY2aG/+dxXwjhuQSuKyIiCck77EMI04G9E6hFREQKJOoVtJWVlWmXsI4s\n1gTZrEs1tYxqarms1pWERBZVteiFzEKxXktEpFyYGSELA7QiIpJ9CnsRkQgo7EVEIqCwFxGJgMJe\nRCQCCnsRkQgo7EVEIqCwFxGJgMJeRCQCCnsRkQgo7EVEIqCwFxGJgMJeRCQCCnsRkQgo7EVEIqCw\nFxGJgMJeRCQCCnsRkQgo7EVEIqCwFxGJgMJeRCQCCnsRkQgo7EVEItAh7QJEJE5PPw0zZkDXrv6x\nySb+66abQp8+YJZ2heVFYS8iRbV8OVx0EYwdC0OGwIIFsHCh/7pgAcyeDT/5Cdx4Y9qVlheFvYgU\nzcyZcMIJsNVWMGWKt+LXNm8eDBgAu+8O//Efxa+xXOXdZ29mvc3s/5nZ22Y23czOT6IwESkv48Z5\niB9/PDz22PqDHmCzzWDMGLj4YnjppeLWWM4shJDfBcx6Aj1DCK+bWRdgCnBUCOG9tR4X8n0tESk9\nK1bAVVfBfffB//4vDBzYsuc9+SScfjq8/LL34cfKzAgh5D2CkXc3TghhDjCn7vPFZvYu0At4r8kn\nikjZCwFOPBGqq73bpkePlj93yBDv2z/ySJg4Ebp0KVydMci7Zb/Gxcy2A3LAN0IIi9f6M7XsRSJz\n883w17/Ciy9Cx46tf34IcMYZMH8+PPIItItwsnhSLfvEbl1dF87DwAVrB72IxGf8eLjpJnjoobYF\nPfj0y9tv90HbK69Mtr7YJDIbx8w64EF/Xwjh8cYeV1VV9fXnlZWVVFZWJvHyIpIx//oXnHQS3H13\n/v3tFRXeqt93Xxg8GL7znURKzKxcLkcul0v8uol045jZvcAXIYRfNfEYdeOIRGDlSjj8cDjwQLj2\n2uSu+1//BdOnw113JXfNUpBUN04Ss3EOBF4EpgOh7uOyEMJTaz1OYS8SgSuu8Bk0zzwD7dsnd905\nc2C33eCTT6Bz5+Sum3VZmo0zEUjwn1REStW4cd51M3VqskEP0LMnHHAAjB4NP/5xsteOQYRj2yJS\nCDNn+jYHo0a1bopla5xyCtx7b2GuXe4SnXrZ5AupG0ekbIUAhx7qc+N//evCvU5NDfTqBW+/DVtv\nXbjXyZLMTb0UkXjdf79vZvarRqdoJKNTJzj2WHjwwcK+TjlSy15E8jJ/vm9aNmaMT48stBdegHPP\nhTffjGMbZLXsRSQTLr0UjjuuOEEPPs9+4UJ4443ivF650BbHAsCiRfDhh/DPf8Lixb4j4aabQvfu\nqz+vqEi7Ssmal1/2fenfead4r9munQ/U3ncf7LVX8V631KkbJ0KffeZbzE6c6OH+wQce8NtvDzvs\nABtvDF9+6R/z5/uv1dU+IDZokC+YOeww2GKLtP8mkqYVK2CffeCSS3y1bDH94x9QWQmzZkGHMm+y\nZmaevZSGTz6BRx+Fhx/2VYjf/z5897vQt68HfM+eTfd/hgDvvQd//7sPjp19tj9v0CA/jKJYP8JL\ndvzxjz7F8sQTi//au+zi2zA8+6zPAJLmqWVfxlau9DnPI0Z4UB95pPetHn542zemqldbC5Mn+yrJ\nu+6C7bbzKXdDhsS5M2FsZs2Cb33Lu3F22imdGkaMgAkTfI/8cpaZ7RJa/EIK+6JZtcq3la2q8j73\nSy/1Vnyh+txra/31brrJzxe9+GI4+WT18ZezY4+FPfeEoUPTq+GLL/wn048/9sPKy5XCXtaxapX3\nxQ8d6gc9XHONt+KLNT0tBO/muekmH7C78ko480y19MvN2LE+n/7NN2HDDdOt5Zhj4Ac/8JW75Uph\nL2sYPx7OO88Hq665xrtT0pyDPHUqnHOO1zNypG9gJaVv+XLYYw+47Tb43vfSrsYbN7feCgXYETgz\nNM9eAF8+/otf+CDZFVd4P/oRR6S/2GTvvb0/9cQT4aCDvEtp2bJ0a5L83XabD45mIejB3+tvveVd\nOdI0hX0JmzDB+03nzvUfqY87Lv2Qb6h9e2/dT5sGr7/uc6InTEi7KmmruXPhd7/zowazomNHn4I5\nfnzalWSfwr4E1dTAL3/pUx5vugkeeAA22yztqhrXu7dvS3v99d7Sv+gin6MtpWXoUN9aeNdd065k\nTQMGwKuvpl1F9insS8y0ad5CnjPHW/NHH512RS13zDFe81tv+cDx55+nXZG01Ftv+dGAV12VdiXr\n2m8/hX1LKOxLyIMP+hTKq6/2ucWbb552Ra3XvbsfcDFwIPTvD6+8knZF0pwQ/CfJK6/0f7+s2Wcf\nb0QsX552JdmmsC8BK1b43PUrroDnniv+0vSktW/vZ5MOH+7T5u64I+2KpCljx8Ls2b5qOou6dIEd\nd/TAl8Yp7DNu/nyfcTBtms+06dcv7YqSc9RRPmA7bBiccQZ89VXaFcnali+HCy/0w7432CDtaho3\nYIB+SmyOwj7Dpk/3PWf69YOnnsr2IGxb7byzf5NWV8PgwbBgQdoVSUMjRvgq1cGD066kaRqkbZ7C\nPqOeeMKPebv6ap/qVs47+3XpAg89BN/8Jhx8sO/KKen74gufQXXLLWlX0rz99lPLvjlaQZtB997r\nm4o99hgccEDa1RRPCHDddb6x2jPPeD+spOecc3yri+HD066keStWQLduvrtrt25pV5MsbXFcpm6+\n2beOff75+LYYMPNB6B49/DSiceN8Z0UpvunTfXO7d99Nu5KW6dDBV21PnuzTemVd6sbJiBC8NX/n\nnX6oSGxB39BZZ63ee+X559OuJj4h+BYcV11VWuNE6rdvmsI+A2pr4bTTfMn3+PGwzTZpV5S+Y4/1\nfvwf/Qj+9re0q4nL44/7or2f/jTtSlpHi6uapj77lC1dCscfv3oP+s6d064oW157zU/VGjnSp2pK\nYS1bBrvvDn/6U+l1h8yc6YH/2WfZ2iMqX9r1sgwsWQL/9m/Qtau3phT06+rf3/vuzzrL75EU1rBh\nvoVxqQU9wLbbehfUrFlpV5JNGqBNyeLF3mLdfnv4y198VamsX33gf//7/rVa+IUxZ45vrPfyy2lX\n0jZmq6dgbrtt2tVkj1r2KVi40Bep7LyzD8gq6JunFn7hXXaZjx2ldaZsEjRI27hEwt7M/mJm/zIz\n7U7RjOpqn2XSrx/8+c86sq81FPiFM2UKPPmkT30tZVpc1bikouYuICNn12TX/PneFzpggC9DV9C3\nXsPAHzs27WrKQwhwwQW+OV3XrmlXk5999/UjMXVewroSiZsQwgTgyySuVa7mz4dBg3w7gGHDymu2\nQLH17w9jxvgh0889l3Y1pW/UKJ8scNppaVeSv27d/LCcd95Ju5LsUduyCL780oN+0CAfAFPQ52/A\nAHj4YT/5auLEtKspXdXVvqvliBHlM3akHTDXr6izcaqqqr7+vLKyksrKymK+fCqqq/3AkUMOgRtv\nVNAn6aCD4P77/QSsp57y5fLSOpdd5mcKfPvbaVeSnPpB2jPPTLuStsnlcuRyucSvm9iiKjPrA4wJ\nIax3x/UYF1UtWOBBf8AB8Ic/KOgLZfRoX+353HM+R1xaZtIkX6n89tuw6aZpV5OcyZPh9NPL5zCT\nLG6EZnUfwurplfvuq6AvtKOP9kPYv/c9yOV8/3VpWm2tnzx1883lFfTgM93ef9/XsnTpknY12ZHU\n1MsHgZeAnc3sYzMrg6Getlu0CIYM8YPBhw9X0BfDySdDVZWPi3z8cdrVZN+tt/ruoqV+xOX6dOzo\ngT9lStqVZEsiLfsQwslJXKcc1K+M3WMPH/RS0BfPGWf4rJJBg+DFF6Fnz7QryqaZM+GGG7wbp1zf\nn/X99gcfnHYl2aHZOAlassSDfqedfCMpzaMvvgsugFNP9fUM8+alXU32hADnnedbGJdzd5cWV61L\ncZSQmhrf1GyHHeCOOxT0abr8cv9Pd/BgHzuR1UaPhhkz4OKL066ksLRtwrq0xXECamp8+lrv3trr\nJivqW7BvvunTMjt1Srui9C1a5NsXP/CAT1stZyH4wSvvvFP63Xna4jgjli71XRi33lpBnyVmfrzj\nDjv4PPxly9KuKH0XXuhTgcs96MH//XffHd57L+1KskNhn4f6oO/RA+6+W0GfNe3awf/8D2y8sa+0\nra1Nu6L0PPKIr0MYNiztSoqnb1+fgilOYd9GS5d6i3GzzeCeexT0WdWhAzz4ICxfDv/+73FukPXJ\nJ/Dzn/t92HjjtKspnr594YMP0q4iOxT2bVB/wlT37nDffR4okl0VFd6yXbgQTjklrsBfudL/zuef\n7zNUYrLjjmrZN6Swb6X6lbF9+ijoS8mGG8Jjj/nuo6eeGk/g//73fr7xJZekXUnxqWW/Js3GaYUv\nv/Sg32cfuO02Ta8sRUuXwpFHwpZbln/326uv+k+gU6bANtukXU3xzZ8P223ne1SV8uIxzcYpsi++\ngEMPhQMP1MEjpWyjjfyUq88+8/3bV65Mu6LCWLTIt5AYMSLOoAfvZu3QAebOTbuSbFBktcCcOVBZ\nCUccAbfcUtqtBPE592PG+MDl6aeXZ+Cff75vFXD88WlXki515aymsG/Ghx/6N82JJ8J11ynoy0V9\n4M+a5f+25TQP/+67/UCXW29Nu5L0aZB2NYV9E6ZOhYEDfSVmqR/ELOvq3NnPs121yrdXWLQo7Yry\n98wz8JvfeFeVtvdVy74hhX0jnn7aB2OHD4dzz027GimUDTeEhx6C7beHww7zsZlSNXWqryV49FHY\nbbe0q8kGLaxaTWG/Hvfc49PzHn3UT/KR8ta+PYwc6TtlDhxYmvvhf/SR78/0pz/5JAJx6sZZTbPE\nGwgBfvc7/8bP5dQ6iomZj8lssYUH/tNPl86//7x5/lPoJZeocbI2deOsprCvs3y5z2CYNAleesk3\nNpP4/OIXvgXGIYf4YeaDBqVdUdPq1w0ceaSPLcmattzS79GCBdC1a9rVpEvdOMCnn/o39+zZfsKR\ngj5up5wCo0Z5V9511/kAbhatXAk//rGv5r7hhrSrySYz78pR615hzwsvQP/+fmbs44/DJpukXZFk\nQWUlTJ4MTzzhB5pXV6dd0ZoWLfIum8WL4a67tMivKRqkddG+RULwBVI/+pHPS77iCn3DyJp69YLn\nn/eZOv37wxtvpF2R+/BD+Pa3/VCOsWP9gG1pnAZpXZTxtmgRnHCC/6j+yit+oIPI+lRU+OKka6/1\n/vu77/aGQlpeeAEOOADOPttn3lRUpFdLqdAgrYsu7CdM8FZa9+4wfrz3d4o056STfIbWzTf7zJc0\nTkAaOdIbKfff72s/tJq7ZdSN46IJ+4UL4Zxz/Jvl+uvhz3/2BTUiLbXHHjBtmof9wIFw0UXFOdB8\n2TKfafOHP3hjJeszhLJGA7QuirAfNw6+8Q346it4+2344Q/TrkhK1QYbwC9/6e+jefNg113h3nsL\nM2OnthbuuAN23tlnir38Muy0U/KvU+569/aV0TU1aVeSrrLez37uXLjgAu+XHznSl8OLJGnSJG91\nt28PP/2pNyTyPfpvxQp44AG45ho/MP3qq31AVtput93gr3/1Rl+p0X72Taiuhqoq/wfeemuYPl1B\nL4Wx//7emLj4Yhg92veOP+kkn7LZ2gPOFy3yc2L32APuvNM/nn1WQZ8EdeWU2Qra+fNh2DD47//2\nfUL0Y68UQ7t23qL/4Q+9a+ehh+C3v/XDUY46ylvnPXvCVlv5rz17+njRG2/4KVJTpsBrr/l2y/vt\n5weOHHaYBmCTpEHaMunGmTfPB69uv90XwFx2mf9PLpKmDz7w8aLZs/0AnDlz/ISsOXP80Pp+/fyI\ny/qP3XbzMQFJ3vDh8O673hAsNUl14yTSsjezwcAwvFvoLyGEG5O4blOWLPHDJ0aN8oUvJ5zgraPt\nty/0K4u0zI47+n5Lkr6+fT0vYpZ32JtZO+A24DDgU2CymT0eQkh8JvJXX8GTT3rAP/WU92XWr4Dt\n1i3pVxORcqFunGRa9gOAGSGEmQBmNgo4Csgr7Jct84HV+j7NKVP8x7D99/dj5EaMgM03T6B6ESl7\nffp4d9ry5fGuOk4i7HsBsxp8/Qn+H0CzJk70rpe5c30e7Ny5qz9mzvTB1b339v7M006DPff0s0NF\nRFqjosL3OqrPlRgVdTZOVVXV159XVlby0UeVzJjhB0b06wc9evjnW2zhfe8bbVTM6kSknNXvkZP1\nsM/lcuRyucSvm/dsHDPbH6gKIQyu+/oSIKw9SJvGoioRkXo/+5mvYSi1M6WztKhqMtDXzPqYWQVw\nIvC3BK4rIpKY2Adp8w77EMJK4FzgGeBtYFQI4d18rysikqTYV9Em0mcfQngK2CWJa4mIFELsLfuy\nWEErItKcmho/TH7xYt+4rlRkqc9eRCTzOnXyQ4tmz067knQo7EUkGjF35SjsRSQaMZ9Hq7AXkWjs\nuKNa9iIiZU8texGRCMTcstfUSxGJRnW1H0C+aFHpnASmqZciIq3UrZsfCfn552lXUnwKexGJSq9e\n8OmnaVdRfAp7EYnKVlv5WcCxUdiLSFQU9iIiEVDYi4hEQGEvIhIBhb2ISAQU9iIiEYg17LWCVkSi\nUlPj+9ovXVoaq2i1glZEpA06dYKOHX3rhJgo7EUkOjF25SjsRSQ6CnsRkQgo7EVEIqCwFxGJgMJe\nRCQCCnsRkQgo7EVEIqCwFxGJgMK+lczsODN7y8xWmtneSRUlIlJIXbtCbS0sWZJ2JcWTb8t+OnAM\n8EICtYiIFIVZfK37vMI+hPCPEMIMoAS2ExIRWU1hLyISgdjCvkNzDzCzZ4EtG/4WEIDLQwhjWvNi\nVVVVX39eWVlJZWVla54uIpKYrIZ9Lpcjl8slft1E9rM3s+eBC0MIU5t4jPazF5HMuP56WLgQbrgh\n7UqalsX97NVvLyIlI6st+0LJd+rl0WY2C9gfGGtmTyZTlohIYcUW9s322TclhDAaGJ1QLSIiRRNb\n2Gs2johEKbaw14HjIhKlVatgww1h8WKoqEi7msZlcYBWRKRktGsHPXrAnDlpV1IcCnsRiVZMXTkK\nexGJlsJeRCQCCnsRkQgo7EVEIqCwFxGJgMJeRCQCCnsRkQjEFPZaQSsi0aqthU6d4KuvoH37tKtZ\nP62gFRHJ0wYbwKabwty5aVdSeAp7EYlaLF05CnsRiVrPngp7EZGyp5a9iEgEFPYiIhFQ2IuIREBh\nLyISAYW9iEgEYgl7raAVkajV1ED37rB0KVje61STpxW0IiIJ6NQJOnaE6uq0Kykshb2IRC+GrhyF\nvYhET2EvIhIBhb2ISAQU9iIiEVDYN8PMfm9m75rZ62b2iJltklRhIiLForBv3jPAHiGEvYAZwKX5\nlyQiUlwK+2aEEP4eQlhV9+UkoHf+JYmIFJfCvnV+AjyZ4PVERIoihrDv0NwDzOxZYMuGvwUE4PIQ\nwpi6x1wO1IYQHixIlSIiBdS1qx8+vmQJdO6cdjWF0WzYhxAOb+rPzew/gSOAQ5u7VlVV1defV1ZW\nUllZ2dxTREQKzgw+/BA22ijtSiCXy5HL5RK/bl4boZnZYOAW4KAQwrxmHquN0EREWimpjdDyDfsZ\nQAVQH/STQgg/b+SxCnsRkVbKRNi36oUU9iIiraYtjkVEpMUU9iIiEVDYi4hEQGEvIhIBhb2ISAQU\n9iIiEVDYi4hEQGEvIhIBhb2ISAQU9iIiEVDYi4hEQGEvIhIBhb2ISAQU9iIiEVDYi4hEQGEvIhIB\nhb2ISAQU9iIiEVDYi4hEQGEvIhIBhb2ISAQU9iIiEVDYi4hEQGEvIhIBhb2ISAQU9iIiEVDYi4hE\nQGEvIhKBvMLezK4xszfM7HUz+7uZ9U6qMBERSU6+LfvfhxD2DCHsBTwOVOVfUvHkcrm0S1hHFmuC\nbNalmlpGNbVcVutKQl5hH0JY3ODLzsAX+ZVTXFn8h81iTZDNulRTy6imlstqXUnokO8FzOy3wKlA\nDbBf3hWJiEjimm3Zm9mzZvZmg4/pdb/+ACCEcEUIYVvgLmBYoQsWEZHWsxBCMhcy2wZ4IoTwzUb+\nPJkXEhGJTAjB8r1GXt04ZtY3hPB+3ZdHA6839tgkihURkbbJq2VvZg8DOwMrgX8CPwshfJ5QbSIi\nkpDEunFERCS78l5Ba2aDzew9M/s/M/tNI4/5o5nNqFt8tVdrnlvEur7V4Pc/qlssNs3MXi1WTWa2\ni5m9ZGZfmdmvWvv3SaGmtO7TyXWv+4aZTTCzfi19bko1FeQ+tbCuIxu89mtmdmhLn5tSTam8pxo8\nbl8zqzWzY1v73CLX1Pr7FEJo8wf+n8X7QB9gA7zPfte1HjMEGFf3+X7ApJY+N4266r7+J7BpErW0\nsqbNgX2Aa4Fftea5xa4p5fu0P9C17vPBhX5P5VNToe5TK+rq1ODzbwLvZ+BerbemNN9TDR73HDAW\nODbt+9RYTW29T/m27AcAM0IIM0MItcAo4Ki1HnMUcC9ACOEVoKuZbdnC56ZRF4CR/L5BzdYUQvgi\nhDAFWNHa56ZQE6R3nyaFEBbUfTkJ6NXS56ZQExTmPrW0rpoGX3Zh9cLHNO9VYzVBSu+pOucBDwOf\nt+G5xawJ2nCf8r2pvYBZDb7+hDXf5E09piXPLWZdsxs8JgDPmtlkMzuziDUV4rmFvG4W7tMZwJNt\nfG4xaoLC3KcW12VmR5vZu8ATwPmteW6Ra4KU3lNmtjVwdAjhdjxIW/zcFGqCNtynvFfQtkEpTME8\nMITwmZltgd/Qd0MIE9IuKoNSvU9mdghwGjCwWK/ZnEZqSvU+hRBGA6PN7DvAfcAuxXrtxjSoaSBr\n1pTWvRoGJDpumIC1a2qYna2+T/mG/Wxg2wZf9677vbUfs816HlPRguemURchhM/qfp1rZo/hP3Ll\n+4ZrSU2FeG7BrpvmfaobAB0JDA4hfNma5xa5pkLdpxbX1aCO8WbWwcw2a+1zC1jThPqaQgjzUnxP\n9QdGmZnh41RDzGxFa/8+RaipNoTwtzbdpzwHGdqzepChAh9k2G2txxzB6oHQ/Vk9mNbsc1OqqxPQ\npe7zzsBE4LvFqKnBY4cCF7bluUWsKbX7hH+TzAD2b+vfp4g1FeQ+taKuHRt8vjfwQQbuVWM1pf69\nV/f4u1g9QJv69956amrTfUriDTcY+EfdG/2Sut87GzirwWNuq/uLvQHs3dRzk/poa13A9nU3fhow\nPcm6mqsJ2BLvx6sG5gMfN/hHLci9amtNKd+nO4B5wNS613+10O+pttZUyPvUwrp+DbxVV9d4YN8M\n3Ku1a+qf9vfeWo+9kzVnvqRynxqrqa33SYuqREQioGMJRUQioLAXEYmAwl5EJAIKexGRCCjsRUQi\noLAXEYmAwl5EJAIKexGRCPx/wHkUWfd1mcMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbbe7aa3f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "for Mval in xss.keys()[0:1]:\n",
    "    plt.plot(xss[Mval],yss[Mval])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.379652839700475\n",
      "4 0.365714285715105\n"
     ]
    }
   ],
   "source": [
    "fx = lambda rval: exp(-rval**2)\n",
    "print 0,mpmath.quad(fx, [0,Rcval])\n",
    "for Mval, gr2 in grfits.iteritems():\n",
    "    #fx = lambda rval : exp(-rval**2)*gr2.subs(r, rval)\n",
    "    fx = lambda rval : gr2.subs(r,rval)\n",
    "    ival = mpmath.quad(fx, [0,Rcval])\n",
    "    #ival = integrate(exp(-r**2)*gr2/r**2,(r,0,Rcval))\n",
    "    print Mval,ival"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}