{ "cells": [ { "cell_type": "markdown", "id": "aaaf8e92", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "d0c85e40", "metadata": {}, "source": [ "\n", "< [3.7 Algorithms Homework 1](https://ndcbe.github.io/CBE60499/03.07-Algorithms1.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [3.9 Algorithms Homework 3](https://ndcbe.github.io/CBE60499/03.09-Algorithms3.html) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[3.8 Algorithms Homework 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8-Algorithms-Homework-2)", "section": "3.8 Algorithms Homework 2" } }, "source": [ "# 3.8 Algorithms Homework 2" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[3.8 Algorithms Homework 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8-Algorithms-Homework-2)", "section": "3.8 Algorithms Homework 2" } }, "source": [ "**Important for Spring 2021**: We are skipping this homework assignment. The Python solutions will be posted to Sakai. Please use this as exam study material.\n", "\n", "Homework solutions are copyright Alex Dowling (2021). You MAY NOT share them with anyone (post online, etc.) without written permission." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "nbpages": { "level": 1, "link": "[3.8 Algorithms Homework 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8-Algorithms-Homework-2)", "section": "3.8 Algorithms Homework 2" } }, "outputs": [], "source": [ "## Load libraries\n", "import matplotlib as plat\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy import linalg\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib import cm\n", "import sympy as sym\n", "# from sympy import symbols, diff" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.8.1 Finite Difference Approximations](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1-Finite-Difference-Approximations)", "section": "3.8.1 Finite Difference Approximations" } }, "source": [ "## 3.8.1 Finite Difference Approximations\n", "\n", "**Main Idea**: Explore the accuracy of the finite difference approximation for $\\nabla f(x)$ and $\\nabla^2 f(x)$ from Example 2.19." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.1 Finite difference order](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.1-Finite-difference-order)", "section": "3.8.1.1 Finite difference order" }, "tags": [ "gradescope" ] }, "source": [ "### 3.8.1.1 Finite difference order" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.1 Finite difference order](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.1-Finite-difference-order)", "section": "3.8.1.1 Finite difference order" } }, "source": [ "Repeat the analysis from class to show the backward and central finite difference truncations errors are $\\mathcal{O}(\\epsilon)$ and $\\mathcal{O}(\\epsilon^2)$, respectively. We discussed these error orders graphically. Please use a Taylor series for your analysis." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.2 Provided Codes](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2-Provided-Codes)", "section": "3.8.1.2 Provided Codes" } }, "source": [ "### 3.8.1.2 Provided Codes" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.2 Provided Codes](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2-Provided-Codes)", "section": "3.8.1.2 Provided Codes" } }, "source": [ "Please review the following code. You do not need to turn anything in for this section." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.1 Finite Difference Code](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.1-Finite-Difference-Code)", "section": "3.8.1.2.1 Finite Difference Code" } }, "source": [ "#### 3.8.1.2.1 Finite Difference Code\n", "\n", "The code below has been adapted from the finite difference examples presented in class. Notice the second input is a function." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.1 Finite Difference Code](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.1-Finite-Difference-Code)", "section": "3.8.1.2.1 Finite Difference Code" } }, "outputs": [], "source": [ "## Define Python function\n", "def my_f(x,verbose=False):\n", " ''' Evaluate function given above at point x\n", "\n", " Inputs:\n", " x - vector with 2 elements\n", " \n", " Outputs:\n", " f - function value (scalar)\n", " '''\n", " # Constants\n", " a = np.array([0.3, 0.6, 0.2])\n", " b = np.array([5, 26, 3])\n", " c = np.array([40, 1, 10])\n", " \n", " # Intermediates. Recall Python indicies start at 0\n", " u = x[0] - 0.8\n", " s = np.sqrt(1-u)\n", " s2 = np.sqrt(1+u)\n", " v = x[1] -(a[0] + a[1]*u**2*s-a[2]*u)\n", " alpha = -b[0] + b[1]*u**2*s2+ b[2]*u \n", " beta = c[0]*v**2*(1-c[1]*v)/(1+c[2]*u**2)\n", " \n", " if verbose:\n", " print(\"##### my_f at x = \",x, \"#####\")\n", " print(\"u = \",u)\n", " print(\"sqrt(1-u) = \",s)\n", " print(\"sqrt(1+u) = \",s2)\n", " print(\"v = \",v)\n", " print(\"alpha = \",alpha)\n", " print(\"beta = \",beta)\n", " print(\"f(x) = \",)\n", " print(\"##### Done. #####\\n\")\n", " \n", " return alpha*np.exp(-beta)\n", "\n", "## Calculate gradient with central finite difference\n", "def my_grad_approx(x,f,eps1,verbose=False):\n", " '''\n", " Calculate gradient of function my_f using central difference formula\n", " \n", " Inputs:\n", " x - point for which to evaluate gradient\n", " f - function to consider\n", " eps1 - perturbation size\n", " \n", " Outputs:\n", " grad - gradient (vector)\n", " '''\n", " \n", " n = len(x)\n", " grad = np.zeros(n)\n", " \n", " if(verbose):\n", " print(\"***** my_grad_approx at x = \",x,\"*****\")\n", " \n", " for i in range(0,n):\n", " \n", " # Create vector of zeros except eps in position i\n", " e = np.zeros(n)\n", " e[i] = eps1\n", " \n", " # Finite difference formula\n", " my_f_plus = f(x + e)\n", " my_f_minus = f(x - e)\n", " \n", " # Diagnostics\n", " if(verbose):\n", " print(\"e[\",i,\"] = \",e)\n", " print(\"f(x + e[\",i,\"]) = \",my_f_plus)\n", " print(\"f(x - e[\",i,\"]) = \",my_f_minus)\n", " \n", " \n", " grad[i] = (my_f_plus - my_f_minus)/(2*eps1)\n", " \n", " if(verbose):\n", " print(\"***** Done. ***** \\n\")\n", " \n", " return grad\n", "\n", "## Calculate gradient using central finite difference and my_hes_approx\n", "def my_hes_approx(x,grad,eps2):\n", " '''\n", " Calculate gradient of function my_f using central difference formula and my_grad\n", " \n", " Inputs:\n", " x - point for which to evaluate gradient\n", " grad - function to calculate the gradient\n", " eps2 - perturbation size (for Hessian NOT gradient approximation)\n", " \n", " Outputs:\n", " H - Hessian (matrix)\n", " '''\n", " \n", " n = len(x)\n", " H = np.zeros([n,n])\n", " \n", " for i in range(0,n):\n", " # Create vector of zeros except eps in position i\n", " e = np.zeros(n)\n", " e[i] = eps2\n", " \n", " # Evaluate gradient twice\n", " grad_plus = grad(x + e)\n", " grad_minus = grad(x - e)\n", " \n", " # Notice we are building the Hessian by column (or row)\n", " H[:,i] = (grad_plus - grad_minus)/(2*eps2)\n", "\n", " return H" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.1 Finite Difference Code](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.1-Finite-Difference-Code)", "section": "3.8.1.2.1 Finite Difference Code" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "xt = [0 0] \n", "\n", "f(xt) = \n", " 1.6212212164426274e-06 \n", "\n", "grad(xt) = [1.49885593e-04 4.20936251e-05] \n", "\n", "hes(xt) = \n", " [[0.00082024 0.00387044]\n", " [0.00387044 0.00102412]] \n", "\n" ] } ], "source": [ "### Test the functions from above\n", "\n", "## Define test point\n", "xt = np.array([0,0])\n", "print(\"xt = \",xt,\"\\n\")\n", "\n", "print(\"f(xt) = \\n\",my_f([0,0]),\"\\n\")\n", "\n", "## Compute gradient\n", "g = my_grad_approx(xt,my_f,1E-6)\n", "print(\"grad(xt) = \",g,\"\\n\")\n", "\n", "## Compute Hessian\n", "# Step 1: Create a Lambda (anonymous) function\n", "calc_grad = lambda x : my_grad_approx(x,my_f,1E-6)\n", "\n", "# Step 2: Calculate Hessian approximation\n", "H = my_hes_approx(xt,calc_grad,1E-6)\n", "print(\"hes(xt) = \\n\",H,\"\\n\")\n" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.2 Analytic Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.2-Analytic-Gradient)", "section": "3.8.1.2.2 Analytic Gradient" } }, "source": [ "#### 3.8.1.2.2 Analytic Gradient\n", "\n", "It turns out that calculating the analytic gradient with Mathematic quickly becomes a mess. Instead, the following code uses the symbolic computing capabilities in SymPy." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.2 Analytic Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.2-Analytic-Gradient)", "section": "3.8.1.2.2 Analytic Gradient" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The exact gradient is \n", " [1.49885593e-04 4.20936251e-05]\n" ] } ], "source": [ "'''\n", "Encoding the exact gradient with the long expression above is very time-consuming. This is a trick of calculating the \n", "symbolic derivative and converting it to an analytic function to be evaluated at a point. \n", "'''\n", "\n", "# Define function to use with symbolic computing framework\n", "def f(x1,x2):\n", " a = np.array([0.3, 0.6, 0.2])\n", " b = np.array([5, 26, 3])\n", " b1 = 5;\n", " c = np.array([40, 1, 10])\n", " u = x1 - 0.8\n", " v = x2 -(a[0] + a[1]*u**2*(1-u)**0.5-a[2]*u)\n", " alpha = b[1]*u**2*(1+u)**0.5 + b[2]*u -b[0]\n", " beta = c[0]*v**2*(1-c[1]*v)/(1+c[2]*u**2)\n", " return alpha*sym.exp(-1*beta)\n", "\n", "# Define function to use later\n", "def my_grad_exact(x):\n", " x1, x2 = sym.symbols('x1 x2')\n", " DerivativeOfF1 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x1));\n", " DerivativeOfF2 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x2));\n", " #DerivativeOfF2 = sym.lambdify((x1,x2),gradf2(x1,x2));\n", " #F = sym.lambdify((x1,x2),f(x1,x2));\n", " return np.array([DerivativeOfF1(x[0],x[1]),DerivativeOfF2(x[0],x[1])])\n", " \n", "x = np.array([0,0])\n", "print(\"The exact gradient is \\n\",my_grad_exact(x))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.3 Analytic Hessian](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.3-Analytic-Hessian)", "section": "3.8.1.2.3 Analytic Hessian" } }, "source": [ "#### 3.8.1.2.3 Analytic Hessian\n", "\n", "The code below assembles the analytic Hessian using the symbolic computing framework in NumPy." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "nbpages": { "level": 4, "link": "[3.8.1.2.3 Analytic Hessian](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.2.3-Analytic-Hessian)", "section": "3.8.1.2.3 Analytic Hessian" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The exact Hessian is \n", " [[0.00082025 0.00387044]\n", " [0.00387044 0.00102412]]\n" ] } ], "source": [ "def f(x1,x2):\n", " \n", " a = np.array([0.3, 0.6, 0.2])\n", " b = np.array([5, 26, 3])\n", " b1 = 5;\n", " c = np.array([40, 1, 10])\n", " u = x1 - 0.8\n", " v = x2 -(a[0] + a[1]*u**2*(1-u)**0.5-a[2]*u)\n", " alpha = b[1]*u**2*(1+u)**0.5 + b[2]*u -b[0]\n", " beta = c[0]*v**2*(1-c[1]*v)/(1+c[2]*u**2)\n", " return alpha*sym.exp(-1*beta)\n", "\n", "\n", "def my_hes_exact(x):\n", " x1, x2 = sym.symbols('x1 x2')\n", " HessianOfF11 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x1,x1));\n", " HessianOfF12 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x1,x2));\n", " HessianOfF21 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x2,x1));\n", " HessianOfF22 = sym.lambdify((x1,x2),sym.diff(f(x1,x2),x2,x2));\n", " #DerivativeOfF2 = sym.lambdify((x1,x2),gradf2(x1,x2));\n", " #F = sym.lambdify((x1,x2),f(x1,x2));\n", " return np.array([[HessianOfF11(x[0],x[1]),HessianOfF12(x[0],x[1])],[HessianOfF21(x[0],x[1]),HessianOfF22(x[0],x[1])]])\n", " \n", "x = np.array([0,0])\n", "print(\"The exact Hessian is \\n\",my_hes_exact(x)) " ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.3 Gradient Finite Difference Comparison](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.3-Gradient-Finite-Difference-Comparison)", "section": "3.8.1.3 Gradient Finite Difference Comparison" } }, "source": [ "### 3.8.1.3 Gradient Finite Difference Comparison\n", "\n", "Repeat the analysis procedure from the finite difference class notebook to find the value of $\\epsilon_1$ that gives the smallest approximation error. Some tips:\n", "1. Write a `for` loop to iterate over many values of $\\epsilon_1$\n", "2. Use $|| \\nabla f(x;\\epsilon_1)_{approx} - \\nabla f(x)_{exact} ||$ to measure the error. Your choice on which norm(s) to use. Please label your plot with the norm(s) you used.\n", "3. Make a log-log plot" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.3 Gradient Finite Difference Comparison](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.3-Gradient-Finite-Difference-Comparison)", "section": "3.8.1.3 Gradient Finite Difference Comparison" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEPCAYAAACQmrmQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VfX9x/HXJyEQIOwRkLCXTEECiBOcOFrrFlEBERx19Fdbq7X+bGuttrZ1LwQZtoJKXUWqtWoElBGWDNmgEGRHRoBAxuf3R4I/DAkkl3vvuUnez8eDh7kn9577ztcLn5zvOZ/zNXdHREQkEuKCDiAiIhWXioyIiESMioyIiESMioyIiESMioyIiESMioyIiESMioyIiESMioyIiESMioyIiERMlaADBKVhw4beqlWroGMct71791KzZs2gY8QUjUnxNC5H0pgUr6RxmTdv3nZ3b1SWfVXaItOqVSvmzp0bdIzjlpaWRv/+/YOOEVM0JsXTuBxJY1K8ksbFzL4p6740XSYiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiUkkczM1nzbasqL6nioyISCXx0HtL+fEzM9i6Oztq76kiIyJSCfx91jdMnLOeG09tRePaiVF7XxUZEZEKbs66TH773lL6d2zEL87vGNX3VpEREanAvt25n9v/MY/m9Wvw1LU9iY+zqL5/pb13mYhIRZedk8ctr84jOyefSSN7Uad6QtQzqMiIiFRA7s79by1m8cZdvHxjKu0a1wokh6bLREQqoDEz1vH2go38/LwOnNc5ObAcKjIiIhXMjFXb+ePUZVzQJZk7BrQLNIuKjIhIBbJ+xz7umDifdo2T+OvVPYiL8on+olRkREQqiL0Hchn56lzy852Xb0wlqVrwp92DTyAiIsfN3fnl5C9ZuWUP44b1oWWD2FhWWkcyIiIVwJP/XcXUxZu578ITObNDo6DjfE9FRkSknHtnwUae+ngVV/VKYcQZbYKO8wMqMiIi5Vj615ncO3kRfVvX55HLumEW7In+oipEkTGz/mY23cxeNLP+QecREYmGb3bs5ZZX59GsXnVeuqEXVavE3j/pgScys1fMbKuZLSmyfaCZrTCz1WZ23zF240AWkAhkRCqriEis2LU/h5vGpZPvzitDe1O3RtWgIxUrFq4uGwc8C0w4tMHM4oHngPMoKBrpZvYeEA88WuT1NwHT3f0zM0sG/gYMjkJuEZFA5OTlc/s/5rE+cx+vDu9L64axcSVZcczdg86AmbUCprh718LH/YDfuvsFhY/vB3D3ogWm6H6qAq+5+5UlfH8kMBIgOTm516RJk8L1IwQmKyuLpKSkoGPEFI1J8TQuRyqPY+LujF16kGkZudzcrSqnNwv/TS9LGpcBAwbMc/fUsuwrFo5kitMM2HDY4wygb0lPNrPLgQuAuhQcFRXL3UcBowBSU1O9f//+4cgaqLS0NCrCzxFOGpPiaVyOVB7HZNS0NUzLWM5PB7TllxecGJH3COe4xGqRKe7yiBIPudz9LeCtyMUREQneh0s38+i/l3Nxt6bcc150Fx8LVeAn/kuQATQ/7HEK8G1AWUREArc4Yxc/m7SQ7il1+evVJwV+T7LSitUikw60N7PWhedZrgXeCziTiEggvt25n+Hj06lfsyov39iLxIT4oCOVWuBFxswmAjOBjmaWYWbD3T0XuAP4EFgGvOHuS4PMKSIShD3ZBZcq7z+Yx5ihqTSulRh0pDIJ/JyMuw8qYftUYGqU44iIxIyCS5Xns3prFuOG9eHEJrWDjlRmgRcZERE5krvz4DtLmL5qO3++ojunt28YdKSQBD5dJiIiR3o+bQ2T0jdw59ntuLp382O/IEapyIiIxJh3F27k8Q9XcGmPE/j5eR2CjnNcVGRERGJI+teZ/PLNRfRpVZ8/X9k95u6qXFYqMiIiMWLttixGTJhLSv3qjLqxF9WqlJ9LlUuiIiMiEgN2ZB1g2Lh04s0YN7RPzN5Vuax0dZmISMCyc/IYMWEum3dlM3HkKbRoUCPoSGGjIiMiEqD8fOeeN75kwYadPH/dyZzcol7QkcJK02UiIgF67IPlvL94E7++sBMXdmsadJywK3ORMbPzzOxlM+tR+Hhk+GOJiFR8E2Z+zahpaxnSryU3n9E66DgREcp02e3AMOA3ZlYf6BHeSCIiFd9/v9rCb99byrmdkvnfH3Up95cqlySU6bJt7r7T3X8BnA/0DnMmEZEKbVHGTu6cuICuzerw9KAexJeT2/aHIpQi8/6hL9z9PmBC+OKIiFRsGzL3cdO4uTRIqsqYIb2pUbViX391zJ/OzFoU2bSgyLZ3Cx/vdPfdYU0nIlKB7NqXw7Bx6RzMzWPSyL40qlUt6EgRV5oSOp6CpY+PdjznwDh0VCMiUqwDuXmMfHUu63fsY8LwPrRrXCvoSFFxzCLj7gOiEUREpKLKz3funbyI2esyeeraHpzSpkHQkaJGfTIiIhH2149W8O7Cb/nlBR25tEezoONEVSjnZEqiczIiIkVMnLOe5z5dw6A+zbm9f9ug40SdzsmIiETIpyu28pt3lnBWh0b8/tKuFbYX5mh0TkZEJAKWbNzFHf+YT8fkWjw3+GQS4ivn2YmQfmoz+/lhX3cMXxwRkfJv48793DQunTrVExg7rDdJ1Sp2L8zRlOknN7O6wBNARzPLBhYBwym4zYyISKW3a38Ow8bOYf/BPCbfdirJtRODjhSoMh3JFN5OZhgwEZgNtAc+jkQwEZHy5mBuPrf9fR5rt+3lxRt60bFJ5eiFOZpQJwlbuPs8dx8LdAtnIBGR8sjdue+tRXyxZgd/uqI7p7VrGHSkmBBqkYkzszPMLA6oPF1FIiIleOK/q3hr/kb+59wOXNErJeg4MSPUInMvcBIwGngvfHFERMqfN+Zu4OmPV3FVrxTuOqdd0HFiSkiXPLh7HvBsmLOIiJQ701Zu49dvLeaM9g354+XdKmUvzNGEVGTM7LdAX2AjsMDdnwtnKBGR8mDZpt3c/o/5tGucxPOVuBfmaEIdkbrALOARQH0yIlLpbNq1n2Fj00mqVoWxw3pTKzEh6EgxKdQikwnEA1sLvxYRqTT2ZOcwbGw6WQdyGTusN03rVA86UswKtQ31ESAZeBpYEr44oTGzM4DBFPw8nd391IAjiUgFlZOXz+3/mM/qrVmMHdabTk1rBx0ppoV6JPMSBXddHg7MO54AZvaKmW01syVFtg80sxVmttrM7jvaPtx9urvfCkyh4IaeIiJh5+78+q3FTF+1nT9e3o0z2jcKOlLMC/VI5iFgjJnlAguBaceRYRwFV6p9fwdnM4sHngPOAzKAdDN7j4IpukeLvP4md99a+PV1wM3HkUVEpERPf7yaN+dlcPc57bk6tXnQccoFc/eyv8jsFWA90AZ4wN03HFcIs1bAFHfvWvi4H/Bbd7+g8PH9AO5etMAcvo8WwIPuPuIozxkJjARITk7uNWnSpOOJHROysrJISkoKOkZM0ZgUT+NypLKMyYyNOYxefJDTTqjCzd2qVuhLlUsalwEDBsxz99Sy7CvUI5l73X27mdUEniL8Rw/NgMMLVwYFl0wfzXBg7NGe4O6jgFEAqamp3r9//+OIGBvS0tKoCD9HOGlMiqdxOVJpx+Tz1dsZ9585nNauAWOH9qFqlYp9qXI4PyuhFpnfmVlNdx9qZm+HJckPFfcrwlEPudz9oQjkEJFKbvnm3dz66jzaNkrihet7VfgCE26hjtZBYG3h12eEKcvhMoDDJzxTgG8j8D4iIiXavCubYWPTqVEtnrHDelNbvTBlFmqR2QfUMbMEoEUY8xySDrQ3s9ZmVhW4Ft0jTUSiKOtALjeNS2dPdi5jh/bhhLrqhQlFqEXmIWANBVeAvXY8AcxsIjCTgoXQMsxsuLvnAncAHwLLgDfcfenxvI+ISGkd6oVZsWUPzw0+mc4nqBcmVGU+J1N4sj/b3Z8PRwB3H1TC9qnA1HC8h4hIabk7D76zhGkrt/GnK7pxVgf1whyPYx7JmFmcmV1nZu+b2VZgObDZzJaa2eNm1j7yMUVEouO5T1czKX0Dd53djmt6R+JsQOVSmumyT4G2wP1AE3dv7u6NKDjhPwt4zMyuj2BGEZGoeHtBBn/5z0ou79mM/zmvQ9BxKoTSTJed6+45RTe6eybwT+CfhRcAiIiUW1+s3s69kxfRr00DHruie4VutoymYx7JHCowZvaklTDqxRUhEZHyYuWWPdzy93m0bliTF29QL0w4lWUks4D3Ck/8Y2bnm9nnkYklIhIdW3ZnM/SVOVRPiGfssD7Uqa6JmXAq9dVl7v4bM7sOSDOzA8Be4Kh3RxYRiWX7c52bxqWza38Or9/Sj2bqhQm7UhcZMzsHGEFBcWkKDHf3FZEKJiISSbl5+Ty/8ADLM/czekgqXZvVCTpShVSW6bIHKLjLcX/gSuB1Mzs7IqlERCLI3Xnw3SUs3p7HH37SlQEdGwcdqcIqy3TZ2Yd9vdjMLqTg6jKtQiki5crzaWuYOGcDl7RJYFAf9cJEUmmaMUu6omwTcM7RniMiEmveWbCRxz9cwWU9m3FFe53kj7RSNWOa2Z2Fi4J9r/DGlf3MbDwwJCLpRETC6Is12/nl5C85pU19/qRemKgozXTZQOAmYKKZtQZ2AokULIX8H+AJd18YuYgiIsdv5ZY93PLqPFo1qMlL16eqFyZKSlNk/uTud5vZOCAHaAjsd/edEU0mIhImW3cXrAuTmFCwLkydGpomi5bSlPJzCv873d1z3H2TCoyIlBd7D+QybFw63+07yNihvUmpVyPoSJVKaYrMB2Y2E2hiZjeZWS8zS4x0MBGR45Wbl89PX5vP8s0F68KoFyb6jjld5u6/MLM2QBrQGvgx0MXMDgJL3P2ayEYUESm7Q70waSu28cfLuqkXJiCl7ZMZBpzj7qsObTCzJKBrRFKJiBynQ70wt/dvy3V91QsTlNJeXnE38L6ZfWhmV5pZFXfPcvdZkQwnIhKKdxcW9MJc2uMEfnF+x6DjVGqlLTJb3L0D8BhwKbCmcFVM/d8TkZgya+0OfvnmIvq2rs+fr+xOXJx6YYJU2iLjAO7+qbvfAHQD9gFLzKxzpMKJiJTFqi17GDlhLi0a1GDUDalUqxIfdKRKr9T3LgMwszjgEmA40BF4EFgXgVwiImWydXc2Q8emU7VKPGOHqhcmVpS2yCSZ2WPAtcAXwFPu/knkYomIlN7eA7ncND6dzL0HeeOWfjSvr16YWFHaIpMJbAJOdvfMCOYRESmT3Lx87py4gK++3c3oIal0S1EvTCwpVZFxd12qLCIxx9156L2lfLJ8K49c1pWzT0wOOpIUoTvEiUi59eJna/nH7PXcelZbBvdtGXQcKYaKjIiUS+8u3MifPljOj046gXsvUDdFrCpzkTGzH0UiiIhIac0u7IXp07o+f7lKvTCxLJQjmUfCnkJEpJRWb93DiAlzSalfnVE39FIvTIwLpcjoVwYRCcS2PQcKe2HiGD+sD3VrVA06khxDmZoxC3nYU4iIHMO+g7kMH5/OjqyDvH7LKeqFKScqxIl/M+tsZm+Y2QtmdmXQeUQkvHLz8rnztQUs2biLZwb1pHtK3aAjSSkFXmTM7BUz22pmS4psH2hmK8xstZndd4zdXAg84+63ATdGLKyIRJ2789t/LeXj5Vv53Y+7cG5n9cKUJ6FMl20Jc4ZxwLPAhEMbzCweeA44D8gA0s3sPSAeeLTI628CXgUeMrMfAw3CnE9EAjRq2lr+Pms9t5zVhhv6tQo6jpSRuQd/isXMWgFTDt1ZwMz6Ab919wsKH98P4O5FC0zR/cQDb7n7pSV8fyQwEiA5ObnXpEmTwvUjBCYrK4ukpKSgY8QUjUnxyuO4zNmUy/NfHqBPk3huPakacRbe647K45hEQ0njMmDAgHnunlqWfYVyJBMNzYANhz3OAPqW9OTCIvVroCbweEnPc/dRwCiA1NRU79+///EnDVhaWhoV4ecIJ41J8crbuMxZl8noj2bTu1U9JgzvS2JC+C9VLm9jEi3hHJeQioyZxbt7XlgSlPAWxWwr8ZDL3b+m8AhFRMq/Nduyvu+FefnG1IgUGImOUE/8jzKzGgBmdmYY8xySATQ/7HEK8G0E3kdEYsz2rAMMHTuHhHhj3FD1wpR3oU6X/S8wxsxygYXAtPBFAiAdaG9mrYGNFKxjc12Y30NEYsy+g7kMH5fOtj0HeH1kP1o0UC9MeRfqkczDwAoKprDeOJ4AZjYRmAl0NLMMMxvu7rnAHcCHwDLgDXdfejzvIyKxLS/fuWviQhZv3MUzg07mpObqhakIQj2Sudfdt5tZTeAp4OZQA7j7oBK2TwWmhrpfESk/3J3f/Wsp/122hd9f2oXz1AtTYYRaZO4ys07AXo5yNZeISGmMnr6OCTO/YeSZbbhRvTAVSqjTZfXd/SoKrui6PYx5RKSSeX/RJh6ZuoyLuzXlvoEnBh1HwizUInPAzHpScE6mZhjziEglMvfrTP7njYWktqzHX68+SevCVEChFpnHgXMpaGx8PXxxRKSyWLsti5snzCWlrnphKrJQz8nc6O5/DmsSEak0Cnph0ok3Y+yw3tSrqV6YiirUInOpme0DPnL3FeEMJCIV2/6DeQwfP5ete7KZOOIUWjbQjHtFFup02eXAauAyMxsdxjwiUoHl5Tt3TVrAooydPH1tT3q2qBd0JImwkI5k3H2LmaW7+wfhDiQiFZO78/CUr/joqy089KPOnN+lSdCRJApCvUHmP4GtZlYbGO3un4Y3lohUNGNmrGPcF19z8+mtGXZa66DjSJSEOl223N1vc/fBgJY7FpGjmrq4oBfmwq5N+PVFnYKOI1EU6on/gWaWCXwJ5IYxj4hUMPO+yeRnry+kZ/O6PHFND/XCVDKhHskMBLYBpwKNzGx8+CKJSEWxbvtebh4/lxPqJDJ6SG/1wlRCoR7JvEBBkakDvKxzMiJS1I7CdWHMjHHD+lBfvTCVUqhHMisKz8lch87JiEgRh3phNu/KZvSQVFo1VC9MZaVzMhHg7php3lkqp7x852evL+DLjJ28MLgXJ6sXplI7nnMyy9A5mSNs3LmfM/78KW/Nzwg6ikggHnl/GR8u3cKDF3dmYFf1wlR2oRaZnwJ3Ai2Bz919SPgilW/vLNhIxnf7+dU/F5H+dWbQcUSi6pUZ63jl83XcdFprbjpdvTASepGpC8wC/gB0DF+c8m/Kok10blqblHo1uOXVeWzI3Bd0JJGo+GDJJh5+/ysu6JLMAxerF0YKhFpkMoF4YGvh1wKs2ZbFsk27ubJXCqOHpJKbl8/N4+eSdaDsp632HcwlOycvAilFwm/eN99x96SF9Ghel6eu7Um8emGkUEhFxt1/D7wIPA3sCmuicmzKl5swg4u7N6VtoySeG3wyq7dlcffEBeTle6n388Wa7Zz1eBoXPDlNR0IS877evpcRE+bSpE4io7UujBRR5iJjZjXNLN7dv3X34e7+RCSClUfvL/6W3q3qk1w7EYAz2jfioR915uPlW/nzB8uP+fr8fOeZj1dx/ejZ1KpWhZ37crjihS9Yvnl3pKOLhCRz70GGjp2DuzNuWB8aJFULOpLEmGMWGTOLM7PrzOx9M9sKLAc2m9lSM3vczNpHPmbsW7llDyu3ZHFJ96Y/2H5jv1bccEpLXpq2ljfnbijx9TuyDjBk7Bz++tFKfnTSCfzrztN589Z+xJlx9YszmauLCCTGZOfkcfP4dDbtymb0kN60Vi+MFKM0RzKfAm2B+4Em7t7c3RsBZ1Bw8v8xM7s+ghnLhSlffkucwYVdmx7xvf/9UWdOa9eAX7+9uNgrzuasy+Sip6cze10mf7ysG09e04Oa1arQIbkWk2/rR8Okalw/ZjafLt8ajR9F5Jjy8p2fTVrIgg07efKaHvRqqV4YKV5pisy57v6wuy9y9/xDG909093/6e5XAK9HLmLsc3emLNrEKW0a0KjWkdMFCfFxPH9dL5oXueIsP995IW0Ng16eRfWEeN6+/VSu69viB42cKfVq8Mat/WjXOIkRE+byzoKNUfu5REryx6nL+GDpZn5zcWcu7HbkL1YihxyzyLh7TjieU5F9tWk3a7fv5ZLuJ5T4nDo1En5wxdmGzH0MH5/Onz5YzsAuTfjXnafT5YQ6xb62YVI1Jo44hd6t6vOz1xcy9vN1R82Tm5dP+teZPPrvZfz0tfkhXd0mUpJXZqxjzIx1DDutFcPVCyPHEOptZb5nZr9y9z+FI0x5NWXRJuLj7JjdzW0aJfH84F4MGTuH/n9JI96M31/ahRtOaXnM29DUSkxg7LDe/GzSQn73r6/I3HuQn5/X4fvv79qfw7SV2/h42RbSVm5j574cqsQZuflO20ZJP3iuSKg+WLL5+16Y31zcOeg4Ug6UuciY2RuHPwR6AJW2yLg77y/axKltG5TqLrOnt2/Io5d14x+zv+Hhn3Sle0rdUr9XYkI8zw0+mQfeXswzn6xme9YBfFcOL62cRfrXmeTmO/VrVuXsExtzbqdkTm/fkPvfWszL09Zyfd8WNC686k0kFPPXf8fdkxbQo3ldnrxGvTBSOqEcyex295sPPTCzF8KYp9xZvHEX6zP3cceAdqV+zdW9m3N17+YhvV98nPHo5d2oV7MqL6StAeDEJgcZeWYbzunUmB7N6/3gL/+9F3TkP0s388R/V/Ho5d1Cek+Rb3YUrAtzqBemelX1wkjphFJkHiny+IFwBCmvpizaREK8cUGX6N0I0Mz41cATubBrE1Ysms9VF51Z4nNbNqjJ4L4tmTDza4af3op2jWtFLadUDAW9MOnqhZGQlLoZ08xmArj7usLHtcysp7tHvYHDzNqY2Rgzm3y0bZF2aKrsjPaNqFMjIVpv+73uKXVpVOPY/wvvOqc9NatW4bF/H7shVORw2Tl5jJgwl4079zN6SKp6YaTMytLxXw3AzP4G4O57gOfL+oZm9oqZbTWzJUW2DzSzFWa22szuO9o+3H2tuw8/1rZIm79+Jxt37ufiGL+Es37Nqtw2oC3/XbaVWWt3BB1Hyon8fOd/Xl/I/PXf8dQ1PejVsn7QkaQcKkuRMTNrDFxv/38pVPUQ3nMcBevRHL7jeOA54EKgMzDIzDqbWTczm1LkT+MQ3jMi3l+0iarxcZzXJTnoKMd002mtaVonkUenLsO99PdRk8rrj1OX8e8lm3ngok7qhZGQlaXI3A/MAF4DnjCz28v4egDcfRpH3rm5D7C68GjkIDAJuNTdF7v7JUX+xETbe36+M3XxJs7q2IjaidGfKiurxIR4fn5eB77M2MX7izcFHUdi3LjP1zF6xjqGnqpeGDk+pT7x7+4fAB0AzKwfcBUQrumpZsDhN/bKAPqW9GQza0DBBQg9zex+d3+0uG3FvG4kMBIgOTmZtLS0kAOvyMxj8+5sflLFj2s/xysrK6vU79/AnZQk43dvLyRx+wqqVNBLUMsyJpVJacdl/pZcnllwgJ6N4zmz1lY++2xb5MMFRJ+V4oVzXEJqxnT3mcDMsCQoUNy/diXO6bj7DuDWY20r5nWjgFEAqamp3r9//zIHPeSTd5dQrcoG7ryiPzWrHXdPa8jS0tIoy88R32wbQ16Zw/qqrSrsyoVlHZPKojTjsmD9d4z6eBbdm9fltRGnVPhLlfVZKV44xyXURcvCLQM4vHEkBfg2oCzHlJfvTF28mbNPbBxogQnFme0bcnq7hjzzySp27a/UdwOSIg71wjSulciYIeqFkfCIlSKTDrQ3s9ZmVhW4Fngv4Ewlmr12B9uzDhz1XmWxysy478IT+W5fDi9+tiboOBIjvivshclzZ9yw3jRUL4yESVn6ZJoUedzUzMr8STSziRRMtXU0swwzG+7uucAdwIfAMuANd19a1n1Hy5TFm6hRNZ6zT4yZC93KpGuzOlzWsxmvzFjHtzv3Bx1HApadk8fNhb0wL9+YSptGSUFHkgqkLEcyY4o8fhVYbmZ/Kcsbuvsgd2/q7gnunuLuYwq3T3X3Du7e1t2L3lUgZuTm5fPBks2c0ym5XE8n3HN+Bxz420crg44iAcrPd37+RkEvzBNX96B3K/XCSHiVusi4+8VFHp8LtAHGhjtULPtizQ4y9x48YgXM8ialXg2GndqKf87P4KtvtbxzZfXov5cxdXFBL8zF5fwzLbGpLNNlTx7WhAmAF4jZaa1ImLLoW5KqVeGsDo2CjnLcbu/fjtqJCfzmncXs1Zozlc64z9fx8nT1wkhklWW6LAt4z8xqApjZ+Wb2eWRixab8fOfTFds4r3MyiQnld6rskDo1EvjDT7qycMNOrh8zm137dLVZZfHh0s38bspXnNc5mQcv6XzM9YxEQlWWZszfmNl1QJqZHQD2Ake9x1hFExdnfHzPWRXqt/4fnXQCCfFx3DVxAdeMmsmE4X1oXEvrzlRk89d/x10TF3BSSl2evlbrwkhklWW67BxgBAXFpRFwl7tPj1SwWFU7MYGmdUK5ZVvsGti1Ca8M7c36zH1c/eJMMr7bF3QkiZDD14VRL4xEQ1mmyx4AHnT3/sCVwOtmdnZEUknUnd6+Ia8O70vm3oNc9eJMVm/NCjqShNmeg/79ujBjh/bWujASFWW5uuxsd59R+PViCu6Y/IdIBZPo69WyHq/f0o+cPOfql2ayZOOuoCNJmGTn5PHU/Ozv14VRL4xEyzGLTNEryg5x903AOUd7jpQ/nZrWZvKt/aieEM+gUbOYsy7qa9JJmB1aF2bNznytCyNRV5ojmU/M7E4za3H4xsLbv/Qzs/HAkIikk0C0aliTybf1o3HtatwwZjafroiJ1RUkRIfWhbn2xKpaF0airjRXl60C8oC3zawpsBNIBOKB/wBPuPvCyEWUIDStU503bunHkLFzGDF+Lhd0aUJKveqk1KtOs3rVSalXg2Z1q5e7G4RWNoevC3NWLf2yINFXmn8hTnX3kWZ2M9CCgivL9rv7zshGk6A1SKrGayNO4X/fWcLCDTv56KstHMzL/8Fz6tZIIKVedZrXq0G7xkl0SK5Fxya1aN2wJgnxsXL/1crpP4W9MOcX9sJMn1Zx14WR2FWaIvOhmc0EkoEbgS+BStXlX5nVTkzgyWvg1SplAAAPUUlEQVR7AgVz+9uyDpDx3X427txPxnf72Fj49YrNe/hw6WbyC1cBSog32jRMokOTWnRoXPDfni3qqgcnShZu2MldkxbQPaUuT6kXRgJ0zCLj7veYWRsgDWgN/BjoYmYHgSXufk1kI0qsiIszkmsnklw7kV4t6x3x/eycPNZu28vKLXtYsWUPq7bsYeGG7/jXlwVLAyUmxPHgJZ25rk8LdZhH0IbMfdw8Pp1GtaqpF0YCV6oJdXdfa2bnuvv3t+w1sySga8SSSbmTmBBP5xNq0/mE2j/YvvdALiu27OGJj1bywNtLSFuxjT9d0Z36NasGlLTi2rnvIEPGziEnz3l9WB+tCyOBK9WkuZk9DBw4fJu7Z7n7rIikkgqlZrUqnNyiHuOH9eE3F3fisxXbGPjkNGas2h50tArlQG4eI1+dR0ZmwbowbdULIzGgtGdm7wY+MrMPzewqM9MlRVJmcXHGzWe04e2fnkrt6glcP2Y2j7z/FQdy84KOVu7l5zu/eHMRc9Zl8perT6JPa/XCSGwobZHZ4u4dgMcoOCez2sweN7OOkYsmFVWXE+rwrztO5/pTWvDy9HVc/vwXuo3NcXr8Pyv415ffcu/Ajvz4pPK3LLhUXKU9InEAd/8U+NTMagP3AEvM7CR3/ypSAaViql41nj/8pBtndWjMr/65iEuemc6Dl3Tmku4nsGV3Npt2ZbN513427zrA5t37Cx9ns3t/Dk8P6kmqVnD83muz1/NC2hoG9WnBbWe1DTqOyA+UadrLzOKAS4DhQEfgQWBdBHJJJXFe52ROSjmDe978kgfeXsIDby854jkNk6rRtE4izevXYHHGLu5580s+uPtMXTUFTF+1jQffXcJZHRrx8KVddNWexJzSFpkkM3sMuBb4AnjK3T+JXCypTBrXTmT8sD68s3Aj27MO0KROdZrWSaRJ4eXSVav8/6zuzDU7GPTyLP7ynxU8eEnnAFMHb+WWPdz+9/m0b5zEs9f1pIqaXyUGlbbIZAKbgJPdXXdMlLCLizMuPznlmM/r17YBN5zSklc+X8eFXZtU2mmzbXsOMGxsOolV4xkztDe1EhOCjiRSrNLchbkFcBHwNgVHNC1K+FP7GLsSCYv7LjyRZnWr88vJi8jOqXxXpmXn5DHy1bns2HuAMUNSaVa3Yi2iJxVLaY5kxlNw4v9ok70OjAMmhCGTyFHVrFaFP1/RnetGz+av/1nBAxdXnmmz/Hznnje/ZOGGnbwwuBfdU+oGHUnkqEpzW5kB0QgiUhantmvI4L4tGD1jHQO7Ni32NjcV0d8+Wsn7izZx/4UnMrBrk6DjiBxTqabLSvlH02USVfdf1IkT6lTnl5O/rBTTZpPnZfDsp6u5tndzRp7ZJug4IqVS2umyY9F0mURdUrUqPHZFN24YM4cn/ruS+y/sFHSkiJm5Zgf3v7WI09o14OGfdNWlylJuaLpMyrUz2jdiUJ/mvDxtLQO7NKFni2NPm2Xn5BFn9oNLo2PZ2m1Z3Pr3ebSoX4PnB/fSOj1SrugeZFLu/fqigptu/nLyIqbceXqxz8nPd2at28GbczP495JNJCbEc1WvFK7r25LWDWtGOXHpZe49yE3j0qkSZ4wd2oc61XWpspQvKjJS7tVKTODRK7oz5JU5PPXxKvoeti5axnf7+Oe8jUyev4ENmfuplViFy09OYde+HMZ+/jUvT1/HGe0bMrhvS87t1DimGhqzDuQybOwcvt2VzcQRfWnRoEbQkUTKTEVGKoSzOjTimtTmvPTZGur2TmTXwo28OTeDz9cULCdwWtuG/OL8jlzQpQmJCQW3o9m6O5tJ6RuYOGc9t/59Hk1qJ3Jtn+YM6tOC5NrBruCZnZPHiPFzWfLtbl66vhe9WlbOplMp/8plkSlcqfMBoI67X1m4rRMFSxI0BD529xcCjCgBeOCSTny2chuPzsmGOQtpXr86PzunA1f0akZKvSOPAhrXTuSuc9pze/+2fLJ8K3+fvZ4n/7uKZz5Zzfmdk/n1RZ1oXj/6Rw+5efnc8doCZq7dwZPX9ODczslRzyASLlEvMmb2CgU32dzq7l0P2z4QeAqIB0a7+2Ml7cPd1wLDzWzyYduWAbcW3sTz5Ujll9hVOzGBZ6/ryXPvpzNiYC9Oad2AuFKsbV8lPo7zuzTh/C5N+GbHXl6bvZ7XZq/nsue/YNyw3nRtVicK6Qvk5zv3Tl7Ef5dt4feXduEnPZtF7b1FIiGICehxwMDDN5hZPPAccCHQGRhkZp3NrJuZTSnyp3FJOzazHwMzgI8jF19iWWqr+gzpUo1T2zYsVYEpqmWDmtx/USfe/umpVI03rh01i89XR2cFT3fn91O+4q0FG7nnvA7c2K9VVN5XJJLM3aP/pmatgCmHjmTMrB/wW3e/oPDx/QDu/ugx9jP50HRZke3vu/vFxWwfCYwESE5O7jVp0qTj/EmCl5WVRVKSltk9XLjG5LvsfP46N5tNe50R3atxStPIHvi/veog767J4YJWVbi2Y9Ww98Los3IkjUnxShqXAQMGzHP31LLsK1bOyTQDNhz2OAPoW9KTzawB8AjQ08zud/dHzaw/cDlQDZha3OvcfRQwCiA1NdX79+8flvBBSktLoyL8HOEUzjE5+8wcRkyYy4tfZtK4eVtuOr11WPZb1JgZ63h3zVdc1SuFP1/ZPSLNlvqsHEljUrxwjkusFJni/kaVeIjl7juAW4tsSwPSwppKKr06NRKYMLwPP5u0kN9P+Yote7K5b+CJYS0Cb87dwMNTvmJglyY8enk3dfNLhRIrTQEZQPPDHqcA3waUReQHEhPieW7wyQzu24KXPlvLPW9+SU5eflj2/cGSzfzqn4s4vV1DnhrUI6b6dETCIVaOZNKB9mbWGthIwQqc1wUbSeT/xccZf/hJV5JrJ/K3j1ayI+sgzw8+mZrVQv8rtHlXNndPWsBJzevy0g29qFZFy0lLxRP1X5vMbCIwE+hoZhlmNtzdc4E7gA+BZcAb7r402tlEjsbMuOuc9jx2eTemr9rG4NGzOZAb+t2fX/xsDXn5ztPX9jyuYiUSy6L+yXb3QSVsn0oJJ+xFYsm1fVpQs1oV7py4gNdmr2fYaWW/GGDrnmwmzlnPZT2bBdLwKRItmgAWCcEl3ZtyatsGPPvJarIO5Jb59S9PW0tOXj4/HdAuAulEYoeKjEgIzIxfDTyRHXsP8vK0tWV67Y6sA/x91nou7dGMVjF8B2iRcFCREQnRSc3rclG3JoyevpbtWQdK/brRM9aRnZunoxipFFRkRI7DPed3JDs3n2c/WV2q5+/cd5AJX3zNxd2a0q6xOs2l4lORETkObRslcXVqCv+Y/Q0bMvcd8/mvzFjH3oN53Hl2+yikEwmeiozIcbr7nA7EmfHERyuP+rxd+wsWShvYpQkdm9SKUjqRYKnIiBynJnUSGXpaK95euJFlm3aX+LzxX3zNngO53HmOzsVI5aEiIxIGt5/VjlrVqvD4hyuK/f6e7BzGzFjHuZ0a0+WE6K1PIxI0FRmRMKhTI4Hb+rfjk+VbmbMu84jvvzrrG3btz9G5GKl0VGREwmToqa1Irl2NP32wnMPXadp3MJfR09dxVodGnNS8boAJRaJPRUYkTKpXjefuczow75vv+HjZ1u+3/2PWejL3HuSuc3QUI5WPioxIGF2dmkKbhjX584fLyct39h/M46VpazmtXQN6tawXdDyRqFOREQmjKvFx3HN+R1ZuyeLtBRuZOGc927MOcJfOxUglpfuLi4TZRd2a0D2lDk98tJLc/Hz6tq5P3zYNgo4lEggdyYiE2aGbZ27cuZ8tuw/oXIxUajqSEYmA09o15NxOjdmfk8epbXUUI5WXioxIhIy6IRUoOLIRqaxUZEQiJC5OxUVE52RERCRiVGRERCRiVGRERCRiVGRERCRiVGRERCRiVGRERCRiVGRERCRi7PB1LyoTM9sGfBN0jjBoCGwPOkSM0ZgUT+NyJI1J8Uoal5bu3qgsO6q0RaaiMLO57p4adI5YojEpnsblSBqT4oVzXDRdJiIiEaMiIyIiEaMiU/6NCjpADNKYFE/jciSNSfHCNi46JyMiIhGjIxkREYkYFRkREYkYFRkREYkYFRkREYkYFZkKysw6m9kbZvaCmV0ZdJ5YYWZnmNmLZjbazL4IOk+sMLP+Zja9cGz6B50nFphZp8LxmGxmtwWdJ1aYWRszG2Nmk0vzfBWZGGRmr5jZVjNbUmT7QDNbYWarzey+Y+zmQuAZd78NuDFiYaMoHOPi7tPd/VZgCjA+knmjJUyfFweygEQgI1JZoyVMn5VlhZ+Vq4EKcVeAMI3LWncfXur31CXMscfMzqTgL/wEd+9auC0eWAmcR8E/AunAICAeeLTILm4q/O9DwD7gVHc/LQrRIyoc4+LuWwtf9wZws7vvjlL8iAnT52W7u+ebWTLwN3cfHK38kRCuz4qZ/Ri4D3jW3V+LVv5ICfPfocnufsxZkirhiy/h4u7TzKxVkc19gNXuvhbAzCYBl7r7o8AlJezqp4UfoLcilTWawjUuZtYC2FURCgyE9fMC8B1QLRI5oylcY+Lu7wHvmdn7QLkvMmH+rJSKikz50QzYcNjjDKBvSU8u/CD9GqgJPB7JYAEr07gUGg6MjVii2FDWz8vlwAVAXeDZyEYLTFnHpD9wOQVFd2pEkwWrrOPSAHgE6Glm9xcWoxKpyJQfVsy2Euc63f1rYGTE0sSOMo0LgLs/FKEssaSsn5e3qCBHvEdR1jFJA9IiFSaGlHVcdgC3lnbnOvFffmQAzQ97nAJ8G1CWWKJxKZ7G5Ugak+JFdFxUZMqPdKC9mbU2s6rAtcB7AWeKBRqX4mlcjqQxKV5Ex0VFJgaZ2URgJtDRzDLMbLi75wJ3AB8Cy4A33H1pkDmjTeNSPI3LkTQmxQtiXHQJs4iIRIyOZEREJGJUZEREJGJUZEREJGJUZEREJGJUZEREJGJUZEREJGJUZEREJGJUZEREJGJUZEQCVrhg1MLCP7PNTH8vpcJQx79IwMxsFXCGu28OOotIuOk3JpHgTQUWm9mTUPY11EVimdaTEQmQmZ1KwXoeTQtvVEjhCoXDVWSkItCRjEiwrgJWunuuFagddCCRcFKREQnWROAWM1sEzALaB5xHJKx04l8kxhy2hvp5wOhjraEuEstUZEREJGI0XSYiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhGjIiMiIhHzf/8RaIuHGk7PAAAAAElFTkSuQmCC\n", "text/plain": [ "

" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.4 Hessian Finite Difference using Approximate Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.4-Hessian-Finite-Difference-using-Approximate-Gradient)", "section": "3.8.1.4 Hessian Finite Difference using Approximate Gradient" } }, "source": [ "### 3.8.1.4 Hessian Finite Difference using Approximate Gradient\n", "\n", "Repeat the analysis from above. Use `my_grad_approx` and the best value for $\\epsilon_1$ you previously found. What value of $\\epsilon_2$ gives the lowest Hessian approximation error? Note: $\\epsilon_1$ is used in the gradient approximation and $\\epsilon_2$ is used in the Hessian approximation." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.4 Hessian Finite Difference using Approximate Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.4-Hessian-Finite-Difference-using-Approximate-Gradient)", "section": "3.8.1.4 Hessian Finite Difference using Approximate Gradient" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEdCAYAAAAvj0GNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VOX1wPHvycaSBAiBhCUQ9jUsQgTBLagoVFFq1Yraal2wWrva/upa92JrW2vV1rog2laQWhdQFK0SRQVBVAxhJ2whLAkQIIGELOf3RwYbQpaZzJ3cWc7neeYhczP33jMvk5y8u6gqxhhjjFOi3A7AGGNMeLHEYowxxlGWWIwxxjjKEosxxhhHWWIxxhjjKEssxhhjHGWJxRhjjKMssRjjMhG5RUQ+F5FyEZnlwPU6ishrIlIqIltF5Io6379cRNZ4vr9JRE73957G1BbjdgDGGAqAB4HzgDYOXO9J4CiQCowE3hKRlaqaKyITgd8B3wWWAV0duJ8xx7EaizG1iEh7EZktIkUickhEvhKRgP6cqOqrqvo6sLeBmLqJyH9EpFBENovITxq6lojEA98B7lbVElX9GJgHfM/zkvuA+1V1qapWq+oOVd3h8FsyEc4SizHH+x1QBaQD7YGrVLXamxNF5E0RKW7g8WZzgvEktfnASqA7cDbwMxE5r4FTBgBVqrq+1rGVwFARiQYygc4islFE8kXkCRFxopZkzDcssRhzvApgA3DY8xf9Kk8tZpmIlIhIRkMnquoFqtqhgccFzYznZKCzqt6vqkdVNQ94Bri8gdcnAAfqHDsAJFLTNBYLXAKcTk0z2UnAXc2MzZh6WWIx5nhrgV8CpSJyo+fYYeB84BUX4kkHutWu/QB3AKkicqUn2ZWIyNue15cA7epcox1wCDjief64qu5U1SLgT8C3WuB9mAhinffGeIjIhcAtwChV3XDsuKpWAIUi0tT5b1NTE6jPYlWd3IywtgObVbV/A9//V53n64EYEelf6z2MAHJVdb+I5AO2pLkJKKuxGPM/GcAOYDeAiPQUkSRvT1bVyaqa0MCjwaQiIjEi0hqIBqJFpLWIHPujbxlwUER+LSJtRCRaRDJE5OQGYigFXgXuF5F4ETkVuAj4h+clzwM/FpEUz3v7GdCs/h9jGmKJxZj/eZ6aYbrbReQANb+gW6JWfxc1zVS3AVd5vr4LQFWrgCnU9IdsBoqAZ6kZWNCQm6kZtrwHmA3cpKq5nu89ACynpmazBvgSeMjZt2MindhGX8Z4xzN58Q+qusrtWIwJZlZjMcYLIrIAOBd4RkSucTkcY4Ka1ViMMcY4ymosxhhjHGWJxRhjjKMssRhjjHGUJRZjjDGOiqiZ9506ddJevXq5HYbfSktLiY+PdzuMoGPlciIrk/pZuZyosTJZsWJFkap29vZaEZFYRGQKMKVfv358/vnnbofjt+zsbLKystwOI+hYuZzIyqR+Vi4naqxMRGSrL9eKiKYwVZ2vqtPbt29ssrIxxhgnRERiEZEpIvL0gQN1VxM3xhjjtIhILFZjMcaYlhMRicVqLMYY03IiIrFYjcUYY1pORCQWY4wxLSciEos1hRljItm6XYdoyQWHIyKxWFOYMSZSbd1bypQnPuaJDza22D0jIrEYY0wkUlXunZdLXHQUl53co8XuGxGJxZrCjDGRaGHubhatK+Rn5/QntV3rFrtvRCQWawozxkSaw0cruX9+LoO6JHLN+F4teu+ISCzGGBNp/vL+RgoOlPHg1Axiolv2V70lFmOMCTMbdh/i2cV5XDI6jcxeHVv8/pZYjDEmjKgqd7+xivhWMdw+eZArMUREYrHOe2NMpJi3soClefv41XkDSU5o5UoMEZFYrPPeGBMJDpZV8MCbaxiR1p5pY3q6FkdEbPRljDGR4E/vrmdvaTkzr8kkOkpciyMiaizGGBPuVu04wItLtnDV2HSGp3VwNRZLLMYYE+Kqq2s67JPaxvHLcwe6HU7oJhYRGSwiT4nIKyJyk9vxGGOMW+Z+vp0vtxVz+7cG075trNvhuJNYRGSmiOwRkVV1jk8SkXUislFEbmvsGqq6RlV/CFwGZAYyXmOMCVb7So/y8DtrGdOrI98Z1d3tcAD3aiyzgEm1D4hINPAkMBkYAkwTkSEiMkxE3qzzSPGccyHwMfB+y4ZvjDHB4ffvrOVQWSX3Tx2KiHsd9rVJS67Rf9yNRXoBb6pqhuf5OOBeVT3P8/x2AFWd4cW13lLV8xv43nRgOkBqauroOXPmOBK/m0pKSkhISHA7jKBj5XIiK5P6hUu5rN9fxW8/K2NSrxguH+TfnJXGymTChAkrVNXrlqFgGm7cHdhe63k+MLahF4tIFnAx0ApY0NDrVPVpEdkJTElMTBydlZXlSLBuys7OJhzeh9OsXE5kZVK/cCiXo5XV/PbxxXRr35o//uBM4lv59+vcyTIJpsRSXx2uweqUqmYD2d5cWFXnA/MzMzNvaFZkxhgTZJ79OI/1u0t45vuZficVpwXTqLB8oPZONGlAgRMXtiVdjDHhZPu+w/zl/Q2cOySViUNS3Q7nBMGUWJYD/UWkt4jEAZcD81yOyRhjgsqxRSajRbj3wqFuh1Mvt4YbzwaWAANFJF9ErlPVSuAWYCGwBpirqrlO3M/WCjPGhIsFObvIXlfIL84dSLcObdwOp16uNMyp6rQGji+gkY745hKRKcCUfv36OX1pY4xpMQfLKrhvfi5Du7Xj6nHpbofToGBqCgsYq7EYY8LBHxeuo7CknN9+e1iL7wrpi+CNzEHWeW+MCXUrtxfz4tKtfP+UdEb0cHeRyaZERGKxGosxJpRVVlVzx2s5dE5oxa3nub/IZFMiIrEYY0woe2HJVnILDnLPlKG0a+3+IpNNiYjEYk1hxphQVVB8hD++u46sgZ351rAubofjlWYlFhGZKCLPiMhIz/PpzoblLGsKM8aEqnvn5VKtygMXZQTNIpNNae5w45uBHwB3iUhHYKRzIRljjAF4b/Vu3l29m19PGkSPjm3dDsdrzW0KK1TVYlX9JXAucLKDMRljTMQrLa/knjdWMSA1getP7+12OD5pbmJ569gXqnob8KIz4QSG9bEYY0LNn95bT8GBMmZcPIzYIJ6zUp8moxWRnnUfwJd1nr8hIu0CH27zWB+LMSaUrNpxgOc/2cyVY3syOr2j2+H4zJs+lhe8eI1SsytkUNdcjDEm2FVWVXP7qzkkJ7Ti/yYNcjucZmkysajqhJYIxBhjDLy4ZCs5Ow7wxBUn0b5N8M9ZqU+TicXT1OWNYlU96Gc8AWGLUBpjQkHtOSvnD+vqdjjNFhFNYbaDpDEmFNwzL5eqEJuzUh9rCjPGmCCwMHcX763eze2TQ2vOSn1CawybMcaEoZLySu55I5dBXRK59rTQmrNSn2YnFhH5Ra2vg3+5TWOMCVJ/WLiO3YdCc85KfXx+ByLSQUSeBy4RkZtF5DTgNudD8yqWeBFZISIXuHF/Y4zx18rtxbywZAvfOyWdk3omuR2OI3xOLJ6lXH4AzAY+A/oD7/tyDRGZKSJ7RGRVneOTRGSdiGwUEW+S1a+Bub7c2xhjgkXtfVZ+GQL7rHjLnzpXT1VdoarPA8N8PHcWMKn2ARGJBp4EJgNDgGkiMkREhonIm3UeKSJyDrAa2O3He/DajuIjLXEbY0wEmfXpFnILDnLfhaGxz4q3mru6MUCUiJwOfAIk+3Kiqn4kIr3qHB4DbFTVPAARmQNcpKozgBOaukRkAhBPTRI6IiILVLXa53fhhWcX5/Hoe+t5+cZxZHS3ZWGMMf7bUXyEP723nrMHpTApIzT2WfGWqGrzTqypYdwEjAJeV9V5Pp7fC3hTVTM8zy8BJqnq9Z7n3wPGquotTVznGqBIVd9s4PvTgekAqampo+fMmeNLmADsL6vmwaVlVFTDXae0JqWtu51rJSUlJCQkuBpDMLJyOZGVSf3cLhdV5c9flLNmXxW/Pa0Nndq432HfWJlMmDBhhapmenutZtdYVLUKeKK559ejvtlATWY9VZ3VxPefFpGdwJTExMTRWVlZzQpu2KhDfOdvS/hrrvCfm8aRnNCqWddxQnZ2Ns19H+HMyuVEVib1c7tcFuTsZGXhF9x1/mAuOb2Pa3HU5mSZ+DPc+F4ReVtEnhWRHzkQSz7Qo9bzNKDAges6srpxv5REZl6Tyc4DZVw7azml5ZVOhGaMiTAHyyq4d14uQ7u145rxvdwOJyD8qX91AJYCDwFODGdYDvQXkd4iEgdcDvjUvNYQp/ZjGZ3ekSeuGEXOjgP86KUvqKgKSJeOMSaMPfLOOopKyplx8TBiwmDOSn38eVf7gGhgj+drr4nIbGAJMFBE8kXkOlWtBG4BFgJrgLmqmutHfAExcUgqD04dRva6Qm5/NYfm9lEZYyLPiq37+ednW7l6fC+Gp3VwO5yA8WdU2ENAKvAXYFUTrz2Oqk5r4PgCYIEfMTV0P0cXobxibE92Hyzjsfc3kNquFb86LzT3TDDGtJyKqmrueDWHLu1ac+u54TNnpT7+1Fj+Ts1S+dcBKxyKJyACsTXxz87pz7QxPXhy0SZeXLLFsesaY8LTM4vzWLf7EPdflEFCK3/+pg9+/iSWe4DnROQfwMkOxRMQgdiaWER44KIMzhmcyj3zclmQs9Oxaxtjwsu2vYd57L8bOG9oKhOHpLodTsD5k1geANZRMyQ4qJdVCUSNBSAmOorHp53EST068LM5X7E0b6/f1yyvrHIgMmNMsFBV7nw9h9joKO67MMPtcFqEP4nl/1T1XmomSd7jTDiBEYgayzFt4qJ57uqT6dGxDVfPXMb8lc0bIV1RVc1v3ljFsHveZWHuLoejNMa4Zd7KAhZvKOJX5w2kS/vWbofTIvxJLPeJyCxVLQVecyqgUJQUH8fcG8cxPK09P579JX96dx3V1d6PFis8VM6Vz3zGi0u20jE+jh/P/pLPHKj9GGPcVXz4KPfPX82IHh246pR0t8NpMf4klqNAnufr0x2IJWAC1RRWW3JCK/51/SlclpnGXz7YyM3/+oLDR5ueRLlyezEXPvExX+8o5rHLR7Lgp6eTltSG61/8nDU7DwYsXmNM4D389lqKj1Qw49vDiI4K3a2GfeVPYjkMtBeRWKCnQ/EERCCbwmqLi4nid98Zzt0XDOHd1bu45G9LGl0V+ZUV+Vz69yVEifCfm8Zz0cjudIyP4x/XjSU+LoarZy5j+77DAY3ZGBMYyzbvY87y7Vx/Wm+GdGvndjgtyt9RYZuoWer+JWfCCX0iwnWn9ea5a05m+77DXPTEJ6zYuv+411RUVXPvvFx++e+VZKYnMf/HpzG02/+SXvcObXjxujGUV1bz/ZnLKCopb+m3YYzxQ3llFbe/+jVpSW346Tn93Q6nxTVnB8l4EYlW1UpV/auqTm9oZeFg0RJNYXVNGJjCaz8aT3yraKY9vZT/rMgHoKiknKue/YxZn27h+tN68+K1Y+gYH3fC+QNSj61NdoRrZy2nxNYmMyZkPJWdx6bCUh6cmkHbuPCes1KfJhOLiESJyBUi8paI7AHWAjtFJFdEHhGRoE/HLdUUVle/lERev/lUMnslceu/V3L7qzlc+PjHfLW9mD9/dyR3XTCk0bWCRqd35MkrRpFbcJAf/mMFRyttbTJjgt2mwhKeXLSRKSO6kTUwxe1wXOFNjWUR0Be4Heiiqj1UNYWaDvulwMMiclUAYwxpSfFxvHDtGK46pSezl21DPP0pU0/q7tX5Zw9O5eGLh/HxxiJu/fdKn0abGWNalqpy52s5tI6N4u4LBrsdjmu8qaOdo6oVdQ+q6j7gP8B/PB34pgGx0VE8OHUY3xrWlcFd2pFUT9NXYy7N7MHe0qM8/PZakuPjODPRkosxweiVFfkszdvHb789jJTEyJizUp8mE8uxpCIifwZ+rvUs51tf4jEnGt+3U7PPvfGMPhQdKufZjzezp3cs406ronVstIPRGWP8sbeknIcWrCEzPYnLT+7R9AlhzJfO+xJgnojEA4jIuSLySWDCcpYbnfdOExHu+NZgLstMY8HmCs74/SKe/2QzZRW2BIwxweCht9ZQWl7JjIuHERVBc1bq43ViUdW7gNlAtoh8DNwK3BaowJzkVue906KihN9fMoJfn9ya3p3iuW/+akswxgSBjzcU8eqXO/jhmX3pn5rodjiu8zqxiMjZwA1AKdAZ+ImqLg5UYKZhg5OjefnGccy+4RR6eRLMmY8sYpYlGGNaXFlFFXe+nkPvTvH8aEI/t8MJCr40hd0J3K2qWcAlwMsiclZAojJeGdc3mbk3juOlG8aS3jGeez0J5oVPt9gqyca0kMc/2MDWvYd5aGqG9Xt6eD1zR1XPqvV1johMpmZU2PhABGa8N75vJ8b1SWbJpr08+t/13DMvl7dydvL8NScTH+YbChnjpnW7DvH3D/P4zqg0xvdr/uCccOPNBMl6e6FUdSdwdmOvCTQRyRKRxSLylIhkuRFDsBARxvfrxNwbx/Gny0bw+ZZ9/GDWcq8WwjTG+K66WrnjtRwSW8dw5/mRO2elPl5NkBSRH4vIcQtNikgcME5EXgCu9vXGIjJTRPaIyKo6xyeJyDoR2SgiTQ0OUGpGq7UG8n2NIRyJCBePSuPR746sSS7PW3IxJhBeWraNFVv3c+f5Q+pdlimSeZNYJgFVwGwRKRCR1SKSB2wApgGPquqsZtx7lufa3xCRaGoWtZwMDAGmicgQERkmIm/WeaQAi1V1MvBr4L5mxBC2LhrZnUe/O5LlW/ZxrdVcjHHUnoNl/O6dtYzvm8x3Rnm3ikYkkXrmOzb84poZ9p2AI6pa7PfNRXoBb6pqhuf5OOBeVT3P8/x2AFWd0cR14oCXVPWSer43HZgOkJqaOnrOnDn+hu26kpISEhISvHrtkoJKnv66nIEdo/j5qNa0ignf8fW+lEuksDKpn7/l8tevyvhiTxUPntqGLvH+LBIfPBorkwkTJqxQ1Uxvr+VTz65nhv1OX87xUXdge63n+cDYhl4sIhcD5wEdgCfqe42qPg08DZCZmalZWVlOxeqa7OxsvH0fWcDgwTv4xdyvmJXXhpnXnEybuPAcueJLuUQKK5P6+VMui9buYdmu5dw6cQCXnx30a/B6zcnPSqOpVkQmisgzIjLS83y6I3dt5Jb1HGuwSqWqr6rqjar6XVXNbvCiYTDz3h9TT+rOny4byWeb93LtrOUcOWpDkY1pjsNHK7nr9VX0T0ngxjP7uh1O0GqqDncz8CvgKs+clZEBjicfqL3IThpQEOB7RoSpJ3Xnj5eNsORijB8efW89O4qP8NuLhxEXEx5NYIHQVMkUqmqxqv4SOBc4OcDxLAf6i0hvT7/J5cA8fy8aLku6+OvbJ6V9k1yue2G5zdI3xgerdhxg5idbmDamJyf36uh2OEGtqcTylog85vn6fuBFp24sIrOBJcBAEckXketUtRK4BVgIrAHmqmquA/eK6Kaw2r59Uhq/v2QEn27ay+xl29wOx5iQUOWZs5LUNo7bJg1yO5yg12hiUdU38EyCpGZo7+NO3VhVp6lqV1WNVdU0VX3Oc3yBqg5Q1b6q+pBD97IaSy2XjE4jMz2JZxdvpqLKdqU0pikvfLqFr/MP8JspQ2jf1rafaoo3jYTviMgSoIuIXCsio0UkpHawsRrLiX54Zl92FB/hra8DOcjPmNBXUHyEP767jjMHdGbK8K5uhxMSmkwsnv6VK6mZJNkbuBvI8ex5/3KA43OE1VhOdNagFPqnJPDUh5vwZS6TMZFEVfnNG7lUqfLg1AxcWr0q5Hg1rEFV86jZovhuVZ2qqv2pmV/yaECjc4jVWE4UFSVMP6MPa3cdInt9odvhGBOUFubu4r9rdvPzcwbQo2Nbt8MJGV4lFhF5ACivfUxVS1R1aUCicpjVWOp30cjudG3fmqeyN7kdijFB52BZBffMy2Vw13Zce1pvt8MJKd4OxP4p8J6ILBSRS0XE1mIPA3ExUVx3Wm8+27yPL7ftdzscY4LKHxauY8+hcmZcPIzYaJuz4gtvS2u3qg4AHgYuBDaKyCMiMjBwoTnHmsIadvmYnrRrHcPfP8xzOxRjgsYX2/bzj6VbuXpcL0b26OB2OCHH28SiAKq6SFW/BwwHDgOrRGRIoIJzijWFNSyhVQzfG5fOwtW72FRY4nY4xriuoqqaO17NITWxNbeeO8DtcEKST/U7EYkSkQuBfwDfpWaE2OZABGZazjXjexMbHcUzH1mtxZhnF29m7a5D3HvhUBJb25yV5vA2sSSIyMNAHjXLrDymqoNU9WFVPRK48JxhTWGN65zYiktHp/HqFzvYc7DM7XCMcc22vYd57P31TBySyqSMLm6HE7K8TSz7qFkuf5SqXqGqHwQwJsdZU1jTbji9D5XV1cz8ZIvboRjjClXlztdziBbhvguHuh1OSPN2HkuGqj6mqvsCHZBxR69O8UzO6Mq/lm7lYFmF2+EY0+LmrSxg8YYifnXeQLp1aON2OCHNxtCZb/zwzL4cKq/kpc9scUoTWYoPH+X++asZ0aMD3xvXy+1wQp7PiUVEpgQiEOO+YWntObVfMjM/3kx5pS2pbyLHjAVrKT5SwYxvDyM6ypZt8VdzaiyOrDjckqzz3ns/PLMvew6V8/qXO9wOxZgW8VneXl7+fDvXn9abId3auR1OWGhOYgm5dG6d9947rV8nhnZrx98/yqO62hanNOGtvLKKO17LIS2pDT89J3z2r3dbcxKL/bYJYyLCjWf2Ja+wlHdX73Y7HGMC6qnsPDYVlvLg1AzaxtlKVU6xzntzgm9ldKFHxzY8/PYanvpwE59uKuKQjRQzYWZTYQlPLtrIlBHdyBqY4nY4YSVkU7SIRAEPAO2Az1X1BZdDChsx0VHcf2EGv5m3ioffXguACPTtnMDwtPaM7NGB4WkdGNw1kVYx0S5Ha4zvVJU7X8uhdWwUv7kg6FelCjnNSSx+t4+IyEzgAmCPqmbUOj4JeAyIBp5V1YcbucxFQHdqJm/m+xuTOd6EQSksHnQW+0qP8nV+MV/nH2Dl9mI+Wl/Eq1/UdOzHRgv9UxIZ2CWRAamJDOqSyIAuiXRr39o2RDJB7eMdlSzN28eMi4fRObGV2+GEHZ8Ti6pOdOC+s4AngBePHRCRaOBJYCI1iWK5iMyjJsnMqHP+tcBAYImq/l1EXgHedyAuU0fH+DiyBqZ801Sgquw8UMbK7cWszD/Amp0HWZq3l9dqjSJLbBXDgFrJZmyfjgxMTbRkY4LC3pJy5qw7ysm9kvhuZg+3wwlLrjSFqepHItKrzuExwEbPbpWIyBzgIlWdQU3t5jgikg8c9Ty1SRctRETo1qEN3Tq0YfKw/+3/feBIBet3H2LdLs9j9yEW5Oxk9rKayZYpia04rX8nzhzQmVP7daJTgv2VaNzx4FtrKKuEGRcPI8rmrASENHe/cxGJVtVm/0L3JJY3jzWFicglwCRVvd7z/HvAWFW9pYHz2wKPU7N8/1pVfbKB100HpgOkpqaOnjNnTnNDDholJSUkJCS4HUaTVJW9ZcrqvVWsKqoid28VpZ4xAOntohiaHE1Gp2j6J0UR68APeKiUS0uyMjneqqIq/vB5GZN6KJcPtXKprbHPyoQJE1aoaqa31/KnxvK0iPxYVQ+LyBmq+pEf14L658c0mPVU9TBwXVMXVdWnRWQnMCUxMXF0VlZW8yMMEtnZ2YTi+6iqVnILDrB4QxEfrS/k3a37WbC5gvTktiz82em0jvVvIEColksgWZn8T1lFFff8+SN6d4rn4kFq5VKHk58Vf4Yb/wZ4TkT+AZzsQCz5QO0GzzSgwIHr2gTJIBEdJQxP68CPJvTj5RvH8dU953Lb5EFs3XuYVTtsVQQTWH95fwNb9x7moakZxEVbE1gg+ZNYHgDWUVOrmOtALMuB/iLSW0TiqNn3ZZ4D17UlXYJUQqsYLh2dBsCKrftdjsaEs7W7DvL0R3l8Z1Qa4/t1cjucsOdPYvk/Vb0XuAm4x5cTRWQ2sAQYKCL5InKdqlYCtwALgTXAXFXN9SM+EwKSE1rRu1M8n1tiMQFSXa3c8WoOia1juPP8wW6HExGa3ceiqkWef0tF5EYfz53WwPEFwILmxtTI/eYD8zMzM29w+trGf6PTk1i0dg+qakOSjeNeWraNL7YV88dLR9AxPs7tcCJCs2ssIvKAiPxbRGYB/ZwLyXnWFBbcRqcnsbf0KFv2HnY7FBNm9hws43fvrOXUfslcPKq72+FEDH+awpJU9VJqhvL+xKF4AsI674Pb6PQkwPpZjPPum7+a8spqHpw6zGrDLcifxFIuIidR03kf71A8AWE1luDWr3MC7VrHWGIxjvpg7W7eytnJT87qR+9OQf0rKuz4k1geAc4BngZediacwLAaS3CLihJGpSexYus+t0MxYaK0vJK7X89lQGoC08/o63Y4EcefCZLfV9XfOxaJiWijeyaRva6QA0cqaN8m1u1wTIh79L317Cg+wis/HEdcjO0O0tL8KfGLROQWERnoWDQBYk1hwW90r5p+li+2WXOY8c+qHQeY+clmrhjbk8xeHd0OJyL5k1guAYqBqSLyrEPxBIQ1hQW/kT06EB0lfGH9LMYPlVXV3P5qDskJrfj1pEFuhxOx/GkKexwoBNoDzzgTjolUbeNiGNK1nXXgG7+8sGQrOTsO8MQVJ1mTqov8qbGsU9WbVPUKamovxvhldHoSX20vprKq2qfzDhyp4NSHP2D5rsoARWZCwY7iI/zx3XVMGNiZ82tt6WBanj+JZZKI3Coi5wBB/RNtfSyhYVR6EoePVrF21yGfzlu4ahc7io/waUFQfwxNAKkq97yxClW4/6IMm7PiMr8SCzVNYeOBziIStHvOWx9LaMhs5kTJ+V/XLIK9em8V5ZW251skemfVLv67Zg8/n9ifHh3buh1OxPMnsfwNGAcMAp5R1audCclEqm4d2tC1fWufFqQsPFTOJxuLGJ7WnvIq+CzP5sJEmoNlFdwzL5chXdtx7am93Q7HYH0sJsiMSk/yaWTYgpydVCs8cFEGsVHwwdo9AYzOBKNH3llHUUk5My4eRky0zVkJBtbHYoJKZnoSO4qPsPPAEa9eP39lAQNTExnRowODk6M1yGd0AAAd5klEQVRZtK5mlWQTGVZs3c8/P9vK98f1YkSPDm6HYzz87WNZA5yK9bEYh/iyIOWO4iN8vnU/F47sBsCIztFs3XuYvKLSgMZogkNFVTV3vJpDl3at+eV5QT9PO6I0O7GoaqGqLlDV+1T1CutjMU4Y3LUdbWKjvUosb66s6bS/YHjN0NIRnaMBWGTNYRHhmcV5rNt9iPsuHEpCK3+m5Bmn+bMfy70i8raIPCsiP3IyKBO5YqOjGNGjvVeJZd7KAkb06EB6cs3KtZ3aRDEgNcH6WSLA1r2lPPbfDZw3NJVzh3ZxOxxThz9NYR2ApcBDgNVDjWNGpyeRW3CQw0cb7rrbVFhCbsFBpgw/fiLchEEpLNu8j0NlFYEO07hEVbnr9VXERkdx34UZbodj6uFPYtkHRAN7PF+3OBE5XUSe8tSaPnUjBuO8zPSOVFUrK7c3PNhi/soCRGDKiG7HHT9rYAqV1crHG4oCHaZxyRtfFbB4QxG/Om8gXdq3djscUw9/+ljuB54C/gL4PNxKRGaKyB4RWVXn+CQRWSciG0XktiZiWKyqPwTeBIJ28IDxzUk9a0b3NLTSsaoyb2UBY3t3JLXd8b9YRqcn0a51jDWHhaniw0d54M3VjOzRgatOSXc7HNMAn3u8RCQeKFPVKlUtAK5r5r1nAU8AL9a6djTwJDARyAeWi8g8ampGM+qcf62qHvvtcQVwfTPjMEGmQ9s4+qUkNNjPsnrnQfIKS7nutBMnw8VER3HGgM4sWldIdbUSFWVLe4ST3y5YQ/GRCv558TCi7f82aDVZYxGRKBG5QkTeEpE9wFpgp4jkisgjItK/OTdW1Y84sQltDLBRVfNU9SgwB7hIVXNU9YI6jz2e+HoCB1T1YHPiMMEpMz2JFVv3U1194pyUeSsLiIkSJmfUv9DgWYNSKCopZ1WBzVsKJ5/l7WXu5/lcf3pvBndt53Y4phHe1FgWAf8FbgdWqWo1gIh0BCYAD4vIa6r6Twfi6Q5sr/U8HxjbxDnXAc839E0RmQ5MB0hNTSU7O9vPEN1XUlISFu+jMfFHKjhwpII5CxbRLeF/f/+oKq98doQhHaP4evnx3WrHyiXmqCLAzHeWMbVfXAtHHlzC5bNSUa3c/ckROrURRsXtIjt7t1/XC5dycZKTZeJNYjlHVU8YYqOq+4D/AP8REac2PqivbtvoNGpVvaeJ7z8tIjuBKYmJiaOzsrL8CC84ZGdnEw7vozE9C0t4btWHRKX0I2tMz2+Or9i6j70Ll3DnhRlkjUo77pza5fL8xk/YXKZkZZ3WkmEHnXD5rDz23w3sKl3PrB+cTNbAFL+vFy7l4iQny6TJprD6kko9fuFALFBTQ+lR63kaUODvRW3mfejp3SmepLaxJ/SzzPuqgFYxUUwcktro+WcNTGFl/gEKD5UHMkzTAjYVlvDkoo1MGdHNkaRiAq9Zo8JEZG6tx79xruN8OdBfRHqLSBxwOTDP34vaWmGhR0QYnZ7EilojwyqrqnkrZydnDUohsXXjleQJg2p+AWWvs9FhoUxVufO1HFrHRnH3BYPdDsd4qbnDjQ+q6mWex6XU9MH4RERmA0uAgSKSLyLXqWolcAuwkJp1yOaqam4zYzQhbnR6R/IKS9lXehSApXn7KCo5esLclfoM7daOlMRWLLLEEtJeWZHP0rx93DZ5MCmJNmclVDQ3sTxU5/mdvl5AVaepaldVjVXVNFV9znN8gaoOUNW+qlr3Ps1iTWGh6diClMeW0Z+/soCEVjGcNajp5hARYcLAFBavL6LCx62OTXDYV3qU3y5YQ2Z6Epef3KPpE0zQaFZiUdXNdZ4H9e5K1hQWmoantSc2Wvh8637KK6t4e9VOzh2SSuvYaK/OnzAohUPllSzfEtQfT9OAB99aTUl5Jb+9eJjNRwoxXicWEVlS53miiJzkfEjOsxpLaGodG83Qbu35Yut+Fq8v4mBZpVfNYMec1r8TsdHi6GrHT3ywgZXbix27nqnfpxuLePWLHdx4Rl8GpCa6HY7xkS81llYAIvInAFU9BPw1EEE5zWosoWt0ehIr84t59ct8OrSN5bT+nbw+N6FVDGN7Jzu2vMvK7cX84d313PX6KttMLIDKKqq48/VV9Epuyy1n9XM7HNMMviQWEZEU4CoROVYvbROAmBxnNZbQlZmeRHllNQtydjE5oyuxPm49O2FQCpsKS9m297Dfsfxz6VYAcnYc4CNb5DJgnly0kc1FpTw4dZjXzZ4muPjyU3o78DHwEvCoiNzs4/nG+GyUpwMf4EIfmsGOOdbR/8Fa/2ZqHzhcwfyvC7h0dBpd27fmyQ82+nU9U78Nuw/x1Ieb+PZJ3X2qnZrg4nViUNV3PKO1fga8DPSj+QtQtihrCgtdqe1ak5bUhpTEVozp3dHn83t3iqd3p3gWrSv0K45XvsinrKKaa07txfQz+rBsyz6WbbZBAU6qrlbueC2H+FYx3HW+zVkJZc0dFbZEVX+hqsudDigQrCkstN0zZSgPf6f5q9lOGJjCkry9jW4c1hhV5V+fbeWknh0Y2q09l5/ck+T4OJ5YZLUWJ839fDvLt+znjm8NJjmhldvhGD9YU5YJehOHpHLWoMaXcGnMWYNSOFpZzacb9zbr/CWb9pJXWMpVY2v2/2gTF811p/fmo/WFfJ1vI8ScUHionBlvr2Vs745cOjqt6RNMUPNluHGXOs+7ikhI/FlhTWGRbUzvjsTHRfNBM2fh//OzrXRoG8v5tbZB/t4p6bRrHcOTVmtxxANvrubI0Soe+vYw/jc2yIQqX2osz9V5/g9grYj8wcF4AsKawiJbXEwUp/XvxKK1e3weJrznYBnv5u7m0tFpx41QSmwdyzXje7Ewdzfrdh1yOuSIkr1uD/NWFnDzhL70S0lwOxzjAF8678+v8/wcoA+N7IViTLCYOKQLOw+UsSBnl0/nvbx8O5XVyhVjT9wG9wen9qZtXDR/zbZaS3MdOVrF3W+som/neG7K6ut2OMYhvjSF/Vnq1FG1hi0SaYLe1JHdGJ7WnrvfWEVRiXdL6VdVK7OXbeP0/p3o3Sn+hO8nxcdx5diezF9ZwJaiUqdDjgh/fn892/cd4bffHkarGJuzEi58aQorAeZ59rxHRM4VkU8CE5YxzoqJjuKPl46gpLySO1/L8apJ7IO1eyg4UMaVY3s2+JobTu9DTHQUT324yclwI8LqgoM8u3gz383swdg+yW6HYxzkS1PYXcBsIFtEPgZuBW4LVGBOss57A9A/NZFbJw5gYe5u5q1sev+4fy7dSmq7VpwzuOERaSntWvPdzB7854t8CoqPOBluWKuqVm5/LYektrHc/q1BbodjHOZLU9jZwA1AKdAZ+ImqLg5UYE6yzntzzPWn92FUzw785o1cdh8sa/B12/Ye5qMNhVx+ck9imlhG5sYz+6AKT3+U53S4YeufS7eycnsxd18whA5t49wOxzjMl6awO4G7VTULuAR4WUTOCkhUxgRIdJTwh0tHUF5ZxR2vNtwk9q9lW4kSYdqYhpvBjklLasvUk7ozZ/k2r/tvItmuA2U8snAdp/fv1Kxlekzw86Up7CxV/djzdQ4wGXgwUIEZEyh9Oifwf+cN4v21e3hlRf4J3y+vrOLfn+dzzuAUurT3btfCm7L6Ul5ZzXMfb276xRHunnmrqKyu5qGpNmclXDWZWOqOBDtGVXcCZzf2GmOC1TXjezGmd0fun7/6hL6Rt3N2sa/0KFedcuIQ44b07ZzAt4Z15R9LtnLgcIXT4YaNd3N3sTB3Nz89ewA9k9u6HY4JEG9qLItE5MciclybgIjEAeNE5AXg6oBE1wgR6Ski80RkpoiExCACEzyiooQ/XDKCKlV+/Z+vj2sS++fSraQnt+XUvr6trvujrH6UlFfywpItzgYbJkrKK7lnXi6DuiRy/em93Q7HBJA3iWUSUAXMFpECEVktInnABmAa8KiqzvLlpp5ksEdEVtU5PklE1onIRi+SxQDgLVW9Fhjiy/2NAeiZ3JbbvzWYxRuKmL1sOwBrdx3k8637uXJsT5+3wx3SrR1nD0ph5ieb+XLbftsMrI4/LFzHroNlzLh4mM/76pjQ0uT/rqqWqepfVfVUIJ2a5q9Rqpquqjeo6lfNuO8sahLWN0QkGniSmr6bIcA0ERkiIsNE5M06jxTgS+ByEfkAWNSMGIzhyjE9ObVfMg+9tZrt+w7zz6VbiYuJ4tLRPZp1vZ9PHEBllfLtv37K5McW88KnWzhwxJrGVm4v5oUlW/jeKemc1DOpydeb0CZu/VUlIr2AN1U1w/N8HHCvqp7neX47gKrOaOD8XwLLVPUjEXlFVS9p4HXTgekAqampo+fMmeP0W2lxJSUlJCTYmkp1Nbdc9h6p5s6Pj9CzXRTbDlYzKjWG6cObv77qkUplaUEl2fmVbD1YTVwUnNwlhqweMfTrENWiHdbB8FmpqlbuX1rGgXJlxultaBPjfpdsMJRLsGmsTCZMmLBCVTO9vVZMY98UkYnAZcCTqvqViExX1ad9itZ73YHttZ7nA2Mbef07wL0icgWwpaEXeeJ9GiAzM1OzsrL8DtRt2dnZhMP7cJo/5VLZaRu//k8OALdeNIbR6f79VT3Z829O/gFeWraNeV/t4JOCMvqnJDBtTE8uHtW9ReZvBMNn5flPNrP14Gr+euUoJg/r2vQJLSAYyiXYOFkmjSYW4GbgB8BdItIRGOnIXetX358xDVanVHUVNfNpmr6wyBRgSr9+/ZoZmgl3l2X24KP1RewrPcqonh0cu+6wtPbMSBvGXecPZv7KAmYv28b9b67mr9kbeemGUxiQmujYvYLRnoNl/PHd9ZwxoDOTM7o0fYIJC031sRSqarGq/hI4Fzg5gLHkA7UbttOAptfdMMYBIsKTV47in9ePDUhTVXyrGC4f05M3bjmN1390KlEiXPHMUjbsDu8l9x94aw1Hq6q5/8KhNmclgjSVWN4Skcc8X98PvBjAWJYD/UWkt2co8+XAPCcubEu6GG81d/tjX4zs0YGXbjgFEWHaM5+xcU9JwO/phsUbCpm/soCbs/rSq57VoU34ajSxqOobeCZBAotV9XEnbiois4ElwEARyReR61S1ErgFWAisAeY6tSS/LUJpgk2/lARm31DThTjtmaVsKgyv5FJeWcVv3silV3Jbfnim7bMSabwZTP6OiCwBuojItSIyWkS8W+eiAao6TVW7qmqsqqap6nOe4wtUdYCq9lXVh/y5R537WY3FBJ1+KYnMvmEsqsq0p5eSF0bJ5e8f5rG5qJT7L8o4budNExm8mcfyS+BKaiZJ9gbuBnJEJFdEXg5wfI6wGosJVv1TE3nphlOoqlamPbOUzWGwYdjWvaU8sWgj5w/vyhkDOrsdjnGBV9NfVTUPOEdV71bVqaran5qhwI8GNDqHWI3FBLMBnuRSUVVTcwnl3ShVld+8kUtcdBS/ucAWxIhUXiUWEXkAOG49cFUtUdWlAYnKYVZjMcFuYJdE/nX9WMorq5j2zFK27g3N5PLOql18uL6Qn08cQGo7v1rMTQjzdsGenwLvichCEblURJqa/xJUrMZiQsHgru341/WncKSiimlPL2X7vsNuh+STkvJK7pu/msFd23H1OO9Xhjbhx9vEsltVBwAPAxcCG0XkEREZGLjQnGM1FhMqhnRrx7+uH8uhskoeWbjO7XB88th/17PrYBkPTs1octdNE968/d9XAFVdpKrfA4YDh4FVIhL0DalWYzGhZGi39pw9OIVPNxWFzArJa3YeZOYnW5g2poffy+GY0OfTnxUiEiUiFwL/AL5LzQgx2zLPGIeN79uJopKjbAiByZPV1cpdr6+ifZtY/u+8QW6HY4KAt4klQUQeBvKomRH/mKoOUtWHVfVIE+caY3w0rm8yAJ9uLHI5kqa9/Pl2Vmzdz22TB5EUH/iFNU3w8zax7AN2UrMPyxWq+kEAY3Kc9bGYUNOjY1vSktqwJG+v26E0avu+wzz45mrG9UnmklFpbodjgkSTo7s8WxJ/y/M0QUQa2sSgWFUPOhaZg1R1PjA/MzPzBrdjMcZb4/smszB3N9XV6vNullCzD8rhisD10VRXK7/890pEhEcuHd6sGE148mbY8AvUdN439qlRanaFDOQilcZElPF9OzH383xW7zxIRnffB5489eEm/vrRYc48o5L4Vs7PEJj16RY+27yP339nOGlJbR2/vgldTX7aVHVCSwRijDnesX6WJZv2NiuxzPuqgNIKyF5XyPnDnd1ga1NhCb97Zy1nDUrh0kxrAjPHa7KPRUR6evlo1xIBGxMpUtu1pk/neD7d5HsH/paiUtZ59np5J3eXo3FVVlVz69yVtImL5uGLh9k+K+YE3jaFNSWom8JsB0kTqsb3Tea1L3ZQUVVNrA+TDt9dXZNMhiZH8cGa3ZRVVDm2yvDfP8rjq+3FPD7tJFJs2RZTD29WN57gxeMsVQ3KpAI2QdKErnF9OlF6tIqcHb6NaHw3dzeDu7bjvF6xlB6t4hOHhi2vLjjIn/+7nvOHd2XKiG6OXNOEH1t3wZggdkqfjkBNP4u3Cg+Vs2Lbfs4bmsqQ5GgSW8fwzir/m8OOVlbzi7lf0b5NHA9clOH39Uz4ssRiTBBLTmjFoC6JPiWW99fsRhXOHdKFmCjhnMGpvLdmNxVV1X7F8pf3N7B21yEevngYHW0ipGlEyCYWERkiInNF5G8iconb8RgTKOP6JrN8yz7KK6u8ev3C3F2kJbVhcNdEACZldKH4cAXLNu9rdgxfbtvPX7M3cunoNM4Zktrs65jI4EpiEZGZIrJHRFbVOT5JRNaJyEYRua2Jy0wGHlfVm4DvByxYY1w2vm8nyiur+XJbcZOvLSmv5JONezlvaJdvRmud0b8zbWKjm90cduRoFbfOXUnX9m24e0rQrzlrgoBbNZZZwKTaB0QkGniSmoQxBJjmqZUME5E36zxSqFkI83IReQRIbuH4jWkxY3p3JEq862f5cF0hR6uqObdWraJNXDRZAzuzMHcX1dW+z8R/ZOE68opK+f0lw2nXOtbn803kcSWxqOpH1Kw/VtsYYKOq5qnqUWAOcJGq5qjqBXUeezyPHwG3AcG/Up8xzdS+TSwZ3dt7lVgW5u6iY3wcmb06Hnd8UkYX9hwq58vt+3269+qCg8z8ZDPfH5fOqf06+XSuiVzBtBNkd2B7ref5wNiGXiwivYA7gHjgkUZeNx2YDpCamkp2drb/kbqspKQkLN6H08K5XNJij7JwSwUL319Eq+j6JyRWVivv5R4mMzWGxR99CPyvTOIqlRiBZ95ezqFBrby+72NflNE2Bsa2LQyrsg3nz0pzOVkmwZRY6vtpabDerqpb8CSMxqjq0yKyE5iSmJg4Oisrq9kBBovs7GzC4X04LZzLRboVsmDmMtr2zOD0/p3rfc1H6ws5UrmMq88eSZanKax2mZy+fRmr9pRw5plnejVb/uv8Yr585xNunTiA88/u79h7CQbh/FlpLifLJJhGheUDPWo9TwMKnLiwTZA0oS4zPYmYKOHTRprDFubuom1cNKf1r7/JanJGV/L3HyG3wLtFyB99bz0d2sZyzam9mhOyiWDBlFiWA/1FpLeIxFGzodg8Jy5s+7GYUBffKoaRPTo02M9SXa28t3o3Zw7o3ODSLecMSSU6SrwaHfbFtv0sWlfIjWf0JdE67I2P3BpuPBtYAgwUkXwRuU5VK4FbgIXAGmCuqua6EZ8xwWh832S+zi/mYFnFCd9bmV/MnkPlnDu04TkmHePjGNu7o1eLUj763nqS4+P4/rh0v2I2kcmtUWHTVLWrqsaqapqqPuc5vkBVB6hqX1V9yMH7WVOYCXmn9E2mWmF5PRMdF+buJiZKOGtg45MXJ2V0YeOeEjbuOdTga5Zt3sfiDUXclNU3IPu4mPAXTE1hAWNNYSYcjOqZRFxMVL3NYe+u3sUpfZJp37bxZqtzh3QBaLQ57NH31tM5sRVXjrXaimmeiEgsVmMx4aB1bDSZ6UkndOBv3FNCXmFpo81gx3Rp35pRPTs02Bz26aYiluTt5easvrSJc2aZfRN5IiKxWI3FhItxfZJZvfMg+0uPfnNsoSdJTPRyDa9JGV1YteMg2/cdPu64qvLoe+vp0q4108b0dC5oE3EiIrFYjcWEi/H9alYv+mzz/2ot767ezYi09nRt38ara0waWrNN8cI6tZbFG4pYvmU/Pzqrn2ObgpnIFBGJxZhwMTytA23jor9pDtt1oIyV24s5d2gXr6/RM7ktQ7q2O66fRVX503vr6d6hDZfZHvbGTxGRWKwpzISL2OgoTu7V8ZvE8t6a3QDHLTrpjUkZXVixbT97DpYBkL2ukK+2F/Pjs/rRKsZqK8Y/EZFYrCnMhJPxfZPZuKeEPYfKeDd3F306xdMvJcGna0zO6IIqLFy9+5vaSs+ObfnOaKutGP9FRGIxJpyM61vTz7IwdzdLNu1l4tBUr9b+qq1fSgJ9OsezcNUu3lu9m5wdB/jJ2f2JjbZfCcZ/NvvJmBAztFt7ElvH8Nh/N1BZrd/MTfGFiDBpaBf+/lEeOw8coXeneKaO7BaAaE0kiog/T6yPxYST6CjhlD7JFJWU0zmxFSf16NCs60zO6EpVtbKpsJSfnt2fGKutGIdExCfJ+lhMuBnXp6Y5bOKQVKKifGsGOyajezvSktrQLyWBKSOstmKcY01hxoSgswal8Id31/Htk7o3+xoiwgvXjqFVTBTRzUxOxtTHEosxIahXp3hy7zvP5077uvp29m00mTHeiIimMOtjMeHI36RiTKBERGKxPhZjjGk5EZFYjDHGtBxLLMYYYxxlicUYY4yjLLEYY4xxlCUWY4wxjrLEYowxxlGiqm7H0GJEpBDY6nYcDugEFLkdRBCycjmRlUn9rFxO1FiZpKtqZ28vFFGJJVyIyOeqmul2HMHGyuVEVib1s3I5kZNlYk1hxhhjHGWJxRhjjKMssYSmp90OIEhZuZzIyqR+Vi4ncqxMrI/FGGOMo6zGYowxxlGWWIwxxjjKEosxxhhHWWIJMyIyRETmisjfROQSt+MJBiJyuog8JSLPisinbscTLEQkS0QWe8omy+14goWIDPaUySsicpPb8QQDEekjIs+JyCvevN4SSxARkZkiskdEVtU5PklE1onIRhG5rYnLTAYeV9WbgO8HLNgW4kSZqOpiVf0h8CbwQiDjbSkOfVYUKAFaA/mBirUlOfR5WeP5vFwGhPwkSofKJE9Vr/P6njYqLHiIyBnU/KC/qKoZnmPRwHpgIjU//MuBaUA0MKPOJa71/HsPcBgYr6qntkDoAeNEmajqHs95c4HrVfVgC4UfMA59VopUtVpEUoE/qeqVLRV/oDj1eRGRC4HbgCdU9aWWij8QHP4ZekVVm2wJiXEufOMvVf1IRHrVOTwG2KiqeQAiMge4SFVnABc0cKkfeT44rwYq1pbiVJmISE/gQDgkFXD0swKwH2gViDhbmlPloqrzgHki8hYQ0onF4c+KVyyxBL/uwPZaz/OBsQ292PMBugOIBx4JZGAu8qlMPK4Dng9YRMHB18/KxcB5QAfgicCG5ipfyyULuJiaZLsgoJG5x9cySQYeAk4Skds9CahBlliCn9RzrMH2S1XdAkwPWDTBwacyAVDVewIUSzDx9bPyKmFQq/WCr+WSDWQHKpgg4WuZ7AV+6O3FrfM++OUDPWo9TwMKXIolWFiZ1M/KpX5WLicKaJlYYgl+y4H+ItJbROKAy4F5LsfkNiuT+lm51M/K5UQBLRNLLEFERGYDS4CBIpIvItepaiVwC7AQWAPMVdVcN+NsSVYm9bNyqZ+Vy4ncKBMbbmyMMcZRVmMxxhjjKEssxhhjHGWJxRhjjKMssRhjjHGUJRZjjDGOssRijDHGUZZYjDHGOMoSizHGGEdZYjHGJZ6Nlr7yPD4TEft5NGHBZt4b4xIR2QCcrqq73I7FGCfZX0jGuGcBkCMifwYQkaki8oyIvCEi57ocmzHNZvuxGOMCERlPzZ4YXT0LAqKqrwOvi0gS8AfgXRdDNKbZrMZijDsuBdaraqXUaFfre3cBT7oUlzF+sz4WY1wgImOA56jZte8IcDPwBfAw8J6q/tfF8IzxiyUWY4KEiPwEuJqaTZi+UtWnXA7JmGaxxGKMMcZR1sdijDHGUZZYjDHGOMoSizHGGEdZYjHGGOMoSyzGGGMcZYnFGGOMoyyxGGOMcZQlFmOMMY6yxGKMMcZR/w/4F/r3h87kSQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.5 Hessian Finite Difference using Exact Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.5-Hessian-Finite-Difference-using-Exact-Gradient)", "section": "3.8.1.5 Hessian Finite Difference using Exact Gradient" } }, "source": [ "### 3.8.1.5 Hessian Finite Difference using Exact Gradient\n", "\n", "Repeat the analysis from above using `my_grad_exact`. What value of $\\epsilon_2$ gives the lower Hessian approximation error?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.5 Hessian Finite Difference using Exact Gradient](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.5-Hessian-Finite-Difference-using-Exact-Gradient)", "section": "3.8.1.5 Hessian Finite Difference using Exact Gradient" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.6 Final Answers](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.6-Final-Answers)", "section": "3.8.1.6 Final Answers" } }, "source": [ "### 3.8.1.6 Final Answers\n", "\n", "**Record your final answers below**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.6 Final Answers](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.6-Final-Answers)", "section": "3.8.1.6 Final Answers" } }, "source": [ "A. Using $\\epsilon_1 = $ ... gives error $|| \\nabla f(x;\\epsilon_1)_{approx} - \\nabla f(x)_{exact} || = $ ...\n", "\n", "B. Using $\\epsilon_1 = $ ... and $\\epsilon_2 = $ ... gives error $|| \\nabla^2 f(x;\\epsilon_1,\\epsilon_2)_{approx} - \\nabla^2 f(x)_{exact} || = $ ...\n", "\n", "C. Using $\\epsilon_2 = $ gives error $|| \\nabla^2 f(x;\\epsilon_2)_{approx} - \\nabla^2 f(x)_{exact} || = $ ...\n", "\n", "These answers were computed using the ... norm." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.7-Discussion)", "section": "3.8.1.7 Discussion" } }, "source": [ "### 3.8.1.7 Discussion\n", "\n", "What is the benefit of using the *exact gradient* when approximating the Hessian with central finite difference?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.1.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.1.7-Discussion)", "section": "3.8.1.7 Discussion" } }, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.8.2 Analysis of possible optimization solutions](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2-Analysis-of-possible-optimization-solutions)", "section": "3.8.2 Analysis of possible optimization solutions" } }, "source": [ "## 3.8.2 Analysis of possible optimization solutions\n", "\n", "Consider the following optimization problem:\n", "\n", "$$\\min f(x) := \\left[x_1^2 + (x_2 + 1)^2\\right] \\cdot \\left[x_1^2 + (x_2 - 1)^2\\right]$$\n", "\n", "Our optimization algorithm terminates at the following points:\n", "1. $x = [0,0]^T$\n", "2. $x = [0,1]^T$\n", "3. $x = [0,-1]^T$\n", "4. $x = [1,1]^T$\n", "\n", "Classify each point.\n", "\n", "You may solve this problem entirely on paper, entirely in Python, or some combination. Please record your answer below. (If you solve on paper, you can typeset the justification in 1 or 2 sentences.)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.8.2 Analysis of possible optimization solutions](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2-Analysis-of-possible-optimization-solutions)", "section": "3.8.2 Analysis of possible optimization solutions" } }, "source": [ "**Suggested solution approach: define systematic analysis routine**\n", "\n", "Create a function that:\n", "1. Evaluates $f(x)$, $\\nabla f(x)$, and $\\nabla^2 f(x)$ for a given $x$\n", "2. Calculates eigenvalues of $\\nabla^2 f(x)$\n", "\n", "We then reuse this function to analyze each point." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 2, "link": "[3.8.2 Analysis of possible optimization solutions](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2-Analysis-of-possible-optimization-solutions)", "section": "3.8.2 Analysis of possible optimization solutions" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.1 Point 1](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.1-Point-1)", "section": "3.8.2.1 Point 1" } }, "source": [ "### 3.8.2.1 Point 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.1 Point 1](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.1-Point-1)", "section": "3.8.2.1 Point 1" }, "scrolled": true }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.1 Point 1](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.1-Point-1)", "section": "3.8.2.1 Point 1" } }, "source": [ "**Answer**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.2 Point 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.2-Point-2)", "section": "3.8.2.2 Point 2" } }, "source": [ "### 3.8.2.2 Point 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.2 Point 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.2-Point-2)", "section": "3.8.2.2 Point 2" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.2 Point 2](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.2-Point-2)", "section": "3.8.2.2 Point 2" } }, "source": [ "**Answer**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.3 Point 3](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.3-Point-3)", "section": "3.8.2.3 Point 3" } }, "source": [ "### 3.8.2.3 Point 3" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.3 Point 3](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.3-Point-3)", "section": "3.8.2.3 Point 3" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.3 Point 3](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.3-Point-3)", "section": "3.8.2.3 Point 3" } }, "source": [ "**Answer**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.4 Point 4](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.4-Point-4)", "section": "3.8.2.4 Point 4" } }, "source": [ "### 3.8.2.4 Point 4" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.4 Point 4](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.4-Point-4)", "section": "3.8.2.4 Point 4" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.4 Point 4](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.4-Point-4)", "section": "3.8.2.4 Point 4" } }, "source": [ "**Answer**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.5 Visualize in 3D](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.5-Visualize-in-3D)", "section": "3.8.2.5 Visualize in 3D" } }, "source": [ "### 3.8.2.5 Visualize in 3D\n", "\n", "Plot $f(x)$ in 3D over the domain $x_1 \\in [-1.1,1.1]$ and $x_2 \\in [-1.5, 1.5]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.5 Visualize in 3D](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.5-Visualize-in-3D)", "section": "3.8.2.5 Visualize in 3D" } }, "outputs": [], "source": [ "# YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.2.6 Convexity](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.2.6-Convexity)", "section": "3.8.2.6 Convexity" } }, "source": [ "### 3.8.2.6 Convexity\n", "\n", "\n", "Is $f(x)$ convex? Did you need to make the plot to make this determination? Write a sentence or two to justify your answer.\n", "\n", "**Answer**:" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.8.3 Multivariable Taylor Series](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3-Multivariable-Taylor-Series)", "section": "3.8.3 Multivariable Taylor Series" } }, "source": [ "## 3.8.3 Multivariable Taylor Series\n", "\n", "You will use `my_grad_exact`, `my_grad_approx`, `my_hes_exact`, and `my_hes_approx` to construct Taylor series approximations to an arbitrary twice differentiable continuous functions with inputs $x \\in \\mathbb{R}^2$. We will then consider Example 2.19 and visualize the Taylor series approximation in 3D.\n" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.1 Create a function to plot the first order Taylor series using $\\nabla f$](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.1-Create-a-function-to-plot-the-first-order-Taylor-series-using-$\\nabla-f$)", "section": "3.8.3.1 Create a function to plot the first order Taylor series using $\\nabla f$" } }, "source": [ "### 3.8.3.1 Create a function to plot the first order Taylor series using $\\nabla f$\n", "\n", "Create a general function that:\n", "1. Constructs a Taylor series using $\\nabla f$ centered around a given point\n", "2. Plots the true function and Taylor series approximation\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.1 Create a function to plot the first order Taylor series using $\\nabla f$](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.1-Create-a-function-to-plot-the-first-order-Taylor-series-using-$\\nabla-f$)", "section": "3.8.3.1 Create a function to plot the first order Taylor series using $\\nabla f$" } }, "outputs": [], "source": [ "def taylor1(xc, f, grad, dx):\n", " '''\n", " Constructs a Taylor series using first derivates and visualizes in 3D\n", " \n", " Arguments:\n", " xc - point to center Taylor series\n", " f - function that computes function value. Only has one input (x)\n", " grad - function that computes gradient. Only has one input (x)\n", " dx - list or numpy array. creates 3D plot over xc[0] +/- dx[0], xc[1] +/- dx[1]\n", " \n", " Returns:\n", " none\n", " \n", " Actions:\n", " 3D plot\n", " \n", " '''\n", " \n", " # YOUR SOLUTION HERE" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.2 Taylor Series using `my_grad_approx`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.2-Taylor-Series-using-`my_grad_approx`)", "section": "3.8.3.2 Taylor Series using `my_grad_approx`" } }, "source": [ "### 3.8.3.2 Taylor Series using `my_grad_approx`\n", "\n", "Consider $x_c = [0.7, 0.3]^T$ (center of Taylor series) and $\\Delta x = [0.5, 0.5]^T$ (domain for plot)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.2 Taylor Series using `my_grad_approx`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.2-Taylor-Series-using-`my_grad_approx`)", "section": "3.8.3.2 Taylor Series using `my_grad_approx`" } }, "outputs": [], "source": [ "# Specify epsilon1\n", "calc_grad = lambda x : my_grad_approx(x,my_f,1E-6)\n", "\n", "# Specify dx\n", "dx = [0.5, 0.5]\n", "\n", "# Specify xc\n", "xc = np.array([0.7, 0.3])\n", "\n", "taylor1(xc, my_f, calc_grad, dx)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.3 Taylor Series using `my_grad_exact`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.3-Taylor-Series-using-`my_grad_exact`)", "section": "3.8.3.3 Taylor Series using `my_grad_exact`" } }, "source": [ "### 3.8.3.3 Taylor Series using `my_grad_exact`\n", "\n", "Consider $x_c = [0.7, 0.3]^T$ (center of Taylor series) and $\\Delta x = [0.5, 0.5]^T$ (domain for plot)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.3 Taylor Series using `my_grad_exact`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.3-Taylor-Series-using-`my_grad_exact`)", "section": "3.8.3.3 Taylor Series using `my_grad_exact`" } }, "outputs": [], "source": [ "# Specify epsilon1\n", "calc_grad = lambda x : my_grad_exact(x)\n", "\n", "# Specify dx\n", "dx = [0.5, 0.5]\n", "\n", "# Specify xc\n", "xc = np.array([0.7, 0.3])\n", "\n", "taylor1(xc, my_f, calc_grad, dx)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.4 Create a function to plot the second order Taylor series using $\\nabla f$ and $\\nabla^2 f$](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.4-Create-a-function-to-plot-the-second-order-Taylor-series-using-$\\nabla-f$-and-$\\nabla^2-f$)", "section": "3.8.3.4 Create a function to plot the second order Taylor series using $\\nabla f$ and $\\nabla^2 f$" } }, "source": [ "### 3.8.3.4 Create a function to plot the second order Taylor series using $\\nabla f$ and $\\nabla^2 f$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.4 Create a function to plot the second order Taylor series using $\\nabla f$ and $\\nabla^2 f$](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.4-Create-a-function-to-plot-the-second-order-Taylor-series-using-$\\nabla-f$-and-$\\nabla^2-f$)", "section": "3.8.3.4 Create a function to plot the second order Taylor series using $\\nabla f$ and $\\nabla^2 f$" } }, "outputs": [], "source": [ "def taylor2(xc, f, grad, hes, dx):\n", " '''\n", " Constructs a Taylor series using first derivates and visualizes in 3D\n", " \n", " Inputs:\n", " xc - point to center Taylor series\n", " f - computes function value. Only has one input (x)\n", " grad - computes gradient. Only has one input (x)\n", " hes - computes the Hessian. Only has one input (x)\n", " dx - creates 3D plot over xc[0] +/- dx[0], xc[1] +/- dx[1]\n", " \n", " Outputs:\n", " none\n", " \n", " Creates:\n", " 3D plot\n", " \n", " '''\n", " \n", " ### Evaluates function and gradiant\n", " fval = f(xc)\n", " gval = grad(xc)\n", " Hval = hes(xc)\n", " \n", " ### Creates domain for plotting\n", " x1 = np.arange(xc[0] - dx[0],xc[0] + dx[0],dx[0]/100)\n", " x2 = np.arange(xc[1] - dx[1],xc[1] + dx[1],dx[1]/100)\n", " \n", " ## Create a matrix of all points to sample\n", " X1, X2 = np.meshgrid(x1, x2)\n", " n1 = len(x1)\n", " n2 = len(x2)\n", "\n", " ## Allocate matrix for true function value\n", " F = np.zeros([n2, n1])\n", "\n", " ## Allocate matrix for Taylor series approximation\n", " T = np.zeros([n2, n1]) \n", " \n", " xtemp = np.zeros(2)\n", "\n", " # Evaluate f(x) and Taylor series over grid\n", " for i in range(0,n1):\n", " xtemp[0] = x1[i]\n", " for j in range(0,n2):\n", " xtemp[1] = x2[j]\n", " \n", " # Evaluate f(x)\n", " F[j,i] = my_f(xtemp)\n", "\n", " # Evaluate Taylor series\n", " dx_ = xtemp - xc\n", " \n", " '''\n", " print(\"dx = \",dx)\n", " print(\"gval = \",gval)\n", " print(\"Hval = \",Hval)\n", " '''\n", " \n", " temp = Hval.dot(dx_)\n", " # print(\"Hval * dx = \",temp)\n", " \n", " \n", " # T[j,i] = fval + gval.dot(dx_) + 0.5*(temp).dot(dx_)\n", " T[j,i] = fval + gval.dot(dx_) + 0.5*(dx_.dot(Hval.dot(dx_)))\n", " \n", " # Create 3D figure\n", " fig = plt.figure()\n", " ax = fig.gca(projection='3d')\n", "\n", " # Plot f(x)\n", " surf = ax.plot_surface(X1, X2, F, linewidth=0,cmap=cm.coolwarm,antialiased=True,label=\"f(x)\")\n", "\n", " # Plot Taylor series approximation\n", " surf = ax.plot_surface(X1, X2, T, linewidth=0,cmap=cm.PiYG,antialiased=True,label=\"Taylor series\")\n", " \n", " # Add candidate point\n", " ax.scatter(xc[0],xc[1],fval,s=50,color=\"black\",depthshade=True)\n", "\n", " # Draw vertical line through stationary point to help visualization\n", " # Maximum value in array\n", " fmax = np.amax(F)\n", " fmin = np.amin(F)\n", " ax.plot([xc[0], xc[0]], [xc[1], xc[1]], [fmin,fmax],color=\"black\")\n", " \n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.5 Taylor series using `my_grad_approx` and `my_hes_approx`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.5-Taylor-series-using-`my_grad_approx`-and-`my_hes_approx`)", "section": "3.8.3.5 Taylor series using `my_grad_approx` and `my_hes_approx`" } }, "source": [ "### 3.8.3.5 Taylor series using `my_grad_approx` and `my_hes_approx`\n", "\n", "Consider $x_c = [0.7, 0.3]^T$ (center of Taylor series) and $\\Delta x = [0.2, 0.2]^T$ (domain for plot)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.5 Taylor series using `my_grad_approx` and `my_hes_approx`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.5-Taylor-series-using-`my_grad_approx`-and-`my_hes_approx`)", "section": "3.8.3.5 Taylor series using `my_grad_approx` and `my_hes_approx`" } }, "outputs": [], "source": [ "# Specify epsilon1\n", "calc_grad = lambda x : my_grad_approx(x,my_f,1E-6)\n", "\n", "# Specify epsilon2\n", "calc_hes = lambda x : my_hes_approx(x, calc_grad, 1E-6)\n", "\n", "# Specify dx\n", "dx = [0.2, 0.2]\n", "\n", "# Specify xc\n", "xc = np.array([0.7, 0.3])\n", "\n", "taylor2(xc, my_f, calc_grad, calc_hes, dx)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.6 Taylor series using `my_grad_exact` and `my_hes_exact`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.6-Taylor-series-using-`my_grad_exact`-and-`my_hes_exact`)", "section": "3.8.3.6 Taylor series using `my_grad_exact` and `my_hes_exact`" } }, "source": [ "### 3.8.3.6 Taylor series using `my_grad_exact` and `my_hes_exact`\n", "\n", "Consider $x_c = [0.7, 0.3]^T$ (center of Taylor series) and $\\Delta x = [0.2, 0.2]^T$ (domain for plot)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.6 Taylor series using `my_grad_exact` and `my_hes_exact`](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.6-Taylor-series-using-`my_grad_exact`-and-`my_hes_exact`)", "section": "3.8.3.6 Taylor series using `my_grad_exact` and `my_hes_exact`" } }, "outputs": [], "source": [ "x = np.array([0,0])\n", "# Specify epsilon1\n", "calc_grad = lambda x : my_grad_exact(x)\n", "\n", "# Specify epsilon2\n", "calc_hes = lambda x1 : my_hes_exact(x1)\n", "\n", "# Specify dx\n", "dx = [0.2, 0.2]\n", "\n", "# Specify xc\n", "xc = np.array([0.7, 0.3])\n", "\n", "taylor2(xc, my_f, calc_grad, calc_hes, dx)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.7-Discussion)", "section": "3.8.3.7 Discussion" } }, "source": [ "### 3.8.3.7 Discussion" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.7-Discussion)", "section": "3.8.3.7 Discussion" } }, "source": [ "Write 1 or 2 sentences to describe the the shapes for the first order and second order Taylor series approximations? Why do these shapes make sense?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.7-Discussion)", "section": "3.8.3.7 Discussion" } }, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.7-Discussion)", "section": "3.8.3.7 Discussion" } }, "source": [ "Is there a visible difference in the Taylor series approximations using the finite difference versus exact derivatives? Why does this make sense?" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[3.8.3.7 Discussion](https://ndcbe.github.io/CBE60499/03.08-Algorithms2.html#3.8.3.7-Discussion)", "section": "3.8.3.7 Discussion" } }, "source": [ "**Answer**:" ] }, { "cell_type": "markdown", "id": "f96fba79", "metadata": {}, "source": [ "\n", "< [3.7 Algorithms Homework 1](https://ndcbe.github.io/CBE60499/03.07-Algorithms1.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [3.9 Algorithms Homework 3](https://ndcbe.github.io/CBE60499/03.09-Algorithms3.html) >

\"Open

\"Download\"" ] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }