{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import division\n", "\n", "import numpy as np\n", "import numpy.linalg as la\n", "\n", "import matplotlib.pyplot as pt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "np.random.seed(40)\n", "\n", "# Generate matrix with eigenvalues 1...25\n", "n = 25\n", "eigvals = np.linspace(1., n, n)\n", "eigvecs = np.random.randn(n, n)\n", "print eigvals\n", "\n", "A = la.solve(eigvecs, np.dot(np.diag(eigvals), eigvecs))\n", "print la.eig(A)[0]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.\n", " 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.]\n", "[ 25. 24. 23. 1. 2. 3. 22. 4. 21. 20. 5. 6. 7. 19. 18.\n", " 8. 9. 17. 16. 10. 11. 12. 15. 14. 13.]\n" ] } ], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "# Initialize Arnoldi iteration\n", "\n", "# Set up Q and H\n", "Q = np.zeros((n, n))\n", "H = np.zeros((n, n))\n", "\n", "k = 0\n", "\n", "# Pick a starting vector, normalize it\n", "x0 = np.random.randn(n)\n", "x0 = x0/la.norm(x0)\n", "\n", "# Poke it into the first column of Q\n", "Q[:, k] = x0\n", "\n", "del x0\n", "\n", "# ritz_values = []" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 82 }, { "cell_type": "code", "collapsed": false, "input": [ "# Carry out one iteration of Arnoldi iteration\n", "# Run this cell in-place (Ctrl-Enter) until H is filled.\n", "\n", "print k\n", "\n", "u = np.dot(A, Q[:, k])\n", "\n", "# Carry out Gram-Schmidt on u against Q\n", "for j in range(k+1):\n", " qj = Q[:, j]\n", " H[j,k] = np.dot(qj, u)\n", " u = u - H[j,k]*qj\n", "\n", "if k+1 < n:\n", " H[k+1, k] = la.norm(u)\n", " Q[:, k+1] = u/H[k+1, k]\n", "\n", "k += 1\n", "\n", "pt.spy(H)\n", "\n", "# ritz_values.append(la.eig(H)[0])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "24\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD7CAYAAABOrvnfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADAZJREFUeJzt3E9M2/Ufx/HX9+e4bcZtQmnsTHWjDPlTiGQmJkgXVvQi\nQkgMxhAy8OLNaNxOxnrRcjBmmx7MshmOehnuIOgJnbtwsEs0O2giREJaVKBxzMQx8/Hw02a4Cl3/\nl/fzkXwSKOPbz77Zk893Xz6t55xzAmDC/yo9AQDlQ/CAIQQPGELwgCEEDxhC8IAhZQl+dnZWR48e\nVVNTkyYnJ8vxlAULBoPq6OhQV1eXjh07VunpZDU+Pi6fz6f29vbMY2tra4pGowqFQurv71c6na7g\nDLfKNt9YLKZAIKCuri51dXVpdna2gjO829LSko4fP67W1la1tbXp7Nmzkqr7PG/Lldjt27fd4cOH\n3cLCgrt165YLh8Pu+vXrpX7aggWDQbe6ulrpaWzrq6++ct98841ra2vLPPb666+7yclJ55xz8Xjc\nnT59ulLTu0u2+cZiMffuu+9WcFbbSyaTLpFIOOecu3HjhguFQu769etVfZ63U/IVfn5+XkeOHFEw\nGFRdXZ1GRkb06aeflvppi8JV+Z6knp4e7d+/f8tjly9f1tjYmCRpbGxM09PTlZhaVtnmK1X3eW5s\nbFRnZ6ckae/evWppadHy8nJVn+ftlDz45eVlHTp0KPN5IBDQ8vJyqZ+2YJ7n6cSJE+ru7tb58+cr\nPZ2craysyOfzSZJ8Pp9WVlYqPKOdnTt3TuFwWBMTE1V9aby4uKhEIqEnnniiJs+zVIbgPc8r9VOU\nxNWrV5VIJDQzM6MPPvhAV65cqfSU7pnneVV//l9++WUtLCzo2rVr8vv9eu211yo9paw2NjY0PDys\nM2fOaN++fVu+Vgvn+R8lD/6hhx7S0tJS5vOlpSUFAoFSP23B/H6/JKm+vl5DQ0Oan5+v8Ixy4/P5\nlEqlJEnJZFINDQ0VntH2GhoaMsG89NJLVXmeNzc3NTw8rNHRUQ0ODkqqvfP8j5IH393drR9++EGL\ni4u6deuWPv74Yw0MDJT6aQvy+++/68aNG5Kkmzdv6osvvthyZ7maDQwMaGpqSpI0NTWV+QdarZLJ\nZObjS5cuVd15ds5pYmJCjz32mF555ZXM47V2njPKcWfws88+c6FQyB0+fNi9/fbb5XjKgvz4448u\nHA67cDjsWltbq3bOIyMjzu/3u7q6OhcIBNzFixfd6uqq6+vrc01NTS4ajbr19fVKTzPj3/O9cOGC\nGx0dde3t7a6jo8M999xzLpVKVXqaW1y5csV5nufC4bDr7Ox0nZ2dbmZmpqrP83Y856r4FimAomKn\nHWAIwQOGEDxgCMEDhhA8YEm+t/dnZmZcc3OzO3LkiIvH43d9XRKDwajgyCav4HN5BVyl/7IMhvWR\nTV6X9LX8CjjAsryCr9VXwAHW5RV8rbwyCMBWeQVfq6+AA8zL56bd5uame/TRR93CwoL7448/uGnH\nYFThyGaP8rBnzx69//77evrpp/Xnn39qYmJCLS0t+RwKQBmV7NVy/D8fqKxsabPTDjCE4AFDCB4w\nhOABQwgeMITgAUMIHjCE4AFD8tppl6ud9vSwOQcoL1Z4wBCCBwwheMAQggcMIXjAEIIHDCF4wBCC\nBwwp6cabneTyZjtszgGKhxUeMITgAUMIHjCE4AFDCB4whOABQwgeMITgAUMquvEmF2zOAYqHFR4w\nhOABQwgeMITgAUMIHjCE4AFDCB4whOABQwraeBMMBnX//ffrvvvuU11dnebn54s1r3vC5hwgNwUF\n73me5ubmdODAgWLNB0AJFXxJn8vqCqA6FBS853k6ceKEuru7df78+WLNCUCJFHRJf/XqVfn9fv3y\nyy+KRqM6evSoenp6ijU3AEVWUPB+v1+SVF9fr6GhIc3Pz28JPhaLZT6ORCKKRCKFPB2AQrk83bx5\n0/3222/OOec2Njbck08+6T7//PPM1ws4dElIYjBMjWzyXuFXVlY0NDQkSbp9+7ZefPFF9ff353s4\nAGXg/b36Ff/AnldVd/D5PTysydZf1b/jTbHk8sOHHwrY7dhaCxhC8IAhBA8YQvCAIQQPGELwgCEE\nDxhC8IAhZjbe5ILNOdjtWOEBQwgeMITgAUMIHjCE4AFDCB4whOABQwgeMISNN/eIzTmoZazwgCEE\nDxhC8IAhBA8YQvCAIQQPGELwgCEEDxjCxpsSYHMOqhUrPGAIwQOGEDxgCMEDhhA8YAjBA4YQPGAI\nwQOG7Bj8+Pi4fD6f2tvbM4+tra0pGo0qFAqpv79f6XS6pJPcjZxzOw6g2HYM/uTJk5qdnd3yWDwe\nVzQa1ffff6++vj7F4/GSTRBAEbkcLCwsuLa2tsznzc3NLpVKOeecSyaTrrm5+a7vyfHQ2IYkBiPv\nkU1e/4dfWVmRz+eTJPl8Pq2srORzGABlVvCLZzzP+88XgsRisczHkUhEkUik0KcDUIC8gvf5fEql\nUmpsbFQymVRDQ0PWP3dn8AAqL69L+oGBAU1NTUmSpqamNDg4WNRJASiRnW4cjYyMOL/f7+rq6lwg\nEHAXL150q6urrq+vzzU1NbloNOrW19ez3nBCYVQFN34YtTuy8f7+h1V0nufxu+QC8SYZKES2/njH\nmyqWyw9MfijgXrC1FjCE4AFDCB4whOABQwgeMITgAUMIHjCE4AFD2HhT49icg3vBCg8YQvCAIQQP\nGELwgCEEDxhC8IAhBA8Ywu/hDeB39fgHKzxgCMEDhhA8YAjBA4YQPGAIwQOGEDxgCMEDhrDxBpLY\nnGMFKzxgCMEDhhA8YAjBA4YQPGAIwQOGEDxgCMEDhrDxBjljc07t23GFHx8fl8/nU3t7e+axWCym\nQCCgrq4udXV1aXZ2tqSTBFAcOwZ/8uTJu4L2PE+vvvqqEomEEomEnnnmmZJNEEDx7Bh8T0+P9u/f\nf9fjuVzeAagued+0O3funMLhsCYmJpROp4s5JwAl4rkclurFxUU9++yz+vbbbyVJP//8s+rr6yVJ\nb7zxhpLJpC5cuLD1wJ6nN998M/N5JBJRJBIp4tRRjbhpVz2ypZ1X8Ll8zfM8LvsNIvjqka2/vC7p\nk8lk5uNLly5tuYMPoHrt+Hv4F154QV9++aV+/fVXHTp0SG+99Zbm5uZ07do1eZ6nRx55RB9++GE5\n5gqgQDld0ud1YC7p8R+47C+Pol3SA6hNBA8YQvCAIQQPGELwgCEEDxhC8IAhBA8YwjveoOx455zK\nYYUHDCF4wBCCBwwheMAQggcMIXjAEIIHDCF4wBA23qAqsTmnNFjhAUMIHjCE4AFDCB4whOABQwge\nMITgAUMIHjCEjTeoWWzOuXes8IAhBA8YQvCAIQQPGELwgCEEDxhC8IAhBA8Ysm3wS0tLOn78uFpb\nW9XW1qazZ89KktbW1hSNRhUKhdTf3690Ol2WyQL3yjm347DEc9v8jVOplFKplDo7O7WxsaHHH39c\n09PT+uijj/Tggw/q1KlTmpyc1Pr6uuLx+NYDe565k4natFt342Xrb9sVvrGxUZ2dnZKkvXv3qqWl\nRcvLy7p8+bLGxsYkSWNjY5qeni7BdAEU27Yr/J0WFxfV29ur7777Tg8//LDW19cl/f+nyIEDBzKf\nZw7MCo8awQr/LxsbGxoeHtaZM2e0b9++LV/zPG/XnjBgt9nx1XKbm5saHh7W6OioBgcHJUk+n0+p\nVEqNjY1KJpNqaGjI+r2xWCzzcSQSUSQSKcqkAeRn20t655zGxsZ08OBBvffee5nHT506pYMHD+r0\n6dOKx+NKp9PctEPN2q1XqNn62zb4r7/+Wk899ZQ6OjoyJ+Wdd97RsWPH9Pzzz+unn35SMBjUJ598\nogceeGDrgQkeNYLgi4DgUSssBc873sA8S++cw9ZawBCCBwwheMAQggcMIXjAEIIHDCF4wBCCBwxh\n4w2Qg92yOYcVHjCE4AFDCB4whOABQwgeMITgAUMIHjCE4AFD2HgDFEktbM5hhQcMIXjAEIIHDCF4\nwBCCBwwheMAQggcMIXjAEDbeAGVU6c05rPCAIQQPGELwgCEEDxhS1uDn5ubK+XRFUWtzrrX5SrU3\n51qb750Ifge1Nudam69Ue3OutfneiUt6wBCCByxxJdLb2+skMRiMCoze3t6sXXoul60/AHYFLukB\nQwgeMITgAUMIHjCE4AFD/gKflH2/NE4ELQAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 107 }, { "cell_type": "code", "collapsed": false, "input": [ "# Check that Q^T A Q = H\n", "\n", "la.norm(\n", " np.dot(np.dot(Q.T, A), Q)\n", " - H) / la.norm(A)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 108, "text": [ "6.9230906536862802e-06" ] } ], "prompt_number": 108 }, { "cell_type": "code", "collapsed": false, "input": [ "# Check that Q is orthogonal\n", "\n", "la.norm(\n", " np.dot(Q.T, Q) - np.eye(n))" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 109, "text": [ "1.3525100273749116e-05" ] } ], "prompt_number": 109 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot convergence of Ritz values" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Enable the Ritz value collection above to make this work.\n", "\n", "for i, rv in enumerate(ritz_values):\n", " plot([i] * len(rv), rv, \"x\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stderr", "text": [ "/usr/lib/pymodules/python2.7/numpy/core/numeric.py:460: ComplexWarning: Casting complex values to real discards the imaginary part\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX9wVPXV/886ZvxRbRcdiZQgyfCrQiCbaqVT7TwwCNSf\n0MfqYGvKyPpPmc58fWwjw/B0NtonGgR/l3metv4Y2nnaaDsDOgnQMnYzbHlaBQELJWIlu0ohiRJ2\nEV1oonu+f6x3k002m/2cmz177933ayZD7iXnfj75cHnv+ZzP+ZyPj5mZAAAAuJbzSt0BAAAA9oCQ\nAwCAy4GQAwCAy4GQAwCAy4GQAwCAy4GQAwCAy8kr5OfOnaP58+dTIBCg2bNn09q1a4mIqKmpiaqq\nqqi+vp7q6+tpx44dKp0FAAAwEt9YeeTJZJIuvvhi+vTTT+mGG26gjRs30muvvUaXXnopPfDAA1r9\nBAAAMApjhlYuvvhiIiLq7++nzz77jCZMmEBERNhHBAAAzmBMIU+lUhQIBKiyspIWLlxIc+bMISKi\nZ599lurq6igYDFIikSh6RwEAAORmzNCKxenTp2np0qXU0tJCs2fPpiuuuIKIiH7yk59Qd3c3Pf/8\n80XtKAAAgNycX+gPfulLX6JbbrmF9u7dSwsWLMjcv+++++i2224b8fPTp0+no0ePjksnAQCgXJg2\nbRq9++67RjZ5QysnT57MhE3Onj1LO3fupPr6eurp6cn8zJYtW2ju3LkjbI8ePUrMjC9mCoVCJe+D\nU74wFhgLjEX+L4kDnNcj7+7uppUrV1IqlaJUKkUNDQ20aNEi+v73v08HDhwgn89HNTU19POf/9y4\nYQAAAONDXiGfO3cu7du3b8T9X/3qV0XrEAAAADOws1OBoWsK5Q7GYhCMxSAYC3sUnLVi/GCfj4r0\naAAA8CwS7YRHDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdC\nDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAALgdCDgAA\nLqfoQp5IELW3F7sVAAAoX4oq5IkE0bp1RNdfX8xWAACgvCnq4curVzM1NxP5/cVoAQAAvMe4H758\n7tw5mj9/PgUCAZo9ezatXbuWiIhOnTpFixcvppkzZ9KSJUsokUjktG9shIgDAECxySvkF154IYXD\nYTpw4AD97W9/o3A4TH/+85+ppaWFFi9eTO+88w4tWrSIWlpactpv2JAOrwAAACgeY8bIL774YiIi\n6u/vp88++4wmTJhAr776Kq1cuZKIiFauXElbt27NadvcnI6RQ8wBAKB4jCnkqVSKAoEAVVZW0sKF\nC2nOnDnU29tLlZWVRERUWVlJvb29OW39/rSY7949vp0GAAAwyPlj/cB5551HBw4coNOnT9PSpUsp\nHA5n/b3P5yOfz5fTtqmpKfP9F76wgBYsWGCrswAA4DU6Ojqoo6PD1jOMslZ++tOf0kUXXUTPPfcc\ndXR00JVXXknd3d20cOFCevvtt7MfLFh5BQCAcmfcs1ZOnjyZyUg5e/Ys7dy5k+rr6+n222+nzZs3\nExHR5s2bafny5cIuAwAAsEtej/zgwYO0cuVKSqVSlEqlqKGhgRobG+nUqVN011130fvvv0/V1dX0\n8ssvk39YniE8cgAAMEeinUXdEAQhBwAAM8Y9tAIAAMD5QMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDl\nQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgB\nAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDlQMgBAMDl5BXyY8eO\n0cKFC2nOnDlUW1tLzzzzDBERNTU1UVVVFdXX11N9fT3t2LFDpbMAAABGklfIKyoq6Mknn6S///3v\n9Ne//pU2bdpEnZ2d5PP56IEHHqD9+/fT/v376Vvf+pZWfwHIS19fOw0MJLLuDQwkqK+vfdzbam8n\nSmQ3RYlE+v6oNu+0U+JctlHiXILa3xmjf4LG2vv6KDEwkG0yMEDtfX15m+pr76OBRLbdQGKA+tpH\nt5OMe3t7OyWG/U6JRILa8w0gKY67pCGSj7st2IBly5bxzp07uampiTdu3Jj3Zw0fDTxM28mTHO/v\nz7oX7+/ntpMn89qdbDvJ/fFsu/54P59sG92uvz/OR46s5v7+eM7rUfvY1sbxePbPxONxbmtrG9Um\nHmdevTr9Z67rnDZn47y6bTXHz8ZzXo9nY/H+fl595Ehm7Idfj0Z/vJ+PrD6SGfvh1zltBOMej8d5\n9erVmXEffj26ndK4Sxpi+bhbSLSzYItoNMpXXXUVnzlzhpuamnjq1Kk8b948XrVqVc6Bh5ADC01B\nYR4UkWQyWpCIM9sXlWi0oP/jaZvPRSQajxYm4jYas8Y6mkwaiYk11slosqAxZ7Y37tFotKDxHrRT\nGndJQywfd+YiCvmZM2f4mmuu4S1btjAzc29vL6dSKU6lUrxu3TpetWpVzs6EQqHMVzgcNu4c8A6a\ngsLMnExGORwmTiajhfdRKCrRKDNR+s9CicajTE3E0biBkbCxaDLJFA5zNJk0aioZTXKYwpyMFm4n\nGfdoNMpExFGTAWTFcZc0xIWPezgcztLKogh5f38/L1myhJ988sncnY1Guba2duSD4ZGr0dY20lGI\nx9P389odaRvhlcTPxrntyBiGQrQEReIZZvpoKCrwyIfYwCMfNHOSR55KpbihoYHvv//+rPsnTpzI\nfP/EE0/w3XffPS6d0UAarzXl5Mm2ES9yf3+cT57ML5JasVpmYdxQ+KmhJSjSGHn61zATFcTIh9gg\nRj5o5rQYeSQSYZ/Px3V1dRwIBDgQCPC2bdu4oaGB586dy/PmzeNly5ZxT0/PuHRGA7uDXChSQdGM\n1TILvJTWVuZgMPvlDgbT90drQ1FQpB+gknGXfKaJZ0GCxjQXmSXjLnFa0nZK4y50Wuw6ixLtPD9f\nRssNN9xAqVRqxP2bbrppHPNmdPFXVFBzTQ2ti0apccoU2nDsGDXX1JC/omJc26mo8FNNTTNFo+to\nypRGOnZsA9XUNFNFhT9///x+am5upnXr1lFjYyNt2LCBmpubye8fy46osZGopoYoGk1fF4L/Qj81\nXt9INU/XUPT/Rcl/4RiGS5cS7dxJ9OMfE/3nfxL9138N3h+F3R99lDXG1r/B7o8+olsuv3xUu38+\n9U+a8uAUqvCn7Sr8FTTlwSn0z6f+STVNNYX9ggXy1FNP0YMPPpgZZ7/fTw8++CA99dRT1NTUNK5t\nATDuGEt/gRTx0eOCNF5rimTxh1knVstsI254zz3pBaB77im8MUO0pvjMMo8coZUhNgitDJo5LbRi\nBycLuZ2FCBOki24asVpmm6JiIOR2pppai27MsoU3LHYOscFi56CZkxY77eBUIfdijFyatRL6U4hj\n8VjWvVg8xqE/hfJ1MB0TDwbTL7f1fRE8QwutNDhmWSoc0g+H2CD9cNBMOO4Q8gLwYtaKFJFHLljs\nZHa+Z5j+VeCRW8Ajt9EQwyMHyhi/3FL3n809FKfHyCWfaYiRD/2VECPPBYQciBBP8w2QeChaaXDM\nspmQRMjF6YehEHMsOwzGsVj6/mgmXV0cG/ahGUsmOdTVlbeprlAXJ2PZdslYkrtCo9t1dYU4mczu\nXzIZ466uPP0LhTg27HeKxWIcyvM7pe2Mh0IWRpQ0xPJxtygvIbfhGZoiERS3IJ7mm7ShtC5hB2lI\ny3TmLRZyeORD7OCRj7Axtij0wcUWcmmqhgBp8SanI57mG6K1LmEHqagwm62F2RpzxMiH2CFGnmVj\nbFHogzVCK8JBliAt3qSBW2qtOB3HL3YyI2slyw5ZKxkbY4tCH6wVIxcOsgTJi62B4uTE8Uhj5BYm\nomJn3EWCAo98iB088iwbY4tCHwyPXBXFoXA0mkWzpDMhW4KCGDli5LlsjC0KfTBi5OooTk4cjebB\nEqaIY+SKWSta2UKaWSteL5rlXiFH1koWXvTI7Yy7aaxWkrUieQVFaXDWg5U8ci0h1/TIIeRCkEeu\nh1dj5JpHvUlQm+IPb7DIMfJSbMQqdowcoRUhEHI9FCcn6mgdLKGVR87sjqwVzdIIWlkrWOwUACF3\nPlrph3anmibZQpoHS1iopcEpZq0w66QfamatMCP90PzBEHIRji+aJXD/7Uw1NbOFHJ1HrhgjZ9bx\nyDVj5MzwyEVAyGVoZU9k2jN9uW3GDYsdq7VLsfPIxTFyxaPetGLkmke9IUYuBEIuRzrdlGI83RR6\nKaZTTe1sIY08cjfsptUsViZBLVsIRbMg5HaRLgCZIl54M4wbap3KJEV7JgTsLzAXfSYEjxxCbsdD\n0fLIbdfGLtAj93L1Qy+iNROy8+GpvjaBGHl5onnAgRTRNF/gpbih+iEYRHNtwo7TgqJZQ2yMLQp9\ncJkLObMsr9bxnqGXk9Y9huah13ZmoJIwIjzyYTb5/vL999/nBQsW8OzZs3nOnDn89NNPMzNzX18f\n33jjjTxjxgxevHhxzk9RCHka6SHAANhF89BrzZ2diJHnsMn3l93d3bx//35mZj5z5gzPnDmTDx8+\nzI2Njbx+/XpmZm5paeE1a9aMS2e8htZWcVAGCGdCTi5jKw0jImslh43JDy9btox37tzJs2bN4p6e\nHmZOi/2sWbPGpTNewk45VQBGIN0FwzqHXluhlaEz0LFCK8gjz01RhTwajfJVV13FH330Efv9/sz9\nVCqVdW2nM15CM68WlAmCeK3Wodf9/XHu7AxyZ2eQk8lo5vt8jou0jC12duawKeSHzpw5w1/96ld5\ny5YtzMwjhHvChAk5OxMKhTJf4XDYuHMAmOL5D1CDDArNtE+JkMdiMa6trc2I+fDrfAj11ZFZK+Fw\nOEsriyLk/f39vGTJEn7yyScz92bNmsXd3d3MzHzixAmEVoBj8HRIy1C9pFkrdnZ2moRWmAfFOxKJ\nFCziFqh+OMQm31+mUiluaGjg+++/P+t+Y2Mjt7S0MDPzo48+isVO4Ci0FplVt+jbiJGbol0HPhKJ\nMBFxJBIpuI+m+lrWMfJIJMI+n4/r6uo4EAhwIBDg7du3c19fHy9atKhs0g/dcEIQyEYj7dPpRbPs\noFUHXuKRt7YyB4PZ4x4Mpu+PBk4IEuIlIXfLmZ0gjWbap9rGlBJsxCp2HXhpjBxCnsPG2KLQB3tI\nyJl162IDOaWIkatsFVcujaDxvkuzVpgRWhlhY2xR6IM9JuTMshNTgC7SrBVXHPWmVKzMTh551nOK\nnC2Exc4hNsYWhT7YY0IOj9zbSHYZqsbILRTKB0vzyLUOlkj/nHfSD4cDIS8STo+RO72OlVuqH2oc\nLCHaKj704QoHekjBUW9DzOCROw+nZ60oZqaJkE7xSzHuxT7Qw3YNeGG8ttgeuYXG4cvYop/Dxtii\n0Ad7SMjdgHSqqYUbzuzUPtCj2FkrmjFyZnm2kOmHJ4pm5bAxtij0wRBydUxDedrnR0qm+FprE9pH\nvYlitYZoZq2gjO0QM3jkoBTZE+KFNwPsxAw1soU0D/QQZ60oU+w8chz1lhsIuQfQyp7I2Jq+3IpT\nfGbvZQtpfnjaQWPc7X544qi3ITbGFoU+GEIuRiN7YihGL7fixhSnZwtJ0A5nSXDDuMMjH2ZjbFHo\ngx0q5G5JhSt29oSF5sttSleoi5OxbG8mGUtyV6iwRSMgw4tZWoiRC3GqkGvWaJainT0hmuYLp5sm\nSDzDrq4QJ5PZmQbJZIy7ukJF66cKDt8soLmzU5JMgqwVIU4VcmZ7055io5k9YbukqkKuo2msNpmM\n8euv12bEfPh1qXFDGVun7+yER57Dxtii0Ac7WMiZ9Xa7maKZPSGiBLuPTLNWLPGOxyNFFXG1jSlD\nH65wsIRWHjmyVnIDIS8QJ3vkjsfhdbEHuxThcJg4Hi/ssALJB6j0M81W+qHSUW/ScdfY2WmBrJUh\nNsYWhT7YoULuhhg5SCP1DCUeucQ7bGtLh0yHOmyxWGGfaSJBEXiHmvn7Wjs7meGRj7Axtij0wQ4V\ncrdkrXgJzbMj7cTITb1D6//4W2+lHba33lIo3iQIaWnsqMXOziFmiJEDL6I5C7KbtWLqHcZizLW1\nzJFI+s+xTirTPupNq8aN5s5OtbNScUIQhBxk44Z1CQ2PXHNDUGtvLwc7O7M+QIOdndza25vXTiuP\n3PEL+xByCLkmDk8xzuDUTCFm/Ri5Bq09PbmFvKenxD1zCQitQMg1KUFGoDFO98g1s1Y0cfq4Ox4s\ndkLINVHbo6NcNMvJYCZUJiD9EGiisGte5IYiU6h0wCO3iVc88nvvvZcnTpzItbW1mXuhUIgnT57M\ngUCAA4EAb9++fVw64yW0TxVX3DWv3BiQojkT0n7fVfBSjHzXrl28b9++LCFvamrixx9/fNw74yWk\nebUSShKrVXH/9XB8BoUAzZmQ5vuuln7otaJZ0Wh0hJBv3Lhx3DvjNaQ73UxRj9V60CPXPurNyUjT\nD03fdxTNyo2qkE+dOpXnzZvHq1atyvmyQ8jTmNaecDxuSNUQolU+2OnYOVjC5H1H0azcqAl5b28v\np1IpTqVSvG7dOl61alXOzoRCocxXOBw27pzb0fLIVVF0/0txwIHWgR5OR1I0S/K+o2gWczgcztJK\nNSEv5O/K3SPXjBl6Fe0jxxztkZcg19GkaJad9x1Fs7JRE/ITJ05kvn/iiSf47rvvHpfOeAlPruKX\nAK3Dlx0fI1dO+zQdd+n7rlU0q/VQKwe3BrNi5MGtQW491Dq+DbFDY+QrVqzgSZMmcUVFBVdVVfHz\nzz/PDQ0NPHfuXJ43bx4vW7aMe3Js+S13IQfjh2k5VQmuyFox9A6lgqI1E9IsmtV6sJWDrwwT8leC\n3Howj5Cj1gqEHIwPWh65azCM10qm+JK1CYlHrv3haRxagZBDyIF9tMqpugZhvFZji75b1oSMFju9\nFFqRAiEHdtE6BFgT7cOXNbfoOz1Lq+wXOyVAyEGpcLKgaB4sUYpiZU7dN2Hr0GsUzfIGpchnBvZw\nqqAwCz1DAdrFypy8kzn0pxDH4tnb7WPxGIf+FMrfGDxy7wi5dj4zsIeTPXIL0cYUB+P02kLYoi/E\nS0LOjOwJt+D04k3Meh65Jm6o9mk87l4rmiXBa0LOrJPPDOyhKShqniHISdG36MMj956Qa3nkrtiY\nAjKYeoaahy87HTvvOrboD7Mxtij0wR4Scs0YueZWcYjK+OCx0uxqSN919ZkQsla8gXbWilbxJkzz\n7SN02MzxaNVJybuudrCE9WB45ECKVjlVLy68aSEMoao1Jk0/1N5R69jSwYiRQ8jtoF1O1UupcJqL\nnU4/mcmOoJiuCUmzhRztkSNrBUIuRbucqtc8cqdv0beNQtEsC9MsLdP8/dbWVg4Gg1nvejAY5NbW\nPBUJGXnkOW2MLQp9MIRchGbWilaM3K6HYopkQ5ArsoUUi2ZJs7RMdtRKhTz9s8haybIxtij0wRBy\nx6OVtRJLJrn29dczYj78uhiYbtH34sESzDJBkWZp4ai3YWbIWimAEhx9BeRY4h2JxwsWca3T3C00\n1ibENT8Ui2Z1hbo4Gcv+90nGktwVGn0GhaPehpnBIy8Q1ZQBkMHGB2gkHmcKhzlS4L+RxDO0GyMv\ndgZFLB7j2k21GTEffj2euCFrReuoN8TIhaiEVtSSeEEG4cst8ciZ9c6OZNbLFrLEO/JepGgibheN\ncZeGsyTJJKIwIk4IUoyRY1udPoYfoHZj5Bo1brRj5JH3IkxNxJH3IoUZlCCUaDLukpmQdIFZ4ktA\nyIXAI/c4Bh+gdrJWvFjjRuSRK24IYpaNu+nahGatFYRWhCBGrofTN6ZIcXodeMm424qRK20IsjPu\nJtlCdmdBppNxLHYKQNaKHqqfaYqNOf1kJslQiLNWLBQ2BGlmC0nXJaS+BNIPTR+MPHJV1KJM+ADN\nQjW6J2xMKigmaKYfSn2JsvbI7733Xp44cSLX1tZm7vX19fGNN97IM2bM4MWLF+f8FIWQ64N139Kg\nMu4247XF9sg10w8lvkTZx8h37drF+/btyxLyxsZGXr9+PTMzt7S08Jo1a8alM0AO1n31aWtLp70N\nHfdYrEiTE8UNQVprE46vvS+cfTq2aFY0Gs0S8lmzZnFPTw8zM3d3d/OsWbPGpTMaaJ8srgHWfUtD\nLMZcWzuY0zz8utRoZ60Y988N9W0EONIjZx4p5H6/P/N9KpXKurbTGQ3sDnKheLqcqgeRiIqqR14C\ncEatHO0Yue9zw7zEYjG67bbb6ODBg0RENGHCBIrH45m/v+yyy+jUqVNZNj6fj0KhUOZ6wYIFtGDB\ngrGaUiExMEDrolFqnDKFNhw7Rs01NeSvqBjXNgYGEhSNrqOammaqqPCPuAbOIpFI0Lp166i5uZn8\nfv+I63zEYkQ1NUTRKFF1tUp3i85AYoCi66I0pXEKHdtwjGqaa6jCP77/R7xO7OxZqnn9dYrOn0/V\nF1006s91dHRQR0dH5vqhhx6iAmQ5m0LUPldopbu7m5mZT5w44arQioXmSr5p8SZQGuzU/fDS2oTT\n8/fdgOOyVphHCnljYyO3tLQwM/Ojjz7qusVOO4Nsimk5VVBaTFLhVNcmPHpmpwSnhxIdGSNfsWIF\nT5o0iSsqKriqqopfeOEF7uvr40WLFrky/VArRs7sfI9cqx65WzD1yFUFBSvaGZxea8WxWSsSTDpj\ne7ebAVpZK244ckzrhCA34PiDJZi9GccRglorw2yMLQp9sEFnNGs0a6GZtWIHr53ZKcU1qXAO3fVV\nivcdtVaG2BhbFPpgw864oUazVxHVnwD6ONgj1yxjm/451FrJsjG2KPTBgs4Y12gGtoFH7hJcECM3\nXROShrNQayWHjbFFoQ+GR+54ECN3EU5P1fgc6aHXqLUyiGuF3IsxcjeglbXixbIIbkB73KVZWsU+\nJ5XZ+ycEnWe2fag4vLj/RWr7bhtN9U8lIqKp/qnU9t02enH/iyXumbe5ZeYt5L8we9ei/0I/3TLz\nlnFt5/ovfpHWRaOUGBggosGdtdd/8Yvj2g4RUV9fOw0MJLLuDQwkqK+vfdzbktD+TjslzmX3L3Eu\nQe3vjNG/9naiRLYdJRLp+6OgOe4ffPASHT36Y6qpaaaLLqqmmppmOnr0x/TBBy+NatPe3k7vvfce\nbdiwgaLRKG3YsIHee+89as/zO6XtjIfC+xhLf4EU8dGglAi9FK1NWE5P+xSHsxxexranp5U7O4NZ\n497ZGeSentZRbWKxGNfW1nLs8ypjw69H/Z0EQ4HQihAIuUexseimURaBWTbFlxbNkoStxQvMwoU3\n03GXbtGXnNkZi8WyYuSxWKxoWStY7BQAIfcwgpdbsywCs3zRzSSDwk4iiTjlU+GoN2Z5GVtJSQpp\njFySFYj0Q9MHQ8i9jcHLrVkWgVm+6KZVNEvLI2/t7eVgZ2fWuAc7O7m1t7eg5kzL2Gqd2SktHwyP\nXACEXA/1zDTDl1u6ii+J1dqNkUu8QxOHTTNG3trTk1vIPz8UJh+mHrlk3KV55JIDPRAjFwIh10N1\nr4hiY5JYrZ2t4hoeubiukOIis9a4S3d2Sjxy0biHQiM/HWKx9P08lG3RLGAftd3byu6/xpFjzHox\n8lJswjKN1WqVsbVb36boMyF45BDyUuDQekq20ThyzBVZKwK00g8lHrmdipPIWhlmY2xR6IMh5Kqo\neeTKaHnk2mgUKpN6hpLQinRtwk44Sy1bCFkrQAPVGLkiXj1yTMsjt7NVXPIBqrVFX3UmBI/cO0KO\no69Kg9PHXYKbCpVJQloaRbMkIEYuxEtCrukZuuaAAyDCLcfraXjkmqcyibJWXFQ0C0JeIE7OngBg\nPNGKkdtJPzTVV5FHDiH3npAz62RPMOtNNwHIhVbWihTpmpBxjByhFe8JuXb2hEaNZgDcijRLyzhr\nBYud3hFy7ewJeOQAjI3K4cuShiwzpB86C83sCcTIQU6E28W9iqmjbLvGjZc98qlTp/LcuXM5EAjw\n1772NdudAchaKRWOH3dBpSivHrEnCV2LsoXKJUZeXV3NfX1949YZ4E00qx9KkcyEJEkNttIPLfGO\nRMYu98dyQdFa7JR+eEomJ6JxL5eiWdXV1XxylP+MEHJgoblV3FY/DdcmSlI0KxJJx2sjkYJ+XKv6\noWYZW7VxLxePvKamhgOBAF9zzTX8i1/8wnZngHfROqnGbhqcabaQWvEmZmOP3EKy6Ka1RV+6sI+i\nWcNsjC2GcOLECWZm/uCDD7iuro537dqV1ZlQKJT5CofDdpoCTsFGPQDpKr5J/r6dgyWkoqJy5Jjk\nNAW2JygaW/SZcdRbOBzO0sqSZq00NTXxxo0bBx9c5h655gYJVWxONzXOjrTjGUqn+UX3DAXxWjtT\nfHjkNhpiF3nkn3zyCX/00UfMzPzxxx/zN77xDf7DH/5gqzNewu6RY47G8OWWCkpvay93BjuzYrWd\nwU7ubR37zElTz1Cy8FaSGLkB0kVmxMhtNMQui5F3dXVxXV0d19XV8Zw5c/iRRx6x3RmvIS3r6QoM\npptSQelp7ckp5D2t+c+clIy71sESmkWzpOPeFeriZCw7HJCMJbkrNHrWRVdXiJPJ7BlDMhnjrq7Q\nqDahUIhjw2YZsViMQ2NlhYTMk0lw1JsQCHkaSczQ8Uj3RwswneL39rZyZ2cwyzPs7Axyb29rXjvH\npx8KGnP6wRLwyHMDIXcYTvfINTdJ2MFk0a2nJ7eQ9/TkF3Jmh6cfenxtAjHyQSDkDsINMXLN0p5S\nNE+qYXZ4+qFQVDSyhTI2yFoZNEOtFffjlqwVzYOATbGzIUgiKI5OPxQ01nbyJMc+9wgtzzCWTDrq\nqDd45COBkAMRGgcBS5Bu0fdk+qGgsVgyybWvv55ZeBt+PRqIkdtoiBEjByXAyR65BGlIy/Hph4LG\npB65JGtFMgOVZq2oLTLjhCAIuRtoPdjKwVeCWaISfCXIrQfHXhh0KpohLadnrViYxmolHrlk3KUe\nOYQ8h40eRGfmAAAKo0lEQVSxRaEPhpA7ntZDrRzcOkzItwa59ZB7hRxko5W1Ip0JSWLkCK3ksDG2\nKPTBEHIR2nWxvRZaAYPYFRTTrBVptpAkawWLncNsjC0KfTCEXEQpTghy6mInsIedKb70jFrTbCE7\nxxoi/XCIjbFFoQ+GkIvRPLMTHjkYjjTt09Qjt+O0wCMfZmNsUeiDIeS2kG6SMMENG4KAPpK0T0mM\nXBpGRIw8h42xRaEPLraQK64oax45xmzukUu1VWuLvubZkdKsFccXzRI0pnnEnqRoluZRbyiaJaTo\nQq74aamVjpX+Ncynm8KhkGM43bTroZhgN3ui2OOumUeOollDbOCRy1AJrSjGr7TSsexONxUKEqYx\nXACSjLnmzk5me6lwTtzZySyP1WKLvo2GGDFycxRXlLXSsaQIh8Ic4cutsTHFQlo+WLI24dRaKxkT\nFM0atEHWiuGDy9wjt9CqR67mkducbjrVM0z/KvDILeCR22iI4ZEXjsNj5Mx6HnlrK3MwmD0UwWD6\nfj60ti1rbkxBjHyICWLkgzaIkcso96wVzXrkUiHXOj9Sc2OKJHuCWVbASS17QtiYNHsCR73ZaIiR\nteIptOuRS0MrxtNNxTQ4yeHL8MiHmMAjH7SBRy6j3IW8FEgXO40WgBQFRfPw5fSvghi5BWLkNhpi\nxMiBEDWPXNiYpqAwI2slywRZK4M2yFoxfDCEXA3hDNBejNzBggKPfIgJPPJBG3jkI9m+fTvPmjWL\np0+fzi0tLePSGSBDdYu+9XCHCgpi5ENMECMftEGMfCSffvopT5s2jaPRKPf393NdXR0fPnzYdme8\nSjgcLvhnteuRG2NTUMLhcFEFxU21Vqz3wgu1ViTjPnTMM2NRwLuOE4JGch4JeOONN2j69OlUXV1N\nFRUVtGLFCnrllVdy/mwgQOTz5X+e7yEfXfs/12bdu/Z/riXfQ3kMfT6ir389+97Xvz5mY76ODrph\nz56sezfs2UO+jo5RbTp8HbR34d6se3sX7qUOXx6bDh+9+eaSz79P/9ybby6hjo78/bv11luptraW\nEokEERElEgmqra2lW2+9dVQbn49o+fLse8uXFzbud718V9a9u16+K/+4T5hA9OGHRH5/+trvT19P\nmDC6ye7ddPpf/yJ/RQV1dHSQv6KCTv/rXzRh9+68/ds9YTcNnBugCn8FERFV+Cto4NwA7Z4wut3B\ng7dTV9fDWfe6uh6mgwdvz9vW8uXL6bHHHsu699hjj9Hy4QM7hH//d6Jnn82+9+yz6fujcefv7qTn\n3nyOiAbfi+fefI7u/N2deftH3/0u0W9+k33vN79J3x+Fezs7aeuHH2bd2/rhh3RvZ2fepo6sPkKn\ndp7Kundq5yk6svrIqDbvvruGTp/+S9a906f/Qu++u2ZUm8cee4wOHTpERINjcejQoRH/DsP53/8l\nev/97Hvvv5++Pxp7ju+h0+dOZ/fv3Gnac3zPKBZybrn8cvJXVGTd81dU0C2XXz7ubWUwln5m/t3v\nfsf33Xdf5vrXv/41//CHPxzxqVJXlw6j1tXlf941/30NUxPxNf99Tc7rnMyfn374/Pm5r0fh+jfe\nYAqH+fo33sh5nYs9C/ZwmMK8Z8GenNe52Lt3MYfDxHv3LuZQKJR1nY+bb76ZiYgnT57M0WiUJ0+e\nzETEN99886g2y5alf/Vly3Jfj8adL93J1ER850t35rzOyfe+l374976X+zoHqw4dYgqHedWhQxwK\nhbKu83H4B4c5TGE+/IPDOa9z8fbb/8HhMPHbb/9HzuvRWLt2LRMRr127Nud1Lh5+OP2rP/xw7utc\nbPjzBqYm4g1/3sChUCjrOi+bNqUfvmlT7uscvHj8OFM4zC8eP57zejR6Xu7hMIW55+WenNe5+PDD\nbRwOn8cffrgt53UuIpEIX3DBBRyJRDgUCmVd5+Ott5gvuyz9Z67rXMTiMa7dVJvJJR9+nRPp4pNN\nJLIsEvLf//73BQl5ISJuYYn3BQ9fMLaIW1jifdFFBYm4hSXeXwiHxxRxC0u8w5eGxxRxC0u8V668\noCARt7DE3PrKJ+IWlnhPnFiYiFtY4v3lDV8eW8QtLPG+6qoxRdzCEu8Jq1YVJOIWlnj/39X/N6aI\nW1ji/de/BgoScQtLvOfPnz+miFtY4v1v/za2iFtY4j3929MLE3ELS7y/850xRdzCEu97//73gkTc\nwhLvI/cfGVPELSzxPnr0oTFF3MIS75tuuqkgEbewxPvll8cWcQtLvCPvRcYWcQtpOpgN1IT8L3/5\nCy9dujRz/cgjj4xY8CSaliVG+MIXvvCFr7G/pk2bZqzJvrTomvHpp5/SrFmz6LXXXqMvf/nLdN11\n19Fvf/tbuvrqq00fBQAAwCbni4zOP59+9rOf0dKlS+mzzz6jYDAIEQcAgBIh8sgBAAA4B1H64Vjs\n2LGDvvKVr9CMGTNo/fr1xWjCNVRXV9O8efOovr6errvuulJ3R5VVq1ZRZWUlzZ07N3Pv1KlTtHjx\nYpo5cyYtWbIkk2bpdXKNRVNTE1VVVVF9fT3V19fTjh07SthDPY4dO0YLFy6kOXPmUG1tLT3zzDNE\nVJ7vxmhjYfxuSBY781HIZqFyorq6mvv6+krdjZKwa9cu3rdvH9fW1mbuNTY28vr165mZuaWlhdes\nWVOq7qmSayyampr48ccfL2GvSkN3dzfv37+fmZnPnDnDM2fO5MOHD5fluzHaWJi+G+PukZtsFioX\nuEyjV9/85jdpwrCNQq+++iqtXLmSiIhWrlxJW7duLUXX1Mk1FkTl+W5ceeWVFAgEiIjokksuoauv\nvpqOHz9elu/GaGNBZPZujLuQHz9+nKZMmZK5rqqqynSsHPH5fHTjjTfStddeS7/85S9L3Z2S09vb\nS5WVlUREVFlZSb29vSXuUWl59tlnqa6ujoLBYFmEEoYTi8Vo//79NH/+/LJ/N6yx+PrnO9ZN3o1x\nF3LfWPvCy4zdu3fT/v37afv27bRp0yaKRCKl7pJj8Pl8Zf2+/OAHP6BoNEoHDhygSZMm0Y9+9KNS\nd0mVjz/+mO644w56+umn6dJLL836u3J7Nz7++GP6zne+Q08//TRdcsklxu/GuAv55MmT6dixY5nr\nY8eOUVVV1Xg34xomTZpERERXXHEFffvb36Y33nijxD0qLZWVldTT00NERN3d3TRx4sQS96h0TJw4\nMSNY9913X1m9GwMDA3THHXdQQ0NDpp5Nub4b1ljcc889mbEwfTfGXcivvfZa+sc//kGxWIz6+/vp\npZdeottvz1+wyKskk0k6c+YMERF98skn9Mc//jEra6Ecuf3222nz5s1ERLR58+a8Ram8Tnd3d+b7\nLVu2lM27wcwUDAZp9uzZdP/992ful+O7MdpYGL8bRViI5W3btvHMmTN52rRp/MgjjxSjCVfQ1dXF\ndXV1XFdXx3PmzCm7sVixYgVPmjSJKyoquKqqil944QXu6+vjRYsW8YwZM3jx4sUFHyTgdoaPxfPP\nP88NDQ08d+5cnjdvHi9btox7esauZ+IFIpEI+3w+rqur40AgwIFAgLdv316W70ausdi2bZvxu4EN\nQQAA4HKKsiEIAACAHhByAABwORByAABwORByAABwORByAABwORByAABwORByAABwORByAABwOf8f\nSUMc7mhho+8AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 110 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }