{ "metadata": { "name": "Labs17_19_sympy" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab 17\n", "\n", "Key by Ryan Brunt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1" ] }, { "cell_type": "code", "collapsed": false, "input": [ "\"\"\"\n", "lab17key.py - written by Ryan Brunt\n", "\"\"\"\n", "import numpy as np\n", "\n", "\n", "def deriv(order, accuracy, type, func, vec, h=2e-5):\n", " \"\"\"\n", " Return the numerical derivative of a function\n", "\n", " Parameters\n", " ----------\n", " order : integer\n", " order of derivative. must be 1 or 2\n", "\n", " accuracy : integer\n", " number of terms in approximations. must be 1,2,3 for forward,\n", " or 2,4,6 for centered\n", "\n", " type : string\n", " type of numerical derivative to be calculated. must be 'forward'\n", " or 'centered'\n", "\n", " func : function\n", " function whose derivative you wish to calculate\n", "\n", " vec : array_like, dtype=float\n", " points to evaluate the derivative at\n", "\n", " h : float, optional(default=2e-5)\n", " Step size\n", "\n", " Returns\n", " -------\n", " a : array_like\n", " vector of derivative values calculated at the points in `vec`\n", "\n", " \"\"\"\n", " try:\n", " a = []\n", " q = vec.shape[0]\n", " if type == 'forward':\n", " if order == 1:\n", " if accuracy == 1:\n", " for i in range(q):\n", " p = ((-func(vec[i])+func(vec[i]+h))/h)\n", " a.append(p)\n", " elif accuracy == 2:\n", " for i in range(q):\n", " p = (((-3./2)*func(vec[i])+2*func(vec[i]+h)-(1./2)*func(vec[i]+2*h))/(h))\n", " a.append(p)\n", " elif accuracy ==3:\n", " for i in range(q):\n", " p = (((-11./6)*func(vec[i])+3*func(vec[i]+h)-(3./2)*func(vec[i]+2*h))+(1./3)*func(vec[i]+3*h))/h\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy: must be 1, 2, or 3'\n", " elif order == 2:\n", " if accuracy == 1:\n", " for i in range(q):\n", " p = (func(vec[i])-2*func(vec[i]+h)+func(vec[i]+2*h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==2:\n", " for i in range(q):\n", " p = (2*func(vec[i])-5*func(vec[i]+h)+4*func(vec[i]+2*h)-func(vec[i]+3*h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==3:\n", " for i in range(q):\n", " p = ((35./12)*func(vec[i])-(26./3)*func(vec[i]+h)+(19./2)*func(vec[i]+2*h)-(14./3)*func(vec[i]+3*h)+(11./12)*func(vec[i]+4*h))/(h**2)\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy: must be 1 2 or 3'\n", " else:\n", " print 'invalid order'\n", " elif type == 'centered':\n", " if order == 1:\n", " if accuracy == 2:\n", " for i in range(q):\n", " p = ((-1./2)*func(vec[i]-h)+(1./2)*func(vec[i]+h))/h\n", " a.append(p)\n", " elif accuracy ==4:\n", " for i in range(q):\n", " p = ((1./12)*func(vec[i]-2*h)-(2./3)*func(vec[i]-h)+(2./3)*func(vec[i]+h)-(1./12)*func(vec[i]+2*h))/h\n", " a.append(p)\n", " elif accuracy ==6:\n", " for i in range(q):\n", " p = ((-1./60)*func(vec[i]-3*h)+(3./20)*func(vec[i]-2*h)-(3./4)*func(vec[i]-h)+(3./4)*func(vec[i]+h)-(3./20)*func(vec[i]+2*h)+(1./60)*func(vec[i]+3*h))/h\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy:must be 2,4, or 6'\n", " elif order == 2:\n", " if accuracy == 2:\n", " for i in range(q):\n", " p = (func(vec[i]-h)-2*func(vec[i])+func(vec[i]+h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==4:\n", " for i in range(q):\n", " p = ((-1./12)*func(vec[i]-2*h)+(4./3)*func(vec[i]-h)-(5./2)*func(vec[i])+(4./3)*func(vec[i]+h)-(1./12)*func(vec[i]+2*h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==6:\n", " for i in range(q):\n", " p = ((1./90)*func(vec[i]-3*h)-(3./20)*func(vec[i]-2*h)+(3./2)*func(vec[i]-h)-(49./18)*func(vec[i])+(3./2)*func(vec[i]+h)-(3./20)*func(vec[i]+2*h)+(1./90)*func(vec[i]+3*h))/(h**2)\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy: must be 2, 4, or 6'\n", " else:\n", " print 'invalid order'\n", " elif type == 'backward':\n", " if order == 1:\n", " if accuracy == 1:\n", " for i in range(q):\n", " p = (func(vec[i])-func(vec[i]-h))/h\n", " a.append(p)\n", " elif accuracy ==2:\n", " for i in range(q):\n", " p = ((3./2)*func(vec[i])-2*func(vec[i]-h)+(1./2)*func(vec[i]-2*h))/(h)\n", " a.append(p)\n", " elif accuracy ==3:\n", " for i in range(q):\n", " p = (((11./6)*func(vec[i])-3*func(vec[i]-h)+(3./2)*func(vec[i]-2*h))-(1./3)*func(vec[i]-3*h))/h\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy: must be 1, 2, or 3'\n", " elif order == 2:\n", " if accuracy == 1:\n", " for i in range(q):\n", " p = (func(vec[i])-2*func(vec[i]-h)+func(vec[i]-2*h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==2:\n", " for i in range(q):\n", " p = (2*func(vec[i])-5*func(vec[i]-h)+4*func(vec[i]-2*h)-func(vec[i]-3*h))/(h**2)\n", " a.append(p)\n", " elif accuracy ==3:\n", " for i in range(q):\n", " p = ((35./12)*func(vec[i])-(26./3)*func(vec[i]-h)+(19./2)*func(vec[i]-2*h)-(14./3)*func(vec[i]-3*h)+(11./12)*func(vec[i]-4*h))/(h**2)\n", " a.append(p)\n", " else:\n", " print 'invalid accuracy: must be 1, 2, or 3'\n", " else:\n", " print 'invalid order'\n", " else:\n", " print 'invalid type'\n", " return np.asarray(a)\n", " except AttributeError:\n", " print 'you n00b. you didnt enter all the parameters'\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 2" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import pandas as pd\n", "from itertools import product\n", "\n", "\n", "def fun17_2(x):\n", " \"\"\"\n", " Analytical derivative is -3 * sin(x)\n", " \"\"\"\n", " return 3 * np.cos(x)\n", "\n", "\n", "def dfunc17_2(x, order=1):\n", " \"\"\"derivative of fun17_2\"\"\"\n", " if order == 1:\n", " return -3 * np.sin(x)\n", " elif order == 2:\n", " return -3 * np.cos(x)\n", "\n", "test_x = np.arange(-2, 2, .025)\n", "\n", "orders = [1, 2]\n", "acc = [1, 2, 3]\n", "h_vals = [1, 1e-3, 1e-5, 1e-7, 1e-9]\n", "der_types = pd.MultiIndex.from_tuples(list(product(orders, acc, h_vals)))\n", "results = pd.DataFrame(index=der_types, columns=['centered', 'forward'])\n", "results.columns.names = ['errors_by_type']\n", "results.index.names = ['Order', 'Accuracy', 'Step (h)']\n", "\n", "for ord_i in orders:\n", " actual = dfunc17_2(test_x, ord_i)\n", " for acc_i in acc:\n", " for h_i in h_vals:\n", " cent = deriv(ord_i, 2 * acc_i, 'centered', fun17_2, test_x, h=h_i)\n", " forward = deriv(ord_i, acc_i, 'forward', fun17_2, test_x, h=h_i)\n", " \n", " results.ix[ord_i, acc_i, h_i]['centered'] = np.abs(actual - cent).max()\n", " results.ix[ord_i, acc_i, h_i]['forward'] = np.abs(actual - forward).max()\n", "\n", "results" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
errors_by_typecenteredforward
OrderAccuracyStep (h)
111 0.4755828 1.458758
0.001 4.999959e-07 0.0015
1e-05 7.426371e-11 1.5e-05
1e-07 3.758555e-09 1.525071e-07
1e-09 3.465643e-07 6.126496e-07
21 0.08876399 0.9332589
0.001 6.346035e-13 9.999933e-07
1e-05 5.761036e-11 1.885256e-10
1e-07 4.868778e-09 1.073327e-08
1e-09 5.516318e-07 1.320662e-06
31 0.0176361 0.6691347
0.001 7.482903e-13 7.512357e-10
1e-05 7.683065e-11 2.793474e-10
1e-07 7.228002e-09 1.948056e-08
1e-09 6.662133e-07 2.166962e-06
211 0.2418138 2.768717
0.001 2.503239e-07 0.002999979
1e-05 9.967765e-06 3.802239e-05
1e-07 0.1004527 0.106455
1e-09 890.8112 1333.889
21 0.03049187 2.413696
0.001 1.609021e-09 2.750664e-06
1e-05 1.758813e-05 3.128172e-05
1e-07 0.1721033 0.3668233
1e-09 2107.803 3776.379
31 0.004586739 2.09567
0.001 2.134336e-09 9.640946e-09
1e-05 2.682928e-05 9.166569e-05
1e-07 0.2289471 0.7884678
1e-09 2158.798 8217.208
\n", "
" ], "output_type": "pyout", "prompt_number": 3, "text": [ "errors_by_type centered forward\n", "Order Accuracy Step (h) \n", "1 1 1 0.4755828 1.458758\n", " 0.001 4.999959e-07 0.0015\n", " 1e-05 7.426371e-11 1.5e-05\n", " 1e-07 3.758555e-09 1.525071e-07\n", " 1e-09 3.465643e-07 6.126496e-07\n", " 2 1 0.08876399 0.9332589\n", " 0.001 6.346035e-13 9.999933e-07\n", " 1e-05 5.761036e-11 1.885256e-10\n", " 1e-07 4.868778e-09 1.073327e-08\n", " 1e-09 5.516318e-07 1.320662e-06\n", " 3 1 0.0176361 0.6691347\n", " 0.001 7.482903e-13 7.512357e-10\n", " 1e-05 7.683065e-11 2.793474e-10\n", " 1e-07 7.228002e-09 1.948056e-08\n", " 1e-09 6.662133e-07 2.166962e-06\n", "2 1 1 0.2418138 2.768717\n", " 0.001 2.503239e-07 0.002999979\n", " 1e-05 9.967765e-06 3.802239e-05\n", " 1e-07 0.1004527 0.106455\n", " 1e-09 890.8112 1333.889\n", " 2 1 0.03049187 2.413696\n", " 0.001 1.609021e-09 2.750664e-06\n", " 1e-05 1.758813e-05 3.128172e-05\n", " 1e-07 0.1721033 0.3668233\n", " 1e-09 2107.803 3776.379\n", " 3 1 0.004586739 2.09567\n", " 0.001 2.134336e-09 9.640946e-09\n", " 1e-05 2.682928e-05 9.166569e-05\n", " 1e-07 0.2289471 0.7884678\n", " 1e-09 2158.798 8217.208" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab 19\n", "\n", "Key by Spencer Lyon" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from itertools import product\n", "from time import time\n", "import numpy as np\n", "import pandas as pd" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def jacobian(f, x0, h=1e-8, how='centered', *args, **kwargs):\n", " \"\"\"\n", " Numerically estimate the Jacobian matrix for the function f.\n", " If f operates :math:`\\\\mathbb{R}^n \\\\rightarrow \\\\mathbb{R}^m`, the\n", " return value of this function will be an (n, m) matrix of\n", " partial derivatives.\n", "\n", " Parameters\n", " ----------\n", " f : function\n", " The function for which you would like to calculate the jacobian.\n", " Must accept as an input a numpy array of the same size as\n", " x0 (also can accept other args passed as *args or **kwargs) and\n", " return a numpy array of values.\n", "\n", " x0 : array_like, dtype=float, shape=(n,)\n", " The array of values representing the point where the jacobian\n", " is to be evaluated.\n", "\n", " h : float, optional(default=1e-8)\n", " The step size to be used in calculating the derivative.\n", "\n", " how : str\n", " A string specifying the method of differentiation. Acceptable\n", " values are centered (c), forwards (f), backward (b).\n", "\n", " *args, **kwargs :\n", " Other arguments that should be passed to f for evaluation.\n", "\n", " Returns\n", " -------\n", " jac : array_like, dtype=float, shape=(n, m)\n", " The numpy array of partial derivatives representing the jacobian\n", "\n", " \"\"\"\n", " # Make sure we have an array\n", " x0 = np.asarray(x0)\n", "\n", " # Get n and m\n", " f0 = f(x0, *args, **kwargs)\n", " nx = x0.size\n", " mx = f0.size\n", "\n", " jac = np.zeros((nx, mx))\n", "\n", " ident = np.eye(nx)\n", "\n", " if how.startswith('f'):\n", " for xi in xrange(nx):\n", " ei = ident[:, xi]\n", " jac[xi, :] = (f(x0 + h*ei, *args, **kwargs) - f0) / h\n", "\n", " elif how.startswith('c'):\n", " for xi in xrange(nx):\n", " ei = ident[:, xi]\n", " jac[xi, :] = (f(x0 + h*ei, *args, **kwargs) -\n", " f(x0 - h*ei, *args, **kwargs)) / (2 * h)\n", "\n", " elif how.startswith('b'):\n", " for xi in xrange(nx):\n", " ei = ident[:, xi]\n", " jac[xi, :] = (f0 - f(x0 - h*ei, *args, **kwargs)) / h\n", "\n", " return jac.T" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "def test_fun(x):\n", " x1, x2 = x\n", " return np.array([np.exp(x1)*np.sin(x2) + x2 ** 3,\n", " 3*x2 - np.cos(x1)])\n", "\n", "def real_jac(x):\n", " x1, x2 = x\n", " return np.array([[np.exp(x1) * np.sin(x2), np.exp(x1)*np.cos(x2) + 3*x2**2],\n", " [np.sin(x1), 3]])\n", "\n", "x_grid = np.linspace(-1, 1, 100)\n", "y_grid = np.linspace(-1, 1, 100)\n", "\n", "max_err = 0.\n", "\n", "methods = 'cbf'\n", "h_vals = [1., 1e-4, 1e-5, 1e-6, 1e-7, 1e-10]\n", "\n", "ind = pd.MultiIndex.from_tuples(list(product(list(methods), h_vals)),\n", " names=['method', 'h'])\n", "jac_results = pd.DataFrame(index=ind, columns=['max_error', 'time_delta'])\n", "\n", "print 'Computing max error in my Jacobian for various how and h values'\n", "for der_type in product(methods, h_vals):\n", " how, h = der_type\n", " max_err = 0.\n", " for point in product(x_grid, y_grid):\n", "\n", " numerical = jacobian(test_fun, point, how=how, h=h)\n", " analytical = real_jac(point)\n", " this_err = np.abs(numerical - analytical).max()\n", " max_err = this_err if this_err > max_err else max_err\n", " \n", " # Time methods. Should be same for all points so just use 1\n", " t1 = time()\n", " _ = jacobian(test_fun, point, how=how, h=h)\n", " num_time = time() - t1\n", " \n", " t1 = time()\n", " _ = real_jac(point)\n", " ana_time = time() - t1\n", " \n", " jac_results.ix[how, h]['max_error'] = max_err\n", " jac_results.ix[how, h]['time_delta'] = num_time - ana_time\n", "\n", "jac_results" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Computing max error in my Jacobian for various how and h values\n" ] }, { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
max_errortime_delta
methodh
c1.000000e+00 0.9684898 0.000177145
1.000000e-04 9.667494e-09 0.0001790524
1.000000e-05 1.175886e-10 0.0001740456
1.000000e-06 5.293845e-10 0.0001809597
1.000000e-07 6.391792e-09 0.0001749992
1.000000e-10 4.380994e-06 0.0001809597
b1.000000e+00 3.826186 0.0001118183
1.000000e-04 0.0002845317 0.0001239777
1.000000e-05 2.845231e-05 0.0001139641
1.000000e-06 2.845191e-06 0.0001130104
1.000000e-07 2.870152e-07 0.0001199245
1.000000e-10 9.935029e-06 0.0001139641
f1.000000e+00 3.826186 0.000109911
1.000000e-04 0.0002845317 0.0001111031
1.000000e-05 2.845231e-05 0.0001130104
1.000000e-06 2.845191e-06 0.0001108646
1.000000e-07 2.870152e-07 0.0001130104
1.000000e-10 8.019226e-06 0.0001130104
\n", "
" ], "output_type": "pyout", "prompt_number": 3, "text": [ " max_error time_delta\n", "method h \n", "c 1.000000e+00 0.9684898 0.000177145\n", " 1.000000e-04 9.667494e-09 0.0001790524\n", " 1.000000e-05 1.175886e-10 0.0001740456\n", " 1.000000e-06 5.293845e-10 0.0001809597\n", " 1.000000e-07 6.391792e-09 0.0001749992\n", " 1.000000e-10 4.380994e-06 0.0001809597\n", "b 1.000000e+00 3.826186 0.0001118183\n", " 1.000000e-04 0.0002845317 0.0001239777\n", " 1.000000e-05 2.845231e-05 0.0001139641\n", " 1.000000e-06 2.845191e-06 0.0001130104\n", " 1.000000e-07 2.870152e-07 0.0001199245\n", " 1.000000e-10 9.935029e-06 0.0001139641\n", "f 1.000000e+00 3.826186 0.000109911\n", " 1.000000e-04 0.0002845317 0.0001111031\n", " 1.000000e-05 2.845231e-05 0.0001130104\n", " 1.000000e-06 2.845191e-06 0.0001108646\n", " 1.000000e-07 2.870152e-07 0.0001130104\n", " 1.000000e-10 8.019226e-06 0.0001130104" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 2" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def hessian(f, x0, h=1e-4, *args, **kwargs):\n", " \"\"\"\n", " Numerically estimate the Hessian matrix for the function f.\n", " It is assumed that f goes from :math:`\\mathbb{R}^n \n", " \\rightarrow \\mathbb{R}`.\n", "\n", " Parameters\n", " ----------\n", " f : function\n", " The function for which you would like to calculate the jacobian.\n", " Must accept as an input a numpy array of the same size as\n", " x0 (also can accept other args passed as *args or **kwargs) and\n", " return a numpy array of values.\n", "\n", " x0 : array_like, dtype=float, shape=(n,)\n", " The array of values representing the point where the jacobian\n", " is to be evaluated.\n", "\n", " h : float, optional(default=1e-4)\n", " The step size to be used in calculating the derivative.\n", "\n", " *args, **kwargs :\n", " Other arguments that should be passed to f for evaluation.\n", "\n", " Returns\n", " -------\n", " hes : array_like, dtype=float, shape=(n, n)\n", " The numpy array of mixed second ordre partial derivatives\n", " at the point x0\n", "\n", " \"\"\"\n", " # Make sure we have an array\n", " x0 = np.asarray(x0)\n", "\n", " # Get n\n", " nx = x0.size\n", "\n", " hes = np.zeros((nx, nx))\n", "\n", " ident = np.eye(nx)\n", "\n", " for xi in xrange(nx):\n", " ei = ident[:, xi]\n", " for xj in xrange(nx):\n", " ej = ident[:, xj]\n", " hes[xi, xj] = (f(x0 + h*(ei + ej), *args, **kwargs) -\n", " f(x0 + h*(ei - ej), *args, **kwargs) -\n", " f(x0 + h*(ej - ei), *args, **kwargs) +\n", " f(x0 - h*(ei + ej), *args, **kwargs)) / (4*h**2)\n", "\n", " return hes.T" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "def test_hes(vec):\n", " \"\"\"Rosenbrock banana function\"\"\"\n", " x, y = vec\n", " return (1 - x)**2 + 100*(y - x**2)**2\n", "\n", "\n", "def ana_ban(vec):\n", " x, y = vec\n", " return np.array([[1200*x**2 - 400*y + 2, -400*x],\n", " [-400*x, 200]])\n", "\n", "h_vals = [1., 1e-4, 1e-5, 1e-6, 1e-7, 1e-10]\n", "hes_results = pd.DataFrame(index=h_vals, columns=['max_error', 'time_delta'])\n", "\n", "hes_x = np.linspace(-2, 2, 100)\n", "hes_y = np.linspace(0, 2, 100)\n", "\n", "print 'Computing max error in my Hessian for various how and h values'\n", "for h in h_vals:\n", " max_err = 0.\n", " for point in product(hes_x, hes_y):\n", "\n", " numerical = hessian(test_hes, point, h=h)\n", " analytical = ana_ban(point)\n", " this_err = np.abs(numerical - analytical).max()\n", " max_err = this_err if this_err > max_err else max_err\n", " \n", " # Time methods. Should be same for all points so just use 1\n", " t1 = time()\n", " _ = hessian(test_hes, point, h=h)\n", " num_time = time() - t1\n", " \n", " t1 = time()\n", " _ = ana_ban(point)\n", " ana_time = time() - t1\n", " \n", " hes_results.ix[h]['max_error'] = max_err\n", " hes_results.ix[h]['time_delta'] = num_time - ana_time\n", "\n", "hes_results" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Computing max error in my Hessian for various how and h values\n" ] }, { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
max_errortime_delta
1.000000e+00 800 0.0003399849
1.000000e-04 3.043648e-05 0.0003130436
1.000000e-05 0.0029908 0.0003321171
1.000000e-06 0.2581446 0.0003361702
1.000000e-07 27.96591 0.0003869534
1.000000e-10 2.274182e+07 0.0003211498
\n", "
" ], "output_type": "pyout", "prompt_number": 5, "text": [ " max_error time_delta\n", "1.000000e+00 800 0.0003399849\n", "1.000000e-04 3.043648e-05 0.0003130436\n", "1.000000e-05 0.0029908 0.0003321171\n", "1.000000e-06 0.2581446 0.0003361702\n", "1.000000e-07 27.96591 0.0003869534\n", "1.000000e-10 2.274182e+07 0.0003211498" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "# sympy lab\n", "\n", "Key by Spencer Lyon" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import sympy as sym\n", "import matplotlib.pyplot as plt\n", "%load_ext sympyprinting" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print('7 appears %i times in the first 10000 digits of pi' %\n", " str(sym.pi.n(10000)).count('7'))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "7 appears 970 times in the first 10000 digits of pi\n" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 2" ] }, { "cell_type": "code", "collapsed": false, "input": [ "a, b, c, x, x1, x2, x3, h, y1, y2, y3 = sym.symbols('a, b, c, x, x1, x2, x3, \\\n", " h, y1, y2, y3')\n", "x2 = x1 + h\n", "x3 = x1 + 2 * h\n", "eqns = (sym.Eq(y1, a + b * x1 + c * x1 ** 2),\n", " sym.Eq(y2, a + b * x2 + c * x2 ** 2),\n", " sym.Eq(y3, a + b * x3 + c * x3 ** 2))\n", "soln = sym.solve(eqns, a, b, c)\n", "print('(a, b, c) = ')\n", "soln[a], soln[b], soln[c]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(a, b, c) = \n" ] }, { "latex": [ "$$\\begin{pmatrix}\\frac{2 h^{2} y_{1} + 3 h x_{1} y_{1} - 4 h x_{1} y_{2} + h x_{1} y_{3} + x_{1}^{2} y_{1} - 2 x_{1}^{2} y_{2} + x_{1}^{2} y_{3}}{2 h^{2}}, & \\frac{- 3 h y_{1} + 4 h y_{2} - h y_{3} - 2 x_{1} y_{1} + 4 x_{1} y_{2} - 2 x_{1} y_{3}}{2 h^{2}}, & \\frac{y_{1} - 2 y_{2} + y_{3}}{2 h^{2}}\\end{pmatrix}$$" ], "output_type": "pyout", "prompt_number": 8, "text": [ "(\n", " 2 2 2 2 \n", "2\u22c5h \u22c5y\u2081 + 3\u22c5h\u22c5x\u2081\u22c5y\u2081 - 4\u22c5h\u22c5x\u2081\u22c5y\u2082 + h\u22c5x\u2081\u22c5y\u2083 + x\u2081 \u22c5y\u2081 - 2\u22c5x\u2081 \u22c5y\u2082 + x\u2081 \u22c5y\u2083\n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " 2 \n", " 2\u22c5h ,\n", " \n", "-3\u22c5h\u22c5y\u2081 + 4\u22c5h\u22c5y\u2082 - h\u22c5y\u2083 - 2\u22c5x\u2081\u22c5y\u2081 + 4\u22c5x\u2081\u22c5y\u2082 - 2\u22c5x\u2081\u22c5y\u2083\n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " 2 \n", " 2\u22c5h ,\n", " \n", "y\u2081 - 2\u22c5y\u2082 + y\u2083\n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " 2 \n", " 2\u22c5h )" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 3" ] }, { "cell_type": "code", "collapsed": false, "input": [ "y = a + b * x + c * x ** 2\n", "\n", "soln[x] = (x1 + x2) / 2\n", "y_est1 = sym.simplify(y.subs(soln))\n", "print 'y(x1 + x2) / 2 ->\\n'\n", "y_est1" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "y(x1 + x2) / 2 ->\n", "\n" ] }, { "latex": [ "$$\\frac{3}{8} y_{1} + \\frac{3}{4} y_{2} - \\frac{1}{8} y_{3}$$" ], "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAHUAAAAeCAYAAAASN1ElAAAABHNCSVQICAgIfAhkiAAAA+5JREFU\naIHt2kmIHFUcx/FP4hJNMgpelMRlMl4UNUJwRU0GNzwEIWEMaFzQiASUJK4H9RAPwhgVRUXBeGgX\nYhBXMBqIhxE8xIMYd4gg8SIGPahxwajEw78aa9rqqeru6upu7S80Xa/eq/r/+q3/9/7NkP8cB6Wu\nlyafs3AbfsJXvRDVIoOkewRbsBM/VmHwe1yXXF+BXxMR/c6g6L4RG3EAo1UZPRXzkusJ/KE/K6eR\nQdNdaaOm2YJ7emG4QwZBd9cb9eCG9BJciF/wSDcNl8yg6q6Um/AB5vfI/mL/7nBF6LXuIlQ2/Z6D\nvViUpE9KjE9UYTyDmmI/vN90F6HrjTo7+f4Tn+GbJD0mHI5d3TReAoOquzKuxu24Fa/g0owyo3gc\n23BVQ97NeKckLTXFe3MR3VSnvRmr8ZQYqVtxS5NyoyrW+SQOxXp83JD3vhBbBjXlT1FVae+UjnXO\nziuQ4ny8h/24DLtTefOEB/puC++rkkHRXrnOYzAHC/EXVqTyLhHTyikl2aopd6RWqb0TStHZyrbh\n2+R7FX7GW6m8C8Rx3ectvA+exekZ948XZ7n7M/LWiG1LK5Sl/TQ8jVkF7e7C2oJlKbmOD+R80mzH\n6w33pvBqw70RvCwaqFVqio3UPN3taD9bOF0bsUMECzqlbJ1LcY04T34BF3cq8As8kErPwW/YkLrX\n6eF1rc3n8sjTPh+TqfxVIkCwsAtaZiJP54xBjFYcpTpf46hUehKHmb6APyMatd/I0z6Gu3Bikt6O\nw3FeVQIT8nSOi1mQGDiHpB9Or6mLsBy/42i8iQ8zDG4QjfaE6CFniNjgR+3/htJ4GK8JDzKLPO2f\n4Fz/xGOPS76/LEFb0fotovPTVNmVYgDty3rRg6nrWSLikcdssbg/1yS/yul3mTgyHC9YPk87PC86\nShm0U78017kEd2Az5jZ7eLeITdbZnFHmRdM3xCuEh3pyk3e226iPCfe+KEdinXAmxpuUaVX7GmxS\n3NPNo0j90rrOGYMY65OHN+Fu2VuN78QRFiwQ09TqJsaoLiJxp3AmpjRv1Fa0LxeNSqxloyVoLFK/\nRXTmBjHSa+pW4c5PCOdgZ4bBtTgTD4mRdKU4uuoll+NtsVbNRFHty8Saty0pV6/EPR3qLFK/RXQW\nDmKMCAdjruiZ9+EH7ffQoofXnbIA16bSU4qvqVmMCYejcQ95RAfvpPz6nTGIUV8vVopDgkdTefcL\nb+ulNg1XwfWiN9c37+uEV/mGGGn9Qk/qd7HwONNM4oRuGewSe3Q2UrtFpfVb/9/vXuFo3CAOjMfF\nXq/ZPqrfOBb34iKxHu4zPcLRawa9focMGTJkyJD/JX8DEq9C4wMoYF0AAAAASUVORK5CYII=\n", "prompt_number": 10, "text": [ "\n", "3\u22c5y\u2081 3\u22c5y\u2082 y\u2083\n", "\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 - \u2500\u2500\n", " 8 4 8 " ] } ], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "soln.pop(x)\n", "\n", "soln[x] = (x2 + x3) / 2\n", "y_est2 = sym.simplify(y.subs(soln))\n", "print '\\ny(x1 + x2) / 2 ->\\n'\n", "y_est2" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "y(x1 + x2) / 2 ->\n", "\n" ] }, { "latex": [ "$$- \\frac{1}{8} y_{1} + \\frac{3}{4} y_{2} + \\frac{3}{8} y_{3}$$" ], "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAAIkAAAAeCAYAAAAYRz0yAAAABHNCSVQICAgIfAhkiAAABCdJREFU\neJzt20mIXEUYwPHfxCWaRcGLIXGZjBdBjRBcUZPBDQ/BYIyCxgWdIAHFGLeDeogHcXAhgqKgQdqF\nGMQtYmJAD+NyiLjFHRckXtSgB02MaFTi4es2L8309Hud6uVJ/6HpV13Vr//zdb16X1X10KdPE/ZJ\neK7pWINN+DXhedvNvOrjJNyIbfi2q0b5KJ33UqzELgx21aQ4P+PK6vFF+F10+F6nrN6l7CTHYmr1\neDH+Uo5gl9W7lJ0kyxrc3m2JFmir977tOnHJmIszsQOruuxShFJ698pIMkdrF8A1eB/T0uoUohX3\nXvDOTa90kop8HqdgK2ZXy0eLv2FxW6zyUdHcvaPek9px0hLxNz7D99XykEgAN3fNKB8d9U7VSZbg\nkerxKK5r0G4QD2I9Lq2ruxavJ/LJy3t4XPiuwAgW4Ju6doPK6U3vuTflYeyP5fi4ru4drE30ORVp\nb3ud8qYH3Tt5uzkdb2MnzsNXmbqpIlN/o4M+eSmrNyV0n4HJmIV/cEGm7hyReB2T6LMq0l2NnfSm\nB92zU63j8CgGcgpsxrKcbeHH6vPF+A0bMnVniGXmzwucD57A8eO8foTY09g5Tt2ImC7mpR3elNu9\nMLtyPLJsxEt1r43hhUx5Op4TAWuFinxXYxH3PN4ni0RyJV4TG3BFqWju3o6Yz8PlYh/uaZxdq0iR\nkwzkeGQ5El9mypNFcN+slpfiJlyYyG8iirg3854mhvNVopM8Jr6cWV32zuNOdJhJWI11eFl1L6gb\n6yTf4ZBMeRQH2J1ArRZB7jWaeQ/hVhxVLW/EgTitU4IT0MwdhsXoTYxE+9UqUu3dzBbz9D9xKF7B\nhw3a3iA6wkNie/sE8fuTjxK57A3340UxI6inmfcnONXu33QcXn3+uk2uqWP+aeZ4kbhQt6cUvjdz\nPCB2JfMwSSRXT45TtzdL/JUW3jtfLHUP52g7kXeNp0SnK0pFPvd2xHwubha3yinZN6Rgofh9A/Hl\n7mjQ7hl7LugsFMPg3Yk8amzDHwXaHyxmGl80qC/qPYIfRMCLkte9HTH/APfhXbwl8YbhcjFluwe3\nGX9qBz+JJWKYKYbmJQ3adnKz8BaRzI0ZfyQp4r1AdBLivj+YyLGelDGfcMMwVU6yVmTLi0WytqlB\nu2U4UfTWGbhELA13k/Pxqri3NyKv93yRH6yvtqsFf0s63f9IGfO2bxhOF8neFHHl3IlftH4F1TYL\nd4lANNosTMFMXJEpj8mXk4zHkEj06tcrDmpdryGpYw6XiaWHFXge59Yq8q6uTsQisej1QOa1u0Tm\n/GyC87eTq8TVVVt8ul7MEtaJ0aBXKV3M54iMPMuoWMApG1u0PpJ0ko7GPMX/3WwVSd/VYrNoWKwz\nNJqz9yKH4Q6cJXKK7fbcMe01/g8x79OnT58+ffr0KcS/hD5Noa4qiRsAAAAASUVORK5CYII=\n", "prompt_number": 11, "text": [ "\n", " y\u2081 3\u22c5y\u2082 3\u22c5y\u2083\n", "- \u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\n", " 8 4 8 " ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 4" ] }, { "cell_type": "code", "collapsed": false, "input": [ "t, c, b, g, h = sym.symbols('t, c, b, g, h')\n", "\n", "P = sym.Function('P')\n", "m = sym.symbols('m', real=True)\n", "sol = sym.dsolve(sym.Derivative(P(t), t) - m * (c + b * P(t) - g - h*P(t)), P(t))\n", "\n", "sol" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\operatorname{P}\\left(t\\right) = \\frac{- c + g + e^{C_{1} + b m t - h m t}}{b - h}$$" ], "output_type": "pyout", "png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAAkCAYAAAApQLALAAAABHNCSVQICAgIfAhkiAAABdxJREFU\neJzt3HuMXVUVwOGPaYUCaSivFguhRVJ8NBblJeVVpqa8UoiIZkRIJFgpAZTWUkhEQxNjoiYYiYBo\nKyG8yvBKUARFbWrig4IEA5EEQwyP8KhYEAotEKD8sfZ470zPvWfOnFt7Z2b/kpt7zt5nr73mzj5r\nr732OofMeOQj21uBTtGzvRXI1OIIfB19+BK+iROG0e64grJJuB0nVtThU1gyzGuPxwMV5bftY+II\nhGW6g0swAT9oKrsS141Q3lvp+48V283H3ypcf19F+W37yAN4dDIHp+OYIeVXYdMIZe6E/XG0sOL3\nYSEexHqciV9iGl7DPan/PjyLKfhvSR+n4EWchAPwRN0+8gAenZyIO7ClqaxH/NNbMRsL0vERYjDA\ne7hGuBU34Q9YhtuwL9ZgI76N3yYZ+6e2D+Fp3InJWFTQ7/O4Px1/HJdjprhR7q/YR6YmF3WJzD5c\nMaTsVOySjk8TVqwV5xSUXYpDsKsYxBNxd6o7NtXDtZiFPTFP+KYzlRvD/fDTdPwtHIkZdfvIi7hq\n7NUlMm8X/7tv4HM4WUzDmzBVDNAdKspcg0NxXvocjLWpbi7uTccbhTXfIGaALcIdeLdE/nFNMnZK\nes6s20fVP3Ks0SPu7vfF9DtJw0oUsSJ92nE0zsDD+AT2xvk1ZPZgKd7Bq8KyXlmiwwqxmHupRf1Z\nuKVExqhgvFvg68S09GO8gi/XlHcQVorpfbWwrk/WlPkzoeNPcBf2qCmPMTJ4GZuLuAu0D9Q/jH4c\nLlby+6Xy3+OvTddNE9a5eZY6RljpATbie03n38eNqZyIFqysIfNjOFtY8LPE1NscNitiKj6KXnET\njUtW4zHhe7wpgs+/we9S+ToR8oDpBv8DitgLu28TTUfOcrGSr8KKkvo3xI1BrMpfVj7LtZPZp1qM\nddzR6sc9U2M6vVbEBU8S4Y05+BFuFYuIMzSC4AMciBfEKpNwxpfiQ51SvAM8L27OZk6vKfMlMWjh\nC2Jh9X4Nef80+LfdQYSqxvva5X+0cyGOTd9rC+r6RezwKmFphnKq8NXWp/MtuEEM4h+OQM9twWp8\nEovF6n0XjVXwSLlYuDCP4ysdkPeo8HsvFj76zviVwfHfTAtuE9ZjSkFdDzanTxF3Kd6SvB47dkS7\n7cNlFa59UMOd6JTMTAVeEP5uEUcKK/Bei/r1+G5B+VJ8tr5qXcly4bMS4bPH5CjPNqfVDzwLH8af\nWtQvEdZ3flNZn1jorRMr4d50fmHTNX83PKs0GtlHhNHOxtfEvn8d/zczDFr5wAPpdkMH8GR8R6S3\nrTN4hdyfPotT/QJbuxgbNMJWRVwvtjOrsESxn/7/ZlnT8c3bTYtxRqsBPC99L9QYzBPS9Q+IhIxf\n4+2Ctr0i1lrkH2/W3gc+t0TfTGYQ7Szwy9rvTG3Abum7mePxixZtdsN/KuiXybSlyAeekT5lic1P\niU2MZmaL3aZWbaendplMRyiywAMuw9qStmtEBtHjTWW9ImPoL+l8ivCbn0vncw3eWh3KSny6pN+h\nLNO4YXJ8NGOVGAizS66biKuHlPUL/3eAK0Q0g9g9+nknFMxkBihyIXqFn/qPkrbviihFc9Rggsie\nJ8Jlm8UjJESS9a0jVTSTacd0kajzhLC+b4tpefkw2n5VI5nnYPxZ5KteqnGD7K2R/JMp53yRlTa3\nhoxpIjPuXx3RKJOpwGT8W8xodVikfYL+qCdvdXYn88UM2GqrfricIB6SHLPkAdydLBCuXJ94amTO\nCGT0iJj87iJNdJXiF5pkMh3nSXwmHS8U70qoyuF4RuMRpAtECuyYYiw+UjTamSFCjuvS+T5iEczw\nH5ci3IdVIo8YDsMjHdU0kyng8yKpaYB+W78DYjisxVHpeGeRGjBV9z3aVYvsA3cfr2u8PmmW8H/L\nHqMfyq5iI+qhdH6ysOib8MUO6Ng11A3TZDrP02IRN108V7fI1glTZRwmfN+B1zFNEhtO04R1L3sJ\nSSaTyWQymUwmk+lSPgBz7CJxXnrmlAAAAABJRU5ErkJggg==\n", "prompt_number": 12, "text": [ "\n", " C\u2081 + b\u22c5m\u22c5t - h\u22c5m\u22c5t\n", " -c + g + \u212f \n", "P(t) = \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " b - h " ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 5" ] }, { "cell_type": "code", "collapsed": false, "input": [ "b, c, g, h, m= (-1.2, 2, 0.1, 1.1, 1.1)\n", "\n", "\n", "def p_t(t, p0=1.):\n", " \"\"\"\n", " I used the sympy solution to write this by hand.\n", " I could have done it with sym.solve and soln.subs,\n", " but this was easier.\n", " \"\"\" \n", " # Constant of integration\n", " cc = c - g - (h - b) * p0\n", " \n", " term1 = (c - g) / (h - b)\n", " term2 = cc * np.exp(m * (h - b) * t) / (h - b)\n", " return term1 - term2\n", "\n", "\n", "times = np.linspace(0, 1, 100)\n", "pbar = np.ones_like(times) * (c - g) / (h - b)\n", "plt.plot(times, p_t(times, 0.5), 'r', label=r'P(0) = 0.5')\n", "plt.plot(times, p_t(times, 1.3), 'g', label=r'P(0) = 1.3')\n", "plt.plot(times, pbar, 'b', label='Equilibrium value')\n", "plt.legend(loc=0)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 13, "text": [ "Legend" ] }, { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdYFHf+B/A3CCplwaX3Ir1IEYwt6mo0dk00GsvZYrlf\nvMRoqpfkEnLp7YzG82KM5aKmqicqtkiyamyoKBYUFUGaSO8d5vfH6CIKiuzCLMv79TzfZ2bZ3ZkP\n8yzvHb5TvnqCIAggIiKdoC91AUREpDkMdSIiHcJQJyLSIQx1IiIdwlAnItIhDHUiIh2idqh//PHH\nCAgIQI8ePTBt2jRUVlZqoi4iImoBtUI9OTkZa9asQWxsLM6fP4/a2lr89NNPmqqNiIgekYE6bzYz\nM4OhoSHKysrQqVMnlJWVwdHRUVO1ERHRI1Ir1C0sLPDKK6/AxcUFRkZGGD58OIYOHap6Xk9PT+0C\niYg6opZe7K9W90tiYiK++uorJCcnIyMjAyUlJdi8efN9hbEJePfddyWvQVsatwW3BbfFg5s61Ar1\nU6dOoV+/frC0tISBgQEmTJiAo0ePqlUQERG1nFqh7uvri+PHj6O8vByCIODAgQPw9/fXVG1ERPSI\n1Ar14OBgzJw5E+Hh4QgKCgIALFiwQCOF6RqFQiF1CVqD26Iet0U9bgvN0BPU7cB50ML19NTuHyIi\n6mjUyU61zn5pKQsLC+Tn50uxatJBcrkceXl5UpdBpBUk2VPnHjxpEj9PpGvU+Uzz3i9ERDqEoU5E\npEMY6kREOoSh3gKrV6/GkiVLmvXalStXYunSpa1cERGRiKHeCDc3NxgbG0Mmk8HOzg5z5sxBaWkp\nAKCqqgoffvghXn/9ddXrz549i7CwMJiYmCA8PBxxcXGq5+bPn4/NmzcjOztb43U+aL33mj17Nrp0\n6QKZTAaZTAYzMzMeXCTSQQz1Rujp6WHXrl0oLi5GbGwsTp06hQ8++AAAEBkZCT8/P9jb2wMQQ378\n+PGYOXMmCgoKMGvWLIwfPx7V1dUAgC5dumDkyJH4/vvvNVrjw9bb2O/0xhtvoLi4GMXFxSgqKuIN\n14h0EEP9IRwcHDBixAhcvHgRALBnzx4MGjRI9bxSqURtbS1eeuklGBoa4sUXX4QgCPj9999Vr1Eo\nFIiKitJoXc1Z7724Z06k+xjqTbgTgKmpqdizZw9CQ0MBABcuXICPj4/qdRcvXlTdIuGO4OBg1ZcA\nIN4j50FdI0FBQZDL5Y22F154odH3NGe991q1ahUsLS0RHh6Obdu2Nfk6Imq/JLmitFk00TXQwj1T\nQRDw1FNPwcDAAObm5hgzZgzefPNNAEBBQQFkMpnqtSUlJTA3N2/wfjMzMxQXF6sey2QyFBYWNrm+\nc+fOPXKNzVnv3RYtWoR//etfMDc3x759+/Dss8/Czs4O/fr1e+R1E5H20t5Ql7CrQE9PD5GRkRgy\nZMh9z8nlchQVFakey2SyBo8BoLCwEGZmZqrHxcXF9wWwupqz3rvd+U8DAEaOHInp06dj27ZtDHUi\nHcPul0cUFBSEK1euqB4HBATct6d97tw5BAQEqB5funQJISEhTS4zICBAdVbKvW3hwoVNvudh6yWi\njoeh/ohGjRqFgwcPqh4rFAp06tQJK1asQGVlJVasWAF9ff0Ge/kHDx7EyJEjm1zmxYsXVWel3NtW\nrVrV6Huas967bdmyBSUlJairq8P+/fuxefNmjBs3roVbgYi0ltCKmlp8K69WbW5ubkJ0dHSjz1VV\nVQkuLi5CRkaG6mdnzpwRwsLCBCMjIyEsLEw4e/as6rny8nLByclJyMrK0nidD1rvpk2bhICAANXj\nAQMGCObm5oKZmZkQEhIi/PzzzxqvRyra/nkielTqfKZ5l8YWWLNmDeLj47Fs2bKHvnblypVIS0vD\nJ5980gaVdUzt/fNEdC91PtMMdWr3+HkiXcNb7xIR6YgbBTfUej9DnYhISxy+cRh91vZRaxkMdSIi\nLbA2di0m/jIRG8ZvUGs52nvxERFRB1BdW40l+5bgt+u/4fCcw/Cx8nn4mx6AoU5EJJGs0ixM+nUS\nZJ1liJkXA/Ou6l95zu4XIiIJnMo4hV5remGAywDsmLpDI4EOcE+diKjNrY1di6XRS/HN6G8w0X+i\nRpfNUG+B1atX4/Lly7z4iIgeSUVNBRbtWYTDKYdxaPYh+Fn7aXwd7H5pRHsZzm7BggXw9fVFp06d\n8N///veBr3399dfh4uICMzMzODk54eWXX0ZNTY3GayKixiUXJGPA+gHIK89DzLyYVgl0gKHeqPYw\nnB0AhISEYNWqVejZs+dDh6abO3cu4uPjUVRUhJiYGOzfvx/fffedxmsiovvturILvb/rjamBU/Hr\npF8h6yJ7+JtaiKH+ENo6nB0ALFy4EEOGDEHXrl0f+lofHx+YmpoCEAcB0dfXV30xEVHrqKmrwZvR\nb+L5qOexbfI2vNz35VYfG1jtUC8oKMAzzzwDPz8/+Pv74/jx45qoS3KClg9n1xKffPIJZDIZnJ2d\nMWbMGIwfP15jyyaihtKK0jDkv0NwKuMUTi84jf4u/dtkvWofKH3ppZcwatQobNmyBTU1Naq+Z3Xp\nvaf+t5nwru4OZ9cSS5cuxdKlS3HmzBk89dRTCA8Px4QJE9pk3UQdye6ru/Fc5HNY1HsRlj6+FPp6\nbdcpolaoFxYW4vDhw6qDdHdCUBNaGsia0B6Gs1NHaGgoFi5ciI0bNzLUiTSoqrYKb//+Nn688CN+\nnfQrBrgOaPMa1Ar1pKQkWFtbY86cOYiLi0NYWBiWL18OY2Nj1WsiIiJU8wqFAgqFQp1VSq6x4ey+\n/PLLBq85d+4cXnzxRdXj5gxnl5KS0uhzM2bMaHL0I3VUV1fDxMRE48sl6qiu5V3D1K1TYW9qjzN/\nPQMrY6tmv1epVEKpVGqmEDUG5xBOnjwpGBgYCDExMYIgCMJLL70k/OMf/1A939Ti1Vxtq3vQyEfb\ntm0TnnzySdXjqqoqwdXVVVi+fLlQUVEhLF++XHBzcxOqq6tVr5k/f77w+eefa7zOqqoqoby8XOjX\nr5+wZs0aoby8XKirq7vvdXV1dcI333wj5OfnC3V1dcKJEycEe3t7YevWrRqvSQra/nki3VZXVyd8\nf/Z7weozK+HrE183+jf4qNT5TKv113Dz5k3Bzc1N9fjw4cPC6NGjH1qYtv8Rtpfh7AYNGiTo6ekJ\n+vr6gp6enqCnpyccPHhQEISGw9nV1tYKI0aMECwsLASZTCYEBgYKa9eu1Xg9UtH2zxPprvzyfGHK\nlimC30o/4ezNsw9/QzOp85lWe+SjgQMH4rvvvoO3tzciIiJQXl6OTz/9FIDujnzE4ey0S3v/PFH7\ndDD5IGZun4lxPuPw2dDPYGRopLFlSzqcXVxcHObNm4eqqip4eHhg/fr1qoOCuhrqpF34eaK2VFlT\niX/88Q9sOrcJ3437DqO8Rml8HRyjlDo0fp6orcRlxmHG/2bA08ITq8eshrWJdausR53PNG/oRUT0\nEDV1Nfji6Bf417F/4Ysnv8CMoBmtfmVoSzHUiYge4HLOZczePhumnU1xcv5JuHZzlbqkB+K9X4iI\nGlFbV4svj36JAesHYFbwLPw24zetD3SAe+pERPeJz47Hc5HPoatBV5yYdwLd5d2lLqnZuKdORHRb\ndW01Pjr8EQauH4hZwbPw+6zf21WgA9xTJyICII4ZOm/HPNia2uL0gtPtoqulMdxTb0MpKSmQyWSq\nU5UUCgXWrl0LANi8eTOGDx+ueq2+vj6uX7/e7GXf+35tlJycDH19fdTV1UldCpFKaVUpXtn/Csb8\nMAav9nsVe6fvbbeBDjDUG3X3cHZ32qJFi9RerouLC4qLi1WnQunp6anmp0+fjn379rV42eq+n6gj\n2n11NwL/E4hbJbdw/vnz+EvQX7T2VMXmYvdLI+4MZ9fYrXe1UW1tLTp16iR1GUTtRkZxBhbvXYzY\nm7FYPWY1nvR4UuqSNIZ76o+orq4Or776KqytreHh4YF///vfDboU3NzcEB0drXp9REQEZsyYAeDB\n3Q8bNmzAgAEN770cFRUFDw8PWFtb4/XXX1d122zYsAH9+/fHyy+/DCsrK0RERDR4f2Prubur5+73\ny+VyeHp64ujRo1i/fj1cXFxga2vb5JiqP//8M3r16tXgZ8uWLVONohQVFYXQ0FCYm5vDxcUF7733\nXpPb8kHbCgCOHz+Ofv36QS6XIyQkBAcPHmxyWUTNUVNXgxUnViD4m2B4W3rj/PPndSrQAYZ6k5q6\nRPfbb79FVFQUzp49i1OnTmHLli0N/l27u0vlzuOW2r59O06fPo3Y2FhERkZi3bp1qudiYmLg4eGB\nrKwsvPXWWw9d1r11xcTEIDg4GHl5eZg6dSomT56M2NhYJCYmYtOmTXjhhRdQVlZ233LGjRuHhIQE\nXLt2TfWzH374AdOnTwcAmJqaYtOmTSgsLERUVBT+85//IDIyslk13T2fnp6OMWPG4J133kF+fj6+\n+OILTJw4ETk5OQ/9XYkaczztOHqt6YXtl7fj0OxD+GDIBxq9CZe20NpQ19NTv7WUcHs4u7vHCr2z\nl/vLL79gyZIlcHR0hFwux5tvvvnAezSoc0+SN954A926dYOzszMWL16MH3/8UfWcg4MD/va3v0Ff\nX79ZA0/fy93dHbNmzYKenh4mT56MjIwMvPPOOzA0NMSwYcPQuXPnBsF9h5GREcaPH6+q5erVq0hI\nSMC4ceMAAIMGDUJAQAAAoEePHpgyZUqz97Dv3labNm3CqFGjMGLECADA0KFDER4ejt27dz/y70od\nW3ZpNubtmIeJv0zEa/1eQ/TMaPhZ+0ldVqvR2lAXBPVbS90Zzi4/P1/V5s6dCwC4efMmnJ2dVa91\ncXFR91dt0r3rycjIaPS5lrC1tVXNGxmJeyvW1tYNflZSUtLoe6dNm6YK9R9++AFPP/206ovlxIkT\nGDx4MGxsbNCtWzesXr0aubm5j1zfjRs38Ouvvzb4Yj1y5AgyMzMfeVnUMdXU1WBlzEr4r/KHWRcz\nxC+Mx7Qe09r9gdCH4YHSR2Rvb99g6Ll7h6EzMTFpMPi2OiGUkpICPz8/1byjo6PquQd9MO8MU1dW\nVgZTU1O167jX0KFDkZ2djbi4OPz000/46quvVM9NmzYNixYtwr59+9C5c2csWbKkyS6TxrbVnd/L\nxcUFM2bMwLfffquxuqnj+CPpDyzetxhWxlZQzlIiwCZA6pLajNbuqUutqW6TyZMnY8WKFUhPT0d+\nfj4++eSTBgEbEhKCn376CTU1NTh16hS2bt3a4j2DL774AgUFBUhNTcWKFSvw7LPPNut91tbWcHR0\nxMaNG1FbW4t169YhMTGxRTU0xtDQEJMmTcKrr76K/Px8DBs2TPVcSUkJ5HI5OnfujJiYGPzwww9N\n/v6Nbas7/vKXv2Dnzp3Yv38/amtrUVFRAaVSifT0dI39HqR7rudfx8RfJuK5Hc/hHwP/gQMzDnSo\nQAcY6k0aO3Zsg/PUJ06cCACYP38+hg8fjuDgYISHh2PixIkNvgDef/99JCYmQi6XIyIiQnUA8Y6m\nAu7eg4YAMH78eISFhSE0NBRjxoxRdQE19tp7f7ZmzRp8/vnnsLKyQnx8PPr37//AdT3qF8+0adMQ\nHR2NSZMmQV+//mO0atUqvPPOOzAzM8P7779/3xfR3et50LZycnJCZGQkPvroI9jY2MDFxQVffvkl\nL1yiRhVVFuHv0X/HY2seQ0+7nohfGI9n/J/R+a6WxnCQDDUlJyeje/fuqKmpaRBu1HZ06fNEj6am\nrgZrY9ci4mAERnqOxAdDPoCDzEHqstTGQTKIqEMRBAFRV6PwxoE3YGNig93TdiPUPlTqsrQCQ10D\nOuK/eERSOZVxCq/99hpuldzCZ8M+w2iv0fwbvAu7X6jd4+epY7iaexVv/f4W/kz5ExGKCDwX+hwM\n9HVzv1SdzzQ7gYlIq90svonno55H37V9EWIXgqsvXsWCsAU6G+jq4lYhIq2UW5aLT498irVn1mJO\nyBwkvJAAS2NLqcvSegx1ItIqhRWFWHZ8GVbGrMTkgMk493/n4Gjm+PA3EgCJQl0ul/PABmmMXC6X\nugTSgOLKYqw4sQJfnfgKo7xGIWZ+TLsbSk4bSBLqeXl5UqyWiLRQcWUxVsasxFcnvsKw7sNw5Lkj\n8Lb0lrqsdovdL0QkicKKQnwd8zVWnFiBYR7DoJyl1Om7J7YVhjoRtancslwsP7Ecq06uwkivkTg8\n5zB8rHykLktnMNSJqE1kFGdg2fFlWHdmHSb4TcDxecfhaeEpdVk6R+3z1GtraxEaGoqxY8dqoh4i\n0jHX8q5hwc4FCFwViKraKpz961msGbuGgd5K1N5TX758Ofz9/VFcXKyJeohIR8Skx+Dzo59DmazE\n8+HP48qLV2BlbCV1WTpPrT31tLQ07N69G/PmzeNl2kSEOqEOUVeioNigwKRfJ6G/c39cX3Qd/xz8\nTwZ6G1FrT33JkiX4/PPPUVRU1ORrIiIiVPMKhQIKhUKdVRKRFiqvLsfGcxux7PgydDXoilf6voJn\nA56FYSdDqUtrF5RKJZRKpUaW1eIbeu3atQt79uzBv//9byiVSnz55ZfYuXNnw4XzRktEOi2jOAP/\nOfUffHv6Wzzm+Bhe7vMyFG4KXlyoJknup3706FHs2LEDu3fvRkVFBYqKijBz5kx8//33LV0kEbUT\nJ9NPYkXMCkRdicLUHlNxaPYhnpaoJTRy692DBw/iiy++4J46kQ6rrKnEr/G/4uuYr5FVmoW/9fob\n5obOhdyIt2nQNK0Y+Yj/bhHpppTCFKw+vRprY9eih20PvDXgLYz2Go1O+p2kLo0aIckgGUSk3Wrr\narE/cT++Of0N/kz5EzOCZuD/wv8Pvla+UpfWIaiTnQx1IlLJLMnEujPr8O3pb2FlbIW/hv0V03pM\ng0lnE6lL61C0ovuFiNqn2rpa7EvchzWxa6BMVuIZ/2ewZfIWhDuES10atQD31Ik6qOv517Hh7AZs\nOLsBdqZ2mN9zPqYEToGsi0zq0jo87qkTUbOUVpVi26VtWH92Pc5nncf0HtOxc+pOBNsFS10aaQj3\n1Il0XJ1Qh8M3DuO/cf/F/y7/D/2c+2F28GyM8xmHLgZdpC6PGsEDpUR0n4ScBGw8txGbzm2CaWdT\nzA6Zjek9psNeZi91afQQ7H4hIgDi2Ss/XfgJm89vRmphKqb1mIbtU7Yj2DaY15J0ENxTJ2rn8svz\n8b/L/8OPF37EqYxTGO8zHtN6TMMQ9yEw0Od+W3vE7heiDqa4shg7Enbgl/hfoExWYmj3oZgSMAWj\nvUfD2NBY6vJITQx1og6guLIYUVej8MvFXxCdFI0BLgMwyX8SnvJ9CuZdzaUujzSIoU6kowoqCrDr\nyi5sid+CP5L/wOMuj2Oi30Q87fs0b6SlwxjqRDoksyQTOxJ2YNulbTiaehSD3Qdjgu8EjPMZxyDv\nIBjqRO1cQk4CIhMisf3ydlzKuYQRniMwwXcCRnqNhGlnU6nLozbGUCdqZ2rqanAs9Rh2XtmJyIRI\nlFSVYKz3WDzt+zQUbgpeFNTBMdSJ2oH88nzsS9yHXVd2Yc+1PXA2c8Z43/EY5z0OPe178jxyUmGo\nE2khQRAQdysOe67uwe5ruxGXGYeBrgMx1nssRnuPhpOZk9QlkpZiqBNpidyyXBy4fgB7E/di77W9\nMO1sihGeIzDaazQGuQ6CkaGR1CVSO8BQJ5JIVW0VjqUew2/Xf8P+xP1IyE3AQNeBGOExAsM9h8PT\nwlPqEqkdYqgTtZE6oQ4Xsi4g+no0DiQdwOEbh+Ft6Y1hHsMw3GM4+jn3Q+dOnaUuk9o5hjpRKxEE\nAdfyruGP5D/we9Lv+D3pd5h3NccQ9yEY6j4UQ9yHwNLYUuoySccw1Ik0RBAEJOYnQpmsxMEbB/FH\n0h/Q09PDYLfBGOw2GEPch8C1m6vUZZKOY6gTtVCdUIdL2ZdwOOUwDt04hEM3DgEAFG4KDHIdBIWb\nAp4WnjzdkNoUQ52omSprKhF7MxZ/pvyJwymHcST1CLp17YaBrgMx0GUgBroORHd5d4Y4SYqhTtSE\nrNIsHEs9hqNpR3Ek5QjOZJ6Bt6U3Hnd5HANcBuBxl8fhIHOQukyiBhjqRACqa6tx7tY5HE87juPp\nx3Es9Rhyy3PR27E3+jr1RX+X/ujt2BuyLjKpSyV6IIY6dTiCIOB6/nWczDiJmPQYnEg/gbjMOLjL\n3dHHqQ/6OPZBH6c+8LP2g76evtTlEj0ShjrpNEEQkF6cjtMZp3Hq5imcTD+JUxmnYGRohF4OvdDb\nsTcec3wMYQ5hMOtiJnW5RGpjqJPOEAQBqUWpiL0Zi9M3T4vTjNOoE+oQ7hCuar0cesFeZi91uUSt\nQrJQT01NxcyZM5GVlQU9PT0sWLAAixYt0khhpPtq6mqQkJOAs5lncfbWWZy5eQZnMs+gc6fOCLUL\nRZhDGMLsw9DTvieczZx5Rgp1GJKFemZmJjIzMxESEoKSkhKEhYVh+/bt8PPzU7sw0i05ZTk4f+s8\nzt06h7hbcYi7FYdL2ZfgbO6MELsQBNsGI9QuFKH2obAztZO6XCJJqZOdBuqs2M7ODnZ24h+gqakp\n/Pz8kJGRoQp16njKq8txKecSzt86j/NZt9ut8yitLkUPmx4Isg1Cb8feWBC2AIE2gRzVh0jD1Ar1\nuyUnJ+PMmTPo3bt3g59HRESo5hUKBRQKhaZWSRIqry7H5ZzLiM+OR3xOPC5mXcSFrAtIL06Hl4UX\netj2QKB1IF7o9QKCbIPgYu7C7hOiJiiVSiiVSo0sSyMHSktKSqBQKPD222/jqaeeql84u1/avbzy\nPFzOuYzLOZdxKecSLmVfQnx2PDKKM+Bl6QV/a3/4W/nD39ofgTaB8LTwhGEnQ6nLJmrXJD37pbq6\nGmPGjMHIkSOxePFijRVGbaeqtgrX86/jSu4VJOQkICH3dstJQEVNBfys/eBr5QtfS1/4WfvBz8oP\n3eXdGd5ErUSyUBcEAbNmzYKlpSWWLVum0cJIs6prq3Gj8Aau5l7FtbxruJp3FVfzruJK7hWkF6XD\nxdwF3pbe8LL0go+lj9isfGBvas9uE6I2Jlmo//nnnxg4cCCCgoJUf/gff/wxRowYoXZh9OhKqkpw\nPf86rudfR2JeIhLzxXYt7xrSitLgIHOAl4UXPC084WXhBS9LL3hZeMFd7s6BHYi0CC8+6iAqayqR\nWpSKpPwkJBUkIbkgGUkFSbiefx1J+UkoqSqBu9wdHnIPdJd3h4fcAx4WHvCQe8Ctmxu6GHSR+lcg\nomZgqOuIkqoSpBSmIKUwBTcKbuBG4e1WcAPJBcnILsuGg8wB7t3c4S53h5u5G9zl7ugu7w73bu6w\nM7VjVwmRDmCotwOVNZXIKM5AalEqUgtTxent+ZTCFKQWpaKsugwu5i5wNXdVTV27ucLV3BXucnc4\nyBxgoK+xs1CJSEsx1CUkCAIKKwuRXpSOjOIMpBenI70oXZwWpyOtKA3pRenIK8+DvcwezmbOcDJz\ngrO5M5zNnOFi7gJnM2c4mzvD2tiae9pExFBvDYIgIK88D5klmbhZchM3i2+K09vzGcUZqtZJvxMc\nZY5wNHOEo8wRDjIHOJk5wcnMCY4yRziZOcHGxAad9DtJ/WsRUTvAUG8mQRCQX5GPrNIs3Cq5hVul\nt+qnt+czSzKRWZKJW6W3YGRgBHuZPexN7WFnagcHmQPsTe1VP3M0c4S9qT0HXSAijeqwoX6n6yOn\nLAfZpdnILstGdmk2skqzxPmy2/Ol2bhVegvZpdkwNjSGjYkNbE1tYWtiCxsTG9iZ2sHWxFacmtrC\n3tQetqa26GrQtdVqJyJqik6Eep1Qh4KKAuSW5SKvPA+55bnIKctBbllu/fzt6Z2WW5aLrgZdYW1i\nDWtja1gZW8HGxEb12MbERnx81zxP6yMibafVoR55ORJ55XnIL89HfkW+OH97encrrCiErIsMlkaW\nsDCygKWxJSyNLGFlbFU/vf0zaxNr1c8Z0kSka7Q61EdvHg0LIwvIjeSQd5XDwshCfHx73tJYDPFu\nXbvxdD0iImh5qGvTgVIiovZAnezkMOtERDqEoU5EpEMY6kREOoShTkSkQxjqREQ6hKFORKRDWv3E\ncN50kIio7bR6qPM0dSKiR6POzjC7X4iIdAhDnYhIhzDUiYh0CEOdiEiHtP5tEdeuBUxMAFPThk0m\nE6cmJoAB785IRKQJrZ+mR48CJSViKy4GSkvr5+/8vEsXMeQf1MzM6qcPavyCIKIOTPpb7woCUFZW\nH/LFxfe3oqL66d3zhYUNf15UBHTtKoa7ufn9rVu3hvONNZmMJ9cTkaR4P/U7BEH8T6CwUGx3gv/e\nlp8vTgsKxHb3fFlZfejL5WK7M29hcf/0TrO0BIyN+YVARGpjqGtSTU198DfW8vLun+bmivM1NQ1D\n/u5mZVU/vXteLgf0ebyaiOox1LVFRUV9yDfWcnIazmdni11OcrkY8NbWYrOyAmxsxPm7pzY24pdB\np05S/6ZE1IokDfW9e/di8eLFqK2txbx58/DGG29opLAOo7pa/CLIzq4P+rtbVpbY7swXFIhfAra2\nYsjb2jZsdnb1UxsbHjgmaockC/Xa2lr4+PjgwIEDcHR0RK9evfDjjz/Cz89P7cKoCTU1YvhnZQG3\nbjXeMjPFlpsrfgHY24shb2/fsDk41E+7dpX6NyOi29TJTrV242JiYuDp6Qk3NzcAwJQpUxAZGakK\ndWoFBgZiQNvZPfy1tbXiHn5mJnDzZn1LSACUSnE+I0N83tRUDHcHB8DRsWFzchKblRUPBBNpObVC\nPT09Hc7OzqrHTk5OOHHiRIPXREREqOYVCgUUCoU6q6RH0alT/RdASEjTr6urE/fqMzKA9PT6aWws\nsHMnkJYmPi4pEUPe2bk+6J2dxebiIja5nMFP9IiUSiWUSqVGlqVWqOs144/37lAnLaWvX3+QNji4\n6deVl4sBn5YGpKaK7fx5YPducT4lReweuhPwrq71U1dXwM1N/E+AB3qJGrh3h/e9995r8bLUCnVH\nR0ekpqYH6MTAAAANcUlEQVSqHqempsLJyUmdRZI2MzICvLzE1pSiIjHcU1KAGzfEFhUlTpOTxf8I\nnJzEgHdzA9zd61v37uJBXu7pE7WYWgdKa2pq4OPjg+joaDg4OOCxxx7jgVJ6sIqK+sBPSqpv16+L\n07Ky+oD38KifeniIP+/cWerfgKjVSXpK4549e1SnNM6dOxd///vfNVIYdVBFRWK4JyaKQX/9OnDt\nmvg4LU3svvH0FP9buDP19mbgk07hxUfUMVRXi104iYnA1ati2F+9KrbUVLFbx9tbbD4+4tTXV/wi\nYJcOtSMMdaKqKnEP/8oVsSUk1LfSUjHkfX0BPz9x6u8v7ukbGkpdOdF9GOpED5KfD1y+LAb8pUv1\nLTVV7LP39xdbQIA49fFhVw5JiqFO1BIVFWLXzcWLQHy82C5eFLt4uncXQz4wEOjRQ2zdu/Pma9Qm\nGOpEmlRZKe7Vnz8PXLggTs+fF0/HDAgQz+UPChJbcLB4q2YiDWKoE7WFwkIx3M+dA+LixOn58/UX\nbYWGilfuhoaKV9ny4Cy1EEOdSCp1deLZOGfPiu3MGbFVVQE9e9a3sDDxXHsGPTUDQ51I29y8KYZ7\nbKzYTp8W9/R79gR69QLCw8Xm5sagp/sw1Inag+xsMdxPnwZOnhRbVZUY8r16AY89BvTuLd4Nkzo0\nhjpRe5WeXh/wJ06IU2trMdz79AH69hX763k+fYfCUCfSFbW14pk3x4/Xt+vXxW6bvn2Bfv3EZm0t\ndaXUihjqRLqsqAiIiQGOHhXb8ePi3Sz79wcef1xsXl7sm9chDHWijqS2VrxI6sgR4M8/gcOHxb75\nxx8HBg4UW1AQL5RqxxjqRB3djRvAoUNiwB88KI5h+/jjgEIBDB4s9stzcJJ2g6FORA1lZoohr1TW\nj0c7aBAwZIjYAgLYXaPFGOpE9GCZmWK4//478Mcf4jnzTzwhtqFDxfPlSWsw1Ino0dy4AURHi+3A\nAcDMDBg2TGxDhvB+NhJjqBNRy9XVifew+e03YP9+4Ngx8R42w4eLLSyMB13bGEOdiDSnvFzsj9+3\nD9i7F8jJEcN95EhxamkpdYU6j6FORK3nxg0x3HfvFvvlAwOB0aPFFhTEA66tgKFORG2jslI8ZTIq\nCti1C6ipAcaMAcaOFU+d7NJF6gp1AkOdiNqeIIjDAu7cKbYLF8QDrePGiXvxFhZSV9huMdSJSHpZ\nWeIefGSkeOrkY48BTz8NPPUU4OgodXXtCkOdiLRLaal4Js3//id20/j4ABMmABMnimO90gMx1IlI\ne1VViRc8bd0KbN8uDvU3aZLYPDykrk4rMdSJqH2oqRFPl/z1V2DbNsDJCXj2WWDyZF7VeheGOhG1\nP7W14pk0P/8s7sV7eQFTp4oBb2cndXWSYqgTUftWXS3eruDHH4EdO8Th/aZPF/vhzcykrq7NMdSJ\nSHeUlYln0WzeLPbFjxgBzJghXs3aQYb1Y6gTkW7KzRX73zduBK5dE7tnZs4EQkN1+kpWSUL9tdde\nw65du9C5c2d4eHhg/fr1ML/nzm4MdSLSmGvXxHD//ntAJgPmzBG7aGxspK5M49TJzhbfeu3JJ5/E\nxYsXERcXB29vb3z88cctXRQR0cN5egLvvQckJgLLlwNnzgDe3uIFTnduWUAtD/Vhw4ZB//btOHv3\n7o20tDSNFUVE1CR9ffE+M99/D6SkAKNGAR98ALi6Am+9BSQlSV2hpAw0sZB169Zh6tSpjT4XERGh\nmlcoFFAoFJpYJRGReGbM/Pliu3gR+O478fYEoaHAggXifWg6d5a6yodSKpVQKpUaWdYD+9SHDRuG\nzMzM+37+0UcfYezYsQCADz/8ELGxsdi6dev9C2efOhG1tYoK8cKm1auBK1eAuXPFgHdxkbqyZpPs\n7JcNGzZgzZo1iI6ORteuXTVaGBGR2i5dAr75Bti0CejXD/jb34Ann9T6kZwkCfW9e/filVdewcGD\nB2FlZaXxwoiINKasTLywaeVKoKQEeP554LnngG7dpK6sUZKEupeXF6qqqmBx+57Jffv2xapVqzRW\nGBGRxgmCOAbrypXAnj3iee8vvAD4+0tdWQO8+IiI6FFlZIj97qtXA8HBwOLF4lWrWtA1w1AnImqp\nykrgp5+AZcvE+ZdeAmbNAoyMJCuJoU5EpC5BEAfWXrYMOHFC7HdfuFCSK1YluaKUiEin6OmJFzXt\n2CHeEvjmTXHEpv/7P/EWBe0EQ52I6F6+vmJf++XLgJUV0KePeJ/32FipK3sohjoRUVNsbcVbECQl\nicE+bpx4K+BDh8TuGi3EPnUiouaqrBQvZPrkEzHw335bPGNGw7cB5oFSIqK2VFsr3uf9gw/Es2Te\nflvci9dQuDPUiYikUFcHbN8OvP+++Pjdd4Hx49UOd4Y6EZGUBAHYuROIiBCDPiJCrXBnqBMRaQNB\nEE+JfPddwMAA+Oc/gZEjHzncGepERNqkrk68/e+77wLm5sCHH4rnwDcTQ52ISBvV1oq3IHjnHcDD\nA/joIyA8/KFv4xWlRETaqFMncXDsy5eBCRPEfvZJk8TBO1oJQ52IqLUZGoq3G7h6FejZUxywY+FC\noJGR5dTFUCciaivGxsDf/w4kJABduwIBAeLB1NJSja2CoU5E1NYsLYF//Qs4dUoccs/bG1i7VuyD\nVxMPlBIRSS0mBnjlFaC4GFi2DHpDhvDsFyKidk0QgK1bgddeg15yMkOdiEgnVFRAz8iIoU5EpCt4\nnjoREQFgqBMR6RSGOhGRDmGoExHpEIY6EZEOYagTEekQhjoRkQ5hqLcRpVIpdQlag9uiHrdFPW4L\nzVA71L/88kvo6+sjLy9PE/XoLH5g63Fb1OO2qMdtoRlqhXpqaip+++03uLq6aqoeIiJSg1qh/vLL\nL+Ozzz7TVC1ERKSmFt/7JTIyEkqlEsuWLYO7uztOnz4NCwuLhgt/xBG0iYhI1NJ7vxg86Mlhw4Yh\ns5Hhlj788EN8/PHH2L9//wML4M28iIjaVov21C9cuIAnnngCxsbGAIC0tDQ4OjoiJiYGNjY2Gi+S\niIiaRyO33m2q+4WIiNqWRs5TZ985EZF20Eior1q1Cv369YOXlxc+/fTTRl+zaNEieHl5ITg4GGfO\nnNHEarXS3r174evr2+S22Lx5M4KDgxEUFIT+/fvj3LlzElTZNh62Le44efIkDAwMsG3btjasrm01\nZ1solUqEhoYiMDAQCoWibQtsQw/bFjk5ORgxYgRCQkIQGBiIDRs2tH2RbeC5556Dra0tevTo0eRr\nWpSbgppqamoEDw8PISkpSaiqqhKCg4OF+Pj4Bq+JiooSRo4cKQiCIBw/flzo3bu3uqvVSs3ZFkeP\nHhUKCgoEQRCEPXv2dOhtced1gwcPFkaPHi1s2bJFgkpbX3O2RX5+vuDv7y+kpqYKgiAI2dnZUpTa\n6pqzLd59911h6dKlgiCI28HCwkKorq6WotxWdejQISE2NlYIDAxs9PmW5qbae+oxMTHw9PSEm5sb\nDA0NMWXKFERGRjZ4zY4dOzBr1iwAQO/evVFQUIBbt26pu2qt05xt0bdvX5ibmwMQt0VaWpoUpba6\n5mwLAPj666/xzDPPwNraWoIq20ZztsUPP/yAiRMnwsnJCQBgZWUlRamtrjnbwt7eHkVFRQCAoqIi\nWFpawsDggSfqtUsDBgyAXC5v8vmW5qbaoZ6eng5nZ2fVYycnJ6Snpz/0NboYZs3ZFndbu3YtRo0a\n1Raltbnmfi4iIyPx/PPPA9DdYzPN2RZXr15FXl4eBg8ejPDwcGzcuLGty2wTzdkW8+fPx8WLF+Hg\n4IDg4GAsX768rcvUCi3NTbW//pr7hyjcc5KNLv4BP8rv9Mcff2DdunU4cuRIK1YkneZsi8WLF+OT\nTz5RDbJ772dEVzRnW1RXVyM2NhbR0dEoKytD37590adPH3h5ebVBhW2nOdvio48+QkhICJRKJRIT\nEzFs2DDExcVBJpO1QYXapSW5qXaoOzo6IjU1VfU4NTVV9S9kU6+5c167rmnOtgCAc+fOYf78+di7\nd+8D//1qz5qzLU6fPo0pU6YAEA+O7dmzB4aGhhg3blyb1tramrMtnJ2dYWVlBSMjIxgZGWHgwIGI\ni4vTuVBvzrY4evQo3nrrLQCAh4cH3N3dkZCQgPDw8DatVWotzk11O/urq6uF7t27C0lJSUJlZeVD\nD5QeO3ZMZw8ONmdb3LhxQ/Dw8BCOHTsmUZVtoznb4m6zZ88Wtm7d2oYVtp3mbItLly4JTzzxhFBT\nUyOUlpYKgYGBwsWLFyWquPU0Z1ssWbJEiIiIEARBEDIzMwVHR0chNzdXinJbXVJSUrMOlD5Kbqq9\np25gYICVK1di+PDhqK2txdy5c+Hn54fVq1cDAP76179i1KhR2L17Nzw9PWFiYoL169eru1qt1Jxt\n8c9//hP5+fmqfmRDQ0PExMRIWXaraM626Ciasy18fX0xYsQIBAUFQV9fH/Pnz4e/v7/ElWtec7bF\nm2++iTlz5iA4OBh1dXX47LPPdPLCxqlTp+LgwYPIycmBs7Mz3nvvPVRXVwNQLzc1ckUpERFpB458\nRESkQxjqREQ6hKFORKRDGOpERDqEoU5EpEMY6kREOuT/AXJLVglb8EqlAAAAAElFTkSuQmCC\n" } ], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }