{
 "metadata": {
  "name": "",
  "signature": "sha256:c563b8805cda58c0871b8f3830643023e3e1269b903decbb7ad8c3c86a233497"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "The Matrix Exponential"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook gives a brief introduction to the matrix exponential,\n",
      "an interesting function which can be found within many areas of applied mathematics.\n",
      "After defining the function we will introduce some applications and algorithms before giving links\n",
      "to currently available software and further information."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Definition"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import scipy as sp\n",
      "import scipy.linalg as la\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The matrix exponential is a direct generalisation of the corresponding scalar function.\n",
      "Remember that the Taylor series expansion for $\\exp(x)$ is\n",
      "$$\\exp(x) = \\sum_{k = 0}^\\infty \\frac{x^k}{k!}.$$\n",
      "We can apply this to square matrices using matrix multiplcation:\n",
      "$$\\exp(A) = \\sum_{k = 0}^\\infty \\frac{A^k}{k!}.$$\n",
      "\n",
      "Two of the most useful properties of the matrix exponential (see [11] for more) are that\n",
      "$$\\exp((a + b)A) = \\exp(aA)\\exp(bA),$$\n",
      "and when two matrices $A$ and $B$ satisfy $AB = BA$,\n",
      "$$\\exp(A + B) = \\exp(A)\\exp(B).$$\n",
      "\n",
      "A thorough treatment of the theory behind the matrix exponential\n",
      "(and matrix functions in general) can be found in [11]."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Applications"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The matrix exponential has a wide variety of applications in areas including\n",
      "cancer research [10] (which also uses the derivative of the exponential),\n",
      "nuclear burnup equations [12], and computer graphics [13].\n",
      "We'll cover two here: solving systems of linear differential equations and \n",
      "finding the important nodes in a complex network."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Differential Equations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One major attraction of the matrix exponential is the fact that it gives the explicit solution to the differential equation\n",
      "$$ y'(t) = A y(t),\\qquad y(0) = y_0,$$\n",
      "for $y, y_0 \\in \\mathbb{C}^n$ and $A \\in \\mathbb{C}^{n \\times n}$.\n",
      "The solution is\n",
      "$$ y(t) = \\exp(At) y_0. $$\n",
      "\n",
      "Knowing this we can compute the solution $y(t)$ at any time point directly\n",
      "without using a time-stepping scheme such as the Euler method, for example.\n",
      "Let's do a simple example showing this for the differential equation\n",
      "$$ y'(t) = \\begin{bmatrix}1 & -20\\\\3 & 4\\end{bmatrix} y(t),\\qquad y(0) = \\begin{bmatrix}1\\\\0\\end{bmatrix},\\qquad t \\in [0, 1].$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.integrate as sint\n",
      "\n",
      "A = np.array([[1, -20],\n",
      "              [3, 4]], dtype=np.double)\n",
      "\n",
      "def yprime(y, t):\n",
      "    return sp.dot(A, y)\n",
      "\n",
      "t = np.linspace(0, 1, 51)\n",
      "yzero = np.array([1, 0])\n",
      "\n",
      "# Solve equation using odeint\n",
      "y = sint.odeint(yprime, yzero, t) # Calls the Fortran library odeint\n",
      "plt.plot(y[:, 0], y[:, 1]) # Blue line is y(t) for t in [0, 1]\n",
      "\n",
      "plt.xlim([-4, 9])\n",
      "plt.ylim([-3, 5])\n",
      "\n",
      "# Solve for 10 t values using the exponential\n",
      "tvals = np.linspace(0, 1, 11)\n",
      "\n",
      "for tval in tvals:\n",
      "    yval = sp.dot(la.expm(A*tval), yzero)\n",
      "    plt.plot(yval[0], yval[1], 'rs') # Red squares are evaluated via the matrix exponential\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEACAYAAAB4ayemAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH0lJREFUeJzt3Xl4k2W+xvFvWiiLiICCgCytgGdgQMAFN7BBFuEAMhxl\nnFFB3D06ox4cFWd0CGeOgvs4LqOICoh43GZU8CigkqKCoAJaWQZZqggIsokgUKE5f/xaU6ALkLd5\n8r65P9eVq2kSkl+g3H3yrCAiIiIiIiIiIiIiIiIiIiIiIiIiPhDy4DkKgG3AXuAnoIsHzykiIlVk\nFdDAdREiIukgw6Pn8aLFLiIilfAitGPAO8AnwFUePJ+IiFShJsVfGwILgW4OaxERCbRqHjzHuuKv\n3wH/xAYi3y+5s1WrVrEVK1Z48DIiImnlM6DT/jcm2j1SGziy+PoRQG8gv/QDVqxYQSwWC8Rl5MiR\nzmvQ+wn+e9H7Se1Lst4L0LGs0E20pX0s1rouea7ngekJPqeIiJQj0dBeRRnNdxERqRpeTflLC+Fw\n2HUJngrS+wnSewG9n1Tm+r0kY351rLh/RkREDlIoFIIyMtqL2SOSpiLDhkFBwYF3ZGcTGT8+ydWI\npAeFthy+ggIieXkH3BxJfiUiaUN92iIiPqLQFs9t3AjbtrmuQiSYFNriuTVr4LjjIByGMWNg4ULQ\nWLSINxTa4rmOHeHbb+GWWyzABw+Gpk3hssvgxRdh82bXFYr4lwYi5fBlZ5c96JidzRFHQL9+dgFY\nvhymTYNJk+Cqq+CXv4Q+fexyyimQmZnEukV8TPO0Jel274YPPoC337bLunXQq5cF+LnnQuPGrisU\nca+8edoKbXHum2/iAf7uu5CTE2+Fn3EGVK/uukKR5FNoiy/89BPMnRsP8eXL4Zxz4iHeooXrCkWS\nQ6EtvrRhA0yfbgE+bRo0bBgP8LPPhpo1XVcoUjUU2uJ7RUUwf368Ff7559C1qwV4377QujWEdFqp\nBIRCWwJnyxZ45514iNesGW+Fd+8Odeq4rlDk8Cm0JdBiMfjii3iAz5sHXbrEQ7x9e7XCxV8U2pJW\ntm+HmTMtwN96y6YZlnSj9OwJ9eq5rlCkYgptSVuxGHz5ZbwV/v77tmpz0CC46CJo0sR1hSIHUmiL\nFNu5E2bNsiX1//wnnHYaDB0Kv/oV1K7tujoRU9WhnQl8AnwDDNjvPoW2pKwff4TXX4eJE+Gjj6z1\nPWQI5OZChnbmEYeqOrSHAycDRwLn7XefQlt8Yd06eOEFC/AtW+CSSyzAf/EL15VJOiovtL1oSzQD\n/h0YV9YLiPhFkyYwfLhtJTtlChQW2mrMLl3gkUfgu+9cVyjiTci+DNwN1AX+gLpHJED27LH9UJ57\nDqZOtW6ToUOhf3+oUcN1dRJkVdXS7g9sABaU9eQifletmu08OGkSrF5tfd6PP277g197LcyerQMe\nJLkSDdq7gSHAHqAm1tp+FRha6jGxkSNH/vxNOBwmHA4n+LIibq1eDc8/b/3fhYXW933JJdCqlevK\nxK+i0SjRaPTn70eNGgVVPOUvF3WPSJqJxeDTT6375IUX4IQTrPtk8GCoX991deJnVTkQWZrSWdJK\nKGQn7zz8sB2tdtttth9KTo4F95Qptt2siFe0uEakCmzZAi+/bC3wf/0LfvMba4GffLL2QJGDoxWR\nIo6sWGEDmc89Z6fwDB1q/d/Nm7uuTFKZQlvEsVgM5syxwcuXX4ZOnWwA8/zz4cgjXVcnqUahLZJC\ndu2CN9+0AM/Ls1Prhw61HQh1Mr2AQlskZX33nW1eNXGiHXJ80UVw9dU2E0XSV7Jmj4jIIWrYEH73\nOzu44b33ICvLjlG74AL45BPX1UmqUUtbJAXt2AHjxsEDD8C//RuMGGH7oGjmSfpQ94iIDxUWwuTJ\ncM89Nlg5YoTt+61tY4NPoS3iY0VFtu/36NGwbZst4rn4YutKkWBSaIsEQCxmZ1+OGQNLlthWsldd\npZPng0gDkSIBEApZ3/b06fDaa7bLYE4ORCKwaZPr6iQZFNoiPnXyybZI58MPbd+TNm3gv/7LdiCU\n4FJoi/jcCSfAU09Bfr4tzOnYES6/HJYudV2ZVAWFtkhAHHcc3H8/LF9uXSa5ufAf/2HzvyU4NBAp\nElA7dsAzz1iQt25t0wV79tRcb7/Q7BGRNPXTT3ZAwz33QK1aFt6DBmmPk1Sn0BZJc0VFdijD6NGw\neTPceqvtMqgDilOTQltEAJvrnZdnc72/+MJmnFx9tbaHTTWapy0igPVph8Pw9tvW8p43D44/Hv78\nZ9txUFKbQlskjXXubNvCzp4N69fb5lQ33ghff+26MilPoqFdE5gLLAQWA6MTrkhEkq5NG3jySesu\nqVHDwvzmm+GHH1xXJvtLNLR3Ad2BTsCJxde7JlqUiLjRtCnce68tzNm8Gdq1s1WXGpZKHV4ORNYG\n8oBLsVZ3CQ1EivjU++/DdddZmD/6qLXIJTmqciAyA+seWQ/MZN/AFhEf69YN5s+H3r3hjDNg5EjY\nudN1VemtmgfPUYR1jxwFTAPCQLT0AyKRyM/Xw+Ew4XDYg5cVkWSoXt36ty+80KYHdugAjzwCffu6\nrixYotEo0Wi00sd5PU/7TmAncH+p29Q9IhIgb79tZ1p26gQPPQTNm7uuKJiqqnvkGKBe8fVaQC9g\nQYLPKSIprE8fm2XSoYPNMrn/flsqL8mRaEu7AzABC/8M4Dngvv0eo5a2SEAtX26t7m++gb//3frA\nxRtaxi4iVSIWg3/8A266CXr0sCmDjRq5rsr/tIxdRKpEKATnnw+LF8Mxx0D79vDEE7B3r+vKgkkt\nbRHxVH4+/Od/QmGhdZmcfLLrivxJLW0RSYoOHWDWLFuU06+f9Xlv3eq6quBQaIuI5zIyYNgw6zLZ\nswfatoVJk7Qc3gvqHhGRKjd3rnWZ1K0Ljz9ue5pIxdQ9IiLOnHYafPwxXHCBHTg8YoSdYSmHTqEt\nIkmRmWn92/n5Nq+7XTt47TV1mRwqdY+IiBMzZ9pgZYcOdmp8nTquK0ot6h4RkZTSvTssXGhnU551\nFhQUuK7IHxTaIuJMjRowbhxccQWcfjocxCZ3aU+hLSJOhUJwww02JfDCC21BjpRPfdoikjKWL4eB\nA23jqb/9DbKyXFfkjvq0RSTltW4Nc+bA2rXQsyd8953rilKPQltEUkrdujYV8Oyz4dRTbbBS4tQ9\nIiIp68UXbW7344/D4MGuq0ku7actIr40fz4MGgSXXgqRiO1rkg4U2iLiW+vX257dDRvCxIk2tzvo\nNBApIr517LHw3nt2yMKZZ8LKla4rckehLSK+kJUFY8fCNddYcL/3nuuK3Eg0tJsDM4FFwBfADQlX\nJCJSjlDIBiYnT4bf/hYefTT9NpxKtE+7cfFlIVAH+BT4FbCk1GPUpy0inlu5Es47D844Ax57LHgL\ncZI1EPka8AjwbqnbFNoiUiV++AGGDoU1Hw6jZ+uCA4M7O5vI+PEuSktYeaFdzcPXyAY6A3M9fE4R\nkXIdeSS8+ipcllPA3XPyDrg/kvySqpxXoV0HeAW4Edi+/52RSOTn6+FwmHA47NHLiki6y8iAnBzg\na9eVJCYajRI9iG0OvegeqQ5MBd4C/lrG/eoeEZEqFQmHieSV0dLOzSXi0/1eq2qedgh4GlhM2YEt\nIiIeSjS0zwIuAboDC4ovfRItSkTEC3v2uK7Ae4n2aX+AFuiIiGvZ2fsMOsZisGwZrC7IZu9eO1Q4\nKLT3iIgEUmEh9OkDnTrBgw+6rubQae8REUkrWVk2HfDNN+HJJ11X4x0v52mLiKSU+vVh6lTo2hVa\ntbLTcPxOLW0RCbQ2beCll+Cii2DJksofn+oU2iISeLm5cO+90L8/bNzouprEaCBSRNLG7bfDBx/A\nO+9AjRquq6mYTq4RkbRXVGRnTR5xBEyYYFu9pirNHhGRtJeRAc89B4sXw+jRrqs5PJo9IiJppXZt\neOMNOP10G6T02ynv6h4RkbS0YAH07m3zuLt0cV3NgdQ9IiJSSufO8PTTMGgQfO2jbV3VPSIiaeu8\n8+DLL2HAAPjoI6hVy3VFlVNLW0TS2vDh0LIlPPWU60oOjvq0RSTtffopDBwIK1akzvxt9WmLiJTj\n5JPhxBPh2WddV1I5tbRFRIA5c+C3v7U+7urVXVejlraISIXOOANat4ZJk1xXUjG1tEVEiuXlwZVX\n2m6A1RzPrVNLW0SkEmefDY0bw4svuq6kfF6E9jPAeiDfg+cSEXEmFII774S77rLNpVKRF6H9LDqB\nXUQColcvOPJIO6osFXkR2u8DWzx4HhER50pa2//zP3aqe6pRn7aIyH769bNtXKdMcV3JgZIyPhqJ\nRH6+Hg6HCYfDyXhZEZHDEgrBHXfAX/5i+5Ik47CEaDRKNBqtvDaPXi8bmAJ0KOM+TfkTEd8pKrJV\nkvffD30cjNppyp+IyCHIyIA//cla26nU7vQitF8AZgMnAKuByzx4ThER5379a/juOziIXouk0YpI\nEZEKjB0LM2bAyy8n93XVPSIichhOPRWWLXNdRZxa2iIiFdi6FVq0gO+/T84skhLltbR13FiCIsOG\nQUHBgXdkZxMZPz7J1cSlal0iflOvHmRmwqZNcMwxrqtRaCeuoIBIXt4BN0dKXd+5E5Yvh7VrYfNm\n+8fftGnf61u2wJ49NkpdVGRfS1/PzIS6deGoow681KsHTZvakUktW9r3B1OXiBycnBxYtUqhHWjL\nltkeBsuWwYYN9o/erBkcfTQ0aGBfjz/e+ssaNID69W3j9VDILhkZ8euhkAX6tm32Ea30ZfNmWLnS\nVm599ZVdMjKg6x7XfwMiwVES2qee6roShfZh27MHPvsMVq8u+/7ateHmm+GEE6z1m5mZnLpiMWu1\n39ULmH/g/XPnwuDB0KkTnHQSdO1qm+OISPmOP95COxUotA/B55/D1Kkwa5YdTdSsGXTcWfZjmzVz\ntYrKWu7lBXH79tB5kP3CufdeC/AOHeCcc6B7dzjzTPuFIyJxOTn2fyYVaMpfJb76CkaPtmAbMMAm\n2l97rZ3avGiRtaT95Igj4KKL4J57YOZM2LgR7r7bulQiEWjUCHr0gCeesG4dEYl3j6QCtbTLsGmT\nTaR//nlYuhQuuAD+/ndrhWbs/2suO7vswb3s7Cqvs0IHWVfNmtbC7t7dvt++Pb6QYMQIOOUUWxU2\naBA0bFjFNYukqFQKbc3TLmX9ehgzBiZMgN694eKL4dxzISvLdWVu7NwJb78NL71kX/v1gxtugC5d\nXFcmkly7dtmsrB07kjc+Vd48bYU2NgPjvvtsueoll8Dtt9s5cRK3dSs8/TQ8+qj93dxwAyx5axgZ\nXxcc+GDNBZcAOu44G8tq0SI5r6fFNWXYtg0eeggeeQTOPx8WLoTmzV1XlZrq1bPZMDfdZIOxDz4I\nR3xcwP/t1FxwSQ8lXSTJCu3ypO1A5Ecf2eDiypUwbx48+aQC+2BkZsLAgbbrWatWrqsRSZ5U6ddO\nu5Z2LGYt67vugnHjbEaIHLpQyBYIiaQLhbYD27bBlVfadL05c2zCvHhv927XFYh4b9cuqFPHdRVp\n1D1SsgS1fn348EMFdlWaP9/2WhEJklmz4OyzXVeRJi3tHTusH/aaa2D4cNfVBEg5c8Eb7M6mRw/I\ny3M/XV3EC9u3Q34+nHaa60oCGtqltyWNxWDxYmiXCd9/lg2Md1dYwFQ0re9vf4P+/WHBAtsIS8TP\nPvoIOneGWrVcVxLQ0C53W9KvHNSSpn7/e1uQ89e/wi23uK5GJDF5eZCb67oK40Wfdh9gKfAlcJsH\nzycBEArFZ+kUFrquRiQxqdKfDYm3tDOBR4GewBrgY+ANYEmCzys+V9JF1W0v3HaaHdYAaLWk+M6u\nXfDpp7b3UCpINLS7AMuBguLv/xcYiEJbSndRLYzfHHFSjMjhmzcP2rVLnX3nE+0eOQ4ofQzAN8W3\niYgEQip1jUDiLe2D2gkqEon8fD0cDhMOhxN82YrFWmZzbpYtU99nMrzmn4nIIcrLsw3Sqlo0GiUa\njVb6uERDew1QeseO5lhrex+lQzsZel89ntcWwrRPkvqyIhIwP/1k0/1eeqnqX2v/Bu2oUaPKfFyi\n3SOfAG2AbCALuBAbiHTq++/tdHIRkUTMn28bo9Wv77qSuERb2nuA3wHTsJkkT5MCg5CZmVBU5LqK\nNJeqJ/qIHIK8vNTqzwZvFte8VXxJGRkZCm3XNK1P/G7vXnjlFfjjH11Xsq9AbhiVkWF/4SIih+ux\nx6BGDTjvPNeV7CuQy9irVbMJ8SIih2PVKvjv/4bZs8s4zNuxQJ4RuX07NGsGy5ZBo0ZJfWkR8blY\nDHr1ssO9b73VXR3lnRGZYr9DvFGnjn2kmTzZdSUi4jdPP20z0FJ1G+dAtrQB3n0X/vAH2xpURORg\nrFkDnTrBe+/Z4jyX0qqlDRAOw8aN8PnnrisRET+IxeDaa+H6690HdkUCG9qZmTBkCEyc6LoSEfGD\nF16ws1NSbYrf/gLbPQI2Atyli02Qb9fOSQki4gMbNsCJJ8KUKXaWbCpIu+4RsCPvx4yB3/xGUwBF\npHw33ABDh6ZOYFck0C1te3EL7UaN7CQVEZHS58hu3AgrVsApp0Bmq9Q5pKO8lnYgF9eUFgrBk0/a\noZy9eqXe6iYRcaCsc2Q/gEimm3IORaC7R0rUqwfPPw9XX21TekRE/CotQhvsfLff/966SnbscF2N\niLjk58Om0ya0AUaMgDZtoEcP2LTJdTUi4sLy5f5edJdWoZ2ZaUtUw2Ho1g1Wr670j4hIgMybZ//3\nmzev/LGpKvADkfsLhWwaYKNGcNZZMG0atG3ruioRqWpTp8Jll8Ezz8Cnr2YTKet0Kx8c0hH4KX8V\nee45uPvqYfRtW0DduvvdmZ06U39EJDFjx8LIkfDaa3Daaa6rOThpO+WvIkOGwLx7C3hwQd4B90WS\nX46IeCwWs7CePBlmzbIxLb9L69AGOPpo1xWISFX46Seb5rtokR1mEJS99RMZiBwMLAL2Aid5U46I\nSOJ++AEGDLA9RWbODE5gQ2KhnQ8MAmZ5VEtK+fprnTMp4kfffmszxJo3h9dfhyOOcF2RtxIJ7aXA\nMq8KSTWbN9vskqVLXVciIgdr4UJbSDdwoA0+VgtgB3AA39Ihys4uc9CxY8tszj3d5nTeeqsdPZTp\ng30JRNLRpk1w553w6qvwwANwySWuK6o6lYX2DKBxGbf/EZhysC8SiUR+vh4OhwmHwwf7R6tcZdP6\n+vSBK66wH4Znn9WcbpFUsmePtagjEfj1r2HJEmjQwHVVhycajRKNRit9nBfztGcCNwPzy7k/Zedp\nH6yiInjiCftNfsEFcMcd/l5RJRIE0ajtg3300fDww3aIQZBU9SEIyVik40xGBlx3HSxbZr/FO3Wy\nzafWrXNdmUj6+eora1UPGwZ//rMdwhu0wK5IIqE9CFgNnA68CbzlSUUp7OijYfRo+wiWlQXt28PN\nN9u0IhGpWjt3wqhRcNJJ8MtfwuLF9sk3FOgm44ESCe1/As2BWli/d19PKvKBRo1ssCM/37Z4bNvW\ndhDUzoEi3ovF4JVX7P/ZokUwf76tcqxd23VlbqT13iNeWb0a7rrLfrCuu85mmtSr57oqEf/Lz4cb\nb7QjwR5+GLp3d11R8qTlwb7J0ry5DVR+/DF8843tbxCJ6JQckcO1ebONG/XoAeefb63rdArsiii0\nPZSTY9s+zp4N69dDhw52JuXUqVpdKXIwNm6E++6zrpC9e2386Prrg7lI5nCpe6QK7dgBL75o80jX\nrLH53pdfDi1auK5MJHXEYpCXZwdwv/WWrWYcPhw6dnRdmVvldY8otJPk88/hqadsi8jTT7fdx/r1\nUwtC0temTTBhgjVqMjPhmmtsu+T69V1XlhoU2inixx9twHLsWFi1ylreV1zhiwMzRBIWi8H771ur\n+s03rfvw6qttn590m7pXGYV2Clq0yFrfkybBKafAlVdC377B25VMZPNmmDjRGiuxmLWqhw7175Lz\nZFBop7Bdu2xvk/HjYe5cyM21fr0BA+DYY11XJ3J4YjH48EML6jfegP79rVXdrZta1QdDoe0TW7fa\nYMzrr9uhw7/4hQX4wIF2XT/skuq2bLHzV8eOtdNjrrkGLr1Up0QdKoW2DxUW2qY4b7xhIV6zZjzA\nzzxTW8VK6ti+3WaAvPyy/az27Wut6txcNTQOl0Lb52IxWLDA/kO8/rpNIezXzwK8d2/1g0tyFRXZ\ngQPTpsH06fDJJ3DqqdYFMmQINGzoukL/U2gHTEEBTJliAT5vnrVouneHrl2hc2eoXt11hRI0a9fC\njBkW1DNmWHfHuedaoyE3F+rUcV1hsCi0A2zLFvuPlJdnAz8rV1qr56yz7HLGGdoLRQ7dzp02PW/6\ndPv5WrMGeva0kO7VC1q2dF1hsCm008jWrTBnjgX4Bx/YR9fjj7dW+Fln2dcWLeJ9jZFhw6zpvr/s\n7EpP9pHgiMXgiy/iIT1njq1KLGlNn3KKxlGSSaGdxgoLrf/xgw/iQV69ejzAC8aHeWB+3gF/LpKb\nS+Qgjj8S/9qwAd55x4J6+nSoVSse0t27w1FHua4wfZUX2lpEnQaysqBLF7sMH24tqhUr4gG+ZknZ\nf06bXAXH7t2wdKltdZqfb9sq5OfbrI9w2IL6zjuhVSvXlUpl1NIWIuEwkbwDW9rnhHJZ1TJKu3bs\nc2nbFurWdVBomjqU7qtYzPZ3LwnlkoBescJ2oTzxRNt9skMHu96ypabkpSq1tOWQdesGT46zY50W\nL7az+B591FpsDRpwQJi3a6fNfqpEQUGZv1Tv2GOflkoHdH6+nehSEsp9+sCtt9rCrJo1HdQunksk\ntO8D+gOFwArgMuB7L4qS1BAK2YEObdrYfPASRUV2uGpJmM+eDeOKw71OnXiA5+RA06bQpIl9bdpU\n88kPRSwGP/xgm4yVZc4cmDE83noePNi+HnNMcuuU5EoktKcDtwFFwBjgdmCEF0VJkmVnEynn9rJk\nZFgg5+TYAp8SsZid3FMS5l9/bXuprF0bv2Rl7RviJZfStzVpEuzz/4qKbFvSdev2vXz77YG3hULQ\nrZyxha5dYdSBDXAJuERCe0ap63OB8xOsRRzxalpfKGRHrzVvbgNb+4vF4Pvv9w3xtWtti9oPP4x/\nv26dfZTfP9SPOcbCvHZtm+VQ1vXS39esab9gvBSL2X4au3fveyks3Pf7LVsODOCSy4YNNibQpIld\nGje2rzk5tj1Bye1Nmtgnl0gYKCOc1Rednrzq074ceMGj55KACoVskU+9etZ9Up5YzEKvdLCXtER3\n7rTLjz/GL6W/L319924L7vJCvXZtm/q4f+BWFMiFhTZXuUaNii/16sWDt3Nn24ujdEhnZSXv712C\npbLQngE0LuP2PwJTiq//CevXnuxhXZLGQiEb6GzQANq3P/znKSqybW8rCvfCQgvQykK45JKV5X3r\nvVKH2H0lwZboB6xhwFVAD2BXOY+JjRw58udvwuEw4XA4wZcVEQmWaDRKtNRitlGjRoHHKyL7AA8A\nucDGCh6nedoiIoeoKpaxfwlkAZuLv58DXFfG4xTaIiKHSHuPiIj4SHmhnewhFRERSYBCW0TERxTa\nIiI+otAWEfERhbaIiI8otEVEfEShLSLiIwptEREfUWiLiPiIQltExEcU2iIiPqLQFhHxEYW2iIiP\nKLRFRHxEoS0i4iMKbRERH1Foi4j4iEJbRMRHFNoiIj6SSGj/BfgMWAi8CzT3pCIRESlXIqF9L9AR\n6AS8Boz0pKIUFo1GXZfgqSC9nyC9F9D7SWWu30siof1Dqet1gI0J1pLyXP9jeS1I7ydI7wX0flKZ\n6/dSLcE/fxcwBPgROD3xckREpCKVtbRnAPllXAYU3/8noAUwHnioakoUEZESIY+epwXwf0D7Mu5b\nDrTy6HVERNLFZ9iY4T4S6R5pA3xZfH0gsKCcx7VO4DVERMQjr2BdJQuBV4FGbssREREREUlTNwNF\nQAPXhSToPmAJ1t/0D+Aot+Uclj7AUqx76zbHtSSqOTATWAR8AdzgthxPZGLdjVNcF+KBetin8iXA\nYvw/y+x27GctH5gM1HBbTtVpDrwNrML/od2L+KybMcUXP8nEBoezgepY91ZblwUlqDHxwZo6wL/w\n9/sBGA48D7zhuhAPTAAuL75eDX82ckpkAyuJB/WLwKXJLiJZe488CNyapNeqajOwTwwAc4FmDms5\nHF2w0C4AfgL+FxtI9qtvsV88ANuxFl1Td+UkrBnw78A4vJvd5cpRQDfgmeLv9wDfuysnYduw/zO1\nsV9AtYE1yS4iGaE9EPgG+DwJr5Vsl2NTHf3kOGB1qe+/Kb4tCLKBztgvU796CLiFeMPAz3KA74Bn\ngfnAU1jQ+dVm4AHga2AtsBV4J9lFeBXa5S3COQ/rAyq9L4kfWg+VLSoCW1hUiPVr+UnMdQFVpA7W\nd3oj1uL2o/7ABqw/2w//TypTDTgJeLz46w5ghNOKEtMKuAlrHDTFfuYudllQVWgPrMf6sldhHy0K\n8P/0wGHAh0BNx3UcjtOx8YUSt+P/wcjqwDTsP5Sf3Y19CloFrMNCbqLTihLTGHsvJboCUx3V4oUL\nsW6rEkOAxxzVkjRBGIjsg40eH+O6kMNUDViBtRay8P9AZAgLtqBto5BLMGaPzAJOKL4eAe5xV0rC\nOmIzlGphP3cTgOudVpQEK/F/aH8JfIV9hF2AffTzm77YLIvlWEvbz7pi/b8Lif+b9HFakTdyCcbs\nkY7Ax/h7imxptxKf8jcB+5QnIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIsnw/zAxi5dqB+EQAAAA\nAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3a8f710>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For our relatively simple ODE both methods perform well and,\n",
      "since the matrix exponential is more expensive to compute,\n",
      "we would probably opt to use <b>scipy.integrate.odeint</b>.\n",
      "For solving more complex differential equations the exponential forms the basis of a class of methods called\n",
      "<b>exponential integrators</b> [3].\n",
      "Some more detail on these methods can be found in [7], for example.\n",
      "Note that in practice we would calculate $\\exp(A)b$ without first computing $\\exp(A)$ as this is cheaper.\n",
      "We could compute $\\exp(A)b$ directly by using <b>scipy.sparse.linalg.expm_multiply</b>, which implements the recent algorithm from [3]."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Important Nodes in a Network"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Another recent application of the exponential is in finding the most influential nodes of a network.\n",
      "A network is a set of <b>nodes</b> connected by <b>edges</b>.\n",
      "We can represent a network by its <b>adjacency matrix</b>:\n",
      "we set $A_{ij} = 1$ if node $i$ is connected to node $j$ and $A_{ij} = 0$ otherwise.\n",
      "For example take the following network."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import networkx as nx\n",
      "\n",
      "Adj = np.array([[0, 1, 1, 0, 0, 0], \n",
      "                [1, 0, 1, 1, 1, 0], \n",
      "                [1, 1, 0, 0, 0, 0], \n",
      "                [0, 1, 0, 0, 0, 1], \n",
      "                [0, 1, 0, 0, 0, 0],\n",
      "                [0, 0, 0, 1, 0, 0]])\n",
      "\n",
      "G = nx.from_numpy_matrix(Adj)\n",
      "nx.draw(G, node_color='y', node_size=1000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFBCAYAAAA2bKVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFXe//F3QjpBkCIIUhQk9BZaqAGkl9C7C6wIIj0o\noquui6urIoSSUKQqHQMoTRSp0kPoCAFCSRDpNQlpM/P7g0f2x4qamczck0k+r+vielYezvl+uS7N\nZ8495z7HzWKxWBARERHDuDu7ARERkZxG4SsiImIwha+IiIjBFL4iIiIGU/iKiIgYTOErIiJiMIWv\niIiIwRS+IiIiBlP4ioiIGEzhKyIiYjCFr4iIiMEUviIiIgZT+IqIiBhM4SsiImIwha+IiIjBFL4i\nIiIGU/iKiIgYTOErIiJiMIWviIiIwRS+IiIiBlP4ioiIGEzhKyIiYjCFr4iIiMEUviIiIgZT+IqI\niBhM4SsiImIwha+IiIjBFL4iIiIGU/iKiIgYTOErIiJiMIWviIiIwRS+IiIiBlP4ioiIGEzhKyIi\nYjCFr4iIiMEUviIiIgZT+IqIiBhM4SsiImIwha+IiIjBFL4iIiIGU/iKiIgYTOErIiJiMIWviIiI\nwRS+IiIiBlP4ioiIGMzD2Q2IGM1isXDp0iWio6O5cuUKqampeHl5UbRoUQIDAylatChubm7OblNE\nsjGFr+QIFouF7du3ExExka1bt+Pmlk5AgCcFC6bh4WEmPd2d69c9OXUqlVy5vGjatAnDhr1B/fr1\nFcQiYnduFovF4uwmRBzFYrGwbNkyPvjgLdLTb9OuXSINGlgoWBCelKkWC1y7Bj/95Ma6dX7kzl2I\nDz+cSOfOnY1vXkSyLYWvZFtXrlxh0KB+/PzzLl5/PZHq1Z8cuH/EYoGoKJg+3Y9atZoxY8Y8ChYs\n6LiGRSTH0IYryZZ2795NlSoB5Mu3lenTE6lRw7rghYd/vnZtmDUrCXf376lU6UUOHDjgmIZFJEfR\nyleyne3bt9OlS1vGjk2kdm37zbtzJ0ye7M/69T9Sp04d+00sIjmOwleylRMnTtC4cV3efjuBwED7\nz79nD0yalIfdu6N58cUX7V9ARHIEha9kG2lpadSqVZmXXjpNu3aO+9d65Up3Dhyows6dB8iVK5fD\n6ohI9qXvfCXb+OSTj/H1jadtW8d+nuzUyUxy8hmmTAlzaB0Ryb608pVs4cqVK5Qr9zxffJHMM884\nvt4vv8DQob6cP/8LTz/9tOMLiki2opWvZAtz5nxBo0YYErwAxYpBnTpuLFgw35iCIpKtaOUrLi89\nPZ1SpYrwwQc3KVPGuLrHjsHkyc9y5swl3N31OVZEMk7HS4rL279/P/7+qTYF70cfwcGDkJwMefNC\nmzbQt2/GxlaqBCbTfY4ePUq1atWsLy4iOZY+rovLO3DgAOXLp9k0tndvWLoU1q+HTz6BVatg//6M\njXVzg/LlLURHR9tUW0RyLoWvuLz9+3dQpkyyTWOffx68vP77z7lyQb58GR9funQiUVG7bKotIjmX\nHjuLyzt16gR169o+PiwMfvgBUlNh5EgoWzbjY59/Hr799ojtxUUkR9LKV1xeUtIDfHxsHz96NGzY\nABMnwty5cPJkxsf6+EBSUpLtxUUkR1L4istzd8/8fbtublCtGgQHw5YtGR9nNkOuXPrPSESso58a\n4vKeeuop7t+3z1zp6Vi1ik5IgDx5nrJPcRHJMRS+4vKqVq3J2bPWj7tz5+Eq98EDMJke7nLevh3q\n18/4HGfPulGtWia+cBaRHEnhKy6vdu0GxMbmtnqcmxusWQPdu0NICMyfD2+/DeXKZXyO2Fh/atcO\nsrq2iORs2u0sLq9OnTq8844Fk+nhq0IZlTcvTJ5se930dDh+PJ3a9rw0WERyBK18xeVVqFCBZ54p\nSlSUsXV37oTy5ctTqlQpYwuLiMtT+IpLO3XqFK+99hpnzvxCZKSxD3LWrvVn+PC3DK0pItmDwldc\njsViYfPmzbRt25bGjRtTuHBhTpw4waVL/hw7ZkwPBw7AtWs+dOzY0ZiCIpKt6DtfcRkpKSksXbqU\nsLAw0tPTGT16NJGRkfj6+gIwY8Y8Ro3qy6xZSZk6dOOvJCbCpEl+zJ+/CK///2xKEZEM0pWCkuVd\nv36dmTNnMn36dKpWrcro0aNp0aIFbm6/P1yjR4+OpKd/x/DhqQ7pxWKBzz7zomjR7sydu9AhNUQk\n+9NjZ8myfv75ZwYNGkTZsmWJi4vjxx9/ZOPGjbRs2fKJwQsPV78nTjzL0qX2f6hjscC8ebnYsSOd\n5s3b2n1+Eck5FL6SpVgsFn744QdatWpF06ZNee6554iJiWH27NlUrFjxL8fnz5+fzZt3sXlzEebN\n88Rksk9fJhPMnOlFdHQJVq5cz+jRo5k5c6Z9JheRHEePnSVLSE5OZvHixYSFheHm5kZoaCi9evXC\nx8Yvb69du0anTq1JTIzhjTcSee4523u7cAEmTMhN4cJViIxcR/78+Tl79iytW7eme/fu/Pvf//7D\nlbhIdnf58mXOnz9PcnIyHh4e5MuXj/Lly2s/xF9Q+IpTXb16lRkzZjBjxgxq1qzJ6NGjadasmV3C\nzGw2M3XqZMaPf5f27VNp185EoULW9AZr13rw3XfefPTRBAYPfu2xvq5fv0779u0pW7Ysc+bM0Q8b\nyRESExNZvHgxq1cv5tChI6SkJFOihA9eXg8vGrlzx8zly8mUK1eKBg2aMmjQUCpXruzstrMcha84\nxfHjxwkLC2PVqlX06NGDkSNHUr58eYfUio2N5fPPP2bp0qXUqOFOUFAiZcvCc889fiKWyQRxcXD6\nNOzenZujRy307fsyb7zxNiVLlnzi3ElJSfTq1YukpCRWrlzJU0/pkgXJnq5cucJHH33AokULqVzZ\nnSZNEihfHgoXfnhU6//vwQM4exaio3Px3XdelC4dwLhx/6JDhw7OaT4LUviKYcxmM99//z1hYWEc\nP36coUOHMnjwYAoWLGhI/fv377No0SI2b15LdPRBrl+/RaFCPnh4uJGWZuH69WSKFClAjRqBtGzZ\nkV69epE791+fGZ2ens6IESPYtWsXGzZsoFixYgb8bUSMYbFYWLJkCaNGDaFp02Q6d06jcOGMj09P\nh127YMGC3AQGNmHGjHkUsuYRVDal8BWHe/DgAQsXLmTy5Ml4eXkRGhpKjx498Pb2dmpft2/f5tq1\na6SkpODt7U2RIkXImzevTXNZLBY+/fRTZsyYwYYNGzK0OUwkq0tOTqZv324cOrSVsWMTCQiwfa6U\nFFiwwIstW3xZvvwbgoOD7danK1L4isNcuXKFiIgIZs2aRZ06dQgNDSU4ODhbb05atGgRoaGhfP31\n1zRu3NjZ7YjYLCkpiTZtmuLpeZQ333yAvbY0REfDxx/7sXjxSlq1amWfSV2QXjUSuzty5Aj9+/en\nfPny3Lx5k59++om1a9fSpEmTbB28AH379mXp0qV069aN5cuXO7sdEZuYTCY6d26Dj88Rxo2zX/AC\nBAbC+PFJ9OnThd27d9tvYhej4yXFLsxmMxs2bCAsLIxTp04xbNgwzp49S4ECBZzdmuGaNWvGjz/+\nSNu2bbl06RKhoaHZ/kOHZC+TJk3g2rUDfPZZslXXdGZUxYowZkwSvXp14sSJWPz9/e1fJIvTY2fJ\nlKSkJL788ksmT56Mv78/oaGhdOvWTa/dAPHx8bRu3ZpmzZoxadIkcjnip5iInZ06dYp69WoQHv6A\nokUdW2vCBB+KF+/FzJnzHFsoC1L4ik0uX75MeHg4s2fPpn79+oSGhtKwYUOt8P7HnTt36NSpEwUK\nFGDhwoWPLoEQyaqaNg2icuX9dOpkdnithAQYONCPdeu2U7NmTYfXy0r0na9Y5eDBg7z88stUqlSJ\n+/fvs2fPHr755hsaNWqk4H2CfPnysXHjRjw9PWnevDk3b950dksif+jnn3/m2LEjtGvn+OAF8PeH\nTp2SmTLlM0PqZSUKX/lLZrOZNWvWEBwcTEhICJUrVyY2NpZp06ZRpkwZZ7eX5Xl7e7N48WLq1atH\n/fr1OX/+vLNbEnmiiIjJtGmThqencTVbtTKzZs3aHPfBVOErfyghIYHw8HACAgL48MMPGTx4MOfO\nnWPs2LE8/fTTzm7Ppbi7u/PZZ58xdOhQGjRoQHR0tLNbEnmM2Wxm8eLFtG6dbvMcly5Bixbw8ccZ\nH5M3L9Sp405kZKTNdV2Rwld+59KlS4wbN45SpUqxdetW5s+fz/79++nVqxeeRn4kzoaGDx9OeHg4\nrVu3ZuPGjc5uR+SR2NhYfH2hSBHb55gyBWw5JbZChST27t1ue2EXpPCVRw4cOEDv3r2pUqUKycnJ\n7N+/n5UrV9KgQQN9n2tHnTp14ptvvqF///7Mm5fzdnlK1hQdHU1AgO078rdsefgdbo0a1o8tWxYO\nHNhnc21XpPd8DWCxWDh69Ch79+5l//6dHDlygHv37mEymfH19eHFFwOoXbsxtWrVomHDhoYeu2gy\nmVizZg2TJk0iLi6OESNGMGPGDJuPWZSMqVevHtu3b6d169bEx8fz/vvv6wOOONXx48coVSrBprGJ\nibBgAUyaBOvWWT++dGmIibmIxWLJMf8dKHwdKDExkaVLlzJt2qfcvPkrVataKFMmiX79IE8ecHeH\n5GS4ePEChw9vYdkyXy5ftvDKK4MZMmTYH96kYw/3799n/vz5TJkyhUKFChEaGkrnzp3x8NC/EkYJ\nCAhg9+7dtGvXjvj4eGbMmKHH+uI09+/fJndu2948nTcP2rSBggV/f8NRRnh7PxyXmprq9DPfjaKf\ntA6yevVqhgz5O2XLptOnTwI1az4M2ycpWxaaN08D0oiPh3XrplK9ejivvPIaH374H5svlH+SuLg4\npk2bxrx582jWrBmLFi0iKCjIbvOLdYoUKcK2bdvo3r07ISEhrFixIkee9iPOZ+uRD2fPwsGDMHv2\nb/PYWv/hpq+cQt/52tnt27fp3j2EUaP68o9/3GH8+ARq1/7j4P1fxYvDkCGpzJ2bzIEDs6hSpSxR\nUVGZ7mvfvn307NmT6tWrYzabiY6OZsWKFQreLMDf3581a9ZQtGhRGjduzJUrV5zdkuRAuXM/RVKS\n9eMOH4YrV6BHD+jSBVasgB07YPDgjM+RlgZms8WuC42sTuFrR5cvX6ZevRqkp29k1qwkKle2fa6n\nn4Z//vMBPXvG07p1MBs2bLB6jvT0dCIjI6lfvz49e/akbt26nD9/nokTJ1KqVCnbmxO78/DwYPbs\n2YSEhFCvXj1iYmKc3ZLkMBUrViYuzvqnLu3bw5IlMGfOw9Vvhw5Qty58ZsW5GefPw4svFssx3/eC\nHjvbzfXr1wkODqJx48v06WP7e3L/q2lTKFIkiZdf7srSpd/QokWLvxxz79495s6dy5QpUyhWrBih\noaGEhITo+9wszs3Njffff5/ixYvTuHFjVq1aRb169ZzdluQQgYGBvP229c+Mvb0f/vqNry94eT18\nfzejYmIgMLC21bVdmVa+dmA2m+nWrT21a/9q1+D9TYUK8MEHD+jVqzMXLlz4wz934cIFQkNDef75\n59m3bx/Lly9n165ddOnSRcHrQgYMGMCXX35Jx44dWbVqlbPbkRyibNmy3L1rIrMHTfXrB++8Y92Y\nkyd9qVMnZ91/rfC1gxkzIrh16zgDBqQ5rEblytC1azIDBvT63caIPXv20K1bNwIDA8mVKxeHDh1i\n2bJl1KlTx2H9iGO1bNmSjRs3Mnz4cKZNm+bsdiQHyJUrF127dmbjRmNv30pIgF27zHTt2tXQus6m\n8M2kixcv8v7743jjjUSH3Hv5/+ve3cT168f44otZpKens3z5curWrUvfvn1p2LAhFy5cYMKECZQo\nUcKxjYghatSowa5du4iIiGDs2LE5aieoOMewYWNYt84bk8m4mj/84EaLFs159tlnjSuaBehKwUwa\nOfJ1btyYzauv2v9x85OcOgXvvpsbP7/8lCpVitDQUNq3b6+7YrOxW7du0aFDB4oXL86CBQtyzHuQ\n4hw1apSnefNTtG7t+FoPHsCrr/qxdOlGGjZs6PiCWYhWvpmQmJjIV199Sfv2xgQvQLlykC9fCqNH\nj2bHjh107NhRwZvN5c+fn02bNpGWlkbLli25c+eOs1uSbCg5OZkPP/yQc+d+ZeZMD27ccHzNefO8\naNSoVY4LXlD4ZsqKFSuoVMnN6oPIR42Cli0fngjTps3DDQrW6No1ne+//8a6QeLSfH19Wb58OdWq\nVaNBgwbEx8c7uyXJJiwWC2vXrqVixYocPnyYI0eOMGLEm4SF+eHIbzoOH4adO/0ID5/tuCJZmB47\nZ0KfPl0oWnQVbdtaN270aGje/GHw2uLePejTx5u7d5Nwz+jpHZJtTJo0ibCwMNavX0+VKlWc3Y64\nsLNnzzJy5EhiY2OZOnXqo1cZU1NTady4Ds899zOvv55q05GRf+b8eRg71pfFizP2+mR2pJ/cmRAd\nHUXZsraNzcxHnqeegnz5PDl9+rTtk4jLCg0N5fPPP+ell15i8+bNzm5HXFBiYiLvvvsudevWJTg4\nmKNHjz4Wgl5eXnz33VbOnn2BadPsuwErJgbeesuXqVPn5NjgBYWvzZKSkrh48VdsPShqzhzo2BGG\nD3/4+MVaZcu66UL2HKxHjx5ERkbSu3dvFi1a5Ox2xEVYLBYiIyOpUKEC586d48iRI7z55pt4eXn9\n7s/my5ePbdv2cvt2dcaMyU1mv+kwmWDp0ly8805uZs9eSq9evTM3oYvTyQs2unHjBk8/7Y2np/Wb\nrQYNglKlwNMTNm+Gf/zj4bFsRYtmfI4CBVK4evWq1bUl+2jUqBFbtmyhTZs2xMfHM27cuBx1PJ9Y\n5+TJkwwfPpyrV6/y1Vdf0bjxXx9qkTdvXjZv3kV4+FRGjvwHXbsm0769mTx5Ml7XYoEjR2DOnNwU\nKlSZ6OilOt4WrXxtlpqaioeHbT/oypd/eASbh8fDjVeVKsE+K++R9vAwkZqaalN9yT4qVqzInj17\nWL58OUOHDsVk5Aua4hLu3bvHG2+8QaNGjejQoQOHDh3KUPD+xt3dnREjRhEVdZR79zrQt68Pn3/u\nw6FDD+/xfRKzGeLiYPVqNwYO9GfGjOKMHj2ZLVt2K3j/j1a+NvLy8iItzXl71dLScul9TwGgaNGi\n7Nixgy5dutC5c2eWLl2Kn5+fs9sSJ7NYLCxZsoSxY8fSokULjh8/TuHChW2er3Tp0ixbtppr164x\ne/YslixZyokTZylc2JvnnnPHx8eCyeTGnTtunD79gPz581KvXj3mzBlJcHCwnsr8D+12tlFycjJP\nP52Hb79N5wlfl/yhhAT4+WeoVg1y5YItW2DSpIffARcrlvF5PvjgKV5/fTbdu3e3vnnJllJTU3n1\n1VeJiYlh7dq1FCpUyNktiZMcOXKEYcOGkZSURHh4uMOuDk1PT+fkyZPExsaSnJyMp6cn+fLlo1q1\nahQoUMAhNbMLhW8mVKr0AkOHnqd8+YyPuXsXxo17+EjG3R1KloQBAyAw0LravXr5sXPnUUqXLm3d\nQMnWLBYL7733HsuXL2fjxo369yOHuX37Nu+//z4rVqxg/PjxDBw4UIfwZFF67JwJNWvWISbGuvDN\nmxdmzMhc3du34cEDCy+88ELmJpJsx83NjX//+98UL16cBg0a8O2331K7ds66qi0nMpvNzJ8/n3/8\n4x907NiRn3/+WSvPLE4brjKhVasQdu+2/vLpzNqxA5o0aaTvUOQPDR48mC+++IK2bduydu1aZ7cj\nDhQVFUVQUBBz5sxh/fr1zJw5U8HrAhS+mdCpUycuXHAnLs64mhYLrFvnz/DhY40rKi6pffv2rF+/\nnkGDBjFr1ixntyN2duPGDQYNGkSHDh0YMmQIu3btItDa76/EaRS+meDt7c0rrwzm22+t2HGVSQ8P\n5MhLkyZNDKsprqt27dr89NNPfP7557z77ru/uwtaXI/JZGL69OlUqFABPz8/Tp48Sf/+/XXUrIvR\nhqtMunr1KpUqvci//32fgADH1kpNhdde8+Pjj+fSs2dPxxaTbOX69eu0b9+egIAAZs+e/cQTjSTr\n27VrF8OGDeOpp55i2rRpOtvbhemjUiYVLlyYyZNnMGFCbhx95sWXX3pRpUojevTo4dhCku0UKlSI\nLVu2cOfOHdq2bcu9e/ec3ZJY4cqVK/Tr148ePXowduxYtm3bpuB1cQpfO+jduzeVKjVk6lSfTF2Y\n8Gd274Yff/Rl1qwF2mglNvHz82PlypW8+OKLNGrUiMuXLzu7JfkLaWlphIWFUalSJYoUKcLJkyfp\n1auXfgZkAwpfO3Bzc2PRoq+5erUM4eHedr8Dc+9emDQpN+vWbcrUCTUiHh4eRERE0LNnT+rVq8eJ\nEyec3ZL8ga1bt1KtWjW+++47du7cyaeffkoeaw5VlixN3/na0Z07d2jdugk+PjGMGvWAvHkzN5/Z\nDKtWubN8eW7Wrv2BunXr2qdREWDRokWMGTOGFStWWHXWrzjWpUuXGDNmDHv37iUsLIxOnTpppZsN\naeVrR/ny5WPLlt1UrTqAgQN92bHD9rl++QVCQ3MTHV2F3bujFbxid3379mXJkiV069aNFStWOLud\nHC8lJYVPPvmEatWqERAQwMmTJ+ncubOCN5vSytdBdu3aRf/+PfHyuku7dvdp0gR8fP58zG9Xb61f\n70dUlIV//vNDRowYpePhxKGOHj1K27ZtGT16NKNHj9YPeyfYuHEjI0aMoFy5coSFhelY0BxA4etA\nJpOJ77//nqlTP2Xv3n2UL+9N6dKJlCljIk+eh2c7JyfDxYsQG+vPyZMW/PyeZtiwN+nXrx95M/vc\nWiSD4uPjad26NS+99BITJ07UBz6DnD9/ntGjR3P8+HGmTJlC27Ztnd2SGETha5ArV64QFRVFVNQ+\nDh7czf37dzGZTPj4+FK+fFVq1apHYGAgFSpU0MpDnOLOnTt06tSJAgUKsHDhQnx9fZ3dUrb14MED\nPvvsM6ZOnUpoaChjxozB568ejUm2ovAVkUdSUlLo378/8fHxfPvttzoj2M4sFgtr1qxh1KhR1KxZ\nk4kTJ1KiRAlntyVOoPAVkceYzWbGjRvHmjVr2LhxI6VKlXJ2S9nCmTNnGDFiBBcuXGDq1Kk0b97c\n2S2JE2m3s4g8xt3dnc8++4yhQ4dSv359Dh486OyWXFpiYiLvvPMOQUFBNGvWjCNHjih4ReErIk82\nfPhwwsPDadWqFd9//72z23E5FouFFStWUL58eeLi4jh69ChvvPGGztUWQI+dReQv7N69m86dO/Of\n//yHAQMGOLsdl3DixAmGDx/OjRs3CA8Pp1GjRs5uSbIYrXxF5E/Vq1eP7du38+GHHzJ+/HhdS/gn\n7t27R2hoKMHBwXTq1ImDBw8qeOWJFL4i8pcCAgLYvXs3a9asYdCgQaSnpzu7pSzFYrGwcOFCypUr\nx927dx+tfD08PJzdmmRReuwsIhmWkJBA9+7dcXNzY/ny5fj7+zu7Jac7fPgww4YNIzk5mYiICOrU\nqePslsQFaOUrIhnm7+/Pt99+y7PPPktwcDBXr151dktOc+vWLYYOHUrLli3529/+xr59+xS8kmEK\nXxGxiqenJ7Nnz6ZDhw4EBQURExPj7JYMZTabmTNnDhUqVMBisfDzzz8zaNAgHckpVtEXEiJiNTc3\nN95//32KFy9O48aNWbVqFfXq1XN2Ww63f/9+hg0bhoeHBxs2bKBGjRrObklclFa+ImKzAQMGsGDB\nAjp27Mjq1aud3Y7DXL9+nYEDBxISEsLQoUPZuXOnglcyReErIpnSqlUrNm7cyLBhwwgPD3d2O3aV\nnp5OREQEFStWJE+ePJw6dYp+/frh7q4fnZI5euwsIplWo0YNdu3aRatWrYiPj+c///mPywfUzp07\nGTZsGPny5WPLli1UqlTJ2S1JNqJXjUTEbm7dukWHDh0oUaIE8+fPx9vb29ktWe3XX39l7NixbNu2\njc8///zRq1Ui9uTaH01FJEvJnz8/mzZtIjU1lVatWnHnzh2b57p8+TI//fQTmzZtYsuWLRw6dIiU\nlBQ7dvu4tLQ0Jk6cSOXKlSlWrBgnT56kR48eCl5xCK18RcTuTCYTY8aMYfPmzWzYsIHixYv/5ZjE\nxESWLFnCypULOXToKKmpKZQs6Y23N5jNcOeOmUuXknnxxRI0aBDMq68OpXr16nbpd/PmzQwfPpzi\nxYszdepUAgIC7DKvyB9R+IqIw0yaNImwsDDWr19PlSpVnvhnrl69ykcffcDChV9RubI7TZokUL48\nFC4M/7voTEmB2FiIjs7Fhg3eFC/+PGPHfkCXLl1sWqHGx8czZswY9u/fz+TJkwkJCdFKVwyh8BUR\nh1q+fDnDhw9n2bJlNG3a9NHv/3bl3vDhgwgOfkDnzmkUKZLxeU0m2LMHFizITfny9Zg9+yuKZHCC\nlJQUJk6cyMSJExk+fDhjx47Fz8/P2r+aiM0UviLicDt27KBbt25MmjSJPn36kJycTL9+PTlw4Efe\neCOR8uVtnzs1FRYu9GTjRh+WLFn5lxfVb9iwgZEjR1KhQgXCwsJ44YUXbC8uYiOFr4gY4sSJE7Rp\n04aBAweydet3uLkd5q23HmCvu+WPHoV//cuX+fOX0aFDh9/9/8+dO8eoUaM4efIkU6ZMoU2bNvYp\nLGIDha+IGObSpUtUr16eChUe8P77Jux9HHJMDLzzji8rV24gODgYgKSkJD799FPCw8N54403CA0N\ndclXoCR70SEbImKYlSu/5tlnzbz3nv2DFyAgAMaNe0CfPl04fvws27ZtY/To0dSuXZtDhw5RokQJ\n+xcVsYFWviJiiLNnz1K7dhWmTHlABt48ypSwMC8OH85P7txPM23aNJo1a+bYgiJWUviKiCFatWrM\nCy/spHt3s8NrJSVB//6efP31Jho3buzweiLW0glXIuJwZ86cISpqPx07Oj54Afz8oEePdMLDPzek\nnoi1FL4i4nDTp0+lZUuT3XY2Z0TLlhZ++GETV65cMa6oSAYpfEXEoSwWC4sWfUnbtmlWjVu9GgYP\nhhYt4NORuWlsAAAZIklEQVRPra/r7w/167uxYsUK6weLOJjCV0QcKi4uDkinaFHrxhUsCC+/DK1b\n2167YsVk9u7davsEIg6i8BURh4qOjqZcOc/fndP8Vxo2hAYNIG9e22sHBEB09AHbJxBxEIWviDjU\n8ePHKVkywebxmXkfo1QpOH/+Mmlp1j3yFnE0ha+IONT9+3fJndv2Xc6ZuWTIwwO8vHKRnJxs+yQi\nDqDwFRGHyuxRAvY4icBsNuYVJ5GMUviKiEP5+z9FUpLty9fMrHxNJkhJMeHr62v7JCIOoPAVEYeq\nWLEicXH+Vo8zmR5eF2gyPf6/rXHxIpQoURgvI18wFskAXawgIg4VGBhITIyVqQksXAhfffXff/7x\nR+jX7+GvjDp9GmrWrGl1bRFHU/iKiEM9//zzpKS4ceUKFCmS8XH9+z/8lRk//+z96GpBkaxEj51F\nxKHc3Nzo1asnGzca+1n/wQPYsQO6detmaF2RjFD4iojDDR06mg0bPDHyddtNm6BRo0YUd/T9hSI2\nUPiKiMNVrFiRChUq8913mdi6bIWUFFi1KjcjR75lSD0Ra+k+XxExxIkTJ2jYsBbTpz+w6rtfW8yc\n6UV6egsiI9c6tpCIjbTyFRFDVKxYkTFj3iYszM/qV4ascewYbN3qy4wZ8xxXRCSTFL4iYpi33nob\nX9/KTJ3qbZeTq/7XhQswfrwv8+YtplChQvYvIGInCl8RMYyHhwdr127i11/LMmmSj11XwGfOwNix\nvkyaNJO2bdvab2IRB1D4ioih8uTJw5Ytu0lNrc3Ikbm5cCFz85lMsHSpOyNGuDFhwgxefvlvdulT\nxJG04UpEnMJisTBr1kz+8Y83CQlJJiTEZNXdvRYLHD0Kc+fmJl++CpQpU4nr16/z7bff4u6udYVk\nbQpfEXGqCxcu8O67bxIZuZLgYC9eeimFgADIk+f3f9ZshkuXIDrajfXrc+Puno8333yPV14ZSHp6\nOsHBwYSEhPDWW3rFSLI2HS8pIk5VqlQpevbsx88/x9K0aQ8iIxdx/PgZ8uf3pEQJd7y8LJjNbty7\n58bp0w8oUCAvQUH1mD17JMHBwbj937VHXl5eLF++nFq1alGnTh0dKylZmla+IuJ0rVq1onfv3vzt\nbw+/rzWZTMTExBAbG8uDBw/w8PDg6aefpmrVquTPn/9P5/rhhx8YMGAABw4c4NlnnzWifRGrKXxF\nxKlOnz5NgwYNiIuLw8fHxy5z/utf/2Lr1q38+OOPeHjoAZ9kPdqVICJONX36dAYOHGi34AV49913\n8fb25r333rPbnCL2pJWviDhNQkICJUuW5NChQ5QoUcKuc1+/fp3AwEAiIiJo3769XecWySytfEXE\naRYvXkyjRo3sHrwAhQoVYvny5QwcOJBz587ZfX6RzFD4iohTWCwWwsPDGTZsmMNqBAUF8c4779Ct\nWzeSk5MdVkfEWgpfEXGKHTt2kJ6eTtOmTR1aZ8SIEbzwwguMGjXKoXVErKHwFRGn+G3V+9t7uo7i\n5ubG3Llz2bp1KwsXLnRoLZGM0oYrETHcL7/8QuXKlbl48SJ5nnSUlQMcO3aMpk2bsnXrVipVqmRI\nTZE/opWviBhu1qxZ9O7d27DgBahcuTITJ06ka9eu3L9/37C6Ik+ila+IGColJYWSJUuydetWypcv\nb3j9QYMGcffuXZYtW+bwR94if0QrXxEx1MqVK6lUqZJTghdg6tSpnDlzhoiICKfUFwGtfEXEYPXq\n1WPs2LF07NjRaT3ExsYSFBTE2rVrqVOnjtP6kJxLK18RMczBgwf55ZdfaNeunVP7KF26NLNnz6Z7\n9+7cvHnTqb1IzqTwFRHDREREMGTIkCxx2UFISAg9evSgb9++mM1mZ7cjOYweO4uIIW7evEmZMmU4\nffo0hQoVcnY7AKSlpdGsWTNatGjBu+++6+x2JAfRyldEDDFv3jw6dOiQZYIXwNPTk2XLljF9+nR+\n/PFHZ7cjOYhWviLicCaTiTJlyrBixQpq1arl7HZ+Z+vWrfTu3ZsDBw5QrFgxZ7cjOYBWviLicBs2\nbOCZZ57JksEL0KRJE4YPH06PHj1IS0tzdjuSAyh8RcThHH17kT2MGzeOvHnz8vbbbzu7FckB9NhZ\nRBzq9OnTNGzYkIsXL+Lj4+Psdv7UrVu3qFGjBpMmTaJz587ObkeyMa18RcShpk+fziuvvJLlgxcg\nf/78fP3117z22mucOXPG2e1INqaVr4g4TEJCAiVLluTQoUOUKFHC2e1k2PTp05k1axZ79+7F19fX\n2e1INqSVr4g4zKJFi2jcuLFLBS/AkCFDqFixYpb/nlpcl8JXRBzCYrG4xEarJ3Fzc+OLL75gz549\nzJs3z9ntSDbk/DPeRCRb2rFjB2azmSZNmji7FZv4+/sTGRlJ48aNCQwMpGrVqs5uSbIRrXxFxCHC\nw8MZOnSoS9+ZW6FCBaZMmULXrl25e/eus9uRbEQbrkTE7i5dukSVKlW4ePEiefLkcXY7mfb6669z\n9epVIiMjXfrDhGQdWvmKiN3NmjWLPn36ZIvgBQgLCyMuLo7Jkyc7uxXJJrTyFRG7SklJoWTJkmzb\nto1y5co5ux27uXDhAnXq1GHVqlXUr1/f2e2Ii9PKV0TsKjIyksqVK2er4AUoVaoUc+fOpWfPnly7\nds3Z7YiLU/iKiF1FRES45OtFGdGuXTtefvll+vTpg8lkcnY74sIUviJiN9HR0fzyyy+0a9fO2a04\nzPjx40lPT2f8+PHObkVcmMJXROwmIiKCIUOGkCtXLme34jAeHh4sXbqUOXPmsHHjRme3Iy5KG65E\nxC5u3rxJmTJlOH36NIUKFXJ2Ow63Y8cOunfvzv79+13u+ExxPq18RcQu5s6dS0hISI4IXoBGjRoR\nGhpK9+7dSU1NdXY74mK08hWRTDOZTJQpU4avv/6amjVrOrsdw1gsFjp16kTJkiWZMmWKs9sRF6KV\nr4hk2oYNG3jmmWdyVPDCwwsYFixYwLp161ixYoWz2xEXovAVkUxz1duL7CFfvnx8/fXXDB06lFOn\nTjm7HXEReuwsIpkSExNDo0aNuHjxIj4+Ps5ux2lmz57NlClT2LdvH7lz53Z2O5LFKXxFJFNGjhyJ\nv78/H330kbNbcSqLxUL//v0xm8189dVXuoBB/pTCV0RslpCQQMmSJTl8+DDFixd3djtOl5iYSN26\ndRk+fDiDBg1ydjuShXk4uwERcV2LFi0iODhYwft/cufOTWRkJA0aNCAwMJDAwEBntyRZlDZciYhN\nLBYL4eHhDB061NmtZCkBAQFERETQrVs3bt++7ex2JItS+IqITbZv347ZbKZJkybObiXL6d69O+3b\nt6dfv36YzWZntyNZkMJXRGzy2+tF2lj0ZBMmTODatWt8/vnnzm5FsiBtuBIRq8XHx1O1alUuXrxI\nnjx5nN1OlhUXF0ft2rVZvnw5jRs3dnY7koVo5SsiVps1axZ9+/ZV8P6FEiVKsGDBAnr37s2VK1ec\n3Y5kIVr5iohVUlJSKFmyJNu2baNcuXLObscl/POf/2THjh1s2rQJD4/fv2SSlJREdHQ00dHRREX9\nxJUrv5CWloq3tw/PPfc8tWs3oGbNmlStWhUvLy8n/A3E3hS+ImKVxYsXs2DBAjZt2uTsVlyGyWSi\nVatW1KpVi48//vjR7x87doyIiDCWLVtG8eKelCmTQpkyKRQqBLlyQVoaXLkCsbF+nD7twa1bFv7+\n90EMGTKMUqVKOe8vJJmm8BURqwQFBTFu3DhCQkKc3YpLuX79OjVq1GDGjBnUqFGDwYP7sX//Ltq2\nTaN163QychNjXBysW+fFpk3udO7cjbCwcJ566inHNy92p/AVkQw7cOAAXbp04dy5c+TKlcvZ7bic\nnTt30qZNGzw8TLRvn0bfvml4elo/T0ICfPGFD4cO+TNv3hKaN29u/2bFobThSkQyLCIigtdff13B\nawOLxcLq1SvImzeFTz5JYsAA24IXwN8fQkOTGTHiBn37hhAREW7fZsXhtPIVkQy5efMmZcqU4cyZ\nMxQsWNDZ7bgUi8VCaOgIfvhhHv/5TxL2fFJ8+TK8+aYf77zzKUOH5sxrHV2RVr4ikiFz584lJCRE\nwWuDWbNmsm7dfLsHL0DRojBhQhLjx7+lTXAuRCtfEflLJpOJ0qVLExkZSc2aNZ3djku5cOECgYGV\nmDQpkZIlHVdn/36YNq0gJ07EahOWC9DKV0T+0vr16ylSpIiC10oP7/jtSdeuyQ4NXoDataFatQRC\nQ4c7tpDYhcJXRP7Sb+c4i3W2b99OXNxxunc3GVJv8OBkVq5cwcWLFw2pJ7ZT+IrIn4qJieHIkSN0\n69bN2a24nGnTJtC+fRJGbQ7394eXXjIzc2aEMQXFZgpfEflT06dP59VXX8Xb29vZrbiUK1eusHnz\nFpo3t25bzb178N570Lo19OwJmzdbV7d9+1Tmzp1FWlqadQPFUL8/ZFRE5P/cv3+fhQsXcuTIEWe3\n4nK2bdtG9eqe+PsnWzVuyhTw8oLVq+HMGXj7bShdGjJ6mmSJEpA3r4WjR48SGBhofeNiCK18ReQP\nLVq0iCZNmlC8eHFnt+JyoqL2Urp0glVjHjyAn36Cv/8dfHygcmWoXx+sfYOobFkT0dHR1g0SQyl8\nReSJLBaLNlplQlTUT5Qta90j50uXHl6oUKzYf3+vdGm4cMG62mXKJLFv30/WDRJDKXxF5Im2bdsG\nQHBwsFP7cFW//HKZIkWsG/PgAfj5Pf57fn6QlGTdPEWKwKVL560bJIZS+IrIE0VERDBs2DDc3Nyc\n3YpLSklJw9qrd319fx+0iYm/D+S/4ukJKSnWfdcsxlL4isjvxMfHs2XLFvr27evsVlyWl5cH1m44\nfu45MJngl1/++3uxsRnfbPWbtDTw8tLu9KxM4SsivzNr1iz69u1Lnjx5nN2Ky3rmmULcuGHdGF9f\naNgQ5s2D5GQ4dgx274YWLayb58YNKFKk2F//QXEavWokIo9JSUlh9uzZbN++3dmtuKSkpCSio6Px\n8spNTAxUr27d+FGj4LPPoFMnyJsXQkOx+mjKs2d9aN68kXWDxFAKXxF5zNdff03VqlUpV66cs1vJ\n8iwWC2fPnmXv3r2Pfp06dYpKlSqRN29eYmK8gRSr5syTBz78MHN9nT3rxbhxesc3K1P4ishjwsPD\neeedd5zdRpZ09+5doqKiHgtbPz8/goKCqFu3Ln379qV69er4+Phw4cIFatQoT0oKGHk42LVrcPly\nOtWqVTOuqFhNVwqKyCMHDhyga9euxMbGksuoA4mzKLPZzMmTJ9mzZ8+joH0YqDWoW7cudevWpU6d\nOhQr9sffrbZo0YDAwF20bGlc3/Pne5Anz98JD59lXFGxmla+ItlcSkoKR48eJTo6msOH93Pv3m3M\nZjN+fv5UqFCdwMBAatSoQd68eYmIiGDIkCE5Mnhv3LjBvn37HgXt/v37eeaZZ6hbty5BQUG8/vrr\nVK5cGU9PzwzPOXz4W7z9dm9atEjAiDe2UlNhwwZPduwY5fhikila+YpkU1FRUUyb9jmrV39L0aLe\nvPhiGs8//4DfNjAnJ0NcnBdnz/py9uwDatWqTlTUcWJjYylcuLBzm3ewtLQ0jh079iho9+zZw7Vr\n16hdu/Zjq9qCBQtmqo7JZKJGjQq0anWG1q0d/6N27lxP7t9vxjfffOfwWpI5Cl+RbOb48eMMHNiH\nS5fO0q5dMq1amcmX78/HpKTAtm0QGelJaurThIfPoX379ob0a4Rff/31scfHBw8epFSpUo+CNigo\niHLlyjlkxX/06FGCg+sya9YDChWy+/SPnDoF772Xh2PHTlPE2qO1xHAKX5FsIi0tjU8++ZiwsE/5\n+9+Tad3aYtM9socOwaRJfjRq1IZp02aRP39++zfrQMnJyRw6dOixTVEJCQmPQrZu3brUqlWLvHnz\nGtbTBx+8x5o1k/jkkySrT73KiLt3YeRIPz7+eDa9e/e2fwGxO4WvSDaQlJRESEgr7t2LJjQ0icw+\nNX7wAObM8ebQofxs2bKbUtYesWQQi8XCxYsXH3t8fPz4ccqVK/doVVu3bl3KlCnj1GMyTSYTvXp1\n5vLlH3n/ffsG8N27MG5cbkJCXuOTTz6338TiUApfEReXnJxMy5aN8fM7yhtvJNu02v0jq1e7s2pV\nfnbtiqZEiRL2m9hGCQkJHDhw4LFVrZub26MVbVBQEIGBgfhZexiyAVJTU+nbtxtnz/7IuHFJPPNM\n5uc8dw7+/W8/unYdzKefTtQ53C5E4Svi4nr27MSNGxt5+237Bu9vIiNzsWlTMQ4dOmloqJnNZs6c\nOfNY0J4+fZqqVas+tqotXry4y4SOyWTio4/GM3nyBF55JZk2bSw27YJOT4flyz1YtcqbTz8NY+DA\nV+3frDiUwlfEhUVGRvLmm/2YOTPJoQc5fPSRLxUrDmDy5AiH1bhz5w779+9/tDFq37595M2b97Gg\nrVatGt5GnljhIEePHuVvf+tOauovtG+fQNOmGTuI4/59+P57N9at86NMmWrMm7ckSzyREOspfEVc\n1PXr16lU6UXef/8uFSs6ttbdu/Dqq76sWvUDDRo0yPR8JpOJEydOPLaqjY+PJzAw8NHj4zp16mTr\nXbsmk4mNGzcybdqn7N8fRY0aHpQunUDZslCwIHh4PLyd6Ndf4fRpN86d8+fw4TRat27JsGFvUL9+\nfZdZ8cvvKXxFXNS4cW8SEzONkSOtOzvYVps2wfbt1dm586DVY69du/bYARZRUVEULVr0sVVtpUqV\n8PDImef+XLx4kZ9++omoqD1ERe3k+vUbpKam4ePjRdGiRalVqxG1atWhcePGPGOPL4vF6RS+Ii4o\nJSWF554rxMSJ9zHqqWN6OvTu7cvWrVFU/JOldmpqKkeOHHlsVXvz5k3q1Knz2AEWrvYKk4g95cyP\nmSIubuXKlbzwAoYFLzx8DNqmTRrh4ZOYMWPuo9+/dOnSY0F7+PBhXnjhBYKCgnjppZd49913CQgI\nwN1d14eL/EbhK+KC1q5dQaNG960ak5YGYWFw8ODDjTtFi8Krr0Lt2hmfo0mTdN5+exWlS5d/FLYp\nKSmPXvUZP348tWrVIs9vZ1iKyBMpfEVcUHT0AZo1s26MyQSFC8OUKQ//75498K9/wdy5kNF9TcWL\nw61bd4mJiaFTp0589tlnPP/889r4I2Ilha+Ii0lISODSpWtYe+iUjw/06/fffw4Kehi6Z85kPHzd\n3aFChTz06NGDl156yboGROQRfQkj4mLi4+MpXNiHzG4MvnULLl3C6hAvWtTEhQsXMldcJIdT+Iq4\nmOTkZLy9M/eYNz0dPvoIWrV6+CjZGp6eJlJSjHm9SSS7UviKuBhPT0/S020fbzbDxx+DlxeMGGH9\neJPJ3aoL5UXk9/Sdr4iLKViwIDdvpmCxYPW5wBYLTJgAd+7AJ59g01nQd+54UKBAAesHisgjWvmK\nuJgiRYrg5eXD1avWjw0Lg7i4h4+cbb3WLiYmnRo1atg2WEQArXxFXFL16lU4ffqnDO9SBrhyBdat\nexi6Xbr89/fHjCHDry3dugWpqW5Z9n5fEVeh8BVxQcHBrdm9O4pGjZIzPKZIEdiyJXN19+2D+vVr\n671ekUzSY2cRF9S//wC2b4eEBGPrrl+fhyFDxhhbVCQbUviKuKAiRYrQsmULNm40bgUaEwP37vnQ\nqlUrw2qKZFcKXxEX9dZb77N8uQ+3bzu+lskEM2bk5q233iOXLVukReQxCl8RFxUYGEj//oOZNs3P\n4bVWr3bH378cQ4YMdXgtkZxA9/mKuLDk5GSqVg2gQ4d42rZ1zH/KMTHwzjt+7N9/lNKlSzukhkhO\no5WviAvz8fFhzZofWLjwKbZutf/8587Be+/5MnfuYgWviB0pfEVcXEBAAD/8sINZs/KxcqU7ZrN9\n5j1wAMaO9WXq1Hl07NjRPpOKCKDHziLZRmxsLD17dsRiOc+YMYk8+6xt8yQmwhdf+BAdnZsvv1xO\nM2svDhaRv6SVr0g2Ubp0afbuPUz37v9g6FBfIiK8iIvL+Pjbt2HJEndeecWPAgU6c+JErIJXxEG0\n8hXJhuLj45kxI5w5c2ZSvLiFSpUSKVvWTJky4O8P7u7w4AGcPw+nT0NMjD+HDqXTpUtnhg0bo7Ob\nRRxM4SuSjaWmprJx40b27t3F/v07OH78JAkJyZjNZnx9vSlb9nlq1WpArVpBtGvXjqefftrZLYvk\nCApfERERg+k7XxEREYMpfEVERAym8BURETGYwldERMRgCl8RERGDKXxFREQMpvAVERExmMJXRETE\nYApfERERgyl8RUREDKbwFRERMZjCV0RExGAKXxEREYMpfEVERAym8BURETGYwldERMRgCl8RERGD\nKXxFREQMpvAVERExmMJXRETEYApfERERgyl8RUREDKbwFRERMZjCV0RExGAKXxEREYMpfEVERAym\n8BURETGYwldERMRgCl8RERGDKXxFREQMpvAVERExmMJXRETEYApfERERgyl8RUREDKbwFRERMZjC\nV0RExGAKXxEREYMpfEVERAym8BURETGYwldERMRgCl8RERGDKXxFREQMpvAVERExmMJXRETEYApf\nERERg/0/EyJBYWj8WPoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4319f90>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One measure of the importance of each node is its <b>centrality</b> [4], [5], [6].\n",
      "It is well known that the number of walks of length $k$ from node $i$ to node $j$ is given by $A^k_{ij}$.\n",
      "We can add up all these contributions to obtain the centrality of node $i$ as\n",
      "$$c(i) = \\alpha_1 A_{ii} + \\alpha_2 A^2_{ii} + \\alpha_3 A^3_{ii} + \\cdots,$$\n",
      "where the coefficients $\\alpha_k$ remain to be chosen.\n",
      "Typically we assume that shorter walks are more important than longer ones so $\\alpha_k \\ge \\alpha_{k+1}$.\n",
      "There are a number of proposed formulae for $\\alpha_k$ but taking $\\alpha_k = \\frac{1}{k!}$ we obtain $c(i) = \\exp(A)$.\n",
      "\n",
      "Looking at our network above we might expect node 1 to be the most important:\n",
      "it has the highest degree and links nodes 3, 4, and 5 to nodes 0 and 2.\n",
      "Let's work out the centrality of each node and rank them in order of importance."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "centralities = np.diag(la.expm(np.array(Adj, dtype=np.double)))\n",
      "nodeorder = np.argsort(centralities)[::-1]\n",
      "\n",
      "print np.array([nodeorder, centralities[nodeorder]])\n",
      "\n",
      "# Note: This is already built into networkx using the following command\n",
      "# print nx.communicability_centrality_exp(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.          0.          2.          3.          4.          5.        ]\n",
        " [ 4.44723536  2.86427609  2.86427609  2.36018456  1.71615913  1.59432922]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we expected node 1 was the most important, followed by nodes 0 and 2 which have identical centralities since switching the two nodes wouldn't change the appearance of the graph."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Computing the Exponential"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section we will give a basic algorithm to compute the matrix exponential.\n",
      "We will use the Taylor series (from the first section), \n",
      "which converges for every $A$ but the convergence can be <b>very</b> slow when $\\|A\\|$ is large.\n",
      "To speed it up we can use the relationship \n",
      "$$\\exp\\left(\\frac{A}{p}\\right)^{p} = \\exp(A)$$\n",
      "repeatedly with $p = 2$ to make $\\|A/2^s\\|$ small before applying the Taylor series and squaring the result $s$ times.\n",
      "This is the basis of the <b>scaling and squaring</b> algorithm.\n",
      "\n",
      "A number of potential methods for computing the exponential are described in detail within [8], [9]."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sconst = 1\n",
      "\n",
      "def TaylorSS(A):\n",
      "    taylordegree = 15 # Use order 15 Taylor approximation\n",
      "    s = np.ceil(sp.log2(la.norm(A))) + sconst # Find s such that norm(A/2^s) is small.\n",
      "    X = A/(2**s)\n",
      "    eX = np.eye(np.shape(A)[0])\n",
      "    for k in range(taylordegree): # Compute the Taylor series\n",
      "        eX = eX + X/sp.misc.factorial(k+1) \n",
      "        X = sp.dot(X, X)\n",
      "        \n",
      "    for k in range(np.int64(s)):\n",
      "        eX = sp.dot(eX, eX) # Do the squaring phase of the algorithm\n",
      "        \n",
      "    return eX\n",
      "\n",
      "#Let's test it against la.expm\n",
      "A = np.random.randn(4, 4)\n",
      "E1 = TaylorSS(A)\n",
      "Eexact = la.expm(A)\n",
      "print la.norm(E1 - Eexact)/la.norm(Eexact) # Relative error"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0116677518269\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Our algorithm worked but the relative error compared to <b>scipy.linalg.expm</b> is quite poor.\n",
      "If we increase sconst then $\\|X\\|$ will be much smaller and the Taylor series will be more accurate."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "relerrs = np.zeros(15)\n",
      "for sconst in range(15):\n",
      "    E = TaylorSS(A)\n",
      "    relerrs[sconst] = la.norm(E - Eexact)/la.norm(Eexact)\n",
      "    \n",
      "plt.plot(relerrs, color='blue', lw=2)\n",
      "plt.yscale('log')\n",
      "plt.xlabel('sconst', fontsize=12)\n",
      "plt.ylabel('rel. err.', fontsize=12)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<matplotlib.text.Text at 0x4441fd0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEVCAYAAADKN2OaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuczGX/x/HX2nUId4misJkiIZJOWm63LbodKp3P3FLp\nSOeiqBlJ0YkODhVy6udUIdWddNcUKaSEnBWhSESo5DC/Pz4z7cFgZndnv9/vzPv5eHiYa3Zn9pPY\nz17X57quD4iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIgkseOBocBEpwMRERFvUyIREXFQCacD\nEBERb3NDIhkObAQW5nu+NbAUWAF0Cz/XAegPVC226ERExPWaAY3Im0jSgZWADygJzAfq5ntdRWAI\neRONiIikKB95E0kW8H6ucffwLxERcZkMpwM4gGrA2lzjdUDjeN+kZs2aoVWrVhVZUCIiKWIVUCvW\nT3ZDjSSaUFG8yapVqwiFQp795ff7HY8hVeP3cuyK3/lfXo8fqBnP91q3JpL1QGaucSY2K4lbIBAg\nGAwWRUwiIkktGAwSCATifp1bE8mXwIlY7aQUcBXwdkHeKBAIkJ2dXWSBiYgkq+zsbM8mkrHALKA2\nVhfpBOwBugDTgMXAeGCJUwE6xesJ0Mvxezl2UPxO83r88UpzOoAEC/n9frKzs1Puf6yISLyCwSDB\nYJBevXpBHPkh6RNJuHAkIiIxSktLgzjygxuWtkRExMOUSEREpFDSnQ4gwQKRBz6fz7koREQ8IBgM\nMmLECD755BOAXrG+TjUSERHJQzWSfC64AGbMAOUTEZHESPoZSeS2lawsePBBaNcOSiR9+hQRKTjN\nSPJ59FGoWBE+/xwuuQTq1YNhw2DXLqcjExFJDkk/IwmFQuzcacnj2Wfhhx/sA8ceC3ffDbfcAkcc\n4WyQIiJuEu+MJCUSScTu3TBhAvTrBwvD3U8OPxxuvdWSyrHHOhSliIiLKJHkFXXXVigE06ZZQolc\nDFyqFHToAA88ACedVLxBioi4SbyJJCXPkaSlQa1acP310LYt/PorfPstfPUVDBwI8+eDzwfVqxd7\nvCIijtE5kuhiPkeyfLnVUEaMgL/+sueaNYNu3SzZpCX7n5SISFgqLW1dBJwPHA4MA6ZH+Zy4DyRu\n2AAvvACDBsG2bfZc/fq25HXNNVCyZOGCFhFxu1RKJBEVgGeAm6J8rMAn23/7DV55Bfr3hx9/tOcy\nM+Gee6BzZyhfvqDhioi4WyomkmeAMcD8KB8r9BUpf/0Fr78OTz0FS5fac0ceCbffDnfdBUcfXai3\nFxFxHS8eSBwObAQW5nu+NbAUWAF0Cz/XAegPVMX+I/sB/yV6EikSpUpBp05WjJ8yBZo0seJ8nz5w\nwgnw+OOwc2eivrqIiPu5YUbSDNgBjAIahJ9LB5YBLYH1wFzgGvK2270T+E/4Y/OBl6O8d0Iubfzs\nM0sg779v46pV4bHHbBdYerLvgxORpOfVpS0fMJWcRJIF+LFZCUD38O9943zfkN/v/3tQ1C13P/rI\nivBffWXjk0+2syna5SUiXhJpsRvh1Va7PvImksuBVkDn8Lg90BjoGuf7Jvwa+X37YPx4ePhhWL3a\nnjvnHKupnHFGQr+0iEhCeLFGEo1nLn0vUcK2BS9daudQjjwSPv4YzjwTrr0Wvv/e6QhFRBLLrYlk\nPZCZa5wJrCvIGwUCgTxTtkQpXRruvRdWrbLlrtKlYexYqFMH7rsPtmxJeAgiIoUSDAYJBAJxv86t\nS1sZWLG9BfAjMIf9i+2xcKxD4po10LMnjBlj4woVoEcP6NIFypRxJCQRkZh4cWlrLDALqA2sBToB\ne4AuwDRgMTCe+JOIo2rUgNGjYd48aNECtm7NuRByzBirrYiIJAO3zEgSJeT3+4t8t1bcQYTggw8s\nkUSur2/UyAryLVs6FpaISB6R3Vte3bWVKI4tbUWzd6/NUnr2hPXr7bnWrW3L8CmnOBubiEiEV8+R\nJIqrEknE77/D88/Dk0/C9u125qRjR+jdW1fXi4jzlEjycmUiidi0yU7IDxoEe/ZYEf6ee+zqerX/\nFRGnqLFVXoHIg9yNrdyiXDlo08bOm2zYAAsWwMyZMHQoHHaY1VF05YqIFBc1torO1TOS/GbPhvvv\nt2QCULu2LYG1bn3w14mIFCUvbv+VsMaN4dNP7Zbh2rWta2ObNnDRRXbQUUTEjZRIXCYtDdq1s23C\nTz9tDbTefhvq1bMDjbqyXkTcRktbLvfTT9C9O4waZeNq1eCZZ+Cqq3TDsIgkhnZt5eX5RBLx+efQ\ntaudlAdo1gxefBEaNnQ2LhFJPqqRJKmsLJgzx3Z0HX00zJgBp50Gd9wBmzc7HZ2IpDIlEg8pUQJu\nvNGK8HfdZUtbgwZZYX7wYDs5LyJS3JL9lEIg8sCN50gKqkwZ2xJ86aXWB2XJEnj3XdvtdfLJdmGk\niEi8dI4kuqSpkRxIKARvvmk9T374wZ675hrb8VWtmrOxiYg3pVKNpA4wGJgA3OhwLI5JS4PLL7dZ\nid9vs5WxY+26+iefhF27nI5QRJJdMsxISgDjgCujfCzpZyT5rV5ts5O33rJxzZowYACcf762C4tI\nbFJpRgJwIfAulkgE8PlsqWv6dKhb107EX3ihJZLly52OTkSSkRsSyXBgI7Aw3/OtgaXACqBb+LkO\nQH+gang8FWgDdEx8mN7SsiV88w307w+HHw7//S/Ur283C2/f7nR0IpJM3LDY0QzYAYwip2d7Otaz\nvSWwHpjL/j3bmwOXAmXCzw+I8t4pt7QVzc8/w0MPwfDhNj7mGOvOeN11tqVYRCQ3r55s92Gzi0gi\nyQL82KwEoHv4975xvm/I7/f/PXC65a7T5syx0/Fz5ti4aVM7h6LujCKpLdJiN8KrrXZ95E0klwOt\ngM7hcXugMdA1zvfVjCSfffvs3q5u3Wymkp5uyaVXL1sCExFJlmK7vvsnSIkScP31sGwZdOli51AG\nDIA6dWDcOBuLiMTDrYlkPZCZa5wJrCvIGwUCgTxTNjEVKtilj3PnWh+Un36yg4wtW9ppeRFJPcFg\nkEAgEPfr3Lq0lYEV21sAPwJz2L/YHgstbcVg3z4YNsyuq9+yBUqWtLMoPXtaO2ARSS1eXNoaC8wC\nagNrgU7AHqALMA1YDIwn/iQiMSpRAjp3tuWum26C3buhb19rpjV5spa7ROTg3DIjSZSQ3+9P+d1a\n8friC7j9dvj6axu3bWvLYCec4GxcIpJYkd1bXt21lSha2iqgvXvtavoePeC336B0aTuL0q2b3ecl\nIsnLq+dIEkWJpJA2bIAHHoAxY2xcs6bNTtq0cTYuEUkcJZK8lEiKyCef2HLX4sU2vvRS2zacmXnw\n14mI98SbSNTYSmLi88HNN9uhxVmzYMECePllyMiAs86yg40i4m1qbBWdZiQJsG4d3HMPvPGGjevW\nhYED4ZxznI1LRIqGF7f/isdUrw4TJ8K0aXDiidZU69xz7RLIn35yOjoRKW5KJFJg//43LFwIvXvb\nTq7/+z+7auX552HPHqejE5HioqUtKRLffw933QVTp9q4YUPbPpyV5WxcIhI/LW2JI44/Ht5+G6ZM\ngRo1rKlWkyZWoN+yxenoRCSRlEikSLVrZ1uEH37Y7ux69VU46SQYOVJXrYgkKy1tScIsWQK33WZn\nUAD+9S9b7qpXz9m4ROTgdI4kr0Dkgc6RFL+jj4aOHe00/MyZllheeQX++MOWvUqWdDpCEcktVc+R\nlAOCWMJ4N8rHNSNxiS1bbLnr5Zdt7PPZVSsXXOBoWCISRaoV2x/ErpgXl6tYEYYMgc8/tx1dq1fD\nhRfaVStr1zodnYgUhpcTyXlYr5JNTgcisTv7bPjyS3juOShfHiZNspPxzz5rfVBExHvckEiGAxuB\nhfmebw0sBVYA3cLPdQD6A1WB5sDZwLVAZ7y/TJcyMjLsipUlS+Cyy2DnTrj/fjj9dLvHS0S8xQ3f\nfJsBO4BR5LTaTcda7bbE+rfP5cCtdjtis5L3onxMNRIPeO896NLFDjWCdWns2xcqVXI2LpFU5cUa\nyQzg13zPnQWsBFYDu4FxwEUHeP1IoicR8Yi2bWHRImuiVbIkDB1qV63o7ImIN2Q4HcABVMP6t0es\nAxoX5I0CgcDfj9Vy173KloXHH7eLH2+/HYJBuP56GD5cZ09EEi3SYreg3LC0BeADppKztHUZViPp\nHB63xxJJ1zjfV0tbHhQKWUfG++6DTZuspnL//fDII5ZwRCSxvLi0Fc16IHfvvUxsVhK3QCBQqEwr\nxS8tDTp0gGXL4JZbrH983742K3nnHaejE0lewWAwzypOrNw6I8nAiu0tgB+BORy42H4wmpEkgS++\ngFtvtYsgAS6+2K6qP+44Z+MSSVZenJGMBWYBtbG6SCdgD9AFmIadFRlP/ElEkkTk7En//nb2ZPJk\nm53o7ImIO7hlRpIompEkmXXr4O674c03bXzKKXbtytlnOxuXSDLRpY15BSIPdGljcjj8cLjySjjr\nLDu8uGyZ7ezasAGaNoXDDnM6QhHvStVLGw9FM5Ik9vvv1ub3mWestW/lyrb8dc01VrAXkYKJd0aS\n7P/clEhSwKJFVoz/7DMbt2wJgwbBiSc6G5eIV3mx2C5SKPXrw6ef2on4ihXhww+hQQN47DHYtcvp\n6ESSn2okkhTS0uC00+CGG+Dnn2HePDsdP2GCJZrjj3c6QhH3U40kOi1tpahg0Ja7li2zcYcOVkup\nXNnRsEQ8QUtbIkB2th1g7N0bSpeG0aPtIshXX4V9+5yOTiS5aEYiSW/lSrjjDvjgAxs3bWrdGuvX\ndzYuEbfSjEQkn1q14P33YexYqFLFdnc1agTdullTLREpHM1IJKVs3Wp9TwYPtluGa9SAgQPh/POd\njkzEPXSOJC8lEolq9mwrxs+fb+NLL7WLIKtXdzYuETfQFSl5BSIPtP1Xcqte3Vr6HnkkzJwJCxZY\nIb5cOTjjDCihRV9JQU5t/x0M3FbI90gkzUjkkNauhTvvtFuFwc6jDBkCZ57pbFwiTklEsb0EcC5Q\nuoAxJUo21u99MNDc2VDEyzIzYdIkmDLFepx89RU0bgxdu8K2bU5HJ+J+sSSSfcDbQLTLJpycjewD\ntmMJrkDdE0Vya9cOvv3W2vqWKAEvvQR168Ibb1hhXkSii3Xq8h7QG/g8gbHEKw0IAZWB57C+7vlp\naUsK5JtvrM3v7Nk2btvWdnep1CapIFHnSNYA/wVGYAkl8uux+MKLajiwEViY7/nWwFJgBdAt/FwH\noD9QFUsiAFtx37KbeFzDhtbvZPBgOOIIeO8968r41FPqyiiSX6wZ5zVyZgC5XxvCWuMWRjNgBzCK\nnJ7t6VjP9pbAemAu+/dsvwRoBVQABgGfRnlvzUik0DZsgHvugXHjbNyggXVlzMpyNi6RREnEOZJ0\nbBttH+DPAkV1aD5gKjmJJAvwY7MSgO7h3/vG+b5KJFJkpk2D22+H776z8c03Q9++toVYJJnEm0gy\nYvicvVhR3V/AmAqiGrA213gd0LggbxQIBP5+nJ2dTXZ2dmHikhTWqpU10erTx5a4XnnFtgyrK6N4\nXTAYJBgMFvj1sf7Vfw5YBQws8Fc6OB95ZySXYbORzuFxeyyRdI3zfTUjkYRYvNhOxs+YYePzzrOu\njLVqORuXSFFIVLG9MTAAK7rPxM5vzCB6XaIorAcyc40zKeAW30AgUKhMKxJNvXrW82TYMOvKOH26\n3Sb8+OPqyijeFQwG86zixCrWjHP9AZ4PASPj/qr785F3RpKBFdtbAD8Cc9i/2B4LzUgk4TZtsrMn\no0bZuE4dOxnfXMdkxaO8eGnjWOxkeiXgZ+BRbJdYG2wWlA4MA54swHsrkUix+fhjW+5avtzGnTpZ\nLeWoo5yNSyReiUokJYCbgKuBo7GZw7+AY4AJ8YVYrEJ+v19Fdik2u3bZTq4nnoC//oJKlazFb8eO\nKsaL+0WK7r169YIEJJLewL+xGcIQ4AigJjAROC2+UIuVZiTiiOXL4bbb4KOPbNy8uS131anjbFwi\nsUjUjGQd0AjYBPwKHInNUrZgBwLdSolEHBMKwZgxcO+98MsvULIkdO8ODz8MZco4HZ3IgSVq11YJ\n7PR5buWwSxNFJIq0NOjQAZYts94nu3dD7952Mv7DD52OTqToxJpI/oudJYn8HFUCW+6amoigRJJJ\nxYrWNGvGDNs2vHKlnTtp3x5+/tnp6EQKL9YOiZ9ghwKHAWWxSxQ3A3cQ/Xp5twhEHqhDojjtuOOg\nc2coWxY++wy+/trOoVSqBKeeqmK8OK+4OiRWAWpg15f8FOdrnaAaibjSd9/ZvV3Tptm4aVO7CPLk\nk52NSwS8eY4kkZRIxLVCIRg/Hu6+GzZuhIwMeOAB6NnTZi0iTlEiyUuJRFxv61Z46CGbkYRCcMIJ\ndm9Xq1ZORyapSokkLyUS8YzPP7eujAvDLd6uvtpuFj7mGGfjktSTqO2/IpJgWVkwb55dq3LYYdZI\nK3Jv1759TkcncmCakYi40OrVcMcd1uIX4OyzbenrlFMcDUtSRLwzkli3/x7Idgp2mWJxCUQeaPuv\neEmFCnDttXY1/cyZsHSpnUXZuROaNIFSpZyOUJJRcW3/za8Z1pfErTQjEc/bts12cg0caMX4GjWs\nGN+2rdORSbJKpWJ7GvA48A/gS2BUlM9RIpGkMWeOFePnz7fx5ZfD889D1arOxiXJpygTyY1Y46qD\nvTYEDI/1ixWxS4CLgF+A94CPonyOEokklT174IUX4NFHbZnrH/+AJ5+0PijphV2oFgkrykQS5OCJ\nJOKcWL9YEeuG3T78Knad/RVRPkeJRJLSDz9A167w9ts2PvNMeOUVu2pFpLC8uLQ1HDgf647YINfz\nrcnpkDgU6Ad0wPqfPI0lsL+wJDIeuCrKeyuRSFKbPBm6dIH1621Gctdd0KsXlC/vdGTiZYlMJJWw\nb/jHAE8B1cKvXxfHe0TTDLuifhQ5iSQd69neElgPzGX/nu2HAS8Cv4efHxzlvZVIJOlt3w6PPAIv\nvmjnTTIz4aWXoF07pyMTr0pUImkOvIkVtZtiBe5s4D7gwrgijM6HXUkfSSRZgB+blQB0D//eN873\nVSKRlDFvnhXj582z8SWXWD2lenVn4xLviTeRZMT4ec9j/do/xDokAnwBNI4nuDhUw24YjlhX0K8V\nCAT+fqze7ZLMTj8dZs+2bcI9esCkSTB9OvTpY4cbVYyXA4n0ai+oWDNOpL1u7sfpWF2jUoG/eg4f\neWckl2Gzkc7hcXsskXSN8301I5GUtG4d3HmnJROwJPPyy/a7yKEk6q6tJeQsM0W0ABbG+oXitB7I\nzDXOpIC1mEAgUKhMK+JF1avDW2/BlClWM5k3D846y66s364G2XIAwWAwzypOrGLNOGcD72DnNa4A\nRmO1kYuAOXF/1f35yDsjycCK7S2AH8NfI3+xPRaakUjK27ED/H47vLh3ryWZF1+Eiy92OjJxq0TN\nSOYADYFvgdeA74AzKZokMhaYBdTG6iKdgD1AF2AasBjb3htvEhERbCvws8/C3Ll23mTdOivEX3wx\nrF176NeLHEosGScDu5yxAu7uzx5NyO/3q8guErZ3LwweDA8/bEtc5cpB7952uDEj1q03krQiRfde\nvXpBArb/LgDaYLULL9HSlkgU69dbveSNN2zcqJGdjD/jDGfjEndI1DmSB7Htvy9gy0+5vztHu+PK\nLZRIRA7i3Xdta/CaNVCihD1+/HE4/HCnIxMnJSqRrA7/Hu278vGxfjEHKJGIHMLOnXatynPP2dJX\n1apWjL/kEkhzwyVKUuwSVWz3hX8dH+WXq2n7r8jBlStn7X3nzYPGjeHHH+Gyy+yKlTVrnI5OilOi\nt/96lWYkInHYu9dqJd27w2+/Qdmy8NhjdhmkivGpw4u3/yaSEolIAfz0kxXjJ0ywccOGdjK+caIu\nRRJXSdTSloikkGOPhfHj4b33wOeDb76BrCy7sn7bNqejE7dRIhGRA2rTBr79Frp1s0sfBw6EunVh\n4kTrHy8CWtoSkRgtXGjX1H/+uY3btrXE4vM5GpYkgJa2RCQhGjSAmTNhyBA44ghb9qpXz3Z87d7t\ndHTiJM1IRCRuGzbAvffC2LE2btDAivFZWc7GJUUj3hlJsre6CUQe+DT/Fiky5cvbWZOsLPjsM1i2\nDIYPh40boWlTKFPG6QilIILBICNGjOCTTz4B6BXr6zQjEZFC+f13u1bl6adhzx445hgYMACuvFIn\n470qlc6R/BO4DruduB7WSz4/JRKRYrJoEdx6q81QAFq1gkGD4IQTnI1L4pdKiSTiIqAy8GqUjymR\niBSjfftg2DB48EHYutWWuPx+uO8+KFnS6egkVqmYSMYDNwA7o3xMiUTEARs3WvJ4/XUb169vxfgm\nTZyNS2Ljxe2/w4GN7N//vTWwFFgBdAs/1wHoD1QNj48DthE9iYiIQ6pUgTFj4IMPoGZNW/Zq2tTO\nofz6q9PRSVFzw4ykGbADGEVOz/Z0rGd7S6yZ1lyi92wPAO8DXxzgvTUjEXHYH39Anz45500qV7Zi\n/NVXqxjvVl5d2vIBU8lJJFmAH5uVAHQP/943zvcN+f3+vwdquSvinG+/tWL8zJk2Pu88a/tbs6az\ncUlOi92IRLXaTTQfeRPJ5UAroHN43B5oDHSN8301IxFxkX374LXX4IEHbImrTBl45BG4/34oVcrp\n6CTCizWSaPTdXyQJlSgBN94IS5dC+/bw55/Qo4f1jI/MVMR73JpI1gOZucaZwLqCvJE6JIq4T+XK\nMHo0TJ8OtWrB4sXQrBl07gxbtjgdXeryeodEH3mXtjKwYnsL4EdgDtGL7YeipS0Rl/vzT3jiCejb\n14rxRx9t/eOvu07FeKd4cWlrLDALqA2sBToBe4AuwDRgMXZWJN4kIiIeUKaMtfP95hublWzaBB06\nWDF+xQqno5NYJHu+D/n9fu3WEvGIfftgxAgrxm/ZAqVLQ8+eNi5d2unokl9k95ZXd20lipa2RDxo\n0ybbyTVqlI3r1rU+KP/6l7NxpQovLm2JiORx9NEwciT8739w4omwZAk0b247vjZvdjo6yU+JRERc\n69xzYcECu/ixVCnreVKnju340mKDe6ixlYi4WkYGZGfDFVfYnV1LlsCkSTBjhjXWqlTJ6QiThxpb\nRacaiUgSCYVsyev++22Jq3RpO9D44IMqxhclr961lShKJCJJ6JdfbCfXiBE2Pukku6a+eXNHw0oa\nKraLSNI76ii7s+ujj6B2besZn50NN9xgSUaKlxKJiHjWOedYMT4QsGL8a69ZMX7kSBXji5OWtkQk\nKSxbBrfdBh9/bOPsbDt7ctJJjoblSVraEpGUdNJJdu5k5EjbyRUMwimnQK9esGuX09ElN81IRCTp\n/PKL7eR67TUb165txXjdlBSbeGckOkciIkmnbFm46CKroXzxhS17jRwJq1fDP/9pH5f96RxJdJqR\niKS4XbusX3yfPva4UiV45hno2FHX1B9IKp0jqQ68APwKLAf6RfkcJRIRAWD5civGf/SRjZs3t2J8\nnTrOxuVGqVRsbwC8CdwINHI4FhFxudq14cMP7Ubho46CTz6xYrzfb821pOC8PCM5Angba4I1GhgR\n5XM0IxGR/WzeDN26wbBhNj7xRJudnHuus3G5hRdnJMOBjcDCfM+3BpYCK4Bu4ec6AP2BqlgnxZ5Y\nO97ziyVSEUkKlSrB0KHw6afW62TFCmjRAv7zH+uFIvFxw4ykGbADGEVOz/Z0rGd7S2A9MJf9e7af\nAjwKbAK2Aw9GeW/NSETkoP76C55+Gnr3tmJ8xYpWnO/UCUq44UdtB3i12O4DppKTSLIAPzYrAege\n/r1vnO8b8vv9fw/UcldEDmTlSrj9dpg+3cbNmtlyV716zsZVHCItdiO82mrXR95EcjnQCugcHrcH\nGgNd43xfzUhEJGahEIwbB3ffDT//DCVL2sHGHj3gsMOcjq74eLFGEo2++4tIsUtLg2uugaVL4eab\nYfduO3/SoEHOTEX259ZEsh7IzDXOBNYV5I0CgUCeKZuIyKEceaRdqTJzJtSvD6tWwb//DdddBxs3\nOh1d4gSDQQKBQNyvc+vSVgZWbG8B/AjMYf9ieyy0tCUihbJ7Nzz3nF3++McfUKEC9OsHN92UvMV4\nLy5tjQVmAbWBtdi23j1AF2AasBgYT/xJRESk0EqWtDMnixZB69awdSvccosV4xctcjo6d3DLjCRR\nQn6/X7u1RKRIhEIwcSLcdRds2AAZGdY//pFHkuMiyMjuLa/u2koULW2JSJHbutV2cg0ebMnl+ONh\n0CCbsSQDr54jSRQlEhFJmC++sGWuBQtsfNVV0L8/HHuss3EVlhdrJCIinnT22fDll3YyvmxZGD/e\nrlwZMgT27XM6uuKjxlYiIoWQng5NmtjW4JUrYeFCePdd+OADaNwYqlRxOsLYqbFVdFraEpFiEwrB\nm2/CnXfCTz9ZkrnvPnj0UShXzunoYqcaSV5KJCJS7LZtg549YeBASy4+nz1u29bpyGKjRJKXEomI\nOGbOHCvGz59v4yuugAEDoGpVZ+M6FBXbRURc4qyzYO5cePZZK8ZPnGjF+EGDYO9ep6MrOpqRiIgU\ngzVroGtXmDrVxo0b231eDRs6G1c0mpGIiLhQjRowZQq89RZUqwazZ8Ppp8MDD8DOnU5HVziakYiI\nFLPffrNrVV580YrxNWrASy/BBRc4HZmJd0aicyQiIsWsdGlo08Z2cX35pfU/GTvWLoFs2hQOP9yZ\nuFLxHEk9rB3vZuB/wJtRPkczEhFxtT17bGbyyCO2xPWPf8ATT8Btt9k5FCek0vbfe7E+JTOBKcBF\nUT5HiUREPGHtWivGT5li4zPPhFdegVNPLf5YUqnYPhq4GngKqORwLCIihZKZCZMnw6RJUL26bRs+\n4wy7pn7HDqejOzg3JJLhwEZgYb7nWwNLgRVAt/BzHYD+QFVgE9b86iHgl2KJVEQkwS6+GBYvhrvv\ntkL8s89CvXo524bdyA1LW82AHcAoclrtpmOtdlti/dvnsn+r3RrAw0A5YBDWZTE/LW2JiGfNmwc3\n3wxffWWjOlMaAAAItklEQVTjSy+F55+3GUsieXFpawbwa77nzgJWAquB3cA49q+BrAFuAdoTPYmI\niHja6afbeZMBA6B8eTuDUreujXfvdjq6HBlOB3AA1bD+7RHrgMYFeaNAIPD3Y7XcFRGvyciw1r6X\nXWa3Ck+aBPfcA0OH2m6vc84p/NeItNgtKDcsbQH4gKnkLG1dhtVIOofH7bFE0jXO99XSlogklXfe\nsfrJqlU2vvJKeOYZK9YXFS8ubUWzHsj9x5KJzUriFggECpVpRUTc5IIL7ODi44/DYYfBhAlQp46d\nPdm1q3DvHQwG86zixMqtM5IMrNjeAvgROy+Sv9geC81IRCRp/fCDbQ+eONHGtWpZMb6wfU+8OCMZ\nixXLa2N1kU7AHmxr7zRgMTCe+JOIiEhSO+44m5F8+KFtEV65Es4/Hy68MGfpqzi4ZUaSKJqRiEhK\n2L3bLn70+2H7dihVym4Wfuih+Nv86tLGvAKRB7q0UUSSWXo6ZGVBp07wyy929mTGDBg92s6d1KsH\naYdIDal4aWMsNCMRkZQ0axZ06QJff23jc8+17cL16h36tV6skYiISBFr0sTu6xoyBCpWhI8+sm6M\n994L27YV7dfSjEREJMlt3mzX1A8ZYvd3VakC/fpBhw5QIsp0QjWSvAKRB6qRiEiqKls2ZzfXokXW\nSGvyZJg+HRo1gmOPtc9TjSQ6zUhERHIJhWDMGNvRtXGjFeBvvhn69IFK4YYcqpGIiMgBpaXZktby\n5XDffbbb6+WXoXZtGDwY9u4twHsWfZiuohmJiMhBLF5sl0H+7382PvVUmD8/dVrtxkKJRETkEEIh\nu6L+3nvt2pVwatDSloiIxCYtza6pX7LEdnfF/fqiD8lVNCMREYmTtv/mFYg80PZfEZGD0/bf6DQj\nERGJU7Ju/z0eGAqEb92nHDASeAW41qmgRETEO4nke+CmXONLgQnAzUA7RyISERHAO4kkv2pYEyyA\nAhyf8Qavtwj2cvxejh0Uv9O8Hn+8ijuRDAc2AgvzPd8aWAqsALqFn+sA9AeqRnmfdeT0dPdqMjwk\nr/9l9HL8Xo4dFL/TvB5/vIr7m/BrWNLILR14Kfx8Paw3e11gNHAP1rO9IjAEaIQlmreAy4BBwNvF\nEbiIiESXUcxfbwbgy/fcWcBKYHV4PA64iLw92rcAt+Z73Q1FH56IiMTLie2/PmAq0CA8vhxoBXQO\nj9sDjYGuRfC1VgI1i+B9RERSySqgVqyfXNwzkmgSedAj5j8IEREpGDcUqteTUzgn/HidQ7GIiIgH\n+Mi7aysDm0b5gFLAfKzYLiIisp+x2C6sXdg5kE7h59sAy7CaxkNF8HWibSf2ikzgY+BbYBFwp7Ph\nFFg68DVWD/OaCsAb2IaPxcDZzoYTt4ewvz8Lgf8DSjsbziFFOxZQEZgOLAc+wP6fuFW0+J/G/v58\ng+0yPcKBuGJxoCMZAPcB+7D/FyknHUtIPqAk3pvhHAOcGn5cHkuwXoo/4l7gdby5PXskObsCM3Dv\nN4FofMB35CSP8UBHx6KJTTNsa3/ub2ZPAQ+GH3cD+hZ3UHGIFv955JQO+uLe+KPFDvYD7fvYrSIp\nmUiysD+AiO7hX141GWjhdBBxqg58CJyD92YkR2DfiL2qIvbDx5FYEpwKtHQ0otj4yPvNbClQJfz4\nmPDYzXxE/6ke4BJgTPGFEjcf+8c+ETiFGBOJG4rtRS339SlghftqDsVSWD7sp4XZDscRr/7AA9i0\n2GuOBzZhh2e/Al4FyjoaUXy2AM8CP2DLyFuxpO41VbAlF8K/VznI57rdDcB7TgcRh4uw75sLYn1B\nMiaSZLk3vjy2Tn8XsMPhWOJxAfAzVh/xYpuCDOA07NaE04CdeGtGWxO4G/shpCr29+g6JwMqAiG8\n+++6B/AXVqvygrLAw4A/13OH/HecjIkkGbYTlwTexKbDkx2OJV5NsBuZv8c2V5wLjHI0ovisC/+a\nGx6/gSUUrzgDmAVsBvZghd4mjkZUMBuxJS2AY7EfTrzmeqAt3krkNbEfQr7B/g1XB+YBlR2MyRFe\n306chn3j7e90IEWgOd6rkQB8CtQOPw4A/ZwLJW4Nsd1+h2F/l0YCdzgaUWx87F9sj+y47I57i9UR\nPvLG3xrbOXeUI9HEx8eB6zspW2yHot9OXJz+idUW5mPLQ1+z/0WXXtEcb+7aaojNSNy+dfNAHiRn\n++9IbIbrZpFjAX+RcyygIlbb8cL23/zx34AdPVhDzr/hQY5Fd3AHOpIR8R0pnEhERERERERERERE\nREREREREREREREREREQkr3SnAxCRuAXJ6fci4rhkvGtLJNl59QJDERHB7n9aB/yG9cg4F/uB7GHs\nSp7fgC+xy+7ALkyci13nPgfrlxMRBB4DZoZfNw2oFP5YGezSzl+AX8OvrQz0wS5j/APYDrxQ5P+F\nIiKSMCdhfT4it9IeB5yA9V5ZAJwYfr4Bdj9RRSwJXIclm6uxfiFHhj8viN3JVAtLHB8DT4Y/dgt2\nT1kZ7PLFRsA/wh/7mJwOjiKO09KWSOz2Yi1sT8YuQvwBu9TuRqzvxIrw5y3EEsb52OWhr2MXcY7D\nZjHtwp8XwhporQT+BCaQ02b5L2x2cmL4877GZiARXuz1IklKiUQkdiuxplEBrF/GWKx5VCbWuiC/\nqliyyW1N+PmIDbke/4E1ogIYjS11jcN67PTDWiREqE4irqFEIhKfsUAzoAb2zbwfdv12rSifuz78\nebnVCD9/KHuw+snJWJ3lAuA/4Y8piYirKJGIxK42VlwvjfVv+BP7hj8U6I0lkzTgFKw+8l74Nddg\ns4mrgDrAO7ne80BLVOdgtZZ0bElrN7a0BjYbqllE/00iIlKMGgCzsR1Wm7Fi+DHYD2Q9sHrJb+HP\niSxfNcV2cW3Fdm/lbnubv2jeEevOCFaYXwrswJa/BpDzg9/ZWO1lS/h5ERERERERERERERERERER\nEREREREREREREREREREREZEk9v+ltQax0rVtyAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4343990>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This helps a lot but we are still not getting the same accuracy as <b>scipy.linalg.expm</b>.\n",
      "\n",
      "The extra accuracy of <b>scipy.linalg.expm</b> is because it implements the algorithm from [2] which:\n",
      "<ul>\n",
      "<li> uses Pade approximants instead of Taylor series.</li>\n",
      "<li> analyses the backward error of the function.</li>\n",
      "<li> uses cost analysis to perform the least possible number of matrix multiplications.</li>\n",
      "</ul>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Available Software"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is a lot of software available for computing the matrix exponential.\n",
      "\n",
      "In Python:\n",
      "- <a href=\"http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.linalg.expm.html\"><b>scipy.linalg.expm</b></a> - computes $\\exp(A)$ using [2].\n",
      "- <a href=\"http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.sparse.linalg.expm_multiply.html#scipy.sparse.linalg.expm_multiply\"><b>scipy.sparse.linalg.expm_multiply</b></a> - computes $\\exp(A)B$ where $B \\in \\mathbb{C}^{n \\times k}$ (useful for vectors $B$) using [3].\n",
      "- <a href=\"http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.linalg.expm_frechet.html#scipy.linalg.expm_frechet\"><b>scipy.linalg.expm_frechet</b></a> - computes the Frechet derivative of the matrix exponential using [1].\n",
      "- <a href=\"http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.expm_cond.html\"><b>scipy.linalg.expm_cond</b></a> (SciPy 0.14) - computes the condition number of the matrix exponential in the Frobenius norm using [11].\n",
      "\n",
      "In MATLAB:\n",
      "- <a href=\"http://www.mathworks.co.uk/help/matlab/ref/expm.html\"><b>expm</b></a> - built-in function evaluating $\\exp(A)$ (currently using an older algorithm than SciPy).\n",
      "- <a href=\"http://eprints.ma.man.ac.uk/1442/\"><b>expm_new</b></a> - implements algorithms from [2].\n",
      "- <a href=\"http://www.mathworks.co.uk/matlabcentral/fileexchange/29576-matrix-exponential-times-a-vector\"><b>expmv</b></a> - implements algorithms from [3].\n",
      "- <a href=\"http://www.maths.manchester.ac.uk/~higham/mftoolbox/\"><b>expm_frechet_pade</b></a> - computes the derivative of $\\exp(A)$ as explained in [11].\n",
      "- <a href=\"http://www.maths.manchester.ac.uk/~higham/mftoolbox/\"><b>expm_frechet_quad</b></a> - computes the derivative of $\\exp(A)$ as explained in [11].\n",
      "- <a href=\"http://www.maths.manchester.ac.uk/~higham/mftoolbox/\"><b>expm_cond</b></a> - computes the condition number of $\\exp(A)$ in the Frobenius norm as explained in [11].\n",
      "- <a href=\"http://www.maths.uq.edu.au/expokit/\"><b>expokit</b></a> - computes $\\exp(A)$ for small dense $A$ or $\\exp(A)B$ for sparse $A$.\n",
      "\n",
      "NAG Library (Fortran but callable in MATLAB, Python etc.):\n",
      "- <a href=\"http://www.nag.co.uk/numeric/FL/nagdoc_fl24/html/F01/f01conts.html\"><b>F01(E/F)CF</b></a> - Real / Complex matrix exponential using [2].\n",
      "- <a href=\"http://www.nag.co.uk/numeric/FL/nagdoc_fl24/html/F01/f01conts.html\"><b>F01(E/F)DF</b></a> - computes the matrix exponential of a real symmetric / complex Hermitian matrix.\n",
      "- <a href=\"http://www.nag.co.uk/numeric/FL/nagdoc_fl24/html/F01/f01conts.html\"><b>F01(G/H)AF</b></a> - Real / Complex versions of $\\exp(A)B$ using [3].\n",
      "- <a href=\"http://www.nag.co.uk/numeric/FL/nagdoc_fl24/html/F01/f01conts.html\"><b>F01(J/K)AF</b></a> - estimates the condition number of the matrix exponential for real / complex input using [1]."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Further Information"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Some further information not directly referenced here can be found at:\n",
      "- <a href=\"http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/\">Cleve's Corner</a> - discusses the accuracy of computing $\\exp(A)$.\n",
      "- <a href=\"http://en.wikipedia.org/wiki/Matrix_exponential\">Wikipedia</a> -\n",
      "a more detailed introduction to the matrix exponential.\n",
      "- <a href=\"http://eprints.ma.man.ac.uk/2102/\">Software Catalogue</a> - A catalogue of available matrix function software (currently version 1.0)\n",
      "- <a href=\"http://nickhigham.wordpress.com/2014/06/02/videos-of-lectures-from-gene-golub-siam-summer-school-2013/\">Nick Higham's Blog</a> - An short course on matrix functions (video lectures and written notes)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "References"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<a href=\"http://eprints.ma.man.ac.uk/1218/\">[1]</a> - A. H. Al-Mohy and N. J. Higham, Computing the Fr\u00e8chet Derivative of the Matrix Exponential, with an Application to Condition Number Estimation, <i>SIAM J. Matrix Anal. Appl.</i>, 30:1639-1657, 2009.\n",
      "\n",
      "<a href=\"http://eprints.ma.man.ac.uk/1217/\">[2]</a> - A. H. Al-Mohy and N. J. Higham, A New Scaling and Squaring Algorithm for the Matrix Exponential, <i>SIAM J. Matrix Anal. Appl.</i>, 31(3):970-989, 2009.\n",
      "\n",
      "<a href=\"http://eprints.ma.man.ac.uk/1426/\">[3]</a> - A. H. Al-Mohy and N. J. Higham, Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators, <i>SIAM J. Sci. Comput.</i>, 32:488-511, 2011.\n",
      "\n",
      "<a href=\"http://arxiv.org/pdf/0707.0756.pdf\">[4]</a> - Ernesto Estrada and N. Hatano, Communicability in Complex Networks, <i>Phys. Rev. E.</i>, 77(3):036111, 2008.\n",
      "\n",
      "<a href=\"http://epubs.siam.org/doi/abs/10.1137/090761070?journalCode=siread\">[5]</a> - Ernesto Estrada and D. J. Higham, Network Properties Revealed Through Matrix Functions, <i>SIAM Rev.</i> 52:696-714, 2010.\n",
      "\n",
      "<a href=\"http://arxiv.org/pdf/0905.4102.pdf\">[6]</a> - Ernesto Estrada, D. J. Higham, and N. Hatano, Communicability Betweeness in Complex Networks, <i>Physica A: Statistical Mechanics and its Applications</i>, 388:764-774, 2009.\n",
      "\n",
      "<a href=\"https://na.math.kit.edu/download/papers/acta-final.pdf\">[7]</a> - Marlis Hochbruch and Alexander Osterman, Exponential Integrators, <i>Acta Numerica</i>, x:209-286, 2009.\n",
      "\n",
      "<a href=\"http://www.cs.cornell.edu/cv/researchpdf/The%20Exponential%20of%20a%20Matrix.pdf\">[8]</a> - C. B. Moler and C. F. Van-Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, </i>SIAM Rev.</i>, 20:801-836, 1978.\n",
      "\n",
      "<a href=\"http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf\">[9]</a> - C. B. Moler and C. F. Van-Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later, <i>SIAM Rev.</i>, 45:3-49, 2003.\n",
      "\n",
      "<a href=\"http://www.tandfonline.com/doi/abs/10.1080/00207160.2013.765560#.U4XffHVdVNA\">[10]</a> - B. Garc\u00eca-Mora, C. Santamar\u00eca, G. Rubio, and J. L. Pontones, Computing Survival Functions of the Sum of Two Independent Markov Processes. An Application to Bladder Carcinoma Treatment, <i>Int. Journal of Computer Mathematics</i>, 91(2):209-220, 2014.\n",
      "\n",
      "<a href=\"http://www.maths.manchester.ac.uk/~higham/fm/\">[11]</a> - N. J. Higham - Functions of Matrices: Theory and Computation, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008, ISBN: 978-0-898716-46-7, xx+425 pp.\n",
      "\n",
      "<a href=\"http://www.vtt.fi/inf/julkaisut/muut/2010/nse_v164_n2_pp140-150.pdf\">[12]</a> - Maria Pusa and Jaakko Lepp\u01cenen Computing the Matrix Exponential in Burnup Calculations, <i>Nuclear Science and Engineering</i>, 164(2):140-150, 2010.\n",
      "\n",
      "<a href=\"http://www.cc.gatech.edu/~jarek/papers/SAM.pdf\">[13]</a> - Jarek Rossignac and \u00c0lvar Vinucua, Steady Affine Motions and Morphs, <i>ACM Trans. Graph.</i>, 30(5):Article 116, 16 pages, 2011."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}