{
 "metadata": {
  "name": "",
  "signature": "sha256:6aa0c8d61bfbacdedb51051c430cae41ae815ce688a0fd6dc2c0431bb5a395f4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import numpy as np\n",
      "import scipy.linalg as la\n",
      "\n",
      "import matplotlib.pyplot as pt\n",
      "\n",
      "import random"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/numpy/oldnumeric/__init__.py:11: ModuleDeprecationWarning: The oldnumeric module will be dropped in Numpy 1.9\n",
        "  warnings.warn(_msg, ModuleDeprecationWarning)\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_random_sparse_matrix(n, row_fill):\n",
      "    nentries = (n*row_fill) // 2 # because of symmetry\n",
      "    data = np.random.randn(nentries)\n",
      "    rows = np.random.random_integers(0, n-1, nentries)\n",
      "    cols = np.random.random_integers(0, n-1, nentries)\n",
      "    \n",
      "    import scipy.sparse as sps\n",
      "    \n",
      "    coo = sps.coo_matrix((data, (rows, cols)), shape=(n,n))\n",
      "    \n",
      "    # NOTE: Cuthill-McKee applies only to symmetric matrices!\n",
      "    return (100*np.eye(n) + np.array(coo.todense() +  coo.todense().T))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prec = 1e-10\n",
      "\n",
      "np.random.seed(15)\n",
      "random.seed(15)\n",
      "\n",
      "A = make_random_sparse_matrix(200, 3)\n",
      "print \"%d non-zeros\" % len(np.where(np.abs(A)>prec)[0])\n",
      "pt.figure()\n",
      "pt.spy(A, marker=\",\", precision=prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "792 non-zeros\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "<matplotlib.lines.Line2D at 0xb4e9a90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD7CAYAAACBpZo1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGZtJREFUeJztnU9sFdUXx7+j7Q6NxNhH00fyTFuo1VKboDtCDS2GBRU2\nlWK0EdjgUqLoqsVEKQvjAiUxBJNuVNhIWUjTuHhI2DwSKSE2BohFSymNUJsQWDTg/Bba93t9zMy7\nM3P/nHvnfJIm0L43c865537vn7lzr+f7vg+GYTLNE6YNYBjGPCwEDMOwEDAMw0LAMAxYCBiGAQsB\nwzAwIATj4+Noa2tDa2srjhw5ovv20ikUCtiwYQO6urrw6quvAgAWFhbQ29uLdevWYevWrVhcXDRs\npTh79uxBLpdDR0dH+XdR/hw+fBitra1oa2vDxMSECZNjEeTf8PAw8vk8urq60NXVhbNnz5b/Zpt/\nifE18vDhQ7+5udmfnp72l5aW/M7OTn9qakqnCdIpFAr+3bt3V/zugw8+8I8cOeL7vu+PjIz4Bw8e\nNGFaIn7++Wf/l19+8V966aXy78L8+fXXX/3Ozk5/aWnJn56e9pubm/1Hjx4ZsVuUIP+Gh4f9zz//\n/LHP2uhfUrT2CEqlElpaWlAoFFBfX49du3ZhbGxMpwlK8KvWZJ05cwaDg4MAgMHBQZw+fdqEWYnY\ntGkTVq9eveJ3Yf6MjY1hYGAA9fX1KBQKaGlpQalU0m5zHIL8Ax4vQ8BO/5KiVQhmZ2exdu3a8v/z\n+TxmZ2d1miAdz/PQ09ODjRs34vjx4wCA+fl55HI5AEAul8P8/LxJE1MT5s+tW7eQz+fLn7O5PI8e\nPYrOzk7s3bu3PPRxyb9aaBUCz/N03k4LFy5cwKVLl3D27Fl89dVXOH/+/Iq/e57nlN+1/LHR1/37\n92N6ehqTk5NobGzEgQMHQj9ro38iaBWCpqYmzMzMlP8/MzOzQnFtpLGxEQDw3HPPYefOnSiVSsjl\ncrh9+zYAYG5uDg0NDSZNTE2YP9XlefPmTTQ1NRmxMQ0NDQ1lgdu3b1+5+++KfyJoFYKNGzfi2rVr\nuHHjBpaWlnDy5En09fXpNEEqDx48wL179wAA9+/fx8TEBDo6OtDX14fR0VEAwOjoKHbs2GHSzNSE\n+dPX14fvv/8eS0tLmJ6exrVr18pPTmxibm6u/O8ffvih/ETBFf+E0D07+eOPP/rr1q3zm5ub/c8+\n+0z37aXy+++/+52dnX5nZ6f/4osvlv25e/euv2XLFr+1tdXv7e31//77b8OWirNr1y6/sbHRr6+v\n9/P5vP/NN99E+vPpp5/6zc3N/vr16/3x8XGDlotR7d+JEyf8t99+2+/o6PA3bNjgv/HGG/7t27fL\nn7fNv6R4vs+vITNM1uGVhQzDyBcC11YOMkwWkDo0ePToEdavX4+ffvoJTU1NeOWVV/Ddd9/hhRde\nkHULhmEUILVH4OrKQYZxHalC4OLKQYbJAnUyLyay6urll1/G5cuXZd6WYRhBOjs7MTk5+djvpfYI\nRFYOXr58GUNDQ/+95OHD9937WfbP1R/AXf9cL7uwRliqEMRdOejzCgYrGRoybQEjG6lDg7q6Onz5\n5Zd4/fXX8ejRI+zdu7fmEwPPY0GgDpeR+0gVAgDYtm0btm3bFvmZ7u7u8r99P36iUU/MSv9coDrW\nOvwzVca2lJ3s+GhfYux5HsJumcY56uLA0CKr+RJW/0gtMV7uHST9rgjUXienZk9WyKIIREFKCIB0\nYhDF8jWpJQA1e5jHsUGs09pITggANWLAFc4MMstRdk6IXs+G3ElrI0khANT1DBi9LCeojLKUOjnm\nyByBrDpCVggAFgOXoFbpqNkTF9lDXdJCAMQXA13CQbnLmwZKtjDhyBYy8kIAxHNal9LLvA+l1omS\nLYw+rBACgFsqRh26c4tiLlsjBLImnSgWQlqCfHLRT1Xo7gXVup+qsou6LgkhiON42glEF7u+QT65\n6GdWUFV2UdclIQRxHeenCUwQnBPJISEESWAxoENlOZjo1i5jshdkey5aKwSAWTHI2n2jqKyAJrq1\naZAVT9uHYlYLAWBODNIWvOqXqxgxbIinjvwmJQRpKkf1dykuQqok7v4LMqHYq6CIyWFOJbLEivxT\ng2XSOFwtBhQXISVFtn3U/aVCdZxkPbo2FX/yTw1kwQnOqCRtflHOTyuEgHo3XxY2256WLPtOASuE\nQERJq7tdlBMrzDbKLYZqKPquModEry37c2FYIQQiVCcS5XUGFJPeVeLmQNJ5priIXlv258JwRgiC\noCwGlehYkJNFkmw+klWRdloIADvEIGpBDnXbKZPVSp0EkkKgYr9C0xXK9AIi0/7bQFSMXI8fSSFQ\noeSmxcB062T6/jYQFSPX40dOCGpV1rSvILuu7LbC5WIWckJQS3ldXtSRZbhczEJOCHTArY8c0rzf\nwdAik0IQNkSgkMhJbDBld9DaDcZOMikEQPAKRAobW1TaIDqLrcNuCiLJqMMZIZD5CrMJ4u47SG3D\nTRegkAemcEYIZL7CzNBApExklptNYic7X50RgrhUB5LFILn/quImUjF1LLgykRe17pnJk45UENYV\nt1EMTO+7Z6Il1Xlatgr/dFf0WmRWCMKwVQxqQc2ntPbY1I0Pgpr9LAQB2CYGOrvQsqBmzzI2lbtM\nWAhCoJqojFqyWu4sBBFQaR2o2MG4CwtBBCLbnqVZoSj6OROTVdWfYTGKh23xYiEQoHLOIOixY9Dn\nRa9rirjzCkltDYub69g2xGAhEGRZDGwrYNMsxysqbq6JhA5/ZL+TwkIgyLIIyCxk1ypAUuJuz0Y9\nbjoWOSW5R2YOOFFJZctmegFPXKhXnGpU70lhCzr9tEIIqC3/TCsGum3OSsVhkmNECGQeAmnyHDld\nG5J6nn2tOmMXRoTAthYq6mQiHRXU9+2LGWMXVgwNTEOxR8IwMqlL+sVCoYCnn34aTz75JOrr61Eq\nlbCwsIA333wTf/zxBwqFAk6dOoVnnnlGpr0kUfVYkR9XMoCePEjcI/A8D8ViEZcuXUKpVAIAjIyM\noLe3F1evXsWWLVswMjIizVDKpD14VdehqDzPYCc6GoNUQwO/ysIzZ85gcHAQADA4OIjTp0+nubx1\nJJ0z0NXqZ7F3weInRqoeQU9PDzZu3Ijjx48DAObn55HL5QAAuVwO8/Pzcqw0TJxkMvEKM4Xju6kS\nJH62+6SCxHMEFy5cQGNjI/766y/09vaira1txd89z4MXEvHh4eHyv7u7u9Hd3Z3UDC0kOVFX5/ie\nwvHdNmGLTzJyqFgsolgs1r6XX92/T8ChQ4ewatUqHD9+HMViEWvWrMHc3Bxee+01/Pbbbytv6HmP\nDSlchScRs0nS8tFRrmH1L9HQ4MGDB7h37x4A4P79+5iYmEBHRwf6+vowOjoKABgdHcWOHTtSmGw/\ncYYJcYcfDF2STh4bPVcjSY9genoaO3fuBAA8fPgQb731Fj7++GMsLCygv78ff/75Z+jjwyz1CJbh\nFpwWWS6PsPonZWggwxCGYdQjdWhgI6ZnisPuz+8RMBTIjBCY7oRUjhurzy1UP0Gk9vqmoOZXXHtU\nfz4OzgqB6iQxmYQ2TULpxmS5JHnMrPLzcXBCCIIKX3XypzkVKO2iI9ETk10nzX6RzEqcEAKdXWuZ\nuxPJuBaFo9xNXyPudbMsnmGQEwKKhSRzN9/q61L0VxQZIpTmGkn39LNdPFVATghkrw3XGfikLxzZ\n4l8teAu22lC1mZwQBJEmeLW+KzN5o0SsVgula9szlVCyhYmHFUIQB2oz6nEOO5G9V70JotZLMP9C\nMRbOCYHNrZLsvepNEGYPNTtNQjEWzglBNba1UFTtYughM1ecFwLbWqiwFYgycVVsbD8hKS4yc9h5\nIagF5WfhshGZuDSFjjUVFP1OA/cIYiL7DDkV1wi6pux1BqoqApWFUXHssEW8o+AeQQAqKjuFrqYN\ni45Mt7RJejqmbaaGM0KgqlWmgA1iYBIq5SQKxbJ0RghMoDMBq8WAYjLJwFW/KpGVNzxHUEEWEmeZ\nSjGwrRUUxVa/TJ3KHYcoG60XAlsTJylZ89cWbCiXKButFwLKqFrM5NIBKjYQNSRzJTaZFwKVBSl7\nMVPlsCDJNlcuvNhkgqjX0F2JTeaFwFRBBrUytSpqUEKKVm4deyMy9mKdEMhqwVXuKizy/aBKnfSl\nI1e6p7XIip8msE4IZLVqUd32tPfQ3fJmRQyo92hkPN41VY7WCQETTFbWGVBGxpZ2psSOhcAhsrDO\nwDSuCiwLgWNkZZigmrDHhK4KrHVCwEkeTprHi66Tdgs7lQJAoaysEwJXFVkGMsaorkI5HhRss04I\nGHEovEZtC7bGQpbdLAQJsOUcglqLjii0RFQqIIVYBBFnkVkaMiMEMnevUXnOgihx/KE2Z1B9GnSS\n72UFXQKVGSHQsXuN7ESVvU0ZlYpk2zN2EajENimZEQIdyE5U2ef3mRAD2yuIKJRFSgRnhMC1hFPl\nj24xoFhBXMsVGTgjBK4951XpD6Vhgi5E5yWyFpdlnBEClVBs1USotbOza0kvYydrimUtUk5py9KY\nELiWhBQROfBD9otKQXsr6CprV1t6EXFKK2DGhICi8mYR2asRg/ZWoFDWFGygDA8NGG2tpc2tsuuw\nEFSQxUT1vPjbniWFW2W9xClPK4SAwhjTVaq78CpPYXYRyrGKk89WCEEWKyhjB5W5SVkUamGFEDD6\nUHEKc1awucFiIRAgi5WCxSBbsBAIIPOMOZvQJQYm4+VKWUUh4mOkEOzZswe5XA4dHR3l3y0sLKC3\ntxfr1q3D1q1bsbi4WP7b4cOH0draira2NkxMTEgz0jZUdxF1vyvg8tMESt15le+X1CJSCN59912M\nj4+v+N3IyAh6e3tx9epVbNmyBSMjIwCAqakpnDx5ElNTUxgfH8d7772Hf/75R4qRzEpMnJvAqMdk\nnCOFYNOmTVi9evWK3505cwaDg4MAgMHBQZw+fRoAMDY2hoGBAdTX16NQKKClpQWlUkmR2fLQsY7b\nBWyJgS12RmHCh9hzBPPz88jlcgCAXC6H+fl5AMCtW7eQz+fLn8vn85idnZVkpjp0rON2AUqTh0lf\nLkpiv2tvnoaRarLQ8zx4EZGK+tu/f09zd0Y3ulYgitqh43tZaQTq4n4hl8vh9u3bWLNmDebm5tDQ\n0AAAaGpqwszMTPlzN2/eRFNTU+A1hoeHAQBDQ0Cx2I3u7u74ljPGWO4dZKWS2EyxWESxWKz5Oc/3\no4vzxo0b2L59O65cuQIA+PDDD/Hss8/i4MGDGBkZweLiIkZGRjA1NYXdu3ejVCphdnYWPT09uH79\n+mO9As/zUOOWyqGSxFTsSIos+22PQxp0+x5W/yJ7BAMDAzh37hzu3LmDtWvX4pNPPsFHH32E/v5+\nnDhxAoVCAadOnQIAtLe3o7+/H+3t7airq8OxY8dqDg1MQSXpqNiRFFk9A4pxUHnEWWXMqPhes0cg\n/YYSewSy1ZRiy0TRpmpssJH5l7D6Z/XKQp27BpuCok3VUHqiUImJfRYoxkEEq4UgKdQKi5o9SQgS\nA9N+6RLRtLsxxY2TirhmUgiotbLU7EmKzhOEXSJunFTE1WohMN3iyMaEPxROZ3KtHG3EaiFwrcVx\nwR+RRUfVf3PBb9uxWgiY9Mh+6lJ5XcqnMIuQpZ4KWSGgVgjU7EmKSj+C5ghsjlsawUrqt8p4RV2b\nrBBQazWo2ZOUaj907DWQdFbc5gNbKZ74HHVtskIA0GpNXN2pR4fAxRWDuKvuZMXNZrG39sgzESgV\njK4KYwIdAqRymKBbMFyEtBAwetC58MZkZaTUsMjG2rMPGb1QGW/rEIM0189qr4GFICNQ2udQtS1p\nru9yryEKFoIAstoq6IRjTAsWggB0tQpZPoWYyrZnWaFWnFkIYiA7aU28HReXOEuFa30+iDTrDFxB\n11ObKJwTglpBTRN03eNHCgkfd6yfdIPQJOsMXIGCP84JQa2gUgi6KDbZmhbTjxazjnNC4DKuVxQZ\nYiAzRq7HuxIWgpjwgpj4xO326xy+JT0sRSc6co6EENikvFSSwyaS7MCjKydq2SbLjriTrpXIyjny\nbx9y5WKqoZITsuwwucBK5D4khIACNvVKbIPiu/nMSlgI/oNKC+Qiad/NZ0FQj9VCwAmSDfjRonrI\nCoFIwXMr/i8isbK9IqURA9t91wFZIeBKLo5IrFyIZ1Ix4CFGbcgKAVWynEyqfNe1zkBEDLNaviwE\nVeh6pmsjqnyntM7AlvKV7T8LQRWqE8Hz3Dg00zTUliPrRnaeshD8h86VbGkPzWT+RfdyZFsRiREL\nwX9kJSlcg8utNiIxYiGwAJu7sMuYOsGHOlSGidYKgc2FL0LlaT+qWj2dMdRxgo+NOZF0mMiThf9h\nc5cwzmIpU0dgVWPD/oqqVyDGubaOo+RkYq0Q2IyNIqbq1GTZ96HyaNG2MmYhMIiNXVkZ6DjXIKux\nTYrzQkA5IWxrNaKgFmcWg3g4LQQqJ9qi7plFKIpalsUgrt9OC4GJ5KRYIVwh6QtHot9zSTTi5qHT\nQsC4RdoNTlRd3wVYCBxF9KAXl1rBKLLiZ1JYCIiSNnFFD3rRuZZARmVM8woyi0E4LAQKqE64sAQ0\ntae+yvf5VX4/7TWCViAmiYWLgkJKCEQrEHWqkzUseU1tlJHlsTCwsneQ9KxG1yAlBJUBNvHojyIc\ng3AoL0eWhS4bSQlBJaIVwIbCtBkT8RW9Z1KRrOwNUM8fXQ1BpBDs2bMHuVwOHR0d5d8NDw8jn8+j\nq6sLXV1dOHv2bPlvhw8fRmtrK9ra2jAxMSHV0LAC4xZTLS6uxah+44+6GOjA8/3wsJ8/fx6rVq3C\nO++8gytXrgAADh06hKeeegrvv//+is9OTU1h9+7duHjxImZnZ9HT04OrV6/iiSdWao3neYi4pVFc\nG4645o9KTK1C1X/P4PoX2SPYtGkTVq9e/djvgy40NjaGgYEB1NfXo1AooKWlBaVSKYXJchE5hNK1\nSkPFH9Mtruhr37rtpFI+QMI5gqNHj6KzsxN79+7F4uIiAODWrVvI5/Plz+TzeczOzsqxUgIUDqHM\nKqbjyysLaxNbCPbv34/p6WlMTk6isbERBw4cCP2sZ7gpSHN7062Yakz6R2Fhkolrm6SWX3VxL9jQ\n0FD+9759+7B9+3YAQFNTE2ZmZsp/u3nzJpqamgKvMTw8XP53d3c3uru745ohhIzFJzYiMvaU5V+S\ncS6FhUm1rl3tF7X5FlF7isUiisUihoaAimr3+PWiJgsB4MaNG9i+fXt5snBubg6NjY0AgC+++AIX\nL17Et99+W54sLJVK5cnC69evP9YrUDVZSL3gXCUqzknKgFK5UbJFFmH1L7JHMDAwgHPnzuHOnTtY\nu3YtDh06hGKxiMnJSXieh+effx5ff/01AKC9vR39/f1ob29HXV0djh07pnVoILqaj0lOUMWQPfcS\n1iKbYHkCkZJNqqjZI5B+Q4EegamA21LQttjpCnHjTbl8Ej0+NIWKIMbZOZg6tex0dcJLJ5UxjPto\n0ZY8qoSkECTF1Nt81MiSr9XIEsHqNxVdX4HolBBkpQKILI7KKrJzICvLkZ0SgqzAi6PS4cpbizLt\nYCFgMgfldQxxkGkHCwHDJIBKr0AWLAQWwktszRO07ZnNsBBYCJWDUStxpULExcScgYr7WS0EWUw+\nqi9SURk3myCuGKjeoToJVguBjcmnKglkLpjKosCmpVoMonauppi3VguBjahKApnXTXItFg+x3ZEp\nPX6shIUgBnELkGKBq4JiKyeLuMuLa32eYqxYCARIupWZzgKnKjouTKTFHVLJaPV1x42FQACKCl4N\nVRt122UyDtXLkWVdSwcsBIqh2lIz6rGp7FkIFEO1pWbUQ3ViMAgWAoYktlSgWtiyAtEpIaAebEYc\n13pSSXoHOvPZKSGwMXlYvNxB5LEh1Z2OrBGCWgG0tUJV77wcBlX/qNqlm+UVg7LFYPnaqrFGCOLu\n0580eCYT28YNR6japZvlOIjEQ1QMdB7FZ40QxCVp8FQHnVtQmpg495DSCkRnhYAqLrWgVERNhh0m\nyoXS40UWAiYxVESNih1JoGI7C4GjmGhpZM7LUGkpdUDBVxYChVCdeIyLqB8y52XinixkMxQWHbEQ\nKMSVtw/T+KEjual0r9NSOWfAbx/WwHb1VwXVpx2uVFJdVB+8qgvrhEBlgFhkwuEKHY808yUmniZY\nJwQqEUl20QJiUck2Qbkksqdh5cIkftdAAqqCKNoyUmpBWZSCMbGIKOjfUZ/XZaMTQhAULFcm6mSQ\nZudjE+iyy5RYy94DUQZOCIHp1tf0/ZNC1W6qdsmC4t6XTggBw7iO6l4BC4FmqHbHGdqoHiKwEGjG\n5W5vVkVO95yGikVHJISAwoYccd4PZ4IxPUErs3yo7iS0fD/Zi45ICAGFDTlEH+eoJupZMwvR/0n7\nfkKS6y9juhxULDoiIQRhmA64CaKeNcdN9CzGTwemh3cqFh2RFgLTAbcdjp/7yBIDkkIgwzFuDWuj\nM0aex2WiChliQFIIZLRk3BrWRmeMfJ/LRCVpxYCkEDAME580QmtMCLibyDDySVqvjAkBdxMZHWSt\nwUm67RkPDZgyLlaarDY4cecMWAiYMlmtNK4SRwyMCEGxWDRxW2245F9QIlHyT3YvhpJvMhAVAxYC\nBbjkX1AvQaZ/aSuy7F7Ma68V5V6QACJiYN3QIMwhHeNbCmNoCjbEIen5fqb8HBoyc1/V1BID64Qg\nLHF0jG8pjKEp2BAHqofRZpHImPqa2bx5sw+Af/iHfwz8bN68ObBeer7P2sswWce6oQHDMPJhIWAY\nhoWAYRgWAoZhwELAMAyA/wH00+IUObCr8wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x9e20910>"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ainv = la.inv(A)\n",
      "print \"%d non-zeros\" % len(np.where(np.abs(Ainv) > prec)[0])\n",
      "pt.spy(Ainv, marker=\",\", precision=prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "7264 non-zeros\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "<matplotlib.lines.Line2D at 0xb69ced0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD7CAYAAACBpZo1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfU+onsX1/7k/zK6Lry3NNeQKV/PH1BpvA9adGInX4qKp\nbsRQSvDPRpcVka4SFzZx0W5sC0UqZNXaTY2LGkIX15ZuQmniwlAUvKUx/2iNAVEwaJ7fwj63k5Pz\n53Nm5nnf572ZD1ze931m5pwzM2fOv3mic13XddTQ0HBD4/9NW4CGhobpoxmChoaGZggaGhqaIWho\naKBmCBoaGqgZgoaGBpqCITh27Bjt2LGDtm3bRi+//PKk2VfH4uIi3X333bRr1y669957iYjo0qVL\ntLy8TNu3b6eHHnqILl++PGUpcTz55JM0Pz9PO3fuXHtmzefQoUO0bds22rFjBx0/fnwaIocgze/g\nwYO0sLBAu3btol27dtFbb7211jZr88tGN0F88cUX3ZYtW7rV1dXuypUr3dLSUnf69OlJilAdi4uL\n3UcffXTNs+eff757+eWXu67rusOHD3cvvPDCNETLwp///Ofu73//e3fXXXetPdPm8+6773ZLS0vd\nlStXutXV1W7Lli3dl19+ORW5UUjzO3jwYPezn/3sur6zOL9cTDQiOHHiBG3dupUWFxdpw4YN9Pjj\nj9PRo0cnKcIg6Ng7WW+++Sbt37+fiIj2799Pb7zxxjTEysJ9991HN9988zXPtPkcPXqU9u3bRxs2\nbKDFxUXaunUrnThxYuIyRyDNj+j6PSSazfnlYqKG4OzZs3Trrbeu/V5YWKCzZ89OUoTqmJubowcf\nfJDuueceevXVV4mI6OLFizQ/P09ERPPz83Tx4sVpilgMbT7nzp2jhYWFtX6zvJ+vvPIKLS0t0VNP\nPbWW+qyn+XmYqCGYm5ubJLuJ4K9//SudPHmS3nrrLfrlL39Jf/nLX65pn5ubW1fz9uYzi3N95pln\naHV1lU6dOkWbNm2i5557Tu07i/NDMFFDsHnzZjpz5sza7zNnzlxjcWcRmzZtIiKib37zm/Too4/S\niRMnaH5+ni5cuEBEROfPn6eNGzdOU8RiaPPh+/nhhx/S5s2bpyJjCTZu3Lhm4J5++um18H+9zA/B\nRA3BPffcQ++//z7985//pCtXrtDrr79Oe/funaQIVfHZZ5/RJ598QkREn376KR0/fpx27txJe/fu\npSNHjhAR0ZEjR+iRRx6ZppjF0Oazd+9e+t3vfkdXrlyh1dVVev/999duTmYJ58+fX/v+hz/8Ye1G\nYb3MD8Kkq5N//OMfu+3bt3dbtmzpfvrTn06afVV88MEH3dLSUre0tNR9+9vfXpvPRx991O3Zs6fb\ntm1bt7y83H388cdTlhTH448/3m3atKnbsGFDt7Cw0L322mvmfF566aVuy5Yt3R133NEdO3ZsipJj\n4PP7zW9+0/3oRz/qdu7c2d19993dD37wg+7ChQtr/WdtfrmY67r2z5AbGm50tDcLGxoa6huC9fbm\nYEPDjYCqqcGXX35Jd9xxB/3pT3+izZs303e/+1367W9/S9/61rdqsWhoaBgAVSOC9frmYEPDekdV\nQ7Ae3xxsaLgRcFNNYshbV9/5znfonXfeqcm2oaEBxNLSEp06deq651UjAuTNwXfeeYcOHDjw33/k\n0VHXrb+/fn7r9Y9o/c5vve+d5oSrGoLom4Nde4NhJnHgwLQlaKiNqqnBTTfdRL/4xS/oe9/7Hn35\n5Zf01FNPuTcGc3Pr1yBocyud86TXbD3v0ZgwzXWuagiIiB5++GF6+OGHzT67d+9e+9519gLwtrQM\nwZ/XOFzW4eU8NfTz0/oiNDR5kHmisqJr1vfp+6f7l9KQ6KHySuNS3toYj38U/dw4/cj+R9ehb4vI\nnjNXq4Q38VeM5+bmSGNpKT1yUL2xukxyX+8Q5ioeMgeugES+0Yzyr0Wrl60GrxzZvMOqOROrX+Tg\nR+WOGOoeNdb3K3ry+RuVIfiqHT+0kfaayp8DTVkj0ZBHO2fsWJDr4cY0x9zoJ6dPDixDMLp/a5Ae\nFimUsQ5N2q6FQVI75yeNnZuzQytLJok/0t+aq+UtpGhCm5P06cmV057KHF2XUjm8+aUySr8t2aV9\nsPigfSTeFpA+1hqPzhAQXa/EKayNQGApocU3wk8K/TQFqOEFNcNgGT2pX25NA43QojlwFFE5LWOa\n5uz8u8cT0SPPoXl9c/TIkmlUhgCx3N5GaJ4vovS5BwLpl1MU8g65RjM1SNK8kfVGZZS+W/JYNLQ+\n2t5789AMYo16jwY08kONtEa/1h6OwhBwBekVV7KsPMzsgR50TwZLmYbOSbliSgcHOTDRQlfJnFAF\ntAwgEuVofCN7IkVo2u+SyLMkakSMpDQ+knJKGIUh8Dwot5TpgmkbpimKpmDIgloHUnsmySXlo1w2\nSx7PA2tKjRqRiNJ7e5GTE3O61hh+kNDoxOuDpAB8fMQAe7pk8e7nHTHCHkZ3a9DQ0DActPM3ioig\nB+JR0fbSnAlBrfyM04u25fKpQbP2GtSmlcsvGhlF1qFGgbs2RmEItDA5/fRSAe23xksK1bX6A1qX\nQGVA5ZfSGyslscJiLWWKrheHd+OijeeyInJ5Yb92eLVULOWnIZpq9LJH9FRaCylNtOghMB3NGFID\nazOkNm/z8uTyaQ7Bt0EGckCH1IvUINWgGdXxEnrWGJqVF4oQaFZXQ6RfbUubCy0KidKQftdODUrp\nSOPRYlhu9BUttJXOMXrzk2KIdINjaoYgFRyt1Ft0osqNhvkIaikLUiOocZg1pfTCcaRNS1Nqyxnd\nLzQFkMJxNFrkcqJyaelgBKXpw9QMQcQKSt8j4yWFQl5YkehZG4XUMXg+rMlkvReBGFEvFO3HWtdv\nHqQrPi4vsl/IOwTS2kh9vFpGP06TD9EVa41qrGGtdCTioEZRIyin+dXnJGZSKzflyujlo1pbSe5Z\nQ+GGqK3kyjWpXNuiM631QGmO6voQyV2RkEnyHJrHtjy5F32kvNL6hGZpNdklOkT2iz08grDm7M0p\nMh5BZH/Q9YqExhpPrx+yDtraI7JptK094TwRHYvCGrMuIoL/0b7ey5bQyGlH+k3S49XCpKKuUs9a\nc434nCO0PdmGiBwwuWbkv0fQ0NAwHEaVGlhAws30uxdSRn8jcmnjIrSiaU80RETl80JYD5FUDU2l\noqgdJufQyUlp0HZkj0v5jsoQpGEhshhe6JX2sX57xoW3ae0pbU25I+Eml1XKHbUqdfpcm3u63uk4\npIqvyWvdl0voZZBk88ZFjVjEaVi5u0dfk98zmj24HmnjawbW6zo1KFksLZ/zaE4zj7fA5RpaTot+\nqSyRPUiNXWltZ5agrfHMpAYSckNVzVt63rofa33Xxnt39xJ/LXxGoqKIt9bm5/FAwOVA7/QtWTR4\n/SWjjfJAjYX2exJpChI1Rdd4VBGBZMWI9Os1rc2jq/Xx6E3SWyCHtgaPmnfoNWlGeFrPIuNzaQ2N\n0rleO24GIgJ+Ry/lq2lfKx9FvVPfN5rbRqHR1Dx6Pz9JrpqRQ+lc+dpJERPPr6O5vCartaZeP+vd\nDfS5BrSOgCInSoliVBHB8LzHY+UbGqaBmYgIiGRvgHg/NG/y+ud4BdQje2OtOgJC28tVS7xGSe3A\nGy/NFZE/d52nCUS/EVmjUQrnc13bjRQRNDTc6JiZiEBCbU9WqzJeShfxliiNUlmi9KU2r3ah0Yju\nUUlkgMoVQen+RfWght5wjNIQcOXQCnk5m93XCazwX1M0q3iZyolukHUbksrD29BCKLJW0YObQno5\nS7t21fj2Yzhv76ZI04/0WeTAlBienBew+j7pmvFieXEBECyuEo3UEHgHLu3HkSqL1id9zhUK5S3x\nlORCvFWqtNLh4kAPgqSg/adE37qi8movUZS+x4BU/bkOeGvq8YrK5q2ZJUeOHmq0tN8pRmUISkIe\n7SBp7T0QjxqlyelLY6V+mkFJvQfC05MFlQvx7qjni8gW9YhexCbtY7QYmRpFiwcSPXpz43JqehG9\n8h59RGB5cSuET2FZv3QTJWhhpHfYIocROVRaiBn1DtKaRY2ARSunDzoeMarpGCRqkXTDW09PZzQe\n3jsrEXjGROKP9JUwyv9mIRL28n4eHyTPs/pGNteLbJB81AodPfm8CMbj7fHj7RKNksMgpWhe2M/5\nc3hrwQ++5kg0Olw+JC2wohLUOFrrjabJRCP5bxZKYRWiUJai8APlGQ8kdJdk1aB5lGi4i47x8ntN\nmb1IxfOM2gGScl1tTsh+a9ES6gy8CFGjz+XTxqXz9FIUvh+cXpoGonqJtI0+NdCKOukzbUz6uyQ8\nQsJAvtkRekjuKI2XDlI0/JTGSzLmIFI7sIxxjkx8LS1Pi4Af/hrrkyJ3jsjzUllHYQg4rE1AJ8wV\nQvJUGl9LmXK8s/YMHSu1eZ7ce4ZGNtIY7TdCj9ONpF3pGHQ9U+/K6Wj0o205QHjV5OlGdu3NwoaG\nGwejerOwxNuiNHOKe1F5alWHPUSLlQi99YSh5jMJutFaUC2+HFMxBCVhP1pM8nLFNDTmYapVXNOK\nXJaMVqUZAVL/yL0h8CrXSD+Ej/YsUm3XZEFSHvTwWSmEt2/IPKwQPS0Q1jb8XmowihqBtBmRirD2\naVVutTyM3zL0zzR5rMqv1idSYOPPpbWKZFp9f16dluSLFEc1+fpnvBCnAa2WazxRJ5Nz2Ky10ByH\n1EfjHXEYmt5rcnuOcRSGoIekkIhi8TFaf42P5/kQWXlbSfU6lclTGImP5xX5WiC3LYgsXrRmFTCl\nw4x6/shtU9+GRlm5HloyPFEHEL1d0oDIPgpDoFX0pcVLN5H30xTf82ySYngKblXTuew5BkGaZ/pc\n6iuNTfl7MktzSGlqcnp9kDCdy+YdAmk/0f1DIi++/t7aIbc43ppL+o8aCS8i8Hivm1sDrjzRkDnK\niyiPfjo2l05uOjBGlMgmjY3QQ6KhodYZ3XvuUCxji8xlVLcGUXhheOR3TV5RWXgb99heCCq1l6Qf\nCI9IaIzkvRzRNagBJJ1J2y3ZvDSsBEjkgsKLCqYSERB1bsEl0ia11/AMuJWt43EtOYjqefVackfk\nqumpS/ujNIniMiOyaF4+Vw9RfDWn9v8+zEbUQKEbN82QvURhxwzkMJU4nFI5po32P0FtaGioXyNY\nXFyku+++m3bt2kX33nsvERFdunSJlpeXafv27fTQQw/R5cuX8yVm8CrPaN5Us4ag5Ye5sHJzKXeP\n0EN+5yC3dmD1yZUrWsOxxlg6ZdGptR5Wf6luU7qX2YZgbm6OVlZW6OTJk3TixAkiIjp8+DAtLy/T\ne++9R3v27KHDhw+XSbfGS74yS/M47aqO/+ZFk2hwIoWXyDWhZ8gseS0eFs/IlRdCL8ojt7aRGzB6\nV2ratZ5GC7kyjbRpQAuq1rPSIDs7Nbjtttvob3/7G33jG99Ye7Zjxw56++23aX5+ni5cuEC7d++m\nf/zjH9cynGBqMI08bay54Y2KMRQepwG9VlI5NZibm6MHH3yQ7rnnHnr11VeJiOjixYs0Pz9PRETz\n8/N08eLFXPJwWGZd31geTwvvvCuiaMSh8bJksMZE6Ui/0VSKpyY5fC3amizeelly8ef9fqARWER2\n6zmapmjzL00xpEjNjFhzI4Lz58/Tpk2b6N///jctLy/TK6+8Qnv37qWPP/54rc/Xv/51unTpEhN+\njg4cOLD2e/fu3bR79+4cERjdrz5Lq9yWAYlcuVn9orRKeM0ivPnkphylcpSuczoeoRWZp0ZvZWWF\nVlZW1n6/+OKLw90avPjii/S1r32NXn31VVpZWaFbbrmFzp8/Tw888ICbGkgTqKHYOTQQBYxupJa3\nloSriJHJRWQ/IryGkJkb1KEMhLXvpfob1akcHukYUt4jyEoNPvvsM/rkk0+IiOjTTz+l48eP086d\nO2nv3r105MgRIiI6cuQIPfLIIy6tSAEngsiComPSdoS+lCZoRSgrzLYKpZaMOdD2Ay02Rujm0JHG\npZ9DRAnWvlv6i4T3UZ3K6YuMyYoIVldX6dFHHyUioi+++IJ++MMf0k9+8hO6dOkSPfbYY/Svf/2L\nFhcX6fe//z393//937UMKxcLa4fZFg/+HelvjanlJYcsiA0ZfWg0ho4Ic71wTflqIE/29kJRQ8MN\nj5n5R0eRKjJKT/q06OfyM6uySpUYrRBb9BA5StYQBVJNL5Ell05Ov5w9qtUnB6V0R2cIiOwcmEM7\nEDy/5p+cXzoOCf2RnBm9QkLR801llPLR0uKVtJbRfNd6Ll1tab8R+t54K0VIv/N8PZ23dA1pyYvs\nvWawLaeRC2/s1A2BNNESjxH1GnzhPRkjsnjIKWjyZznKUrrOmmxeRDSJqCSKIWTy9MhqjxQaI7J7\nznUmagSeR7OKP0MVdkqLXJLMPIrJlWMSYzV6ROVXjesdUR2tqdMzUyOQgCySpoToVZ2FNCTnobkU\nkltemoea6bj0+guVEUmjtFA/kgLl1h28SCA3ivP6ezwjtCK0EZnSPfNksVKdiB67azq2iMDziqg1\n1KxoDc9U4rn52FpRgDWvKI9peO8cnjWjMm8saji9a0mkDypvDmYmIpBeCpG8btqmeVw+tn9emq8i\nxScPmsfwcm7uSbRDztcMLcDmziGNlrT2CM/a9QsJSPQotWnzRl4OspyBB96/xhr1GF1EMGkM7f1K\nIhgi3WNN2mvPSo4/JjlzZcnVGYz2DEQEqaXj+amVU0ntmtXUIgRENk8mK5+T6KXfuRzcAHipBQJJ\nXi8KsdKKdLy0PlIfTS40ErJk1OT0xiAyeZB0AN17vlboNbRWa/BkFNtu9IigoeFGwkxEBBrQXNMb\n6z33PBcy3uKpjfeinGg/hGdOuzUu6kFzZBh671E60egiB4g8NfmOwhBYE9OKMFo4yT95yC2F1jwk\nQ0I66+04a1w6XpORj7XoS8VVTy5rrTRociBjPXk0OlYaJ+mM1F/j782d64W3BzlphAVPv3KMkdVn\nFIYgPYAl+aFEE/EKSLU3AvRKCq0CR4xOzriSOVuGKAepwUHoovm0Ns67KUD61Y4IcujlrFWKURiC\n9CB4lXPvyqfkd2l4Ky00ukHpn9cPpcuNHXIPjvDmY9DIBrmeS/lYkZwmM7reEY/K9SJy7ZgDJALh\n/WYyIohsnrbxnmJbG8ef9UqKhJKaLJ4MOc/T+/+cvNBTUFShPXktQ51bL9FSQv5buknRYEWJ0q0N\nqmspUGfipTZ83yV5089SjCIi0OBFB0ieWMKTP48eFC8lyJED5S215Ra5kHlohio61/4A5Ch4ZD6I\nc8ihzelE+kjGDl2HUoMwFUOQc0CshUstveaRtI21xngRR67306ILHupxj6S1adAUO8ebeHPzDjBa\n/LJoaHvo1Zg8OhJ49GCF4dLaejys/eOODdVpDzNVI5CeaR7fy93SZ9LG9p8IHSRMy0VqBKRD6oWb\n1sHiMkbrIJGDpeXe3voiITE/EJYxjSCSvmkGVNqrSC1D0/u+nzS25lyJqL1Q1NBwI2G0LxQh3gzx\nhtF8WOOL8soNQaUIB/XQkndFeZbQz+GJoiQXL60ReeP4emhrI8mDpAYeLSsaluiVYOqGAMkHS/p4\nfDkNK8+UrrhSGloIl7ZpYSXv56UimoJo6RGinP280NsEjojS8vXIMaxaqI4aF6lIKfHgNwjSemtp\nnTYPa53TdunmRNK3mSwWepCKY1ZfbrX759JvTg+p3koyoJuAFr4kBeMyp8+sAmlEDk8ub7yWu3vG\nKrIfiAwSbc34WvSlg8fbvCJyhGc6RuIrtSM0ohidIUgPtFY46yG1ax6G00VCwfS55glSaAU56bBy\nL6PR1OYlHTjLI6cyaNGBdKg9Y8znIXlZPjc0dLaeaV6Zyybxk6IuK1RHIgytKC2Ns6LOVJb+z9Id\nBIhDHW2xEPV0NaHxrBmCSRFLjXladLgRsKrOXFlL5fT4ReaB9B1SbzzayHqhNIaYx1d7u07+Byel\nCzQNAzM0JOWZhsKVGAsiORWqhTHvuxSJDSGrZQhGlxp4KF2gSSqDFoqV5nNaCKyF4hLQwpbEz6MZ\nhVU4q4VS5yF9l37nYOi5I3RHYQisxa2x0AidIQoyOQWxKF2pzcuva/KrtT9R1J4jH5/WDCwDWxIB\nlaD2uk/NEFjeUqvcSuN4IcQqavXtUmHIKqQhxS1LxpQ+2j/S1yogavyRAqPWL22ziovaWiPr68li\nFWslIGssFT7RCMuD5QiiUUauEZjJYmFDQ0N9jPbNQgnQdQd49ePxKB2DjMsBGukgNIbmgaZVKP3c\nNK0ksvJ0rSZQ/bb419i7FKM0BFIFXOoTKXhJ4zW+nlxcqaMvj1gGjNP03qWwnlnrZ/GIHlhk/tbV\nJv+N0pPSH/SQIbWUSLrijUV4p0h1TdqXqO57/EZpCKz81lPq9DeyqdGaQG+ArDt4S2EiBkyTWZqb\ntUapUiGyosqlrQOnzw2G9AyBVAuRZELWQxuXjkGu9LTrWm1e6d4jNzuRelkJRmkI+AJIlpB7Zkmx\nJcsrKZG0mdZB955HClnci2kH1HoerWSnXqlGEcw69Jo8WpHR4mPRS+lo62FFEP1vaf8tfbTkLOlj\n9ZfWCzEqptxjLRZqShp9Xgu5qcAQcpTMf+h10vhF+Ur9a8s+5Frk0q5ppCXaNAsvFFmHDfWQWn+L\np5fnWYcP4R3x7lwuHvlo/CJGwFs7Sybehq6zl5NbqVqEBpKepZGDtfcR3ULSFoR2GpFYe2/pbQ6m\nFhGUWsxafREPq33Xxk/aO2upkyWn5HmkPj1NScGtuaNrhtKVeGhtKT1EHpS2B5SHxxMxJrmyaRH5\nKFIDvlncWlsFIERJUSvdw/P+lpIh0OSRZO775Siod/hrhKE11laTGzmwkcMn9dV4efqVIzOn7dHk\nclj8cL4jSw2kxZA++TNpYdJ+WrvFXzvkabtEg8uGIpXX6xPhYcmc/tYUyAp5NX5aCieN0/ggYa02\nN2mfJfpcVk0nUv3x+kuI9EPWXaOnnRkrXTAN2RgigoaGhslg3b1ZqI3RfpeiNj2PRw1+Ho3oOpfI\nxItbJXQmNXbIPY+seTRiysEoUgMNJblWZCwi0ySCGCkfReSMzqVvR1IDLp9Gi3+XQnk0NdBoWDJI\ndHibN1YLpXNSF2QtrX3TUhSLniWTh6kZAmtyWt4p9UOVCuGfU/CL8paUJN38HIMmrYOnQFqObcnn\n0fJyd+u3VTOSxkhyRXRKao/WYdLfkaKiJVNu3UlChMYoU4MIkA0uBWpZo4eY086RNUd5h+AR6Vd7\nbIohwvmcomEuctLiKnxbsbBBAnId1TB7mKli4SSRk/954ydRfON0ahfCeJowaQ+VC01etBBbqxg5\nhsJqBDNnCKJFGY2GVpSzaEaLV9K40n5ozQNV9ty0Jx3rrR16e5H+Rgpp0lgtx0YKnSi8PYjQ1IqU\n/B0QVAaNhzfGNARPPvkkzc/P086dO9eeXbp0iZaXl2n79u300EMP0eXLl9faDh06RNu2baMdO3bQ\n8ePHQ0J6z3poLxJFaWgFK63N4yf1l4yFZIT4OA5eR7CKZIiBsl5U4ry06nY61lufaB3Ho8OLk1HD\nFilkcpo8ZYoUoz15hipsIwVI0xA88cQTdOzYsWueHT58mJaXl+m9996jPXv20OHDh4mI6PTp0/T6\n66/T6dOn6dixY/Tss8/S1atXYSG9Z1EMVXiL9NXmFpUtUoQsXbuh9mMIoLcLtfhI9Gsd0ui42jAN\nwX333Uc333zzNc/efPNN2r9/PxER7d+/n9544w0iIjp69Cjt27ePNmzYQIuLi7R161Y6ceJEdYFL\nQ97o1VjaR6LvRTRW2Iv24894W86VX2ltpCaftF9OyIvwy6kRSOuMyjOJGkHNulC4RnDx4kWan58n\nIqL5+Xm6ePEiERGdO3eOFhYW1votLCzQ2bNniwXkKL3Kyg3H+XMr3PI8SK4n68fyqMJKZySgRoPz\nkMZ5tQlp7XLopM+19MdLByN1nUg/qX9uzaYfG404S4uLN5UMnpubozlDAqvtq/a8UFlDzdBq2iHy\nkLwiRkPqZ+0bkrZY/L18tiQtyknnaoXuQ6ZzKA2r3hM2BPPz83ThwgW65ZZb6Pz587Rx40YiItq8\neTOdOXNmrd+HH35ImzdvFmkcPHiQiIgOHCBaWdlNu3fvjorRMEWkhdWGcWNlZYVWVlb8jp2D1dXV\n7q677lr7/fzzz3eHDx/uuq7rDh061L3wwgtd13Xdu+++2y0tLXWff/5598EHH3S33357d/Xq1evo\nASxhpKQiZJG+nHb659Hpn/PP0n4RWHKm8/Boe3T4p9bf+m7J58kk/daeebJ449A+UV3U9sLShxwe\n/zXdch+LwOOPP95t2rSp27BhQ7ewsNC99tpr3UcffdTt2bOn27ZtW7e8vNx9/PHHa/1feumlbsuW\nLd0dd9zRHTt2TBGKXIVC2vgCeofKUwzkYETktoAceD43yRBptLS1QA1PjswSTc1AeHzQwygZnlKj\nqdFCjROqR7lz1GTG6ckdR/WKcXovnOYzUnetmGfRtZ5FC0lDIJWBrwXv0/frn2lrxMdac0bWSWpL\n++SkDMgYziOVO0JH62v9riUf72utn8czKtP/no34P1WGjctbmFnIZVHDN8Z5ILJPC0OtmWVQcw5o\nVLdRYyHTGvG/NUDuZ7m15fDaI/yQcZE7bGmsd1fM29PDltLTZIrKp/XxZJS8coSfNiZy3RblERmD\njkMOv9QXoV/LmJlrOisRQUNDQzlGFxGg3smyyqgnQb2c5c0QL+7Rl+il/S3ZtT7IXKRnXkRhjUl/\n83lI8vJ10+hodCU5opGHt07SXCKIyGPNxeuvPcvRnRSjSA00pHVTpG/6SXR97sb7apAWLEcWr48k\ns0aLf0Zz8tx8OR3Dx/M1kegjc9Ro89/aHlrpiQdNvuhaWfJ4faNrIz3zdMjb/1H+p8o4vLqA1Mcq\nYnn00AOPyomOi/CM1EQk2tIByKn2c++W0tXmhRTGJMXNNSrImL5PNH+P9NPG5EYf3rjI3o4iItDC\n1P7PO/D8ntD4AAAVT0lEQVTWb+SKCUlBkNAx5yCl43JSlv7QIOFlJL2xZEh5c2XjXltKGyz6nIbU\nR0sNtP3RwmZvz7WxHHwPJLqawbTk0uTjV4+onKNMDaQwj9+lSp5ZMgoWLd63JDzWLGzpJqD8rfaI\nx0sVFzkwlgzeQUqV1jJ4Xiqg9eeHCpGTy5M+T+X2ZObyWBElkv5JfNF1kOT0nGKKUaUGWhHFm5xn\nHDQenIZ3mDQrjygflxENDXM9l5ZTW6lSNKSOHFR+yDxaVmrAn3mGJTV+fIyVr6dri+pGNJ3wdF6j\nx/dTm4/EU0LRvz6sDS3M1CZX4t01WJsesfJeuzcP7dDk5M3pnHj+7XkcjY7Gs8ae5BqfEp6l/ND6\nhSWHlWbVgJl2tfcIGiYJxLs2DIeZeo8gt6Bl8YjQs0L3XLmGnI+V26e/tQIbKqNVAEPpIoU1j39E\nTi+14uMiKeUswV2zFhE0NNw4GF1EkItaHtmy/KXWfxLeIxqljMmj5a7xmOaw3jBzhiC3yMOrwVZB\nkBdVckNZSRbtNiRKL1I05f2HSL2iSKveqJHWKuVR2ay9RAxsjfQBSVmsFDUKb+woagSRDdU2UbpG\nQXLn9Dmae0aNkXVFlfZDFMw61FZerskTqTV48kkGM5Kja3y131K9waJnOQepv/a+hXbLY71zIO0h\nYtgkgx9dT2R9RvEegaWY1j1x7vVdSjM9HMj9svXiBh+LyIcYJuRlkQhSpfWurPj7AFKbRktbB+uO\nX+rv0Y5eY1rRlDQn7f0ChL4nj6Vb3hxQHUDOSisWNowGVsrWUAczUSxEQ2OtT47HRMZFUowoP4su\njzTQ3FbjiYSIHg1EXkQmCdGXaGrui0RHWzdtj1A9snRASqm0dY3q+yhTAy8f4s+1cEkKobV+Gu20\nXQvrPCXQZOLtUm6Y9rP4p5/I23g87ZFqC6jyWrUNLmc6F02ZuRwaDYmfpQtaGqXNU+snpR5IyoKE\n9Rb4Xnn6gACRaRQ1AuSQWYUSJIfOLbSk4Lm1lQNbRkmTWzvc3qZr6xc1Dm4e6Rhv7UAj+bo01jJ2\nXvSoFdm0tv477yvpVE7uzuejzSU1nJp+RZxo2mZGja1G0NBw42B0NQI0v/Esfw4/qQ3NsWrloAiv\n9HMoHtE8E6HJaeXuFaIXCB8vNfG+e7IMDStli4xvEUHDTCOSDzfYGF1EUAO1rfQ0rb7Gv6Zn8uiX\n1E5KaHhjorcJ00COhx4TRmkI0FAu4iV4eJR+WhVmTx5PRusZh6TwvKiEFEalZ1p46BULUbn5d219\npWdoMSw3nUTXAtU7CRFdtAzbEOkaQqulBiMFUqEuoa0ZmMi4Ur65NIh0o+XdrkT5z0pagu/hDZga\nlBQatSudmkVFy3Mhr4WW0E69InI1abVJ3p8/s4pytSIw6YrYKypGUi9UpkhhEolCkIjP20NvTUdh\nCKJhXv+HWn8pFLXCRU7be8EjF949uXZIOKQ+fI4SL8SLWKmDxUNaP2kNrZd0tHlphiw1CJwXyp/r\nDPouRwQRg19qfFI+FkZhCKwDICmb9uKHthjSyxvWIUReALIWniurJadGS3rhR/MQ2oH0XiTRwuuo\nh0YOlgYpTenHeXskGeycSIJ/T3mlBsEzzBY9SWbLGWmy8GcIEKM/ihqBJaTWVtNCp/RyZBlCNkSe\nUiC0I5FXbRrIXpR6bZT3kPuQg4g810aIM/6/RW9oaCjHqIqFSDHF6meNtQosUT5eqOfJF2235IiE\noB6iY9D9kp4hqVupPNG98VKdXLnQcUPsZel8pmIIuEFC867+N68jeLkT33gp5OM0eT8p//RycA9o\nsdAKUb3Cp/Qc/c3l0moTyPp7kOom2gH3ioU5cuSuTS7SWkgPKw2SnqEyIf1GVyxMf0tFt/S5V3zh\n/b0+GvimcSPkGTbehii9VPXmRivtKx0Ca86S3FIBVpKZw8v0pMKn1M7lkvpL65G2IXJKemUdQstZ\nSN81OayxvI/kCCSdQ7JspO+6rhHwRYwUV4jGVRziiMxH6p8q9xDztOjm8kRlLpnTUOsxbfxv7UZU\nI5gUtLDLC5W499LSEA1oWFkSbtYwAh6dSK6ryWdFO1p7NGrTUhYNlldGb4VqIZKuDImZiAiGsPBD\nWP5ST5dDYywebIxR1KTWJrp/pXKVnYcZjgiGCPOQzYrC80TWOOl7dGwEtb2NVQOozRsdn7Mfk9i/\nUuM0RNQyKkOAFKciz610wCr4ePmtBatQVBNe0cnjX3rjgcIrNGpj0k/+PHKQLGMg0dd0ZpJheg1E\njc1MpAYNDQ11MLrUINfCRqx01DN710I1eFgoHZezNtaYnKgjSpePS3/XLLJG5S9FhA5aLI3SktZU\nHTfGiAC5JrIq4dJvj1cJhio89qhFm68rss7S+FzeRPaLPsj+eXSQa8Vcvaq5zzUKhkRx2UYXEUjg\nk0MUIdc7pHQ8yyl5uvSTy4laeCSaSOXTrDrq/aNXiHze3jwt/loxMWoELBq5sqWw5IjonBdFIrSs\nWglfB88AehHBqAwBhyQ4v59GlFPabEvRrEXmh0eig26SJDvnxSMDhF9Kz5LRU+iUn8UXUXqtj7fe\n0vPcynw63lqbHNo5fbx+2nx5NEdk7wFyozMqQxAJ7SOHwluEiIL1spXk81oepxmc3IiDG83+e3S+\nHkoORo1cWjrgHo+cq8VaQG8yLD2rvYejMgQp0GscbbFq5uw8NI8eTO2we7mq9jw6t96YSF4ETSvQ\ntpLUhf/WwnYtrZP6pM84fSktlGh58mrPNFh7qHl/Tw4vNfBgGoInn3yS5ufnaefOnWvPDh48SAsL\nC7Rr1y7atWsXvfXWW2tthw4dom3bttGOHTvo+PHjPncQVijXK3nNgy/x53wk7x39jYSFyPNoGGpF\nGxHa1kFMf0cjO299rLX3nkvjo2lCyR5EULI3Uj+rv2kInnjiCTp27Ng1z+bm5ujHP/4xnTx5kk6e\nPEkPP/wwERGdPn2aXn/9dTp9+jQdO3aMnn32Wbp69SomKYPkEbScPhelVr2Eb8TLInJ6HsNqtzwi\nAs+rDgl0HTx4aUKULpo6WnsbLXaW6rNpCO677z66+eabr3suXT8cPXqU9u3bRxs2bKDFxUXaunUr\nnThxQqXtCY6EwFr6gPzmY6OV3ChSr2PJjYSGqOJquXNJoUyCZ6SjB1ZTbmStNHqesZTWKNfAoRGq\nV9BEeaFt1hyyagSvvPIKLS0t0VNPPUWXL18mIqJz587RwsLCWp+FhQU6e/YsLGT6LCfUj4bnvA0J\nN702CzkhrNUnkppY/dF5o8hJnRAayDiLXjTN4GlgzdC/VLeH6BM2BM888wytrq7SqVOnaNOmTfTc\nc8+pfecGjhMjRahoMUfyDpJnQYuXnnxRROeOFuR4H6uwiIa21jiULiqzxjs3XC+NAJF+SDGwlL9H\n66Yog40bN659f/rpp+n73/8+ERFt3ryZzpw5s9b24Ycf0ubNm0UaBw8eXPu+e/dueuCB3ddVtdMw\nXap4E10fYmu/UQXgvDkfjUf63ApPPQ+nySytjWXdOS8p7NXmqsmnRRlamqGtmURXk02TwdIVa4wH\ni5a3VpbcFh9v/TV66N6trKzQAw+s0IEDzjnoHKyurnZ33XXX2u9z586tff/5z3/e7du3r+u6rnv3\n3Xe7paWl7vPPP+8++OCD7vbbb++uXr16HT2LZd8kddHa0udpmz8zua9Fg/NC5JWeW23ad+nPkxOZ\nhyWrRcOTOZWbf/dks9ZIo12Dh7bmFm0+FtVhlJY1p6i+/9c8iG1mRLBv3z56++236T//+Q/deuut\n9OKLL9LKygqdOnWK5ubm6LbbbqNf//rXRER055130mOPPUZ33nkn3XTTTfSrX/0KTg24l9O8D++v\n9bWee309L9g/92SI8Ef4EenRkORxLfn4c1RWJKqR5ODfLVmtvr0MEm2PRyo/sq+cliUTl8GSW+pj\neXfE63v9EIz6Hx35tOyF14yHNhYN/1AZtRSHKF9ur0+EpxT+av3QdbH6SLQt+T3k7GGufH0buj8a\nLUtWiw8yFtETUv4HJ6N4s9ArWGnPpUVEwcemhwahhSibFilYckven/+WvJXUhnh8PjYa3UQMojY3\njZe1BlLNSJMlZz81PUAiIivqSOWRaEl7lls0jBjEURgCDivkQSy3R5OP4RuA8PQ2xzrgVpjIn0th\nr8Q7J/SMKApqICV4h0daVysN0YqdSGqAyimlE966a9/731qapPW1dN1yBjmY6dSgNo3SsBTpg4aN\nubxL165muD00JiWrlbrUTk+GXueZ+O8R9Ih6qvST07BSC95uHS7+HQlhJRmkEFmiL3lgxKMgIWkk\ncvKMTvo9J2LgMkXoIHP12qJpg+SJNd1DEXUOyD5GMcqIoKGhYRjMVETQ0DB21PTGY8AoDYFVLUbH\nloSptfpqxaWS+Xn9vYIWQrtUvhKU7sEQskppi5T2lYbs0zQu6z41mIXil1aMQgpTY4FWO5mmvEPx\nt/Ylp4gclVO7NcHGzkBqgHg6qUCVWmr+iVybIZace0qLZy64QqV8NF4oT41GtHBm0bUKaeja8vmW\nyIYcEquYrLVJ15ole48Yjkh/lM41bWOLCGp58JpjJAtOFLvCS614OraWnF51H+3PD3WJnMi4Eu84\n6WjP44esl9eHe/sae3Atffn8jcoQ5CqzpLyR8KlUyWspX/RwWmMQeesolmzgtH4RmqUooePNR+s7\nhL5pe55Dh2bhFWOrAKPdyfPF6L/zTyvU7OlEUgOJhzQm8tzaVO29hYgiaGtjyWWtC1dIKURGDpS0\nrsheeGlS/JDI4630KZ07IjuXS5oH/66lIppO5GAUhsDKubzN1BYDVYqIokq8tMOj0eQyIhuohZ4o\nHWldkNyYezpJrhr5ayTt6ts0448YX2StJOcgeWBNdtQJSEYE1YlchyDKNabUoKGhYVjMxK3BNFAa\nUtWiP7QcDZPHLO3pDW8Ihg5OUPotSFp/QOoFY8EoDEFksbQcahILPoZNnWaE4dHM3Zdp7ScCtLag\nQStWazxQ1B4zCkMQrXwjVfpcRUIq5CV00rYSxfJ4DRFhRAq3ETmkKrpXvI3eyOSgJi2uO0ihz5on\ncuPjjUkxCkOQInfxtUoypxu52knHeleFHh2Jf63QEbluQ4FWrYeCZeh5n6Hl9HSKw3NE6PsW/Sey\nFqhsHqZya0DUuQeu/52CvwzDkfuyDL/28a7reHuUJ1cGabyl5B5PrU2SN5UjB5Jh8/ZQo2FdiWlz\nSb9r64ry42MRWhI91GBEDIu3xzjfEd0aRBdKsoz9s7QtsgnoXbLWLkUJqRxe6CopLvc+/I45HWfd\nnUvyaPOp4VklOXu+nlfT+ki/+d55+6RFZWlEIemRx1/ig0SMQ0Qwtdz4un2PAMkzNSuL0CbK24SI\np7FooHJa3sLyskPIM4mxQ8whut41eWoRaj7PEUUEKDwLauVkUSMgjUk9T/qXehRLLi190TbWomeF\niRo8L+V5YYk3b0vHIX2jbSjfHkhqpe2L1+YB4c37p2smzUOLUqw10GCuyXqNCBoaGq7HqCOCqOeX\nLCjqiSyv5cEbi9QfJFrRcR4dhBbK06KpRUySbCgP6xnSVjpGm2P63eqj8UKjPasvHxOd06gjghq5\n2VB9eX5GFMupvdwczdFr5fK5QOXKoRGlI40bej2sPL00b/dqADky2v1GGhFYyuE9z/HuqfJ4PKX8\nK5Kre4cnmnvy8ZwugpyIxqsxlHpvKfeN7A9SZ/Fki+T12u/UWaARKq8BoN6+dlQ0dUPAkS5m7kHL\n4dfT8zYhcqsQadc8i/abf+fKp8mAFMU8A6iNs+D1GTLc1/pa14X9GCkq1JC7Dp5sCEojotEZAuQw\n9pAOT9R4SN52yLvgiHflhztyjYRGPUgO7CG3DpAi5xBFlD/HEHm3QxoNNEKS0s3+d431iuzjqAxB\nJBro+0nj+fcIX8kIRUNQNDREPIFmDLxx6TPJuFjXU97VVSo/kpt64b5knDxj7DkLNAXih1FKPaMp\nnKa/aeSm0U1pWPO0DDqnhWBUhqCHtcnW5Gp5CGuDiMqjglzU5psbAfTIqXEg7V7dqDQVsQycVRyu\nnd6UprE1eU3FEEQKSJqyptbS87KR0NVSQh6xWAfJ8txI2M7nlnoTpHgozSlnbUrDYi/MtXJmLTLT\n9gr13lIEWBJhpLwjcy2pAyAyR3hN/fqwoaFhchjt9WFDQ4OPodPRdWsIhly4SEjmjc2hM3QuOkmU\nyJabBg6FIW9bIrcXOZgJQ5CzAF5eXFIo024r0NpHuqna/b9WTedjkNsJ9B48CqmmobVr4Dc10bVE\n5YvKFaGbe1iRfZGKlaVFXgmjMASWZZeKY1KBx9scTk96ecby9LzNu1nQ5iX1tSrTnCcvlEnKhxqJ\niEey2tM1leSNHGrUyEl9uK5oa4kcWE3PJLrR9zc8vfWMorXn2lzcPrNQLESrwVJf6XcPqQqfg7Ri\njFT2pf5cHn7gOR2rnyYj55EzR4Q+OkaiX7If3lpIPCx+3r5JfHPmgaytpcMRzOT/+7DGIY0gohST\n4l96MCIKNgSiBquURy2HMSSi8iJ6idOaAUPQ0NAQR8wAzsD1oZcbSX37795YNIdC8mrez8sJozy8\n/trcPZml55r8nnzas5x58/7RvUJ5WOtmrTMqu8QHld2SyVp3HnHlyEXUIoKGhnUDJDKYiYigoaEh\nHyX+dWqGIBJqNQyHSezDGPZ6DDJMArnznJohmER2MI3NnzWFm8Q+TDMTjN5WzNr+ceS+PDbTEYFH\nA7n3jtBDMGvlD6/omUujBt0aQA+GZTBm0aFE33JsxcKGhnWM69+XGFGxcGVlBe6LWvJSoFdWGK2V\nYnkiGNJjSbT7/asVVVnXdlr/ksjDGmPNbSjUWDMNaGQwekPgBQ/R8B+lI41F350/cGAFkqmWsiGv\nF0feq/BoW/uHhOLeGvB/B6Lx4d6O/7sRjafEq8cDD6y4vGsjR4cj8iHGYBTXh9qLFtILFfy59nKM\ndqfqvcjRf5f+DYD2m+PFF6/n0cN63VbLtTUPyOXXDlj6D4Kk+VmQ9saSmc9Nkls7vB4/zovvQ8ob\nNY6c5oED17dLfSV6Gl1Jdm2elm5r9LQ+Kbz9HoUhSD0A/843lT/n3sPzJlK7xFfrI/Xj6JXJk8F6\nLq2DNlZbJ68f4lWs+Wp7o423eEj7Z9Hy9kubRzrGWjNJzxD6Wl9E3yJ8LfqejErjZHH//fd3RNT+\n2l/7m8Lf/fffL57Lid8aNDQ0jA+jSA0aGhqmi2YIGhoamiFoaGhohqChoYGaIWhoaCCi/w+SKDjO\nMdWChAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xa9e9d50>"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L = la.cholesky(A)\n",
      "print \"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0])\n",
      "pt.spy(L, marker=\",\", precision=prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1865 non-zeros\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "<matplotlib.lines.Line2D at 0xb882590>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD7CAYAAACBpZo1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG29JREFUeJztnU9sHkf5x78v2LdSESHy1vIb6YXEqQl1XUuBW1SjxEU9\nNLSXUBeBRdJLOVJB21McJIh7gEuhEoqKlBO0F5IcSGRxeEvVi5FoqqoWSlFdSJzEog2RqlYiarq/\nA7/37XoyM/vM7Ozu7O73I0V5ve/s/Nk/35l5nmfm7SRJkoAQ0mo+V3UFCCHVQyEghFAICCEUAkII\nKASEEFAICCGoQAguXLiA6elpTE1N4fnnny+7+OD0+33cf//9mJubwze/+U0AwI0bN7CwsIC9e/fi\noYcews2bNyuupZyjR4+i2+1iZmZmdMzWnpMnT2JqagrT09NYXV2tospO6Nq3vLyMXq+Hubk5zM3N\n4fz586Pv6tY+b5IS+eSTT5Ldu3cnGxsbya1bt5LZ2dlkfX29zCoEp9/vJx988MG2Yz/5yU+S559/\nPkmSJFlZWUmeeeaZKqrmxV/+8pfkb3/7W3LfffeNjpna8/bbbyezs7PJrVu3ko2NjWT37t3J7du3\nK6m3FF37lpeXk1/+8pd3pK1j+3wpdUSwtraGPXv2oN/vY3x8HI8//jjOnj1bZhUKIVFiss6dO4el\npSUAwNLSEs6cOVNFtbw4cOAAduzYse2YqT1nz57F4uIixsfH0e/3sWfPHqytrZVeZxd07QPuvIdA\nPdvnS6lCsLm5iV27do3+7vV62NzcLLMKwel0Ojh06BD279+PU6dOAQC2trbQ7XYBAN1uF1tbW1VW\nMTem9ly9ehW9Xm+Urs7384UXXsDs7CyOHTs2mvo0qX1ZlCoEnU6nzOJK4fXXX8cbb7yB8+fP4ze/\n+Q1ee+21bd93Op1GtTurPXVs61NPPYWNjQ1cvHgRExMTePrpp41p69g+CaUKweTkJC5fvjz6+/Ll\ny9sUt45MTEwAAL785S/jsccew9raGrrdLq5fvw4AuHbtGnbu3FllFXNjao96P69cuYLJyclK6piH\nnTt3jgTuySefHA3/m9I+CaUKwf79+/HOO+/gvffew61bt/Dyyy/j8OHDZVYhKB9//DE+/PBDAMBH\nH32E1dVVzMzM4PDhwzh9+jQA4PTp03j00UerrGZuTO05fPgw/vCHP+DWrVvY2NjAO++8M/Kc1Ilr\n166NPv/xj38ceRSa0j4RZVsn//SnPyV79+5Ndu/enfziF78ou/igvPvuu8ns7GwyOzubfP3rXx+1\n54MPPkgOHjyYTE1NJQsLC8l//vOfimsq5/HHH08mJiaS8fHxpNfrJb/73e+s7fn5z3+e7N69O7n3\n3nuTCxcuVFhzGWr7XnrppeT73/9+MjMzk9x///3Jd77zneT69euj9HVrny+dJOEyZELaDiMLCSHh\nhaBpkYOEtIGgU4Pbt2/j3nvvxZ///GdMTk7iG9/4Bn7/+9/ja1/7WqgiCCEFEHRE0NTIQUKaTlAh\naGLkICFtYCxkZpKoqwceeABvvvlmyGIJIUJmZ2dx8eLFO44HHRFIIgfffPNNHD9+/P8XeSRIkub9\nG7avqf+A5rav6ffO1AkHFQLXyMGEEQy15PjxqmtAQhN0ajA2NoZf//rX+Pa3v43bt2/j2LFjmR6D\nToeCEDu8R80nqBAAwMMPP4yHH37YmmZ+fn70OUncH7TYH8x0+5qAeq3LaF9V97gu9y709Sk9xLjT\n6cBUZJ7GxS4OJC7a+ryY3r+oQoyHowPfcyXEtpw8tvq0hTaKgI2ohADIJwY2hnnG9gDEVh9yJ3UQ\n67x1jE4IgGLEgC9cNYS8j6GfCWl+dXh28tYxSiEAihsZkHIZPqAh7mVQ41hDbASh3pFohQCgGDSJ\n2F662OrjSuipbtRCALiLQVnCEfOQNw8x1YWYCS1k0QsB4NbospQ+ZDkx9U4x1YWURy2EAGBPRYqj\n7Gcrxme5NkIQyugU403Ii65NTWxnUZQ9Csoqr6h7Z8u3NkIwJK8BsYlDX12b6t7ONguZ5N75XB9b\nvrUTAoDehDbgI2RteiZaaSzUQTGIh/R9qGJYO6TKUVDdn8XaCgFQrRi0rVwb6RewqJexqHxDXc+6\nT8VqLQRAdWKQ98YXvbiKyKjD9Szj+a69EAB6MYgxCCmN6/4LIYlxVFE2kmtQ5TQnTSixapTXwIQq\nBjEGIfkSun6xt7cMJNdATRPKdV3V9W+c18AEH3BSJHmfr5ifz0YJAVDvtQYxDNlpBG0njRMCnwjE\n0MEZruUN8zLlWaaNIMZhaxFUbSOQ5h0ynS1NVHsWhi8r7uEYKZ6mbYybl1rsWRiaugQdlRGQ00Z8\nXuomi4CNRgsBUA8xsAXkxF73mGnrS+1D44UAiEMMqg4gqrr9dcB0jVzm9FXRyM1Li6BqMai6d6q6\n/DpgukZJEv/1a+zmpUVQtRjoCFGfJuxHULf6Fgn3IyiB2JQ9RH2asB9B3epbJFUs3GqdEADsfUKR\nZ30HiYtWCoFpihDDg+xTh6rqrfYw7NXrSyuFANBHIMawsUW6DtZIsJLrHYNIkuJorRAMqcKAqCvP\nNs+Xps9TfhZN6+0lI0KGGBdImSHGLjQ9tLSOSO5JW++bb7tbGWLsQiyuxSptBL75FHXdfPYM8EU6\nDSuLrDK5eWmBxCAGVfZuvmVXUefQ98nWhiLaV/aLngWFQCEGMSiC2NqUtz51nw7EVn8KgYa6iUGZ\nQ+hQxFafIXW67yGhEBiI9UElxdLW+04hsBBL7xBLPUhzoRBYkGx7lidCUZquCmOVmoZi5EbdrheF\nQEDaZmDabzDrmCnfqnC1K/jWNaaNWcukblMMCoGQoRjU7QZXTdbGrEDzRKKM9oSON6EQOBDam9C0\nF0CHZETguj1b7NetjCAnnzK4DDkgIcWgrNFFDDszhfzlqbaMyspsJ4XAg7xiYFrcUtRe+1X+ziLL\nrQcUAk/yiIFpHb/tB050y5SLoC0/cFJ1ubFBIchB3g1OXNyHoR9Yug+LpW7Xi0KQE7oPiY66Xa8x\n3xP7/T7uvvtufP7zn8f4+DjW1tZw48YNfPe738U///lP9Pt9vPLKK/jiF78Ysr5RUpRbke5KApTz\nHHiPCDqdDgaDAd544w2sra0BAFZWVrCwsIBLly7h4MGDWFlZCVbRmPH54dU0pvOqmA6Q+CijM8g1\nNVB3Ojl37hyWlpYAAEtLSzhz5kye7GuHrwGxrF6/jaMLip+MXCOCQ4cOYf/+/Th16hQAYGtrC91u\nFwDQ7XaxtbUVppY1wkUMqtxZSGosrHuIsE78Yo38qxJvG8Hrr7+OiYkJ/Pvf/8bCwgKmp6e3fd/p\ndNAxXI3l5eXR5/n5eczPz/tWI0qk4ciheujQUWa6NE0aTZTRlhBlhLANDAYDDAaD7LJCbF564sQJ\n3HXXXTh16hQGgwHuueceXLt2Dd/61rfw97//fXuBkW5eWgQ0IrYT/41Fi7+vQTcv/fjjj/Hhhx8C\nAD766COsrq5iZmYGhw8fxunTpwEAp0+fxqOPPpqjyvWnqGkCRSBufI3Hlf6uhs+IYGNjA4899hgA\n4JNPPsH3vvc9PPfcc7hx4waOHDmCf/3rX0b3YZtGBEPYg8dFm++H6f3j7xoQ0iL4uwYVYxomuvzq\nTRHlk8+o+h7lIW/9KAQlkZ43pm9aEesIVJo6FC7rtw2k98i1PiHT572/FIKSsfm1bT1S+nPeFyC2\nX/UJWa6rcda2sKqoa+FyD0Oms9532giqIVQv3dTennxGyHtMG0FkhNrpyDefEJGFLr2QZLQjzTd3\n71cTm0xZe1AAFIJKGb7ERTyYWcNcaWShbdMUl+hE2/xbrZ80ItN3zmzbAMZ0PFRMgIuIhd5ZypaG\nQlAx6ZsteeDUNKZho5pv+v+s/H2JwbLuupmKZOQgfSFd2m+717p755p/Vr531Ic2gjjQvdC2uWHZ\ntgHaItyJ8ZrRRhA5ummCycOQfsBC98KSXsp2rostwOV76TkxjEqG5JkelA1HBIS0CI4IaoJPkEko\nX7Ov4dLFBpGVR0hc5+wh8glF0facO/LiiCBOinQdxTh3DUFWu5rabhc4IqgZWQ9sUbEDRZwbQ/5V\nUHS0ZMhrRiGIGFuwUJ5diULvaBQqMKpoTG65MsorIp+Q14xCEDk6b0LexSqhX4KqhtuqbSJvtGab\n1xpQCGqA+qK5vnhNnRfn+XEZsh0KQU1Qe7s8vVXML4tPu0JNdfIKrkvZ0rylUzLXzWhVKAQ1Ij1N\ncAkoCuVGK/LcIXnj6/NOm3zzKsrtaxN/V5uHta10HxLSHug+bBB5F9a49oKhDVuhCFFX27V0TSMN\nrBqmldwX6b3miKCluAbHqJZ111tQdDDOMP885cQWUFTWNXM7hyOCRmHaAzErvfpZpapFRy7lmPLM\nMqjatokrglCLjsoYaVEIao70YVN7D9Mw1LQ/gnRYKg1+yRsZaeoNXTcZlRhdQ0dxZsU7qIIabpsy\ny3ecGjQDxtFXj+rN8bkfxU8nODVoNCHjDIgf0umXNI8yoRA0iLQYcHRQDE0VWApBwwi1O3LbMbn2\nmiqwFIIGohMDX3EocuGNazmhFl65ell0f/viEpsQGhoLCSEAaCxsJZIIN/I/6notQtWbQtBgTD5y\n23w3xEIeF0IvcvKdMtiuhXQYX0WIcajFZxSCFqBubuLyWwnp0F8TPi+fWhfXwKW8EZC69LZIRNO1\n07kMpWHOprxM9QsRlGWCQtASJEtWdS+8yzr3MofXujX9WT2sr+XfVVwkedny9Hnh8+5HQGNhy2ja\nwhviBo2FBIB/nIFvby8VgTx1Sv8vmVqY0uVxabrUP6ue0rSuLlHaCMg2pGIgfbHyoBuuSx789Dnp\n/20G0qwQYPW6uAy3XYbzrluVSRZycasy4oVEDNIPoe9KP0k9pPlmlWd7ibNemtCRgyGnQ9JRTh4o\nBC2miIVKEpdbUfjssyA5t2okhsW8BlsKQcsJsWJOl1/ofPMSQx3KwLedFAJSWo9d1+i9NkAhILmH\nlVnRdaHm30UuOiraKOqap68XwPcaUQjICDUC0eU8wB4RF+LBl+RpehFMU5W84da+7sMsJOKcZWh1\n8WRQCMg28vbavhb/IsjrUnOxd7jEZ7i6D0NB9yFxwjfoKM9D63NuUS5N3zrU2SBJISBa1GmC67xa\niku+LsNwW5RgGfaAomwARUEhIEayerg8PWB6bu660EkaDOVSj5gIHdwkEVCrEBw9ehTdbhczMzOj\nYzdu3MDCwgL27t2Lhx56CDdv3hx9d/LkSUxNTWF6ehqrq6v5ak+iYPjShR72po1hPvNlibHQdq7u\ns/q3JDjK5DHxWbGYFW7tO5KRrHi0CsEPf/hDXLhwYduxlZUVLCws4NKlSzh48CBWVlYAAOvr63j5\n5Zexvr6OCxcu4Ec/+hE+/fRTW/akJhQ59w1tVyjCuCYVDskLl05ry8uWv62uvliF4MCBA9ixY8e2\nY+fOncPS0hIAYGlpCWfOnAEAnD17FouLixgfH0e/38eePXuwtraWr3YkGkLaCNReLssWIS3Plk+o\nMqTn+bhhbXlJ0knz0+FsI9ja2kK32wUAdLtdbG1tAQCuXr2KXq83Stfr9bC5uemaPYkUn7m8La/0\n56zFTep5tuG/KR/fntXFJWiqiyuSkYLtuE+6XMbCTqeDjuVK2b4j9SNvBGLoepRxXp1dgi6MuZ7Q\n7XZx/fp13HPPPbh27Rp27twJAJicnMTly5dH6a5cuYLJyUltHsvLy6PP8/PzmJ+fd60GqZCiDIgk\nPIPBAIPBIDNd5lZl7733Hh555BG89dZbAICf/vSn+NKXvoRnnnkGKysruHnzJlZWVrC+vo4nnngC\na2tr2NzcxKFDh/CPf/zjjlEBtyprDqoY6MRheMwkHCZLedaSYsneBKY8beeHEjgXz0HRddmep/79\ns44IFhcX8eqrr+L999/Hrl278LOf/QzPPvssjhw5gpdeegn9fh+vvPIKAGDfvn04cuQI9u3bh7Gx\nMbz44oucGjQc9QW3Wb1d5rkh4hfy2gjy+vKlQlakJ8AFbl5KcqPrbaU9sEsZgD1Pae/qOgrxGYFk\nUcSoR/K96f2jEJAgxGgzKKtOoUWvSLiLMSkUdW0CIIuKU+MIXLDFNkhjBfLOXtMjIJ+IQp+2u6y5\nEOfJEQEh7YEjAlIKpvj8rB5aEjUoLVvynTR6L0TZRednOs8lDwoBCYrPYhvpeSHKDllek6AQkEJI\nhySr4bZlRfilBcfHJeniQvQVIdNoKMsbIMnbpU4UAlIYaTFwWTiTlVY6tLa5AyX55Bk1uAzLVWHI\nysO3XrY6UQhIoegWK0l7avWlcF3mm1WmpGdWz5XUNauOtvbndUPahMSWH4WAFI4tys/WO5teZomh\nT2KcVF/IvKv8pOXr6qKepxNQXd7qNbFNDazXhO5DUhZVB9tUXX6RyGMY6D4kFWPr5coqv6nkbRuF\ngJSKJALRZQqgw3WqoCtbctynbr7ps/IwxW+I8+LUgJD2wKkBiQqf3l7a2/n23nk9By7kjWQMXSaF\ngFRC2k2mGrqy/OfphT4qpu8kL5fNS6EuYdbl6zos1wmPWlbR4jeEQkAqxfSCqcd0L416ro9XIH2u\ner5k0xSd31/yUqbdlbq4Aps7UGpHcQnAohCQyjG9XFk+cZt42EKGdel84whM8Q/SqEZTPU1Ckc4/\nK/ZC0rZRmTQWklhosp8/FmgsJNGTFWfgE+qbByf3WwXxESHLpBCQqLAtVJIY7KRGPV2YsS7s2Je8\nIcZqPXXorpUkPkNbD04NSIzEMk0IVQ9bPr5GTj/DKKcGpEbEIAJAuHrY8ilrfwauPiS1RBLcIw0A\nkvj7q5jn27C5//JcG+15nBqQ2IllmtAE6DUgtUX1JviEAWdFGob0VrhE+WXFQ4SqU9Y5FAJSC9Q9\nEE3BNKawXfVvUxhyiGlDlhtUbYutDtIybMFFItHg1IDUiTzTBE4xODUgDUG6oMjXWOgynJbkL0mb\nVbeskUtWmVnxCABHBKSmuPTuatoYRgZ54woYR0AIzBGIumOmhT223ljNy9Y7Z6XRlWvLy7Qk29am\nvMZFjghIrZH2jKZ06eMxjBR06OrlOyKijYA0EunLIFliHKMIALJ9EVzPV6EQkNrjGhEo9d0XFWlo\nmh7opirSIb9PXEEaCgGpPS47A9m+y9pMxDdvG+qIJB0LkbV7U4igoyEUAtIYbIE86ouuBh9JxCTL\nbaczHOry0H1W/zZtUaYeky5morGQtA6duxCQueRM6WI1JLpCYyFpDerIwOZqkw6f/Xz2snQmu4Ap\nEMgn8IkBRaS1lNWL53VBup4zTO9XFkcEpGVkLf7R4WP0y+uCdD1HYh9wbQeFgDQa29oE13h9G6EE\nR+INkJRnsnMY8+PUgJD2wKkBaTUmA2GWYU5NL8k/6xydQdD0vSmvvO5JFQoBaQWmjTxs83s1sMe1\nPMn3pvBh6foJdZrja6+gEJDWYJs323Y9Mp1jwzc0OGu0oJ7jKhgmKASkddii9iQvvuTl81kUpIqS\nZKGURJQkIwwKAWklWUFHOmyBPC4brJpiDrI+q/WVTiEkUAhIa5Fue6ZiGq5L1iuoG5dK1iXoRisu\no5TcU4OjR4+i2+1iZmZmdGx5eRm9Xg9zc3OYm5vD+fPnR9+dPHkSU1NTmJ6exurqqr1kQiJAfTF1\nPXH6/6w0unxM5drKNJWnS5eFZPRgjSN47bXXcNddd+EHP/gB3nrrLQDAiRMn8IUvfAE//vGPt6Vd\nX1/HE088gb/+9a/Y3NzEoUOHcOnSJXzuc9u1hnEEJEZsc3KJ58AnTBiQlanLW50mmNLcecwjjuDA\ngQPYsWPHHcd1GZ09exaLi4sYHx9Hv9/Hnj17sLa2ZsuekGjQDdNdFidJRSBdhlpmOpYha9owPN8m\nEOqx4HEEL7zwAmZnZ3Hs2DHcvHkTAHD16lX0er1Rml6vh83NTZ/sCakE27A/1CBWHabrPtumImrd\ndHU1HbO1wVkInnrqKWxsbODixYuYmJjA008/bUzb8VnBQUiFSNYlZPn2pesFdGl1QUI+0Y6mkYaJ\nMVl1P2Pnzp2jz08++SQeeeQRAMDk5CQuX748+u7KlSuYnJzU5rG8vDz6PD8/j/n5eddqEFIIkriC\nrGMu9gRTz57+nCeycTAY4Pjxwej4iRPGE+xsbGwk99133+jvq1evjj7/6le/ShYXF5MkSZK33347\nmZ2dTf773/8m7777bvLVr341+fTTT+/IT1AkIVEwfFS3D9I/O65Lo37W5WnKT/3f9NlWj6z2mN4/\n64hgcXERr776Kt5//33s2rULJ06cwGAwwMWLF9HpdPCVr3wFv/3tbwEA+/btw5EjR7Bv3z6MjY3h\nxRdf5NSA1BrVGq+zwtvm5mmG55q+k+SV/t7XS2H8/n9KUR50H5K6YQveSVv4h99lhSub0mTlZQpa\nCuE+dLYRENI2bD2wZK6vYhKKrLxcytB9z2XIhOTE5NeXxh3YzpV4CUx5hoJCQIiQPDNaNT4gK19J\nAJNrfYLGERDSZtLz9uHf6mfbC2eyAajnp6cjpqmEdPWkJI6AxkJCWgT3LCQkIE3zjFMICPFA58ZT\nQ3p1pIfpNmOhLbTZJ+Q4C7oPCclBei4vCfSRuCCzAoqk+brAEQEhObEZAHVIFihJFjm5lGPKZwiF\ngJAASMVA3WtAdQOqHgJdNGF6tYF6nlontVwTFAJCApEVXpxOY4oDsMUY6OwSpvNcv6cQEBIQUwSi\nKY2LYdA24rBNL7IiFgEaCwkJjvRF1SFN52qY5OpDQsgIBhQRUjK+vv0qgpUoBIQUhMSVaDqvbCgE\nhBSIzdI/xEcsQp9DISCkBFxGB5LVgj6jBsYREBIBqhiYtivLSu8L3YeERELWnoTDz2nUACTdHopZ\nrsQsOCIgpGRsL3RaEEyxCNKFSC5QCAipAHU3IsAcZWg6P8vuoE4rbKLBqQEhFaGbDtgWIEk8EGr+\n6fS0ERASKbZhvm3hkuSY7juTGHBqQEjF+Gxr7uJBkLgjKQSEVEzWD5qowqDr4bNiDrKMipwaEBIJ\nui3MbWlDQiEgJCJ0IuD7ewYucGpASGSoBkSJZ0GHix2BQkBIhNi2S0+nCbUxCYWAkEjJGvaH3LeA\nQkBIxGRtXioRA8k6BBoLCYkcSaCRZIpg80ZwREBIDSh6+zIKASE1wGWBkc/3nBoQUhN0NgKdS9F0\nLkOMCWkQriHGkmkFhYCQGiIdCUjTUggIqSmudgNODQhpKLb9CqX7GwAUAkJqj2S1IvcjIKQFZIlB\n1roECgEhDYHbmRNCAPj/IAqFgJAGkQ464jJkQlqOajPImjZQCAhpKC4/vFqJEAwGgyqKLQ22r740\nrW22n1BLQyEoALavvjSxbZIfSeXUgJAWwNWHhBAAGQbDpGQefPDBBAD/8R//VfDvwQcf1L6XnSQJ\n/ZsphJC6wakBIYRCQAihEBBCQCEghIBCQAgB8H9t8bZuMwkIewAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb6b1990>"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cholesky is often less bad, but in principle affected the same way."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Reducing the fill-in"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def degree(mat, row):\n",
      "    return len(np.where(mat[row])[0])\n",
      "\n",
      "print degree(A, 3)\n",
      "print degree(A, 4)\n",
      "print degree(A, 5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2\n",
        "4\n",
        "3\n"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def argmin2(iterable, return_value=False):\n",
      "    it = iter(iterable)\n",
      "    try:\n",
      "        current_argmin, current_min = it.next()\n",
      "    except StopIteration:\n",
      "        raise ValueError(\"argmin of empty iterable\")\n",
      "\n",
      "    for arg, item in it:\n",
      "        if item < current_min:\n",
      "            current_argmin = arg\n",
      "            current_min = item\n",
      "\n",
      "    if return_value:\n",
      "        return current_argmin, current_min\n",
      "    else:\n",
      "        return current_argmin\n",
      "\n",
      "def argmin(iterable):\n",
      "    return argmin2(enumerate(iterable))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def cuthill_mckee(mat):\n",
      "    \"\"\"Return a Cuthill-McKee ordering for the given matrix.\n",
      "\n",
      "    See (for example)\n",
      "    Y. Saad, Iterative Methods for Sparse Linear System,\n",
      "    2nd edition, p. 76.\n",
      "    \"\"\"\n",
      "\n",
      "    # this list is called \"old_numbers\" because it maps a\n",
      "    # \"new number to its \"old number\"\n",
      "    old_numbers = []\n",
      "    visited_nodes = set()\n",
      "    levelset = []\n",
      "\n",
      "    all_nodes = set(xrange(len(mat)))\n",
      "\n",
      "    while len(old_numbers) < len(mat):\n",
      "        if not levelset:\n",
      "            unvisited = list(all_nodes - visited_nodes)\n",
      "\n",
      "            if not unvisited:\n",
      "                break\n",
      "\n",
      "            start_node = unvisited[\n",
      "                    argmin(degree(mat, node) for node in unvisited)]\n",
      "            visited_nodes.add(start_node)\n",
      "            old_numbers.append(start_node)\n",
      "            levelset = [start_node]\n",
      "\n",
      "        next_levelset = set()\n",
      "        levelset.sort(key=lambda row: degree(mat, row))\n",
      "        #print levelset\n",
      "        \n",
      "        for node in levelset:\n",
      "            row = mat[node]\n",
      "            neighbors, = np.where(row)\n",
      "            \n",
      "            for neighbor in neighbors:\n",
      "                if neighbor in visited_nodes:\n",
      "                    continue\n",
      "\n",
      "                visited_nodes.add(neighbor)\n",
      "                next_levelset.add(neighbor)\n",
      "                old_numbers.append(neighbor)\n",
      "\n",
      "        levelset = list(next_levelset)\n",
      "\n",
      "    return np.array(old_numbers, dtype=np.intp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cmk = cuthill_mckee(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P = np.eye(len(A))[cmk[::-1]]\n",
      "\n",
      "A_reordered = np.dot(np.dot(P, A), P.T)\n",
      "pt.spy(A_reordered, marker=\",\", precision=prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "<matplotlib.lines.Line2D at 0xc0ef750>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD7CAYAAACBpZo1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFv9JREFUeJztnU1sVNUbxp+r7U6NxNihmSEZ0xZqtdQm6I4whhbDggqb\nSjHaCGxwKVF0ZTFRysK4QEkMwaQbFTZSFtI0LgYJmyGREmJjgFi0lNIItQmBRQPe/0Jn/tNh5s79\nON/z/JIJZebOPe+5c97nvOfb833fByGkoXlMtwGEEP1QCAghFAJCCIWAEAIKASEEFAJCCDQIwcTE\nBDo7O9HR0YHDhw+rTl442WwW69evR29vL1555RUAwOLiIvr7+7F27Vps2bIFS0tLmq0Mz+7du5FK\npdDd3V16Lyg/hw4dQkdHBzo7OzE5OanD5EhUy9/IyAgymQx6e3vR29uLM2fOlD6zLX+x8RXy4MED\nv62tzZ+ZmfGXl5f9np4ef3p6WqUJwslms/6dO3dWvPf+++/7hw8f9n3f90dHR/0DBw7oMC0WP//8\ns//LL7/4L774Yum9Wvn59ddf/Z6eHn95edmfmZnx29ra/IcPH2qxOyzV8jcyMuJ//vnnj1xrY/7i\nojQiKBQKaG9vRzabRXNzM3bu3Inx8XGVJkjBr5iTdfr0aQwPDwMAhoeHcerUKR1mxWLjxo1YtWrV\nivdq5Wd8fBxDQ0Nobm5GNptFe3s7CoWCcpujUC1/wKO/IWBn/uKiVAjm5uawZs2a0v8zmQzm5uZU\nmiAcz/PQ19eHDRs24NixYwCAhYUFpFIpAEAqlcLCwoJOExNTKz83b95EJpMpXWfz73nkyBH09PRg\nz549paaPS/mrh1Ih8DxPZXJKOH/+PC5evIgzZ87gq6++wrlz51Z87nmeU/mulx8b87pv3z7MzMxg\namoKra2t2L9/f81rbcxfGJQKQTqdxuzsbOn/s7OzKxTXRlpbWwEAzz77LHbs2IFCoYBUKoVbt24B\nAObn59HS0qLTxMTUyk/l73njxg2k02ktNiahpaWlJHB79+4thf+u5C8MSoVgw4YNuHr1Kq5fv47l\n5WWcOHECAwMDKk0Qyv3793H37l0AwL179zA5OYnu7m4MDAxgbGwMADA2Nobt27frNDMxtfIzMDCA\n77//HsvLy5iZmcHVq1dLIyc2MT8/X/r7hx9+KI0ouJK/UKjunfzxxx/9tWvX+m1tbf5nn32mOnmh\n/P77735PT4/f09Pjv/DCC6X83Llzx9+8ebPf0dHh9/f3+3///bdmS8Ozc+dOv7W11W9ubvYzmYz/\nzTffBObn008/9dva2vx169b5ExMTGi0PR2X+jh8/7r/11lt+d3e3v379ev/111/3b926VbretvzF\nxfN9LkMmpNHhzEJCiHghcG3mICGNgNCmwcOHD7Fu3Tr89NNPSKfTePnll/Hdd9/h+eefF5UEIUQC\nQiMCV2cOEuI6QoXAxZmDhDQCTSJvFmbW1UsvvYRLly6JTJYQEpKenh5MTU098r7QiCDMzMFLly7h\n448//m+Rh7/iX1dexfy5+nI5fy7nzff9mpWwUCGIOnPQ91f+SwjRg9CmQVNTE7788ku89tprePjw\nIfbs2VN3xKDYmqAYEKIPoUIAAFu3bsXWrVsDr8nlcgD+LwKVeJ7dwlDMn6u4nD+X8xaE8inGnueh\nVpKVAmC7IBBiGrX8z6gpxr6/MkoIEgFHl4UTogUjhKDS+enkhKjFCCGorPnDiEGYJgMFhZBwGCEE\nRcJGBmEdnP0LhITDKCEIGxnY4uCMSIgtGCUE1WCfASHyMV4IAHsigEpstZs0HlYIARAuKjA5cjDZ\nNkKsEYJi7RrkULJq4KROzIlRxHSMFoJaHYWqa9ewTlzLLooAMR3jhMDz6i9EKopBkCDU+zyOXfWg\nwxNbEb7oKCminEm0U3INBHEZ4yKCsPi+umZCZRpJRYAdh8Q0rBWCIknFQEfIz2iCmIb1QgAkE4Mw\noxGEuI41QhBmEVK9a3QMPRJiA9YIQRhHrXcNnZ2Q6lgjBGGR3V+gm6KNNthK7MFaIag3eSeOo0SN\nGJI0ReLCnZ+JDIwVgqROVK/PQIST1hMdOiuxBWOFQER7P0gMRDpp5b1ERQIM/4kqjBUCUcQZWhQR\njYiAEQVRhfNCAERvJtABSaPhlBDUmyegoy3fSOF9I+XVNZwSgjD9CqYuYSZEJ04JQRjomPLgs7WX\nhhMCgCFsHPjM3MZpIQjqE2DBDg/3X3AfJ4QgTiegCasOKUbEFJwQgiTHn+mMDljLElNwQgiAZGcl\nqtjcxBZsPlmKxMcZIUga6ovY3CQqLgkIsRtnhKBIktpL9R6IJta0JtlEoVSHc0JQjSgFSoUYmORs\nhAANIgRRHc/24UWbbS+HgqmOhhCCONhcCG22neiBQhCAjjMTTL0ncRvrhKCykMss9ComHYlIg45P\nkmKdEFQ6joowWGcHYpwDWHQ0DWyNnsi/WCcERVQX9iRikKQA29LeD7IzTP55yKxerBUCHcQVA9cL\ncJJZnVGuIfKwWgh0hIq2Dy3GReUpUY34fHVjtRDoqkUaUQxYY7uN1UIQhIrOPdfEoFZ+uL2b+zgr\nBK6MJpgAHdN9nBUCVbjkJC7lhUSjKe4Xs9ksnnrqKTz++ONobm5GoVDA4uIi3njjDfzxxx/IZrM4\nefIknn76aZH2RiLKFltJtuNSvZWXCVuHqXq2RA2xIwLP85DP53Hx4kUUCgUAwOjoKPr7+3HlyhVs\n3rwZo6OjwgxdmXa466IUvqTLlwG5zYTye9vmVLLsbYRmmSoSNQ38il/49OnTGB4eBgAMDw/j1KlT\nSW4fkK6U2yZGZp+BaXnWZY/NgmgyiSKCvr4+bNiwAceOHQMALCwsIJVKAQBSqRQWFhbEWJkAHT3e\nrKnEQueXT+w+gvPnz6O1tRV//fUX+vv70dnZueJzz/Pg1fCIkZGR0t+5XA65XC6uGXXRUXCKYsBC\nS3STz+eRz+frXuf5lfF9DA4ePIgnnngCx44dQz6fx+rVqzE/P49XX30Vv/3228oEPe+RJkWggTEX\nF5ngiCJtCLpXmHSS2mLC8yTJqeV/sZoG9+/fx927dwEA9+7dw+TkJLq7uzEwMICxsTEAwNjYGLZv\n357A5H/x/XgFsNp3bG4mhDmjIe73k6YfBJtJdhArIpiZmcGOHTsAAA8ePMCbb76Jjz76CIuLixgc\nHMSff/5Zc/gwakTwiMEaa6a4acf5nsq0SONQy/+ENA1EGEIIkY/QpoEObA8xddsvey8F3fkjybBG\nCGwPInSftRj3+YXphExyf2IG1ghBLUQ7FlctRiOKALiUb9ewXghE10SNtGqRy4tJEeuFICnlzqDS\nMYpiUHypTr9oQxAmiBVRQ8MLQbkzVDqGKkcopmtajZlkZ2ViFw0vBEHY0EwojyhUYZpgkeRQCEIg\ne0WhKX0GpHGhEIRAxXr6JFulJ11DYCI6Ip1GJvbqQ5IcE04oMjXMN9UuV2FEYBg214Kqdmgi4qEQ\nKCBKIdY9AzFJ2jJrcUYIcnFeCEyoSeIuo5Zlu4gjymSlTfTgvBCI3BhENTaMJkS1jzW7mTgvBKJw\n6Xg1VRumEHugEAik0rlEOpsMMSiiMuowPcJpVCgEApE5HFi+NkE0Kmt1UZ2hFBSxOCcELCD/x9Rn\nIWI7NRv6T2zCuQlFLrdZo+bN1Gchyi5T82cjzkUEOlBdM8k8vlxkXlhj24NzQqBrmE91etXyKcIO\n0f0a9dC1FwNZiXNCYGO4GKcT0IQZiCKQtReD7c9FNc4JgQxk7wCcBNs6zWyytZGgEIRA1g7A5dcl\nPZbdFgezMWJrBCgEFlKrf8AWMQjChTzYCIVAATJnBVa+3+gbevCchXhQCGogcndhFsrwuBLZ2AaF\noAbl7XabHLlot4nOVM8m1ub6oBA4ionOVM8m0+ZBNBIUAsuIUtObGBXYiuvPkkJgGVFqPFcmHRH5\nUAgkous4tUqi9BlQNBoTCoFEgo5TU00YMUiyPNiUhVf1PouL7t9PNhSCBqKeGCSd3agS1x1TNRQC\ng5ExOcjUoUWiFwqBwciaEyBLDCgw9kIh0EScQ09EIkoMKs9vNAFb+jlMgkKgiWqFVcfx5iqmTzey\ng9kChcAgXD0E1ZRIoR622CkDCoEl8IDReLicN5FQCBKgcsmviNqqlq2cgUgoBP8Rd99AHeFkHIcN\nM1nI5KFFU+1yBQrBf+hy6lqInvgTZds0E4cWTfptXIRCYCg6C77ssxZVIrNJ5RIUAo2YXMBMbSaY\naJMLUAgUIetAEpmoFoMwaYl8ZmHzZqooiiRQCHbv3o1UKoXu7u7Se4uLi+jv78fatWuxZcsWLC0t\nlT47dOgQOjo60NnZicnJSXlWC6B8P0IV5xaY7vS1SOoEJjyjeqMlYbD19wtLoBC88847mJiYWPHe\n6Ogo+vv7ceXKFWzevBmjo6MAgOnpaZw4cQLT09OYmJjAu+++i3/++Uee5Qkp349Q9rkFtmPTqkRT\nbTCdQCHYuHEjVq1ateK906dPY3h4GAAwPDyMU6dOAQDGx8cxNDSE5uZmZLNZtLe3o1AoSDI7GUnD\nvEbc5MOEvOi0QebBsyYQuY9gYWEBqVQKAJBKpbCwsAAAuHnzJjKZTOm6TCaDubk5QWaKJWkNEWUo\nzhVMaCeXP08d6zJcJlFnoed58AJ+kaDPyL/Y9IhMmYGYZCelavdKgisC0RT1C6lUCrdu3cLq1asx\nPz+PlpYWAEA6ncbs7Gzpuhs3biCdTle9x8jISOnvXC6HXC4X1QyjKBbMOAXUxoIUN69FkjqyqqPb\nRQqOLvL5PPL5fN3rPN8Pzur169exbds2XL58GQDwwQcf4JlnnsGBAwcwOjqKpaUljI6OYnp6Grt2\n7UKhUMDc3Bz6+vpw7dq1R6ICz/NQJ0liCaJrZtOKhYk2JaWW/wVGBENDQzh79ixu376NNWvW4JNP\nPsGHH36IwcFBHD9+HNlsFidPngQAdHV1YXBwEF1dXWhqasLRo0fZNKhC1MIlszCKqJlF2Weiw1Xa\n5KIwFKkbEQhP0PKIwOXCEJckTaO4aem4nwu/fS3/48zCkCQ9l8+E4EiWDapEoJhWWMIcYpt0UpEJ\nv6sIGkoIkswiVDXkKLNgyRp+K48ITESGQLl2YGvkUQObseFHU2WjjJ53056vTHtMy2tSGioiIPKp\nFo6bGikkwbU8UQiIdFyrPV2EQmAhJtdG5c0Ek+1MSi1xszXPFALJuHAgZ5IOVpMdw2TbVEMhiECc\nFWimh8WiNwOpvJ/JYmD6b6MSCkEEahUcmwuU6LHxavczZUNUFYJka1mgEGhARYEU7cxJqScGcey1\n1elMhEKggTiTi6I6kYlOEiQGKpYVm7KM2kQoBAZT7hxBhdiWQ1YA+X0GJgqgDVAILMOUgm7KfgI2\npm8iFAKLcCmkdSkvRWzOE4XAInTu2ScattfNgkJgKa6EtybPM6iGzXNGgqAQRMCkAmuCLaLG8eOI\nga45BDY7exAUggiYVAjq2SKq4IusAYOujyoGqjaJNUFwVUAhcIjyjVdEiZZK8YvbTOBwZHIoBDFI\nstORTIrHt5loW1hst99WKAQxSHJeYhBsx/6LqGaCLEFxUagoBAZhuwOLxOQt0l38nSgEFegM+12s\naZJgWn9B8f4u/k4UggpEhP1xC4uLNU0SKicd1XumKrZUl9Us1A2FQAJJCovu2kZ3+tVQeW5Co0Ih\nCIkqB9Fd2E3ciy/suQlRn121+5kohCqgEITEhANKVFOel3qnBstE1oaotXZTioPtvzuFQDC6a3SR\nhI0ObJh0FJWoadj+u1MILCKpA4j6vu45FCLEQHQzw3YoBAJRESLr3IswyfdFO57IZxF2d2rbw/8g\nKAQCMe00YJPQsSFqOfUWT4XZAs7WZx8GCgGxmrDOWe86l508DBQCYhRxJmPZHLKbYjuFgBhFtclY\n9cSB254lh0LgGEWnMdUpothVPkoRJnTnEub4NOk2gIglysQfHe3iKGnGsc/06ciVtpliJ4XAcqI4\ntymFTjZJxMBkEZEJhcByTCm0oicbJXXIuGJgyvNUDfsIyIplvkmWUIt0IlEbk7DPIBwUArJiUU+t\nXvt6mOpwKqYjR7mHqREHhYDUJWyPvanonFot8h4yoRAQI4hb64qYYqwC3enXg0JAlBB2UlBUokwx\nNt0ZdUIhIEooP3NBl0PqnIFo+gG2FAKiFFGjC0mXIJvojDqhEJBATOhxr3YvER2Aug5AMbHjkEJA\nAqlWaKNs2GHa3ILK+5Xb3shRQqAQ7N69G6lUCt3d3aX3RkZGkMlk0Nvbi97eXpw5c6b02aFDh9DR\n0YHOzk5MTk7Ks5poxZQNO0Q4brkYmFhTq8Lz/drZP3fuHJ544gm8/fbbuHz5MgDg4MGDePLJJ/He\ne++tuHZ6ehq7du3ChQsXMDc3h76+Ply5cgWPPbZSazzPQ0CShERC1NoAkWsMTF6vUMv/AiOCjRs3\nYtWqVY+8X+1G4+PjGBoaQnNzM7LZLNrb21EoFBKYTGRj8nLlsIg4lap4Hx2H0JpytkKsPoIjR46g\np6cHe/bswdLSEgDg5s2byGQypWsymQzm5ubEWEmk4OrxXVHQvSRY5NkKSYgsBPv27cPMzAympqbQ\n2tqK/fv317zWs726IUYRZ1MTmWnIQocNkZcht7S0lP7eu3cvtm3bBgBIp9OYnZ0tfXbjxg2k0+mq\n9xgZGSn9ncvlkMvloppBGpDyCUGyNiONkkY5pvYL5PN55PP5utcFdhYCwPXr17Ft27ZSZ+H8/Dxa\nW1sBAF988QUuXLiAb7/9ttRZWCgUSp2F165deyQqYGehG5jQ0y7b+ZLe30RxqOV/gRHB0NAQzp49\ni9u3b2PNmjU4ePAg8vk8pqam4HkennvuOXz99dcAgK6uLgwODqKrqwtNTU04evQomwYOI+u0I9lb\nmUVBxrZnJooDECIiEJ4gIwIrMKHGF4Wumt1Ep481fEjspnLnoSiULxIqv5eN6JqOHGUjWd1QCBym\ncuehpPeohs65CCrTdX2hEoWAJELEXISwh5BWS1slsicd6RQaCgHRjmlOEYSrkQG3MydGYlonWzmy\nbNOZZ0YERBjl/QUu1prl2LJXYlgoBCQU9ToFi0Nl9ToXXSHstmdRn4Mu4aAQkFDU6xR03fFrEaXP\nIMx1up4jhYBYhaoaM0o6YcUgaAhWNxQCooW4fQmqasw4ZybG3fbMhGiKQkCkEXSWoot9CVG3PTMh\nEijC4UMiDdOdPMqR8mGJslApbHoq1ixQCEjDItq5ykdOTFxwFASbBoQIQta2ZyoEhUJASA2C+jjC\nfl+UHbKhEBBSg2KYH3ctRK1hRZM6CYtQCIjzJK3ZqxHlFOaiDWG/a8125oTYRFDNLoIk0UGta1VD\nISDOUR4BqDjENUp0ENUeVaMPFALiHOURQK32varaudpcBRMPXqUQkIYjTjMhbnRR6yQj0zaHpRAQ\nEoKw4hEl0jAlGgAoBMQBkuzWXO0+SYgiFiaJAacYE+sRtYBJx8pGWUe3RYURASECkBWJqIoYKASE\nRKDesuo4VJt0pLozkUJASAQqOw1Fz1bknoWEaCbOLsxR9hQIez8dYsDOQkL+Q9Yy4rD3qycAMmcZ\nUggIMYRKJ1e5uQmFgBCF1Kr1dc9ApBAQIpig9QVxHVp2dEAhIEQwtWp3UfeSAUcNCDGcyrkLMkYV\nGBEQYiDVdjTiqAEhDUa1EYTK/4sUBQoBIRYgu6+AfQSEWEb5ISqioBAQYikixYBCQIjFiBIDCgEh\nliNCDCgEhDhAUjGgEBDiCElGFigEhDhE3KiAQkCIQ1Tb9iwMFAJCHCRqnwGFgBBHiSIGWoQgn8/r\nSFYZzJ+9uJa3sGJAIZAA82cvLuYtjBiwaUBIA1BPDCgEhDQIgfMMfMVs2rTJB8AXX3xpeG3atKmq\nX3q+b8oJ7YQQXbBpQAihEBBCKASEEFAICCGgEBBCAPwPGJgYNsnkv+8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0xbaf7a90>"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L = la.cholesky(A_reordered)\n",
      "print \"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0])\n",
      "pt.spy(L, marker=\",\", precision=prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1283 non-zeros\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 59,
       "text": [
        "<matplotlib.lines.Line2D at 0xc242d10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD7CAYAAACBpZo1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF09JREFUeJztnU9sFsX/x9+Ptjc0EmMfmj4kj2kLtVpqE/RGqKHFcKDC\npVKMNgIXPEoUPVlMlHIwHlASQzDpRYWLlIM0jYcHCZeHREqIjQFi0VJKI9QmBA4NuL/D7/s8bpeZ\n3dnd2X+z71fS9Onu7PxpO+/5zMxnPluwLMsCISTXPJF0BQghyUMhIIRQCAghFAJCCCgEhBBQCAgh\nSEAIJiYm0NHRgfb2dhw5ciTu4rVTLpexYcMG9PT04NVXXwUALC4uor+/H+vWrcPWrVuxtLSUcC3V\n2bNnD4rFIrq6uurX3Npz+PBhtLe3o6OjA5OTk0lU2Rei9o2MjKBUKqGnpwc9PT04e/Zs/V7W2hcY\nK0YePnxotba2WjMzM9by8rLV3d1tTU9Px1kF7ZTLZevu3bsrrn3wwQfWkSNHLMuyrNHRUevgwYNJ\nVC0Qv/zyi/Xrr79aL730Uv2arD2//fab1d3dbS0vL1szMzNWa2ur9ejRo0TqrYqofSMjI9YXX3zx\nWNosti8osVoE1WoVbW1tKJfLaGxsxK5duzA+Ph5nFSLBcvhknTlzBsPDwwCA4eFhnD59OolqBWLT\npk1YvXr1imuy9oyPj2NoaAiNjY0ol8toa2tDtVqNvc5+ELUPePxvCGSzfUGJVQjm5uawdu3a+s+l\nUglzc3NxVkE7hUIBfX192LhxI44fPw4AWFhYQLFYBAAUi0UsLCwkWcXQyNpz69YtlEqleros/z2P\nHj2K7u5u7N27tz71Mal9XsQqBIVCIc7iYuHChQu4dOkSzp49i6+//hrnz59fcb9QKBjVbq/2ZLGt\n+/fvx8zMDKamptDc3IwDBw5I02axfSrEKgQtLS2YnZ2t/zw7O7tCcbNIc3MzAOC5557Dzp07Ua1W\nUSwWcfv2bQDA/Pw8mpqakqxiaGTtcf49b968iZaWlkTqGIampqa6wO3bt69u/pvSPhViFYKNGzfi\n2rVruHHjBpaXl3Hy5EkMDAzEWQWtPHjwAPfu3QMA3L9/H5OTk+jq6sLAwADGxsYAAGNjY9ixY0eS\n1QyNrD0DAwP44YcfsLy8jJmZGVy7dq2+c5Il5ufn659//PHH+o6CKe1TIu7VyZ9++slat26d1dra\nan3++edxF6+VP/74w+ru7ra6u7utF198sd6eu3fvWlu2bLHa29ut/v5+659//km4purs2rXLam5u\nthobG61SqWR9++23ru357LPPrNbWVmv9+vXWxMREgjVXw9m+EydOWG+//bbV1dVlbdiwwXrjjTes\n27dv19NnrX1BKVgWjyETknfoWUgI0S8EpnkOEpIHtE4NHj16hPXr1+Pnn39GS0sLXnnlFXz//fd4\n4YUXdBVBCIkArRaBqZ6DhJiOViEw0XOQkDzQoDMzFa+rl19+GZcvX9ZZLCFEke7ubkxNTT12XatF\noOI5ePnyZXzyySf/O+RhrfhuyletfaZ+mdw+k9tmWZZ0ENYqBH49By1r5XdCSDJonRo0NDTgq6++\nwuuvv45Hjx5h7969njsGtdkExYCQ5NAqBACwbds2bNu2zTVNb29v/bNl/ScGNQqFbAuDvX0mYnL7\nTG6bG7G7GBcKBciKdApA1gWBkLQh63+pcjF2WgduImDosXBCEiFVQgCIpwqEkGhJnRAAamKgMmWg\noBCiRiqFAHAXA9UOzvUFQtRIrRAAcjHISgenRUKyQqqFAOCaASFxkHohALJjATjJar1J/siEEABq\nVkGaLYc0142QzAhBbXR161BRjcBhOzEdo0jayYwQ1EhizUC1E8vqRREgaSdzQgD8JwZuguB13y8q\nebHDk6ySSSEAku10nO8T08isEADu0wTL0icWzjl+2HwpJCRtZFoIgPBrBkmY/JxCkLSReSEAwomB\nym4EIaZjhBAAamKQxNajStmEJI0xQgB4d+YkTXJOB0iaMUoIgOjXC5KmVscs1JVkB+OEwD7n99tZ\nVCMiqXbGKDorIz+TKDBOCGp4dZS4nY0ISTPGCgHg7WfgNy/nZ9norCtEOwWGxIXRQgAE21oM2wF1\nme00/0lcGC8EgLcYOO+xA5K8kQshAPROE/yQp1X+PLTRVHIjBECyR5hpZZA0kyshANgho4S/2+yS\nOyEAaMIGgb8zs8mlEITdSXDO+8O+fyHtMNSa+eRSCAD/pw79+BGInhFhilCQ7JNbIaiR5HsTOMqS\ntJB7IQDiCW6SFbL8ZikSHArB/9AR3MQvQcrLk18CiQ8KgY24pglhziJE7ZeQptGfYhcfFAIHcYhB\nmjobIQCFQEjWX7ya5brboWDGB4VAQpb/CbNcd5IMFAIX4lwvSHuexGwoBC7oCnXufN7+s44y2PFJ\nWCgECjjXDPy6GDtNdZHpLjPng7yARffUIC2vpKfgRQeFQBG7GPh1MQ7zD5yG+b5KHVQDv4ZJk4bf\nhalQCHwQdDfB9H9gr99JWCEh0UMh8IldDPJy+hDwfktUHFGhTfp9pg0KQQBq//hhTx9miTS8Rcqk\n32faoBAERGQZONExgsU5CoZph856hlk4JcGgEITAaRmI7usoI2nSMsdPw+/CVCgEITHpn9OkthB/\nNAR9sFwu4+mnn8aTTz6JxsZGVKtVLC4u4s0338Sff/6JcrmMU6dO4ZlnntFZ31Qiswp0hPgS5RFV\n6DCZ34NKHdzqpKu+DJkWHYEtgkKhgEqlgkuXLqFarQIARkdH0d/fj6tXr2LLli0YHR3VVtE0I/MX\n0Dk1EHkj6say/Dk7qaaJ6s1PXDPQR6ipgeX4y5w5cwbDw8MAgOHhYZw+fTpM9pkjylOLaRsJk6pP\nHIKYR0JZBH19fdi4cSOOHz8OAFhYWECxWAQAFItFLCws6KllhtApBhzx/h92/ugJvEZw4cIFNDc3\n4++//0Z/fz86OjpW3C8UCihI/pNHRkbqn3t7e9Hb2xu0GqnEazfBTz6EhKFSqaBSqXimK1hO+z4A\nhw4dwqpVq3D8+HFUKhWsWbMG8/PzeO211/D777+vLLBQeGxKYSo6FwvteflZqNORRiW0WhwLeVws\nDI+s/wWaGjx48AD37t0DANy/fx+Tk5Po6urCwMAAxsbGAABjY2PYsWNHiCpnHxWnI5U87N+dn0U/\nu+UTJI1sEdGZRtTGOByNSHgCWQQzMzPYuXMnAODhw4d466238PHHH2NxcRGDg4P466+/pNuHebII\naqiM4EFGeVHaMIFRo0KXxULCI+t/WqYGOipCCIkerVMD4p+kdwCClq/6XKEQ/dSARAeFICZ0hT0L\nW77u5+xTkaDOSCR5KAQxo7qo5kcwwjwbFj8d3W+9aE3EB4UgAURiECbuYNBn4+5ofq0DWhPxQSFI\niJoY2OfWfs796+jEbj4Dussi6YZCkCBRjXhRd1wKg3lQCBJGh+chOyYJC4UgBYjOJsg6t1ecxLCi\nIPJg1PXeBC4WphcKQUpwLiB6dXQvoVDFbcdBVlZYd+kg9SLRQs9CQnIEPQszQpZHwijrnuXfSxag\nEKQMvx6IzsVCHYuHuk1/HeHQaURGC4UgpahGOnK69qocGY670/lZBBVBayB6KAQpJsoYiLrw63QU\nRGSSPqeRBygEKScKMdCZp9d2pk44PYgOCkEG8Oq4QebgInNdlo9o3SHqY81x55V3uH2YIdIYfSgo\nYSMSMaJRMLh9aAAqi252goyYIoci3TDce/oIHM6cJIOfOXmYhbmgz/stI8k8yH/QIsggbnP5KPIO\nGkiFx5izA4Ugg8gWD3W/G9BtTULlIJLMehHFYAhyIIkCow8KQUax763LRmcV8znIWkPYHQO705OK\nWIjwcpyiSPiDQpBxVK0Ale1HP1uUQcWB24fphNuHhsDttJXw9yGG24eGo8tbUOfIHmT+z1E+GSgE\nBmEPiBo2D5V0OtLU0O1bQGvAHxQCw9AZwkz3NqWbSPl1llIpi6hDITAQlfcmeD3v9kyY9xN4rfSr\nHKP2WybxhkJgKGnsCCqvVo+6DCKGQmAwOqcDOsmiQ5HpUw2eNTAYu9OR6kgZd1yBoPEMOPLrhRZB\nDlCNZyDzUlQZfVVH6LgcikwfwXVDiyAniF6iYr8n+iz62asMHWlUcLNyonAmMt0CoUWQI7wsgzD/\n7HF3FNM7ZtxQCHKG7hiIfhf7wixg0tyPDgpBDokqeKkuj0S3Mkg0UAhyShRnE9LSUeN6SatJUAhy\njEgMRDEL3dLoinkgeiZMx8xzpw4CjyETkiN4DJlIUYlcnNX3GjBishoUAuLpgagzPmLYdxnYP+s8\naZl3KASkjnPNQDZP9+tF6Jan7MyBKJ3T8cntSDPFwB9cIyCPodszz2+wkDDlO59lyLKVcI2AKKN7\nRPUbY0Cnh6MOEciDdUEhIELsYhB2oTDoYqTKVqZXmbriJZouBhQCIkV2UEl1QVGU3i3qkdM/QeXF\nKiplylCN4pSHqYWrEOzZswfFYhFdXV31a4uLi+jv78e6deuwdetWLC0t1e8dPnwY7e3t6OjowOTk\nZHS1JrEhCojqx9nHT1rn3N4pDKKAJrIy3X4WleeG6dYA4CEE7777LiYmJlZcGx0dRX9/P65evYot\nW7ZgdHQUADA9PY2TJ09ienoaExMTeO+99/Dvv/9GV3MSG87RWTRau42mQc8XiMpQOXfgVVcddTEN\nVyHYtGkTVq9eveLamTNnMDw8DAAYHh7G6dOnAQDj4+MYGhpCY2MjyuUy2traUK1WI6o2iZs4HXNU\n1xFU8oi63qZYC77XCBYWFlAsFgEAxWIRCwsLAIBbt26hVCrV05VKJczNzWmqJkkaXTsJYYKX+Al5\nrvvEYq4tAi8KhQIKLn8Rt3ske9g9EKPELTxa7bqXV6HO+AhRBXNJE75DlRWLRdy+fRtr1qzB/Pw8\nmpqaAAAtLS2YnZ2tp7t58yZaWlqEeYyMjNQ/9/b2ore31281SIKIdhNEnbN2zX7Py7nIzc3ZOVdX\n2dHwciiK6z0LSVGpVFCpVDzTeXoW3rhxA9u3b8eVK1cAAB9++CGeffZZHDx4EKOjo1haWsLo6Cim\np6exe/duVKtVzM3Noa+vD9evX3/MKqBnoTno9NpLowdgGusUFln/c7UIhoaGcO7cOdy5cwdr167F\np59+io8++giDg4M4ceIEyuUyTp06BQDo7OzE4OAgOjs70dDQgGPHjnFqYDj20V6108jShXk2aDqv\n5/xaF1mGZw1IaPyKgY6yksjPBCHgWQMSGV47CjoW7pyLhCrluS061u6HdSoyxeilEBAt+F3885OH\n/b7fk4xui31+tkRl8ReybiHU4NSAkBzBqQGJBdE+vinmsx3T2kQhIFqhsZdNKAQkEtzm30EWFoPG\nQAgzcrtZNabFS6QQkMjwE89A5R6g7hLsFljFnk7FfdiZn2o9sgSFgESKyDKQnSFwfhelUxEKpwCp\nHpH2qqfJ0x4KAYkcUfQhURrnd6/IRn4PA6meWJTVySt/t+tph0JAYkHF6Ug1spA9T5VyiTcUAhIb\norBn9nuyE4qyub7K3F9kifj1dPQ6Fm0CFAISK7LO7TZ1EJnosg5YExSdgUn8TBGyKgz0LCQkR9Cz\nkKSKsCOn19afjrL9bhdm1RoAKAQkIbzm/6rPhynbTzrTjVgKAUmUMHEMkuicpsYvpBCQxJHtJriZ\n5qqr+F4OSm74dWgKM11JGgoBSQUqsQj8PiNyDXbe90LFTdlen6xaBRQCkhqcloHbdqHXVqLT1VjU\nSVXER8VN2VluFqEQkFShI5qQPR+dUZZVrtMiIEQTsjUDmRjIVvd1rAnUrsvyNcUi8P2CE0LiQMeh\nHr+js2qZYY5RpxVaBCS1qBxXdnvG+dlrtFbN0y2vrFoEdDEmqUcUY0CWroZXurz+C9LFmGQWuwio\nHk12G7lV1xHcTh+aZhFQCEgmcFoEblOAWnr7cyJEloHIqnDzYeCuASExY99N8BO9yM+WpCwvmWD4\njXGQVrhrQDKF2wivczQWdX438cmqJVCDQkAyh10M3EZq2fqAPY2KteC2UOmcrmRVEDg1IJnET7xC\nlTMJKuVk3fx3g9uHJNOojsIqW495+Lfk9iExkiBBRsLkYyoUApJ5RNuHukOLReU3kJbpBoWAZB7n\nYp39u1fwErdgKH6eU3k+zXDXgBiDfaVf1sFl706oXfOKgSDbrbCnySK0CIhR+H0HgpvjkT2giTPI\niUqYNBFptRgoBMQ4nPv+oghFsq1FWQeXeTM6r6W1o3tBISBG4jZNsK8fqMz77fedz/lZnEwzFAJi\nLM7pgajTunkayjwT3Q4hZRUKATEaL3PeTxgyURpV4jgfEQZ6FhKSI+hZSHKN13qA2yKhV75ueTjv\n+fFbiBMKAckFbnN/1c7oFe0oy4YuhYDkBpXow14eg0EDl6b9FCOFgOSOIAFHnD4JNcJYE6LykoJC\nQHKJineg099AZWvRa10hqJBEDYWA5BanGIiCo4qCoLoJiNcZhLQGO3UVgj179qBYLKKrq6t+bWRk\nBKVSCT09Pejp6cHZs2fr9w4fPoz29nZ0dHRgcnIyuloTognZQSJRqDL7PdGUwQ0voUgaVz+C8+fP\nY9WqVXjnnXdw5coVAMChQ4fw1FNP4f3331+Rdnp6Grt378bFixcxNzeHvr4+XL16FU88sVJr6EdA\n0ojOCEWyqMdp+LcP5EewadMmrF69+rHroozGx8cxNDSExsZGlMtltLW1oVqthqgyIfEhmiY4P4vO\nGojSyPIX3Q96ilE3gdYIjh49iu7ubuzduxdLS0sAgFu3bqFUKtXTlEolzM3N6aklITEg2z0QnWJ0\n211Q2aaU/Sy7FjW+hWD//v2YmZnB1NQUmpubceDAAWnaQlqWRAlRxO9IL0qvGuVIlm8S3cZ3hKKm\npqb653379mH79u0AgJaWFszOztbv3bx5Ey0tLcI8RkZG6p97e3vR29vrtxqERIIz7Jlz10A2onsF\nLXGml/2sm0qlgkql4pnO89DRjRs3sH379vpi4fz8PJqbmwEAX375JS5evIjvvvuuvlhYrVbri4XX\nr19/zCrgYiHJCm7bibL0NZzCUfuc9O6BrP+5WgRDQ0M4d+4c7ty5g7Vr1+LQoUOoVCqYmppCoVDA\n888/j2+++QYA0NnZicHBQXR2dqKhoQHHjh3j1IBkGlGkI7eOrBqhyO0YdFLwGDIhHrgFMrGP9LJ7\ntTxk04s44TFkQgIimteLrvk5dKTipRgnFAJCFJB1ZPs1ZzrR/bRCISBEEdEagch/QBbhWHavRpIL\niRQCQnzg3EXwO9Kn1TLgm44I8UmYqER+YiHECS0CQgIge6+BLDai13fZtbigEBASAJVdgyAdOqmp\nA6cGhITA7hvgZfY7nZNEaZKCFgEhIZEdYXZOE2RRjtIQ4pxCQIgGZF6FXsFQnZ9FP8cBhYAQTdiD\nl8q2GZ2WAD0LCTEQ1QjFqvENvPLRBRcLCdGM6HCRarDTpBYMaREQEgE6O3Qc4kAhICQivAKgOtOI\nfpZd0w2FgJCI8Apk4hb6LG64RkBIhDgDk9iviUjqBCKFgJAYUI1ORIuAEMNxeiCKwpnZ7zmvRQmF\ngJAYcbMMeNaAkBzhtojo9YKUqKAQEJIAsmCnXq9MiwoKASEJ4XZqMW64RkBIgri5ITs/RwktAkIS\nxs27MC7rgEJASMKI3q8Y9/SAUwNCUkISAlCDFgEhKcLtxSdRQiEgJGUk8do0Tg0ISSFxvzmZFgEh\nKcUeAzFqKASEpJi4Dh9xakBIymGoMkIIAPfwZjqgEBCSAaK2CigEhGQE1XcmBIFCQEjGiMI6oBAQ\nkkF0uyNTCAjJKDrFgEJASIbRJQYUAkIyjg4xoBAQYgBhxYBCQIghhNlNoBAQYhBBrQIKASEGEfSQ\nEoWAEAPxu2ZAISDEUPyIQSJCUKlUkig2Nti+7GJa21TFgEIQAWxfdjGxbSpiwKkBITnASwwoBITk\nBFc/AytmNm/ebAHgF7/4lcDX5s2bhf2yYFlJvISZEJImODUghFAICCEUAkIIKASEEFAICCEA/g9u\nKnj3e/xq4wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb773910>"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}