{ "cells": [ { "cell_type": "markdown", "id": "e235223c", "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": "682be5ed", "metadata": {}, "source": [ "\n", "< [4.6 NLP Diagnostics with Degeneracy Hunter](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.8 Inertia-Corrected Netwon Method for Equality Constrained NLPs](https://ndcbe.github.io/CBE60499/04.08-Interior-Point2.html) >
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.7 Simple Netwon Method for Equality Constrained NLPs](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7-Simple-Netwon-Method-for-Equality-Constrained-NLPs)",
"section": "4.7 Simple Netwon Method for Equality Constrained NLPs"
}
},
"source": [
"# 4.7 Simple Netwon Method for Equality Constrained NLPs"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.7 Simple Netwon Method for Equality Constrained NLPs](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7-Simple-Netwon-Method-for-Equality-Constrained-NLPs)",
"section": "4.7 Simple Netwon Method for Equality Constrained NLPs"
}
},
"source": [
"**Reference**: Section 5.2 in Biegler (2010)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.7 Simple Netwon Method for Equality Constrained NLPs](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7-Simple-Netwon-Method-for-Equality-Constrained-NLPs)",
"section": "4.7 Simple Netwon Method for Equality Constrained NLPs"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.1 Helper Functions](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.1-Helper-Functions)",
"section": "4.7.1 Helper Functions"
}
},
"source": [
"## 4.7.1 Helper Functions\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.1 Helper Functions](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.1-Helper-Functions)",
"section": "4.7.1 Helper Functions"
}
},
"outputs": [],
"source": [
"# Load required Python libraries.\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy import linalg\n",
"\n",
"## Check is element of array is NaN\n",
"def check_nan(A):\n",
" return np.sum(np.isnan(A))\n",
"\n",
"## Calculate gradient with central finite difference\n",
"## Calculate gradient with central finite difference\n",
"def my_grad_approx(x,f,eps1,verbose=False):\n",
" '''\n",
" Calculate gradient of function 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",
"def my_jac_approx(x,h,eps1,verbose=False):\n",
" '''\n",
" Calculate Jacobian of function h(x) using central difference formula\n",
" \n",
" Inputs:\n",
" x - point for which to evaluate gradient\n",
" h - vector-valued function to consider. h(x): R^n --> R^m\n",
" eps1 - perturbation size\n",
" \n",
" Outputs:\n",
" A - Jacobian (n x m matrix)\n",
" '''\n",
" \n",
" # Check h(x) at x\n",
" h_x0 = h(x)\n",
" \n",
" # Extract dimensions\n",
" n = len(x)\n",
" m = len(h_x0)\n",
" \n",
" # Initialize Jacobian matrix\n",
" A = np.zeros((n,m))\n",
" \n",
" # Calculate Jacobian by row\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_h_plus = h(x + e)\n",
" my_h_minus = h(x - e)\n",
" \n",
" # Diagnostics\n",
" if(verbose):\n",
" print(\"e[\",i,\"] = \",e)\n",
" print(\"h(x + e[\",i,\"]) = \",my_h_plus)\n",
" print(\"h(x - e[\",i,\"]) = \",my_h_minus)\n",
" \n",
" \n",
" A[i,:] = (my_h_plus - my_h_minus)/(2*eps1)\n",
" \n",
" if(verbose):\n",
" print(\"***** Done. ***** \\n\")\n",
" \n",
" return A\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\n",
"\n",
"## Linear algebra calculation\n",
"def xxT(u):\n",
" '''\n",
" Calculates u*u.T to circumvent limitation with SciPy\n",
" \n",
" Arguments:\n",
" u - numpy 1D array\n",
" \n",
" Returns:\n",
" u*u.T\n",
" \n",
" Assume u is a nx1 vector.\n",
" Recall: NumPy does not distinguish between row or column vectors\n",
" \n",
" u.dot(u) returns a scalar. This functon returns an nxn matrix.\n",
" '''\n",
" \n",
" n = len(u)\n",
" A = np.zeros([n,n])\n",
" for i in range(0,n):\n",
" for j in range(0,n):\n",
" A[i,j] = u[i]*u[j]\n",
" \n",
" return A\n",
"\n",
"## Analyze Hessian\n",
"def analyze_hes(B):\n",
" print(B,\"\\n\")\n",
" \n",
" l = linalg.eigvals(B)\n",
" print(\"Eigenvalues: \",l,\"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.2 Algorithm 5.1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.2-Algorithm-5.1)",
"section": "4.7.2 Algorithm 5.1"
}
},
"source": [
"## 4.7.2 Algorithm 5.1"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.2 Algorithm 5.1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.2-Algorithm-5.1)",
"section": "4.7.2 Algorithm 5.1"
}
},
"outputs": [],
"source": [
"def alg51(x0,calc_f,calc_h,eps1=1E-6,eps2=1E-6,max_iter=50,verbose=False):\n",
" '''\n",
" Basic Full Space Newton Method for Equality Constrained NLP\n",
" \n",
" Input:\n",
" x0 - starting point (vector)\n",
" calc_f - function to calculate objective (returns scalar)\n",
" calc_h - function to calculate constraints (returns vector)\n",
" eps1 - tolerance for primal and dual steps\n",
" eps2 - tolerance for gradient of L1\n",
" \n",
" Outputs:\n",
" x - history of steps (primal variables)\n",
" v - history of steps (dual variables)\n",
" f - history of objective evaluations\n",
" h - history of constraint evaluations\n",
" df - history of objective gradients\n",
" dL - history of Lagrange function gradients\n",
" A - history of constraint Jacobians\n",
" W - history of Lagrange Hessians\n",
" \n",
" Notes:\n",
" 1. For simplicity, central finite difference is used \n",
" for all gradient calculations.\n",
" '''\n",
" \n",
" # Declare iteration histories as empty lists\n",
" x = []\n",
" v = []\n",
" f = []\n",
" L = []\n",
" h = []\n",
" df = []\n",
" dL = []\n",
" A = []\n",
" W = []\n",
" \n",
" # Flag for iterations\n",
" flag = True\n",
" \n",
" # Iteration counter\n",
" k = 0\n",
" \n",
" # Copy initial point to primal variable history\n",
" n = len(x0)\n",
" x.append(x0)\n",
" \n",
" # Evaluate objective and constraints at initial point\n",
" f.append(calc_f(x0))\n",
" h.append(calc_h(x0))\n",
" \n",
" # Determine number of equality constraints\n",
" m = len(h[0])\n",
" \n",
" # Initial dual variables with vector of ones\n",
" v.append(np.ones(m))\n",
" \n",
" # Print header for iteration information\n",
" print(\"Iter. \\tf(x) \\t\\t||h(x)|| \\t||grad_L(x)|| \\t||dx|| \\t\\t||dv||\")\n",
" \n",
" while(flag and k < max_iter):\n",
" \n",
" # STEP 1. Construct KKT matrix\n",
" \n",
" if(k > 0):\n",
" # Evaluate objective function\n",
" f.append(calc_f(x[k]))\n",
" \n",
" # Evaluate constraint function\n",
" h.append(calc_h(x[k]))\n",
" \n",
" # Evaluate objective gradient\n",
" df.append(my_grad_approx(x[k],calc_f,1E-6))\n",
" \n",
" # Evaluate constraint Jacobian\n",
" A.append(my_jac_approx(x[k],calc_h,1E-6))\n",
" \n",
" # Evaluate gradient of Lagrange function\n",
" L_func = lambda x_ : calc_f(x_) + (calc_h(x_)).dot(v[k])\n",
" L_grad = lambda x_ : my_grad_approx(x_,L_func,1E-6)\n",
" dL.append(L_grad(x[k]))\n",
" norm_dL = linalg.norm(dL[k])\n",
" \n",
" # Evaluate Hessian of Lagrange function\n",
" W.append(my_hes_approx(x[k],L_grad,1E-6))\n",
" \n",
" if(verbose):\n",
" print(\"*** k =\",k,\" ***\")\n",
" print(\"x_k =\",x[k])\n",
" print(\"v_k =\",v[k])\n",
" print(\"f_k =\",f[k])\n",
" print(\"df_k =\",df[k])\n",
" print(\"h_k =\",h[k])\n",
" print(\"A_k =\\n\",A[k])\n",
" print(\"W_k =\\n\",W[k])\n",
" print(\"\\n\")\n",
" \n",
" # Assemble KKT matrix\n",
" KKT_top = np.concatenate((W[k],A[k]),axis=1)\n",
" KKT_bot = np.concatenate((A[k].T,np.zeros((m,m))),axis=1)\n",
" KKT = np.concatenate((KKT_top,KKT_bot),axis=0)\n",
" \n",
" if(verbose):\n",
" print(\"KKT matrix =\\n\",KKT,\"\\n\")\n",
" \n",
" # Check if KKT matrix is singular\n",
" l, eigvec = linalg.eig(KKT)\n",
" \n",
" if(verbose):\n",
" print(\"KKT matrix eigenvalues:\")\n",
" print(l)\n",
" \n",
" zero_eigenvalues = sum(np.abs(l) <= 1E-8)\n",
" \n",
" if(zero_eigenvalues > 0):\n",
" flag = False\n",
" print(\"KKT matrix is singular. Eigenvalues:\\n\")\n",
" print(l,\"\\n\")\n",
" \n",
" ## STEP 2. Solve linear system.\n",
" \n",
" if(flag):\n",
" b = -np.concatenate((dL[k],h[k]),axis=0)\n",
" z = linalg.solve(KKT,b)\n",
" else:\n",
" z = []\n",
" \n",
" ## STEP 3. Take step\n",
" if(flag):\n",
" dx = z[0:n]\n",
" dv = z[n:(n+m+1)]\n",
" \n",
" x.append(x[k] + dx)\n",
" v.append(v[k] + dv)\n",
" \n",
" norm_dx = linalg.norm(dx)\n",
" norm_dv = linalg.norm(dv)\n",
" \n",
" ## Print iteration information\n",
" print(k,' \\t{0: 1.4e} \\t{1:1.4e} \\t{2:1.4e}'.format(f[k],linalg.norm(h[k]),norm_dL),end='')\n",
" \n",
" if flag:\n",
" print(' \\t{0: 1.4e} \\t{1: 1.4e}'.format(norm_dx,norm_dv),end='\\n')\n",
" else:\n",
" print(' \\t ------- \\t -------',end='\\n')\n",
" \n",
" # Increment counter\n",
" k = k + 1\n",
" \n",
" ## Check convergence criteria\n",
" if(flag):\n",
" flag = norm_dx > eps1 and norm_dv > eps1 and norm_dL > eps2\n",
" \n",
" if (not flag) and k >= max_iter:\n",
" print(\"Reached maximum number of iterations.\")\n",
" \n",
" return x,v,f,h,df,dL,A,W"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.3 Example Problem 1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3-Example-Problem-1)",
"section": "4.7.3 Example Problem 1"
}
},
"source": [
"## 4.7.3 Example Problem 1\n",
"\n",
"Consider:\n",
"\n",
"$$\\begin{align*} \\min_{x} \\quad & x_1^2 + 2 x_2^2 - x_3 \\\\\n",
"\\mathrm{s.t.} \\quad & x_1 + x_2 = 1 \\\\\n",
"& x_1 + x_2 - x_3 = 5\n",
"\\end{align*} $$\n",
"\n",
"### Define Functions"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.3 Example Problem 1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3-Example-Problem-1)",
"section": "4.7.3 Example Problem 1"
}
},
"outputs": [],
"source": [
"def my_f(x):\n",
" return x[0]**2 + 2*x[1]**2 - 1*x[2]\n",
"\n",
"def my_h(x):\n",
" h = np.zeros(2)\n",
" \n",
" h[0] = x[0] + x[1] - 1\n",
" h[1] = x[0] + x[1] - x[2] - 5\n",
" \n",
" return h"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.3.1 Test Finite Difference Approximations](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3.1-Test-Finite-Difference-Approximations)",
"section": "4.7.3.1 Test Finite Difference Approximations"
}
},
"source": [
"### 4.7.3.1 Test Finite Difference Approximations"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.3.1 Test Finite Difference Approximations](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3.1-Test-Finite-Difference-Approximations)",
"section": "4.7.3.1 Test Finite Difference Approximations"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f(x0) = 2 \n",
"\n",
"h(x0) = [ 1. -4.] \n",
"\n",
"df(x0) = [ 2. 4. -1.] \n",
"\n",
"dh(x0) =\n",
" [[ 1. 1.]\n",
" [ 1. 1.]\n",
" [ 0. -1.]] \n",
"\n"
]
}
],
"source": [
"## Define initial point\n",
"x0 = np.array([1,1,1])\n",
"\n",
"## Calculate objective\n",
"print(\"f(x0) =\",my_f(x0),\"\\n\")\n",
"\n",
"## Calculate constraints\n",
"print(\"h(x0) =\",my_h(x0),\"\\n\")\n",
"\n",
"## Calculate objective gradient\n",
"print(\"df(x0) =\",my_grad_approx(x0,my_f,1E-6),\"\\n\")\n",
"\n",
"## Calculate constraint Jacobian\n",
"print(\"dh(x0) =\\n\",my_jac_approx(x0,my_h,1E-6),\"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.3.2 Test Algorithm 5.1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3.2-Test-Algorithm-5.1)",
"section": "4.7.3.2 Test Algorithm 5.1"
}
},
"source": [
"### 4.7.3.2 Test Algorithm 5.1"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.3.2 Test Algorithm 5.1](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.3.2-Test-Algorithm-5.1)",
"section": "4.7.3.2 Test Algorithm 5.1"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter. \tf(x) \t\t||h(x)|| \t||grad_L(x)|| \t||dx|| \t\t||dv||\n",
"0 \t 2.0000e+00 \t4.1231e+00 \t7.4833e+00 \t 5.0553e+00 \t 2.4040e+00\n",
"1 \t 4.6667e+00 \t5.7632e-10 \t8.2956e-04 \t 1.5000e-04 \t 6.2903e-04\n",
"2 \t 4.6667e+00 \t0.0000e+00 \t1.4762e-08 \t 1.9888e-09 \t 8.1407e-09\n",
"\n",
"x* = [ 0.66666667 0.33333333 -4. ]\n",
"\n",
"v* = [-0.33333333 -1. ]\n"
]
}
],
"source": [
"## Run Algorithm 5.1 on test problem\n",
"results = alg51(x0,my_f,my_h)\n",
"\n",
"## Display results\n",
"xstar = results[0][-1]\n",
"print(\"\\nx* =\",xstar)\n",
"\n",
"## Display results\n",
"vstar = results[1][-1]\n",
"print(\"\\nv* =\",vstar)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.4 Example Problem 2](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.4-Example-Problem-2)",
"section": "4.7.4 Example Problem 2"
}
},
"source": [
"## 4.7.4 Example Problem 2\n",
"\n",
"**Can we break Algorithm 5.1?**\n",
"\n",
"Probably. Let us try a model where $\\nabla h(x^k)^T$ is always rank-deficient.\n",
"\n",
"Consider:\n",
"$$\\begin{align}\\min_x \\quad & x_1^2 + 2 x_2^2 \\\\\n",
"\\mathrm{s.t.} \\quad & x_1 + x_2 = 1 \\\\\n",
" & x_1 + x_2 = 1 \\end{align}$$\n",
" \n",
"### Test Algorithm 5.1 with the redundant constraint"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.4 Example Problem 2](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.4-Example-Problem-2)",
"section": "4.7.4 Example Problem 2"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter. \tf(x) \t\t||h(x)|| \t||grad_L(x)|| \t||dx|| \t\t||dv||\n",
"KKT matrix is singular. Eigenvalues:\n",
"\n",
"[-1.05145076+0.j 2.51704031+0.j 4.53361158+0.j 0. +0.j] \n",
"\n",
"0 \t 3.0000e+00 \t1.4142e+00 \t7.2111e+00 \t ------- \t -------\n",
"\n",
"x* = [1 1]\n",
"\n",
"v* = [1. 1.]\n"
]
}
],
"source": [
"def my_f2(x):\n",
" return x[0]**2 + 2*x[1]**2\n",
"\n",
"def my_h2(x):\n",
" h = np.zeros(2)\n",
" h[0] = x[0] + x[1] - 1\n",
" h[1] = h[0]\n",
" return h\n",
"\n",
"x0 = np.array((1,1))\n",
"\n",
"## Run Algorithm 5.1 on test problem\n",
"results = alg51(x0,my_f2,my_h2)\n",
"\n",
"## Display results\n",
"xstar = results[0][-1]\n",
"print(\"\\nx* =\",xstar)\n",
"\n",
"## Display results\n",
"vstar = results[1][-1]\n",
"print(\"\\nv* =\",vstar)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.4.1 Trying Algorithm 5.1 again without the redundant constraint](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.4.1-Trying-Algorithm-5.1-again-without-the-redundant-constraint)",
"section": "4.7.4.1 Trying Algorithm 5.1 again without the redundant constraint"
}
},
"source": [
"### 4.7.4.1 Trying Algorithm 5.1 again without the redundant constraint"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.4.1 Trying Algorithm 5.1 again without the redundant constraint](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.4.1-Trying-Algorithm-5.1-again-without-the-redundant-constraint)",
"section": "4.7.4.1 Trying Algorithm 5.1 again without the redundant constraint"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter. \tf(x) \t\t||h(x)|| \t||grad_L(x)|| \t||dx|| \t\t||dv||\n",
"0 \t 3.0000e+00 \t1.0000e+00 \t5.8310e+00 \t 7.4537e-01 \t 2.3335e+00\n",
"1 \t 6.6667e-01 \t1.3978e-10 \t2.9473e-04 \t 4.5326e-05 \t 1.5285e-04\n",
"2 \t 6.6667e-01 \t0.0000e+00 \t4.5431e-09 \t 1.4393e-09 \t 2.0170e-09\n",
"\n",
"x* = [0.66666667 0.33333333]\n",
"\n",
"v* = [-1.33333333]\n"
]
}
],
"source": [
"def my_h2b(x):\n",
" return (x[0] + x[1] - 1)*np.ones(1)\n",
"\n",
"x0 = np.array((1,1))\n",
"\n",
"## Run Algorithm 5.1 on test problem\n",
"results = alg51(x0,my_f2,my_h2b)\n",
"\n",
"## Display results\n",
"xstar = results[0][-1]\n",
"print(\"\\nx* =\",xstar)\n",
"\n",
"## Display results\n",
"vstar = results[1][-1]\n",
"print(\"\\nv* =\",vstar)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.5 Example Problem 3](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.5-Example-Problem-3)",
"section": "4.7.5 Example Problem 3"
}
},
"source": [
"## 4.7.5 Example Problem 3\n",
"\n",
"**Can we break Algorithm 5.1 in another way?**\n",
"\n",
"Let us try a model where $\\nabla h(x^k)^T$ is full rank but there are multiple local optima.\n",
"\n",
"Consider:\n",
"$$\\begin{align}\\min_x \\quad & x_1^3 - x_2 -x_1 x_2 - x_2^2 \\\\\n",
"\\mathrm{s.t.} \\quad & x_1^2 + x_2^2 = 1\n",
"\\end{align}$$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.7.5 Example Problem 3](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.5-Example-Problem-3)",
"section": "4.7.5 Example Problem 3"
}
},
"outputs": [],
"source": [
"def my_f3(x):\n",
" return x[0]**3 - x[1] - x[0]*x[1] - x[1]**2\n",
"\n",
"def my_h3(x):\n",
" return (x[0]**2 + x[1]**2 - 1)*np.ones(1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.5.1 Visualize](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.5.1-Visualize)",
"section": "4.7.5.1 Visualize"
}
},
"source": [
"### 4.7.5.1 Visualize"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.7.5.1 Visualize](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html#4.7.5.1-Visualize)",
"section": "4.7.5.1 Visualize"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABgfUlEQVR4nO2dd3xUVfr/32cmvfdCQhop9F4E6R1EQJq9oK5ldXV/flfXLa5bbKura9fVtaIroiAgIlUCIr23QAiBdNJ7m8zM+f1xAyIkEJLJ3DvJfb9e9zUzmTv3fHIzuZ97znnO8wgpJTo6Ojo6Oga1Bejo6OjoaAPdEHR0dHR0AN0QdHR0dHQa0Q1BR0dHRwfQDUFHR0dHpxHdEHR0dHR0ABsYghCiqxBikxAiRQhxVAjxaBP7CCHE60KINCHEISHEwLa2q6Ojo6NjW5xscAwz8H9Syn1CCG9grxBivZTy2AX7TAMSGrdhwDuNjzo6Ojo6GqHNPQQpZZ6Ucl/j80ogBYi4aLdZwKdSYQfgJ4QIb2vbOjo6Ojq2wxY9hPMIIWKAAcDOi96KALIueJ3d+LO8Jo5xH3AfgJub26CoqKg2abJYJUUV9QR4ueDsZPspE6vVisGg/akYXadt0XXals6qs8JkptzUQJiHK842Om5qamqRlDK4VR+WUtpkA7yAvcCcJt77Dhh5weuNwKArHTMxMVG2ldp6k3xnyVZ5tqiizcdqik2bNrXLcW2NrtO26DptS2fTmVlcKm/5YLFMevoV+ZvFK2Vpda1NjiullMAe2crruE16CEIIZ2Ap8LmUclkTu2QDXS94HQnk2qLtK+Hm4swD86+1R1M6Ojo6l0VKydL9R3luTTJGYeDFG6Zyfd/uCCHUlgbYYMhIKL/JB0CKlPKVZnZbCTwshFiMMplcLqW8ZLhIR0dHp6OSV17JUyvXs/VUBkNjInlh9hS6+PmoLesX2KKHcC1wO3BYCHGg8Wd/BKIApJTvAquB6UAaUAMstEG7Ojo6OppHSsmy/Ud5fu1mrFbJU9PHcfPgfhgM2ugVXEibDUFKuRW47G/WOK71UFvb0tHR0XEkskvL+cu3G9iWnsmQ6EiemzWJrgF+astqFptGGamNlFIzY3E6OjqdF4vVyqKdB3jth58wCAN/mT6emwb31WSv4EI6hCHUm8y4ujjpZqCjo6M6R3Lz+euqjRzJzWdsYix/vW4CYb7eastqEQ5tCPkllazffoLMs6UE+nkwb0J/Av08z79vtUrNO7KOjk7HoKquntc2bePzXQcJ9HTn5bnTmd470aFuVB3aEF78aCPCIAjx9+JkZhFbD6Qza2wfDqfl0ie+i24GOjo67Y6UklWHT/DS+i0UVlZz85B+/Hb8CHzc3dSWdtU4rCGcyS3hyKk81r79IABb9p1i0Xe72Xc8m+z8MgD+cPdE4ru2bsGejo6OzpU4kV/EP1b/wJ6MHHqFh/DGjdfTL9Jxs/I4rCGcyi4iOjyA8qpafL3cCQ/y4fjpfO6fOwIfT3fe/XorB07k6Iago6Njc0qra3kjeTtf7jmEt5srf79+InMH9MLoAOk3LofDGsLA7pF88f1eFq3ajZ+3Bys3H+aWaYMY3DOq8f2u7DycwbyJ/dUVqnHqG8zklFWQU15BblkF+RVVlNbUUlJTS2lNLdX1JmobzNSaGmiwWLBKiUVKkBInoxFnowEngwF3F2c8XVzwdHXBx80VP3c3/Dzc8fdwJ8Tbs3HzIsTbCxcno9q/to5OqzBbrXy8fR9vJe+gxmTixsF9eWTcCPw8HG94qCkc1hD8fTy4fcYQ1u84QUiANwO6RxJ0wYTyln2nGD80QUWF2qOgsorDOfkczjnLyYJi0gqLySotx6rklwKUBSXnLuR+7m4Eenrg7uKMu7MzLk5GjEIgGjezxYLZasVktlDb0ECNqYGqunpyyiooq6mlvLYOeZEGAYR4exHh50Okvw/RAf5EBfgRE+hHbFCAPU+Hjk6LkVLy/dFUnt+VQmGdiZHdonlyyhjiQwLVlmZTHNYQAMYMimfMoHgA1vyUwt/fW0NxWQ2VNXVUVNdx05TOXYcnp6yCHemZbEvPZNvJdEqT9wPgZDAQE+hH97BgruuTREyAP138vOni60OItxdORtt0ey1WK+W1dRRWVlNQVU1BRRW55ZXklleQU1bBnowcvj10/Bem4e/qTI/MYhJCAkkMCSIxNIiEkCDcnB36q6rjwOw4ncXL63/kcG4+EZ5uvHfrbEYnxKotq11w6P8ys9mCk5MRq1UyaXgSVilZuvEgowbE8bcHpqktz+5IKUk5W8i6YydZn5LGqaISAII8Pejm48mkgf3oGxlGj7AQu1xgjQYDAZ4eBHh6kETTczl1DWaySss4U1zGqcISth89RnltHV/uOUyd2QyAQQjigvzpGR5Cj7AQencJpUd4CF6uLu3+O+h0Xg5m5/Hqxp/YfjqLMB8vnps1Gb+ygg5rBuDghkBjfO/rX2zG1GDhibsmMH1kT5VF2Z/s0nKWHzzGioMpZJWWYxCCIdGRLBjUhxHdoogPDmTz5s2MHa69HpObsxMJIUovYFIP6G6tYezYsVisVrJKy0nNL+J4fiEpeYXsPJ3NykPHAWXoqVtwAL27hNEvMow+EWEkhQbhbNTnJ3TaxpHcfN7evIMfTqQT4OHOk1PGcNPgvrg5O5GcXKi2vGaRUvK33RvadAyHNQQpJU5GA1JKlqw/wCd/vxVQhikcfaa/JVitkuST6Xy+6yDbTmUAMDwuivtHDWVCUjf8Pd1VVtg2jAYDMYH+xAT6M7nnz3NBRVXVHM0t4EhuPodzz7Ll5GmWH1Sqtbo5OdE7IpQBXbswoGs4A7p2wd/Dsc+Djv04nHOWN5N3sPnkaXzcXPnt+BHcPmwAng7QE5VS8vc9G/n4+N42HcdhDcFilTgZBf/5ehu94sJIiArGapUd3gxMZjMrDx3no217OVVUQpiPFw+NvYYb+vciQmOpdNuDIC9PxiTGMiZR6bZLKckpq+Bg9lkOZudxIDuPj7bt5X2rFYBuQQEMjOrCoKgIBkdHEOHn41ArR3XaFyklO09n8f5Pe/jpVAa+7m78dvwIbh3aH283V7XltQgpJS/u38xHKXu4u8dgnm7DsRzWEM5NfH6xZi9v/H4eABLJFRKvOixmi5VlB47yVvIO8iur6BEWzEtzpjG1V0KnHiYRQhDp70ukvy/X9UkCoNbUwJHcfPZl5bIvM5c1R0/y1b4jAIR6ezEoOoKh0REMiYkkLihAN4hOiNliZcPxNP770x6O5OYT5OnBYxOu5Zah/R1uburNw9t458gObknsz1ODJ3Q+QzBbrDgZDXy6ahddw/zpm9gFKTtm70BKyaYT6by8YSunikroFxHGc7MnMyIuSr+QNYO7izNDYiIZEhMJKMNrJwuK2JOZw96MHHafyWb1kRMABHp6MDg6giHRyv4JwYF6ypMOTGVdPV/vO8KinQfILa8gOsCPv82YwOx+PXF1wEi2tw9v5+UDPzInrhfPDJvS5muC450Bfu4dfLxyF889PAMAq5QYO9gFMru0nL999wM/pp0hJtCf1xfMYFKPeN0IrhKDQZAUFkxSWDC3Du2PlJLMknJ2Z2Sz+0w2uzOyWXvsJAC+7m4MiY5gaEwkQ2O6khgSpBtEB+DE2UK+2HOIlYdSqDE1MCQ6kj9OHcO4pDiHvZF8/+hOXty/mVmxPXlpxHUYbHBdcEhDADiSlkdSTAjX9I3pcL0Ds8XKJzv28cam7RgNBv44dQw3D+nXqYeGbIkQguhAP6ID/Zg3sDegmK9iEDnsOpPFhuOnAPB1c2VITCQv3DDV4YYSOju1pgbWp6SxeM8h9mXl4upkZFqvJG4f1p9eXULVltcmPjm+l2f3buK66O68fO0Mm13/HNYQeseH8/oTc4GO1TvIKinjd8vWcDA7j3GJcTx93XiHyaXuyJybh7ihfy8Acssq2NXYe0gvKsHTxbnZz5otVoSg2X/Kc4WbLFYrKw6m8MzqTfx5+jjmDOjVLr9LZ0ZKyaGcsyzdf5TVR05QVW8iOsCP308ezQ39e3WIFBOfHt/L07vWM7lrAq+Ouh4nG94M28QQhBAfAjOAAill7ybeHwusAE43/miZlPLvbW3XuTEnTkfpHaw5msqfV65HIPjX3Glc1ztJHx5SiS5+Pszu35PZ/a+8ruXCld1SSixWiRDKgrpzaT4Aluw9zJd7DhMT5I+pcdGdXrPDNmSWlPHtoeOsPJRCRkkZbk5OTOmVwNwBvRkcFdFhzvGi4/v4y671TIyM543Rs3A22HbUwFY9hI+BN4FPL7PPj1LKGW1tqCP+AzVYLDy/ZjP/232QfhFhvDxvOpH+vmrL0mkBf1y+jhqTidEJsYyMj25M/fHz9/Nc7+CVDVtpsFgZnxRHTlkFccFKDhzd71tPRnEZa4+lsvbYSY7mFSCAITGR/GrkEKb0THCYsNGW8nnqfp7atY6JkfG8PeYGXNphCNkmhiCl3CKEiLHFsa7Ecx+so6rWxAuPXG+P5tqdito6Hlmyih2ns1g4fBCPTbxWnytwIG4c3Ic1R0/y6Y79PL1qI85GAz3DQhgZH8O0XolEB/qxeM8h9mbm8MaN1/P90VROFhYTF+QPoPcArwKrVXI49yybTqSzKTWdE/lFAPSLCOPxSaOY3juJ8A46vLr45EH+tGMt4yO68daY2e1iBmDfOYThQoiDQC7wOynl0as9gNlsYfPeNK7tH2d7dSqQXVrO/Z8vJ7OkjOdnTz4/fq3jOPSLDP9FQZT0whI2paaz83QWRVXVTOwRz4aUND5buAAhBDllFfi6uRHk5XmZo+qco6Cyim2nMtmWnsG2U5kUVddgFIKBURE8OWUMk3vE06WDL8hcdGIfT+1cx5gucbw99gZcje132RZSXpyguJUHUnoIq5qZQ/ABrFLKKiHEdOA1KWWTuamFEPcB9wEEBwcPWrJkyfn3UnMq+GBNGndMjKNXtJ9NdLeVqqoqvLy8rvpz2VW1vHYwDYuU3N8rliT/9r2zaa1Oe+NoOq1SIiXn5wwupKy+gUUnMjhZVs2NCZG4GARHSiqI9fFkbEQwVikv+YyUkkPFFUR4uhHo5tLmHoSjnc+y+gbSyqtILVO2szV1AHg7O9Hd35vegT70DvDBU6U1A/Y+nxsrz7Ko9Az93Px4ODgRZ3Hl+dJx48btlVIObk17djmrUsqKC56vFkK8LYQIklIWNbHve8B7AElJSXLs2LHn39v18UbcXJy45+brcLtM1Ic9SU5O5kKNLeFYXgG//3Qpnu5ufHj7XOKC278OQGt0qoGj65RSnjcJJ6OB/oNL+elUBvuz8tiXmUtpTQ1ZdWZmjx5J/66XllrMKinjwdc/AsDb1ZWk0CASQgKJDwkkPljZAjzdW2wUWj6f5bV1pJwtJCWvgPUZp8mrP0teRSUAHi7ODIqK4NaYSK7tFk330GBNzB3a83wuOr6PRZk72nXO4GLsYghCiDAgX0ophRBDAQNQfDXHkFKyZd8prukboxkzaA3Hzxay8JOv8XBx4ZO75hEV4Ke2JB0bIoT4RQj0uQR9tw7tz+I9h1ifksa4xLhmx7pDfbxYfM9NHM8v5MTZQo7nF7Lq8Akq6+vP7+Pj5kpsoD8xQf5EB/jR1d+PqABfIv18r8os7EVZTR2ZJWWcLi4hraCYtMISUguKyCk7f59IgKsL18THsLDrQPpHdqFneIjN6nI4Iv89totn9vxgVzMA24WdfgGMBYKEENnA04AzgJTyXWAe8KAQwgzUAjfJqxyrSs0opLC0ilEDutlCsiqkFRSz8NOleLi4sGjhfD2SqBNhsVqpazDj6+bKbcP6N7ufi5MT/buG/6L3IKWkoLKakwXFpBcVk15UyumiUnakZ7HiYMovPu/qZCTc14dwH6VcaW1pMRmu+wjwdD9fCc/XzRUvN1e8XV1bfdE1W6xU1Zuoqq9XSq5W11JSU0NRVQ1nK6rIr6gkr7ySrNJyKup+NjNng4GYIH/6RoRx46A+9AgPoWdYCIf37NJsT8bevHNkB//cl8y0qCReGzXTbmYAtosyuvkK77+JEpbaarYeSEcIGNHPMYtT5FdUce9ny3AyCD66c65uBp0Mo8HAXcMHnl9/cDUIIQj18SLUx4uR8dG/eO9cgaGskvLGutiV5JVXkFdRxe6MHAoqKlmXVdDssd2cnHB1dsLNyQk3ZyeMBgNGg8DYOFZtkVaklJitkvoGM3VmM/UNZmoaGpo9po+bK6E+XoT5eNMvMpyu/r7ny6RGBfjpUXSX4c1D2/jXgS3MjOnJKyNn2HTRWUtwmJXK2w+epkdsKAG+HmpLuWpqTA08+MUKKurq+XzhAmIC/dWWpKMSLk62/Ze7sMBQU2zatIkBw4ZTWlN7fiuvraOqvp7KOhPVJhN1DWbqzWbqGsxYrFYsUmKxWhEIDEJgMCjDYOeMw9XZCS9XF3zcXPFydcXP3Y0ATw8CPd0J8PRwiPoBWkNKycsHfuTNw9u4Ia4X/xpxnSoLbh3CEMqrajl66iwLZw1TW8pVY7VKfr9sDcfPFvL2zTPpER6itiSdToQQAj8PN/w83IhFvxHRIlJK/rFnIx+m7OGm+H48e80U1bIvOIQh7DmWhVVKrukTo7aUq+btLTtYfzyNJ6eMYWxix1g/oaOjYxusUvLUznV8nrqfhT0G85fBE1QNCnAIQ9h9NBNPdxd6dgtTW8pVsTUtg7eSdzC7X0/uvGaA2nJ0dHQ0hNlq5ffbV7P01BEe7H0NTwwYo3qEmMMYwsDukQ4VhlZcVcOT36whPjiQp68br/ofWkdHRzvUW8w88uNK1mam8lj/UfymzwhNXCM0f4UtKqsiO7+MAd0j1ZbSYqSU/Hnlesrr6vnX3Gm4O/C6CR0dHdtS02Di3h+WsjYzlaeHTOSRvtdqwgzAAXoIB07kANA/KUJlJS1n6f6jbEpN58kpY0gKC1Zbjo6OjkaoMNVx7w9fs6cwhxeHT2NBQj+1Jf0CzRvCoZO5uLo4kRTtGNE5eeWVPL9mM0NjIrljmD5voKOjo1BQW8WdG5aQVl7E66NmMiOmh9qSLkHzhnD4ZC49Y0NxcnKMxSzPfZ+MVVp5btZkTeRe0dHRUZ+sqjJuW7+YgtpqPhg/n9FdtLnAVtNzCFIqKSv6JHRRW0qL+DHtDOuPp/HA6GH6SmQdHR0ATpYVsWDN55TW1/L5pJs0awag8R5Cg8WK2WKllwOEmzZYLDz3fTLRAX4sHD5QbTk6OjoaYG9BNnf/8DUuRiOLJ99Cz4BQtSVdFk33EExmKwA94rRvCIv3HOJ0cSm/nzLa5ukJdHR0HI9N2ae4df1iAlzdWTb1ds2bAWi9h2C24u/jQYi/tgt8VNbV81byDq6J7co4fTWyjk6n55v0Izz+02q6+wfz8YQFBLk7RoU8TRuCyWyle0yIZmJ0m+ODn/ZQVlvHE5NHa16rjo5O+yGl5N2jO/nnvmSGh0Xx3ti5eLu4qi2rxWjaEMwWK4kaDzetbjCzaPdRpvZMoKeeuE5Hp9NisVr5x56NfHx8L9fH9OBf117XrvWP2wNNq5US4rs2ndZXKyTnFFFtMvHgaMfLxKqjo2Mb6ixm/m/rKr7LOM69PYfwx0HjL6mX7Qho2hBA24ZQa2pgU04hYxJi9RXJOjqdlNK6Wu5LXsrugmz+MHAc9/d23JtDTRuCEBAVpt0c7isOpVDVYObeawerLUVHR0cFsirLuHPjErKrynlz9CxNrj6+GjRtCE4Gg2ZXKFutkk937CfKy53B0Y6TZ0lHR8c2HCrK455NX2OyWPh80k0MCe2qtqQ2Y5N1CEKID4UQBUKII828L4QQrwsh0oQQh4QQLVq55eyk3TG47aczSS8qYXyk9qOgdHR0bMvazFQWrP0cF4ORr6fe1iHMAGy3MO1jYOpl3p8GJDRu9wHvtOSgWq5/8MXugwR4uDMoxE9tKTo6OnZCSskHx3bzQPIykvyDWT79ThL8tDvPebXYZMhISrlFCBFzmV1mAZ9KKSWwQwjhJ4QIl1LmXe64zk7aNITCymo2nUhn4YhBOBssasvRaSVSQmkpZGRAcbHyvKzs0seyMqirA4vl5+2uu+BPfwKjUdnc3MDPD/z9f3688HlgIMTEKK91HBOz1cpnpWfYmJXP1KhE/j3yetydtFXr5MfC/W36vFCu0W2n0RBWSSl7N/HeKuAFKeXWxtcbgd9LKfc0se99KL0IgoODBy1ZssQm+mzJhqwCvj6Vw9NDeuAtzXh5aXslNUBVVVWn1Vlff+lmMimPAC4u4OSkbOcu8Bc/F0LZQHk0m6swGhWdUoLV+rNZmM2XPjeblfaEAFdXZXNx+fn5uc3WdOa/uy2ptpp5p+gkR+rKmeodzgK/KM2FlZZSxSdiIyvGvrpXStmqSBd7TSo3deaadCIp5XvAewBJSUly7Nix7Sirdbz+n8/pFR7CzddNJTk5GS1qvJjOoNNshpQU2LcP9u5VHg8ehOBgSEyE2NhfbnFxyt17a/6vW6NTSqUncvo0pKcrj+e21FQoKYH+/WHQIBg4UHlMSlJMqbV0hr97e5NVWcbCH77iTH0lCwPieHrGArUlXYLZauHxg6/iWtu2uwp7GUI2cOGsSySQa6e2bUpGcRnH8gp4YvJotaV0eioqIDkZNmyAXbvg8GGIjPz5gnrDDTBggHaGaYSAoCBlGzLk0vdLSn42s1Wr4G9/g7NnoV8/GDoUJk6EMWNAwzfSHY6d+Zk8mPwNFilZNOlG6o+fVltSkyzJWkdqZQZP9ljIEv7Z6uPYyxBWAg8LIRYDw4DyK80faJU1x1IBmNYrUWUlnY+GBuXCv369sh08CNdco1woX3xRubv28VFbZesJCFB+l4kTf/5ZWRns3w/btsFLL8GNNyqGN2mSsg0e3LYehE7z/C/1AH/ZuY4obz/eHzeXbr6BJGvQEE5WZvJFxlrGhgxmVHDbqjTaxBCEEF8AY4EgIUQ28DTgDCClfBdYDUwH0oAaYKEt2lWD9Slp9IsII9zXW20pnYKyMli+HL75BjZvVoZ6Jk2Cp5+GkSPBw0Nthe2Lnx+MG6dsf/oTVFXBjz8qhnjvvZCTo7x3ww0waxZ461/LNmO2Wnl27w98lLKHMV3ieH30THxd3NSW1SR1FhMvn1iEv4s3D8bPa/PxbBVldPMV3pfAQ7ZoS03yyis5kpvP/00cqbaUDk1NjTJk8sUX8MMPMH48LFgA778PIZ08f6CXF0ybpmwAeXmwbh18+SU89BBMngw33wzTpyuRTzpXR2ldLQ9vWc5PZzO4u8dg/jhoPE4GbUY7Anx0egVZNfk80+chvJzafnek6ZXKWmNzajoA45P0mge2xmRSLmyFhTB7NgwbplzYPv4YfPVqpM0SHg533qlsJSWwdCm8+abSe5g5E+bMUSbb9ZpNV+Z4aQG/2rSU/JoqXhwxnQXxfdWWdFn2lqSwKvdHZkWMZYB/kk2OqV3r0yCbT54h0s+HuKAAtaV0GPLzlcnT6Gh44QXlDjg1FdauVWL9dTNoOQEB8KtfKb2qo0eVCfXcXGWY7YUXlAgnnaZZnXGcOd8vot5iYcnUWzVvBgClpgrivbpyV+z1NjumbggtxGS2sPN0FqMTYvVUFTZg3z7lrrZ7d2UcfP162LpVCRHt7MNCtiA8HB59VDm/K1fCiRMQH68YxpEmE8x0TsxWK8/v3cSvNy8nyS+Yb6+7k/5BXdSW1SImhg3jlQH/h4vBdovjdENoIQez86hpaGBEXJTaUhwWsxm+/hpGjVKGhXr2hLQ0eO896H3JckYdWzFgAHz0kWIKUVHKPMOECYpRWK1qq1OPkroa7tzwJf85upPbkwby5ZRbCfVwrFl5o7DtJVw3hBay60w2AhgSE6m2FIdDSli2DPr0gZdfhkceURZm/f73SkoHHfsQEgJPPQVnzsDdd8Pf/66YxerVyt+oM3GwKI/rv/uY3QXZvDRiOv8YNhkXPX5XN4SWsi8rl4SQIHzd9dCNq2HTJmWtwN//Dq+8osTTz5+vT3KqiYsL3Hor7N6tzN/87nfKgrdt29RW1v5IKVl0fB/z13wGwNdTb2O+BucLzFYLZ2t/nvSxVYqhK6EbQguwWiUHs/MY0DVcbSkOw759MGWKEu3y6KPK62nTWpcmQqd9EEIZujt8WOkx3Hyzspaho84x1DSYeOynVTy1ax3Xhsew6rqF9A3S5v/0i8c/4fmUD9mYvwsp5fl5y/Y2Bt0QWsDp4hKq6k30i9Tml0dLnD0Lt90GM2YoF5eUFLjlFtBwKHenx2hUIrpOnICxY5X5hXvu6VhRSallhcxa/SnL04/yWP9RfDB+Hv5u7mrLahZnYcTd6MrBslT+c2opp6tyANo9oEX/N20Bh3PyAegTEaqyEu1itSqTw337QteucPIk/PrXyvCEjmPg5gb/7/8pYb/e3tCrFyxa5PjzC0vSDjHzu08ora9l0cSbeKTvtZrLVHoxQwJ7UWupJ8ErioqGar7IXMPWwv08d+wDzlTntltPQR/JbQEpZwtwd3bS1x80Q1qaMuTQ0AAbNyqTxzqOi68vvPqq0tP71a8UU/jvf5UIJUeiusHE07vW8/Wpw4wIi+bVUdcT4u4YmQHHhgwmvSqH/v5J9PNPZGfxYf6XsYbyhqpfDCHZGr2H0AJSzhaSGBqMUR/3+AVWK7zxhjJpPGeOso5AN4OOw+DBysTzuHFKQr0PPnCc3sKR4rPMWPURS08d5jd9RrBo4o0OYwZWqcQCezl58EH6cqI8wujhHUuRqYx47658eHrF+SEkW6P3EK6AlJIT+UVM6ZmgthRNUVSkzA1UVyvRKYl68tcOiZMT/OEPcP31yjzDypXw6afaXUEupeST43t5bu8mAt08+GLyLVwT5lhdG0Pj2oIbIseSmZpHXm0RX2at4/ouoxkW2JvjFWeI9Ypon7bb5agdiOLqGspr64gP1gPmz7Fvn5LPf+BA2LJFN4POQO/eivFHRSm1GY4dU1vRpRTUVrHwh6/46+4NjAyPYfWMux3ODC7E2eBMN69I7t/zDJUNNdwWPZ1E72iu6zKq3drUewhX4FRhCQDdgvX5A1DuDv/v/+Cdd2Be27Pt6jgQLi7KEOEnnyjrFrT0HdiUfYrfbfuOqgYTfx86iduTBnaIFDNTwkZQaqpkQuhQhBBIKW2+OvlCdEO4AmeKSwGICfRXWYm6NDTAY48pSeeSk5UIFJ3OyZ13KnNFc+bAnj3w7LPqFempbjDx/N5NfJa6n+7+wSyefAsJfkHqiGkHPJzcuDN2BkahnOD2NjndEK5ARkkZLkYj4T6OlePEllRUKKmUvb2VimVaKUmpox4DBypmcNNNypqTpUvtX6xoX2EOj21dRUZlKff0GMLjA8fgZnT8S1qD1Yyz4eff45wZ2AN9DuEKZJeWE+Hng8Hg+N3P1lBaqlQo69kTVqzQzUDnZ4KCYM0aJUPtddcp1dzsQb3FzMv7tzB/zWc0WC18MfkWnhoyoUOYQX5dMffs+ht7StSZpNEN4QrklFUQ4efAhXrbQHGxUt93+HB46y19tbHOpTg5KZlU4+Jg6lSlN9meHCk+y6zvPuGNw9uYFduL76937InjC5FS8ubJL6mx1BHlEaaKBpv8iwshpgohTggh0oQQTzbx/lghRLkQ4kDj9hdbtGsP8iuqOmX95IICpXTlxInw73/rOYh0msdoVMqb9u2rpNYuK7N9GyaLhdcObmX26k8pqa/lw/HzeGXkDHw0Wuu4NWzI38m+0uPcFXs9IW7qBLG0uY8lhDACbwGTgGxgtxBipZTy4j7Pj1LKGW1tz56YzBaKqmsI9XGMBS22Ij9fMYM5c5QspboZ6FwJg0HpRf6//6fkQlq/XqngZgsOFOXyxLbVpJYVMSu2J38bOgk/V+3mIWoNRfVlvH/qG/r4xjM9XL2a7bboIQwF0qSU6VJKE7AYmGWD46pOSXUNAMFeniorsR91dUpSuhtugH/8QzcDnZYjhNKbHDlSSXHe0NC249WaG3huzw/M+X4RlaZ6Phg/j9dGzexwZiCl5O20JZilmUcSbz6/ME0NRFuTJAkh5gFTpZT3Nr6+HRgmpXz4gn3GAktRehC5wO+klEebOd59wH0AwcHBg5YsWdImfW0ho7KG5/ee4IHesfQP8mtyn6qqKry8tN+DaKnOM2eUlBRxce2vqSk62vlUG7V0pqWBq6uS6LAlXKzzUG0pn5acochSz1ivEBb4ReFhUH/SuD3O5wmyWWnYxRhrb4bS9lWe48aN2yulHNyqD0sp27QB84H/XvD6duCNi/bxAbwan08HTrbk2ImJiVJNtpw8LZOefkXuzchudp9NmzbZT1AbaInOl16ScsAAKaur219Pc3Sk86kF1NJZViZl9+5S/uc/Ldv/nM6Cmir58OblMvqT5+WE5e/J7XkZ7SeyFdj6fFY2VMtbt/9J/mbvP6XZarbJMYE9spXXc1tYbjZw4X1AJEov4ELTqbjg+WohxNtCiCApZZEN2m83ymvrAPBz71hd1KZYvVqpaLZzp/3jyXU6Hr6+St6jkSOhe3cYPfry+1uk5KOUPbxy4EfqLWb+X7+RPND7Glw7QCjp5fgwfSXlpkqe7nWfXdcbNIctzvZuIEEIEQvkADcBt1y4gxAiDMiXUkohxFCUuQvNl9+oqK0HwNvNVWUl7UtampK4bPnylnfxdXSuREICfPYZ3HijkjU1sply5PsLc/n72cNkZNUwKjyWp4dOJN634+cOO1R2krVntzE3cgIJ3toInW2zIUgpzUKIh4G1gBH4UEp5VAjxQOP77wLzgAeFEGagFripsWujaapNJgC8XDtulRcp4b77lIL3I0aorUanozFpEtx/Pzz0kHLDcWGQQn5NJf/ct5ll6UfwN7rw9pjZTItK6hA5iK5EvcXEGycXE+YWyC3R09SWcx6b9MeklKuB1Rf97N0Lnr8JvGmLtuxJjakBgxC4OXfcbutHH0FlpVL3WEenPfjDH2DAACW9xbx5ypqCj1L28Pqhn2iwWniw9zX0LbMwLbq72lLtxpKs9eTWFvJMn4dwM2rnhrPjXulsQG1DA25OTh32juXsWXjySSVm3En/Jui0E66uysK1BQsk5u7HeTt1M5lVZUyMjOepIROI9vYnOTlZbZl2I6smn6+zNjAuZDAD/JPUlvML9MvAZTCZLbg4qT/R0148+qhS+rJfP7WV6HR0PBJzCf7jRp7cl0N3/2A+nrCAsREqxTariJSSd9K+wtXowj1xs9WWcwm6IVyGBosFZ7Xy+rYza9fC3r3w8cdqK9HpyKSWFfKv/VtYl3WSoCBPrF9M48+P92FkROdMjLWpYDcHy1J5KH4B/i7ay5GmG8JlsFglxg6a5fQf/4BnnoFOEFGrowI5VeW8fugnvjp1GE8nFx7rP4q7ewzmiwpXnntWCXPubFQ2VPN++jd094lharg2Izh0Q7gMEomhA84fbN0KeXnaqXal03E4W1PJW4e3s/jkAQSChd0H83CfEfi7KXced9wBTz8NBw92vqHKT86soqqhlof73KhqeorLoRtCJ+Sf/4THH9cnknVsR05VOf85upPFJw8ikSyI78dDfYbTxfOXwyKurkoCvBdfhM8/V0msCqRWZrAmbxuzIsYQ6xWhtpxm0S8Jl8EgBFbtL5e4Ko4cUSpdffWV2kp0OgIZlaW8c3gHS9MPAzA3rg8P9R1OVy+/Zj9z//1Krqz0dPVyZtkTq7TybtrX+Ll4c6uG1hw0hW4Il8HJYMBitaotw6b8619KdJFbx0kjr6MCB4pyeffIDtZmpuJsMHJzQn/u7zWMCC/fK37Wx0dZDPnvf8Mbb9hBrMqsP7uDE5UZPJZ0Gx5O2p600w3hMjgbjZjMFrVl2Iz6emW16IkTaivRcUQsVisbs9P4MGUPO/Iz8XVx48Hew7mr+yBCPK4uA+g998C118Jrr3XsSnyVDdV8fPpbevnEMT5kiNpyrohuCJfB1clIndmstgybkZys1EYODVVbiY4jUW6q48uTB/n0xD6yq8oJ9/Dmz4PHc1NCP7ycW5fnq1s3pSbzzp1KidaOyqIz31FlruHB+PkOscBVN4TL4ObsTL3ZgsVqxdgBbmNWrFCK3+joXAkpJQeL8/j8xH6+PZNCncXM0NCu/GnQOCZ1TcTJBv8Ps2Yp38mOaghnqnP5Pu8npncZpemJ5AvRDeEyeLg4A1BrasCrA2Q8XblSSVOho9McRbXVLD99lKWnjpBSWoCHkzM3xPXm1sT+9A60beH3mTNh4UJ44QWbHlYTSCl579QyPJ3cuU3jE8kXohvCZTiX5bSq3uTwhlBTo9Q56N558ofptJA6i5nk7FN8c/ooG7PSMEsr/QLD+cewycyO7YW3S/t894cMgbIyOHmyXQ6vKrtKjnCwLJUHus3D29lxSvDqhnAZfBpNoKKunjBfb5XVtI3qahg7Vq+RrKNQZzGzLe8M32UcZ13mSSob6gly82Bhj8EsiO9Lgl9Qu2swGGDMGNi+HaK0UQ7AJjRYzfw3fTmR7iFMC79WbTlXhW4Il8HXXYnNPFc5zZGpq9N7B52dSlM9W/NOsyYzlY3ZaVQ1mPB2dmVKVCIzY3syIizaJnMDV0NSkhL11pEMYVXuFnJrC/lb7wdwMjhWLjTdEC6Dv4cSM1xaU6uykrajG0Lnw2K1kmGq5r2jO9mUc4rd+dmYpRV/V3dmRPdganQSw8OiVC1T2b07fP21UkinI1DRUM0XmWsZ5N+DwQE91ZZz1eiGcBkCPJXiwiXVNSoraTt1dcrdmE7HxSolJ8uK2F2QzfazGfyUd4YyUx2che5+wfyq11DGRXRjYHCE3XsCzXGuh9BRWJy5hlpznSZTW7cE3RAuQ4CnOwIorHJsQ6ipgYYGiIlRW4mOrZBSkl9bxeHisxwuPsuh4jz2FeZQYVLqgIe6ezGhazyBZbXcPW4qYR7anANLTFRqencEztYW813uViaFXUO0Z7jaclqFTQxBCDEVeA2lpvJ/pZQvXPS+aHx/OlAD3CWl3GeLttsTZ6MRfw93Cquq1ZbSJjIzwcUFOmhphw6NlJKC2ipOVZRwuqKEk2VFpJQWcKK0ULn7R8m5Fe8byPSo7gwOjWRISCRRXn4IIUhOTtasGQB4eoKvr3LD4uh8lvEdBmHQVI3kq6XNhiCEMAJvAZOAbGC3EGKllPLYBbtNAxIat2HAO42PmifUx4v8iiq1ZbSJurqOnR7AETFbrVSY6igz1VFWX0txXQ2FtdUU1laRX1tFTlUFOdXl5FRXUGv++Wrp4eRMd/8QpkV3J8kviN6BYfT0D8HDWTt1ea8Wd3dw9ByS6VU5JBfsZW7XCQS5+qktp9XYoocwFEiTUqYDCCEWA7OACw1hFvCplFICO4QQfkKIcCllng3ab1fCfLzJLa9QW0abqK/XDeEcDVYLZ6srKaqroaiumpK6Giob6qluMFHdYKLeYsZktdJgtWC2WrFIK1YpkchfXLQKCwv5MrkUUOpmWKU8v59FSixWK2ZpxWSxUG8xU2+xUGdpoMbcQFVjO80R4OpOhJcv8b6BjOkSR4y3P7E+AcT6+BPu6dPhanS4uoKj55D86PQKPJ3cmRc5UVUdDda2jWbYwhAigKwLXmdz6d1/U/tEAJcYghDiPuA+AP+gMNWLb1urKsgsKm1WR1VVleoar0R1NQQGal8n2O58SikpMNeTbqoi3VRFXkMt+eY6is31NHftcRYCF2HECYGTEBiEwIDAAAjEL9ZwWK1W8vJ+jj4TCAQgAKMQGFE+b0TgIgz4CgNBwhk3J1fcnI24GYx4GpzwMDjhZTDiZXDGz+iMj9EZpwuLp1QD1RWYz1ZwkjNc7RouR/h+/uY3YDJpXyc0fT6zKGSf4ThjrL3Z+9MudYQ1Uuy3vE2ft4UhNHW7cnEHsCX7KD+U8j3gPYDwrnFy7NixbRLXVtJd9rI5dwsDhl1zfl3ChSQnJ6O2xiuxfz9s357M/Plj1ZZyRdpyPustZjblnGJNRipbctMpqVcu2O5OzsT7BnKNdxdivP2J9PIl2N2TIDdPAtw88HZxxdPJ5aoibxzh7w6OofPuu+Gdd7SvEy49n1JKnjj4GoF1vjw65C5cjeoN3ZmtdSw5/WybjmELQ8gGul7wOhLIbcU+l2C2qN+P7Oqv5HfPKi1v0hAcAVdXxx+jvRyFtUqs/ZdpB6kw1RPg6s7YiDiGhHSlX1A4iX7Bmgmz1LkURx7S3F96nGMV6fw6fr6qZgCQUZWsiSGj3UCCECIWyAFuAm65aJ+VwMON8wvDgPKWzB+YLepfxaIC/ADIKC6jdxfHzBsdGQkmkzJO66j/eE1RVFvNO0d28HnqfkxWC9dFd2d+fF9VVtzqtI76eigqcsxyrlJKPs/8nmBXfyaHXaO2HE5VfI+XU5c2HaPNfwYppVkI8TCwFiXs9EMp5VEhxAON778LrEYJOU1DCTtd2JJjmy1WpJSq5hGPauwhnCkuVU1DW/HxUUJOc3Kga9cr7+8I7MrP4sHN31BWX8sNcb15qM9wYn0C1Jalc5WkpSnrYxxxnvxA2QmOV5zh1/HzcTY4q6qluqGAvNo99AtYCCxt9XFs4stSytUoF/0Lf/buBc8l8NDVHtcqoaC0itAA9eKo3V2c6eLrw2kHNgRQSmaeONExDOF/qQf4y851RHn78cXkm0n0C1Zbkk4rOXHCMVfQSyn5LGO1ZnoHpyvXAZI476ltOo7m+9UZuSVqSyAuyJ9ThcVqy2gT5wzB0Xnv6E7+uGMN14bH8M30O3QzcHAc1RDO9Q4WdJ2keu8AIL1yLcFuvfBxiWzTcTRvCGc0YAgJIUGcKizB4sDB0q6ujm8IG7PTeG7vJq6L7s6H4+fh6+KYk/w6P+OohvBl5joCXXyZFKb++trS+nRKTaeI9Z7S5mNp2hAMQpCeo/6deUJIICaLhYySMrWltBpPT/jpJ7VVtJ5KUz1PbFtND/8QXh45o0OUNO3sSKl8JwcNUlvJ1XGsPJ3D5WnMiRyvid7BmaoNCAzEeI1v87E0/V/l7CRIz1bfEJJClWIhqflFKitpPV5ekJEBWVlX3leLfHpiH8V1NbwwfBpuKqZr1rEdKSlKlFH//moruTqWZm/E28mDqeEj1JaClJLTlRsIcx+Iu1Pbgyq0bQhGA6eyi5AqB9HHBwfiZDBwLK9AVR1tZfp0pa6yo1FnMfPBsV2MjYijX5BjZpHUuZQVK5S6yo4UYVRCJTuLj3Bdl1G4GdUvq1tcn0JlQzax3rYpKKFtQ3AyUFVTT0GJusnlXJ2d6BYcQMrZQlV1tJVZsxzTEHblZ1FSX8sdSQPVlqJjQ1auVL6TjsQekYaTMHJ9l9FqSwHgTOUPGHAiymuMTY6neUMAOJmp/oW4V3goR3LzVe+ttIXJk5X6teXlaiu5OrbkpuNiMDI8LFptKTo24uxZOH5cqansKJSZKjlKBuNDh+Lnon5KcSklZ6p+oIvHUFyNPjY5prYNwWhACEjNUH+opk9EKKU1teSUOW7mU29vGD8evvxSbSVXR0pJAT0CQnB3Un8CT8c2fPYZzJih1OlwFFbnbcUsrNwQOU5tKYAyXFRtPku0d9snk8+haUMQAiJD/TihAUPoH6mMXe/P0nzG7svy2GPw0ktgsaitpOUU1FZrusiLztVRXw+vvqp8Fx2FBmsDq/N+IlaG0tVDGylsMqqSERjp6jnSZsfUtCEAJEWHcOKM+oaQEBKEh4sz+7OumJNP04waBcHBsGyZ2kpaTmVDPd7O6k/g6diGzz+HXr1gwAC1lbScn4oOUmqqYJCMV1sKoAwXZVZtJsx9oM2Gi8ABDKF7TCh5RRWUV9Zeeed2xMlooH9kOHszc1TV0VaEgCefhBdecJwMqJ7OLtSYTWrL0LEBFgu8+KLyHXQkvsvdShf3YGIIUVsKAOUNGVQ0ZBHlZdvJbe0bQqzSPUs5k6+yEhgSHUlqfhFlNXVqS2kTM2YoZTU3blRbScvwc3GjqK5GbRk6NmDFCqWGsgOUPjjPmepcjlWkMz382sYySOqTVfUjAF09R9n0uNo3hBjFkVPS1TeEoTGRSGBPRrbaUtqEwaDcof3xj44xl5DoF8zx0gKHjvDSUVKwP/208r1zpLUHa/O24ySMjA8dqraU82RXbyPANRFPZ9v2WDRvCN6ebkSH+3M0/azaUugTEYabkxM7Tjvoct8LuPVWJeHdW2+preTK9A4MpcJUT0ZlmdpSdNrASy9BVJSyGM1RMFkb+KFgNyOC+uHr7KW2HADqLOUU1h0m0tP2K6U1bwgAveO7cORkrup3iC5ORgZFR7A9PVNVHbbAYID334e//x0yNf7rjAyPAeCHnDR1hei0mhMn4N//hnfecazewbaiQ1SZa5gSNlxtKefJrd6BxEqkDaOLzuEQhtA3IZzSylqy8svUlsKIuChOFZWQX6Hu6mlbkJQEv/0tPPigtieYo739SfANYm1mqtpSdFqB1Qr33Qd/+YvSQ3Ak1p3dTqhbIH39EtSWcp7s6m24Gf0Jcu1u82M7hiEkRgBwMFX9CJ+R3ZTVsltPZaisxDY88YTSQ9D6YrWZsT3YmZ/F6Qr106HrXB0ffKCsPXjoqktkqUtBXQmHyk4yMXQoBqGNS6VVWsit2UmEx3BEO2jSxm95BWLCA/DxcuPACfUNITE0iGAvT7amnVFbik1wcYH//hcefVQpZ6hVbkzoh5Mw8NmJ/WpL0bkKDh9WJpHff18p4+pIJBfsQSIZFzJEbSnnKa5Lod5aQYRn+1Rpa5MhCCEChBDrhRAnGx/9m9nvjBDisBDigBBiz1WLNAj6JUZowhCEEIyKj+GnUxmYLY5bMOdChg2Dv/1NmezTap6jEHcvZsT04H8nD1BUW622HJ0WUFSkJK979VXo00dtNVdPcsFeevrEEe4epLaU8+TU7AQE4R7tY1Jt7SE8CWyUUiYAGxtfN8c4KWV/KeXg1jQ0ICmC7PwyCkvVH7sfkxhLRV09+xx81fKFPPCAEht+663aDUV9pN+1mCxm3jjkwJV+OgkNDTBvHixYoHynHI0z1blk1OQxJkRbGXbzanYT6JqEm9G3XY7fVkOYBXzS+PwTYHYbj9csA3so1eH3HVd/DcC1cVE4G438cOKU2lJsymuvQXW10sXXInE+Adyc0J9Fqfs5Uqx+GLJO8zzyiJJM8dln1VbSOrYWHkAguDaov9pSztNgraaw7ijhHq26p24Roi2hnEKIMiml3wWvS6WUlwwbCSFOA6WABP4jpXzvMse8D7gPIDg4eNCSJUsAsFolf/v8EH1j/Zg7Uv00yG8eOkVeTR1P9orC21v7ideqqqrw8rpyHLXZrKQlDg+HwEA7CLuIK+mstpr5Y+5B/Iwu/CWsN0aVYhhbej7VRg2dhYVQUAA9eijhzS1Ba+fzI7EBN1y4Wf4yNYSaOmtdT1AQ9AkhRQtxr28+6mncuHF7WzsSg5TyshuwATjSxDYLKLto39JmjtGl8TEEOAiMvlK7UkoSExPlhfzfy9/IOY/9V2qBr/celklPvyI/XvGd2lJaxKZNm1q879GjUoaHS/npp+2npzlaovO7Myky+pPn5esHt7a/oGa4mvOpJvbW+fbbUnbtKuWpU1f3OS2dz5yaAjl982/k8uxNl7ynps7dhW/KT1NHyQZL7WX3A/bIFlxfm9qu6N9SyolSyt5NbCuAfCFEOEDjY5NpSaWUuY2PBcA3QKvWgA/uFUV2QTl5RerXJBjfvRtOBgN7C0vVlmJzevaEDRvgD3+ADz9UW82lTI/uzuzYXrx6cCv7CtUPNNBReO01JXFdcjLExamtpvVsLzoEwDWB2poJz689QJBbT5wMbu3WRlvnEFYCdzY+vxNYcfEOQghPIYT3uefAZJQexlUzpKeyqmX3UfWX1vp7uDM8Loo9BWWqr6BuD3r2hB9+gL/+VVldqjX+PmwS4Z4+PLR5OQW16gcadHZefBFef93xzQBgR/Fh4jwjCHVTYcy0GRqstRTXHSfUvX+7ttNWQ3gBmCSEOAlManyNEKKLEGJ14z6hwFYhxEFgF/CdlHJNaxqLiwwk0NeTXUe0sShseu9ESupNHMh27KI5zZGYCJs2wT//qdz9aQkfFzfeHXMDZaY6frVpKXXmBrUldVqeeUZZfLZ5M0SrP73XJioaqjlecZqhgb3VlvILiuqOIbEQ4t6vXdtpkyFIKYullBOklAmNjyWNP8+VUk5vfJ4upezXuPWSUrY67kAIwdDeUew+monVqv5d+aTu8TgJwarDJ9SW0m5066b8o7/5Jvzud8qks1boHRjGqyOv51BRHr/d+i0Wa8dYF+IomExK2pMvvlB6BpGRaitqO3tLjmFFMiSgl9pSfkFB7UFAEOLWvkblECuVL2R43xjKKms5roH6CF5urvQN8uX7Iydo0Grwvg2IjoadO5VVp5MnK1EkWmFKVCJ/HjyBNZmpPLVrXYccvtMiubkwZgycPQvbtytRaR2B/aUn8HHyJMFbW0mXCuoO4e8Sh4uxfSMaHc4QrukTQ2SoH6UV2iiYMiw0gJKaWn7qILmNmiMgAFavhmuugcGDYc9VrzdvP+7pOYQHe1/D/1IP8NzeTboptDNbt8KQIUqhpaVLwcd2FRxVRUrJ/rLj9PdPwqiR3EWg5C8qrDtKsHvfdm/Lqd1bsDG+3u4s/dfdass4T68Ab/w93PnmwDHGJjr4bNoVMBrhuecUQ5g+XZlbWLhQbVUKTwwYQ02DifeP7cIsrfxl8ASEI+VZdgCkVOpn/OMf8MknMHWq2opsS1ZtPiWmCvr5Jaot5ReUm87QYK1u9+EicEBD0BpOBgMz+3bnf7sOUlJdQ4Cnh9qS2p05c5RFRzfcoMwvvPoq+Pmpq0kIwV+HTsLJYOSDlN3UW8z8Y+hkjC1dGaVzWYqKlGylx48rQ0SOHknUFIfLTgJoKtU1KBPKAEFuPdu9Lf2/xQbMG9ibBquVbw4cU1uK3ejRQxk28vRUEpetaVXcmG0RQvDnweP5de/h/C/1AL/eslyPPrIB33wDfftC166wY0fHNAOAo+XpBLr4Eu6mnWR2AEV1KTgbvPBx7trubemGYAMSQoIYFNWFJXsPayL6yV54eSlDCB9/rCTHu/NO5U5STYQQPDFwDH8ZMoF1mancsn4xxXXamG9yNPLylOR0TzwBX30F//oXuLurrar9OFaRTg+fWM0NNRbXHyfItXu71D+4GN0QbMRNg/uSUVLGtg5QXvNqmTABjhxRch/16qUYhNrzunf3GMLbY2ZztCSfWas/4USphkKjNI7VCu++q/QKEhLg0CG49lq1VbUvJfXlFNaX0sMnVm0pv8AiGyg1nSLQLcku7emGYCOm9EwgwMOd/+0+oLYUVfDygldeUSKR3nxTiUb64Qd1NU2L7s6XU27BZLEw5/tFegnOKyAlfP89DBwIixYpixKffbZj9wrOcaJSiRJM8olRV8hFlJtOY5UNBLjqhtAiLFYr+08oKbHVDDd0cXJi/qA+bDqRTnapRqvM2IFBg2DXLqVW869+paxb2LtXPT39g7qw8ro7ifcN5P7kZby4bzNmfQHbJWzfrtTDeOwxePppJbS0t7YW67YraVVZGDAQ5xmhtpRfUFKvTHQHuNpnotvhDcFoMPDEv1eQmVeKEEJVU7h5cF8MQvD5roOqadACBgPcfDOkpCiRSNdfDzfeCKkq3aCHeXjz5dRbuTmhP28f2c4dG77Uq641cvQozJ6t/H3uvFNZfHjDDaCxYfR2J60yi64eobgaXdSW8gtK609hFK54O9tnGbhDG8K5EpYDe3Rl8dp9AFhUnNQN8/VmSs8Evt53hOp6k2o6tIKLi5La4ORJ6NcPRoyAu+9WxqTtjZvRieeHT+XFEdPZW5jD9FUf8VPeGfsL0Qh798Jtt8G4cTB6tGLWd98NTp00ED2tKot47/aP4rlaSutP4ecSi0HYpyC1QxuCofE25pZpg1i/4zgATkZ1f6U7hw+ksr6er/e1KqFrh8TTU6nClpqqhCxOm6ZciFassH+5zgXxfVk+/Q58XFy5bf1iXty3mQZrx007ciFmsxItNHKkspakb1/FrB97DNzaL6Oy5ikzVVLWUEmsxoaLAMpMp/FzsV+cr2MbgkExhH6JEUgJB1OV3Phqhn72iwxncHQEH23fi8ncOS40LSUgAP78Zzh9WplfeO45JYrllVeg3I7TLj38Q1g5/U5uSujH20e2M2/NZ6RXlNhPgJ0pKVFWlcfFKSmqf/tbOHVKCSf1bZ/SvA7FmWqlNnqMZxeVlfySeksFtZYi/Fxi7NamQxsC/DxsNHl4dxavUYaNrI3zCGrNJ9w/cihnK6r49lCKKu1rHRcXuOUWJWHe//4Hu3dDbCzcey9UVtqn1+Dh7MLzw6fx1ujZnKko5bpVH7HidMdZWGg2w9q1cOaMkrH22DFYvhx+/BHmzeu8Q0NNkVWjJMqM9ghTWckvqTApIey+uiG0HGNjL2HB5AH8uD8d+HnYSK0FJiPjo+kZHsJ7W3frKZmvwDXXKOmTDx+G7t0hO1tJo/zoo8qq2Pb29OtiurNu5j1cExpFhGfLsrTVWTSUA/wCrFblgv/QQxARoUQLeXgo6SY++UQJJ9W5lOzaAjyMbvi7aCtLX3mDEgrr42K/zKsObwhCCCxWKzFdAugTH05qxs8LkNKzi1RJky2E4IFRQ8koKWP1ET32vSVERCj1Fnr0UPIjBQTAXXcpwxx/+APs26dc8NqDUA9vPpown8Ehl4/kWHrqMM/s2cht6xbzfUVu+4i5SiwWJcz38cchJuZnM9i+XTHUkBAIDVVbpbbJrS2gi3uw5lYoV5iyEBjxdrbfUJbDGwJwPoHZ23+cT1xEAACfr97DYy8v5+VPN/GvT+2/Qmpi93gSQgJ5d8tOvZdwlSQmKne3KSlKHh0plTDWsDDl8YMPINPOC8KXnjrMUzvX4SSMLOwxmOSqAt48tM2+Iho5fRreew/mz1cu+AsXgqursqjs0CFlAr+j5htqD87WFhPurq38RQAVDdl4OYdjEPYb3+swI4lfrT/Ap6t2ISXcft1gjqTl8fQDU0HCix9vZOPOVCYMs19aW4NB8ODoYTz29WrWHD3JdX3ss9KwIyEE9O+vbC+8oJjAhg2wfj08+aTSi5g0SdmGD1cuju3B4eKzPL83maeHTOTGBKWE4fFjx9iRn8l9lmG4GNs3JDAvD7ZtU37v9euhulr5nWfMUDLNRmgvOMZhsEorBfUljAhu39KUraGyIQdvZ/v+cdtkCEKI+cBfgR7AUCllk2VThBBTgdcAI/BfKeULbWn3YnIKyvh28xEWzhxGl2BfPli+g9iIAAYkKUMAcyf2Y+32FLsaAsDUnom8E7KTN5K3M6VnguohsY5OVJQSK3/33crw0cGDikG8/bYyvOTpqayUHjRIGS8fNKjtlbyklLx68EeGhkaeNwOA7IYaPD08bGoGUkJOjjI8tnevsu3bB/X1MGwYTJyoDAn17t35Fo61F2WmSszSQrCrv9pSfoGUksqGHILd7FvKs609hCPAHOA/ze0ghDACbwGTgGxgtxBipZTSZiEdOw9nEODrwZwJyj+s2WLl+Q/XA5B5tpT9x3MIDfBGSmnXcUKDQfDb8dfy0OKVLDtwlAWD+tit7Y6OwQADBijb448rF9PTp3++mL7+uvLo4qKYQ0KCEsl04ebldeV2duRnsiX3NDvnPXz+Z1lVZZRaTCR6dcEq5fn1MC2lokLReuF28qSiXcqfTe3uu5VsslFRugG0F4X1ZQAEufqpquNiTNZKGqxVdp0/gDYagpQyBa4YzTMUSJNSpjfuuxiYBdjMEPokdOHDFTvOvw7y98TP253b/rSI6HB/KmvquW/ucFUmjcYnxTGgazhvbNrO9X264+7ibHcNnQEhlHHzuDglrBKUi2tmJhw4oMTdnzwJ69b9fBH29laMISYGgoLA318p9HPh41c1JxgT2IO6Eg/yjWDBzKb8LHJMtdzk3Y3CAoHFAnV1UFoKZWVNPxYV/dxube0vjSkuThkCGjBAGf7RL/72o7ShAoBAF20tyKhqyAPAy86GIGwRqy+ESAZ+19SQkRBiHjBVSnlv4+vbgWFSyocv3rfx/fuA+wCCg4MHLVmypEUaXvsmhdgwL5K6+rI/rZjoUC+6hXuTVViNv5cLsWFe547fit+weaqqqvC6wq1mWlkV/zpwkpmx4UyPVifWuSU6tYA9dTY0gMmkDMlYLErsvsXyy+f7RQHpopzrahKQEtKcSzjmUkiUu5H+RfGAcgEXQontNxqVrannrq7KZs81APrf/fIc5DTrDPt5wDoVb65c7dBeOmvcjlEY+BlhBQ/h2nB18wjjxo3bK6Uc3Jp2r/jVFEJsAJq6iv1JSrmiBW00dQVu1oWklO8B7wEkJSXJsWPHtqAJ8A7txoYdJ9hyLJcAH29unDmOuIjAFn22LSQnJ3MljWOBA3Ur2Ziexe/nziTQy/5lNluiUwtoTefO/Ex+vXk51QOUr2x6XjnjwuIZVAHTbh2rrrgWoLXz2Rxq6czPNMGZ/UwZPQkXw5V77/bSmVJWRGEhjB42DXengHZv7xxXNAQp5cQ2tpENXJg1KhKweRD3oB5dGZAUSV5ROREhfr94z95zB03x2MSRXP/Wp7yRvJ2/zpigqhadljMsNIqvp97GS/s3E+bhzWP9RzEiLJqdW3/6xX5rM1MZGxGHq7HDBO51CiobqnE1uLTIDOxJjbkAg3DGzWjfyW57fHt3AwlCiFggB7gJuKU9GhKC82ZwoQmobQYAcUEB3DSkH1/sPshNg/vSPSxYbUk6LSTWJ4C3x9zQ7I3FidJC7k9eRqJfEC+OmE7/IG3lxNFpnipzLV5O2qsAVGMuwsMYZPdrV5viIIUQNwghsoHhwHdCiLWNP+8ihFgNIKU0Aw8Da4EUYImU8mjbZDerp8nnWuE3Y4fj4+bKM6s3qVq3Qad1NPedSvIP5qPx86k01TPn+0U8v3eTZtNb6PySGkstHk7aS/Vaay7C3an9h7wvpk2GIKX8RkoZKaV0lVKGSimnNP48V0o5/YL9VkspE6WU3aSUz7ZVtKPi5+HG/00cyZ7MHL49dFxtOTo2ZFxkN9bOvIcb4/vyn6M7uW7VR+wv1EZ6C53mqbOYcDe6qi3jEmotJbgbHcwQtExtXQMbdp7Q3J343AG96RsRxovrtlBRW6e2HB0b4uPixvPDp/HpxBupbWhg7hq9t6B16q0mXA3aqpIGUGcpsfv8AXRgQ1izLYU/vfkdBxprJGgFg0Hw9HXjKamp5dUf1MmFo9O+jO4Sy5qZd7Ogm9JbmPndxxwpPqu2LJ0mMFnNOBu0FQggpZV6SwWuRvuvjeiwhjDt2h74eLnxxfcqVnhvhl5dQrmlcYL5ULZ+oeiI+Li48cKIaXw0YT5l9XXMXv0prx3c2mmqszkKFmnByY7J41qCyVqFxIqb0c/ubXdYQ3BzdWbu+H5s2XeKrPxSteVcwm/HjyDE24s/rViHyawPKXRUxkV0Y93Me5ge3Z1/H9zK/DWfc7oDV2dzNCzSilFo6zJYb1FWT7sa7V+fQVtnwsbMn9QfJ6ORL77fp7aUS/Byc+Vv10/gZGEx727ZpbYcnXbEz9Wd10fP5M3RszhdUcL0VR/xv9QDmpvf6oxo8W9gslYC4GLwtnvbHdoQAv08mTqiO6u2HKG0okZtOZcwNjGOWX178N7W3RzNtX8hHx37MiOmB2tn3sPAoC78cccafrVpKSV12vte6qiLyVoFgIvR/qk8OrQhANw6fTD1DRa+Wn9AbSlN8oepY/H3cOMPy9dhMuvjyx2dMA9vFk26iacGT2BL7mmmffsh285mqC2r02IUBixSWwWsGqzKTYKzwf4pbjq8IcRGBDJqYDe+2nCA2roGteVcgp+HG3+7fiKpBUW8mbxdbTk6dsAgBPf0HMKyabfj6ezCreu+4IW9yfqEswo4GYxYpLbOu7nREJyEbgjtwp0zhlBRVcfyTYfUltIk45O6MXdAL97fuptdZ7LVlqNjJ3oHhrHquru4KaE/7x7dwc1rvyC3ukJtWZ0KZ+FMg1VbQR16D6Gd6ZPQhUE9uvLZ6j3UmbTXSwD449SxRAX48cSy7ymr0ResdRY8nF14fvhUXh81k5TSAqZ/+yHrs06qLavT4GJwwmTV1jXBLOsBcDLYfwV1pzAEgHtuuIaismpWJh9RW0qTeLq68K+50yiuquGpb9drMvpBp/2YGduTb2fcRYSXL7/atJTn927CbNXW2HZHxM3oSp3VpLaMX2BpNASj0A2h3VDSY0fw6apd1Ju01UU8R5+IMH474VrWp6Tx5d7DasvRsTNxPgEsm3Y7tyUO4D9Hd3L7hsUU1FapLatD4250pdasrR65xWoCBALb1etuKZ3GEADunTOcwtJqViRr92K7cPggRnaL5rnvkzmWV6C2HB0742p04plrpvDytdexvzCX61Z9xM78TLVldVg8nNypstSqLeMXSMwYhLMqGZs7lSGc6yV88u0uzc4lGAyCF+dMJcDTnUeXrKJcT4DXKZnbrQ8rpt+Jt7Mrt65bzKIT+/RhxHbAy8mdGnMdVg2FnlqlGYNdStVcSqcyBCEE980dQVFZNcs2ajPiCCDA04NX58/gbHklTyxbg9WqXwg6I0n+wSyffgeju8Ty1M51PLFttZ451cb4OHsikVSbtdNLkFgRKqXT6FSGADCwR1eG9orik293UV2rrcmkC+nfNZw/TB3L5pOneXvLDrXl6KiEj4sb/x0/j0f6XstXpw5z45rPya+pVFtWh8HXWVkNXNagnbkaKSUGlS7Nnc4QAB5cMJKyylo+X71HbSmX5eYhfZnVrwdvJu9g4/FTasvRUQmDEDzWfxTvjZ3DyfIiZn73CYeK8tSW1SHwc1byBZWZtLP+QyIBdSo+dkpD6BkXxoShifzv+70Ul1WrLadZhBD8bcZE+nQJ5fGl33P8bKHaknRUZHJUIkun3Y6TwcD8tZ+z8vQxtSU5PAEuSkbREg0Zgpq0tabyfCHEUSGEVQgx+DL7nRFCHBZCHBBCaOK2/MEF12IyW/jvN9pOF+Hm7MSbN83E282VX3+xgsJK7RqYTvvTwz+EldfdRb/AMB75cSWvHPhRn2xuAwGuShGaYlO5ykp+RiCQqDPJ3dYewhFgDrClBfuOk1L2l1I2axz2pGuoP3PG92VF8mHO5Go7P32ojxdv3zyT0ppaHl68kroGfWKxMxPo5sFnk25mfrc+vH7oJ3679Vt9srmVeBrdcTe6UlSvnZopQhgc0xCklClSyhO2EmNv7pl9Da4uzry5uCV+pi69uoTy0pxpHMo5yxPLvseir2Lt1LgYjbw4YjqPDxjNitPHuG3dYopq9d7j1SKEINjVn4I6DRkCRqRKCffsNYcggXVCiL1CiPvs1OYV8ffx4K6ZQ/lxfzq7j2p/8c/EHvH8fsoY1qWk8ez3yfpQQSdHCMFDfUbw1ujZHCk5y9w1i8io1M6FzVEIcQsgv65YbRnnMQpnrFKdHp+40kVFCLEBCGvirT9JKVc07pMM/E5K2eT8gBCii5QyVwgRAqwHfiOlbPK2vNEw7gMIDg4etGTJkpb+Lq2iwWzllWXHcDYaePSGHhgNVze7X1VVhZeXfQtZfJ2Ww4bsAmbFhjMtuqk/zaWoobM16DpbR1p9Ja8WnkAA/y+4O3Guijat6WwONXVuEAc4SiaPyOsRV4jusYfOMu+NlPtsJCrnGUQr7tnHjRu3t9VD81LKNm9AMjC4hfv+FcU8rrhvYmKitAfJu0/Kobe9LL9cu++qP7tp0ybbC7oCFotV/u7r1TLp6Vfkkj2HWvQZNXS2Bl1n60krK5Ijl74tkz57Sa7LTJVSalNnU6ipc3n2Jjl9829kWX3FFfe1h85DxZ/Kj1OHywZLXas+D+yRrbyWt/uQkRDCUwjhfe45MBllMvqKWCz2GRIZPagbQ3tH897SbZSUa7+kocEgeHbWZEbFx/D0qo18f8Rhp3F0bEg330CWTbuDJP9g7k9expI07a7G1xIR7iEA5NRqI6z7XNprs7R/2pq2hp3eIITIBoYD3wkh1jb+vIsQYnXjbqHAViHEQWAX8J2Uck1Ljl9ZZp8TIoTgsdvGUlPfwNtLfrRLm23FxcnIawtmMLBrFx5ftoYf9IVrOkCwuyf/m3Qz14ZF88S21ayuyFVbkuY5ZwhZNdqoa24UbgBYrA5mCFLKb6SUkVJKVyllqJRySuPPc6WU0xufp0sp+zVuvaSUz7b0+HW1ZjJP28e1YyMCuWnKAL7dcpTDaY7xT+Th4sy7t8yiZ3gIj371HT+mnVFbko4G8HR24YPx85kR04MlZZk8s2cjVj0AoVlC3AJwMTiTVXNWbSkAOBkUQ2iQ9s+vpO2VygKWfLzVbs3dM3s4wf5evPjxRswWxwjr9HJz5f3bbiA+OICHF6/kp1N6wXYdJSz19VEzmegVxn+P7ebxbd/pBXeawSgMdPUIJaNGG+lAnIUn8HMpTXuiaUPw8HRh4/eHyMu2z8IxT3cXfnfHOFIzClm8Zp9d2rQFvu5ufHj7XGIC/XnwfyvYcvK02pJ0NIBBCG71j+ax/qNYeuoIj/y4gnp9AVuTRHuEc6ZaGyMDzkallnKD1f7rSjRtCF7eLhiNBr60Yy9hzKB4Rg3sxnvLtpFTUGa3dtuKv6c7H985j/jgQB5a/K2eDE8HUObHHul7LX8ePJ7VGSe494el1DRoN8uvWsR4dqHEVEFFg/qL+1wMSliryWL/rLaaNgSDUTB19kDWf3uA/Lwyu7QphODxO8ZjNAhe+GijQy3+8vdw56M759IzLJhHl6zSo490znNvz6G8OHwaP509w50bl1DVUK+2JE0R69kFgNPVOSorAReDkoHVZLV/Sm5NGwLAgjtHArD4I/tF/4QGevPrBSPZdSSD739KsVu7tsDX3Y0Pbp9D/8hwHvt6NYt3H1Rbko5GWJDQj9dHzWRfYQ53bFhChUmvxneOWK8IANKr1DcEV+M5Q9B7CJcQEubL1NkDWbdiP2dz7Lcsf+6E/vRN6MK/P9uk6RTZTXFuonl0Qix//e4H3t68w6F6Ojrtx4yYHrw1ejaHi/O4bf1iyuq1UylMTfxdfAh08eVUVZbaUnASHhiEM/UW+2dg1bwhANx8z2gMRsFn7yXbrU2DQfDnX02m3mTm+Y82ONwF1d3FmTdvup5Z/Xrw+qbtLD6ZrSfE0wFganQS/xk7h+Olhdy07n+U1Gl/MaY96ObVlTQNGIIQAleDr24IzREU4sOMeUPYuPqg3dYlAESHB/Dg/JH8uO+Uww0dATgbjTw/awr3jBjE5twiHl78LTWmBrVl6WiA8ZHxfDh+PqcrSrll/Re6KQCJ3lFk1xRQo4H6ym5Gf2ot9k9U6BCGAHDjwlG4ujnz6Ts/2LfdKQPpnxTBy4s2kV/seLVsDQbB45NHc1NCJJtPnub2j5ZQUKmd+rE66jGySwwfjJunm0Ijid7RSCQnK9XvJfT0v5F472l2b9dhDMHP35M5tw7nx43HSD1mv4kfg0Hwl/umYLFY+cf7a7FaHWvo6BxjI4J5++aZnC4q5cb3F+vlOHUAxRT+O26uYgrrvqC4E5tConc0ACmV6q/jifeZTrT3OLu36zCGADD3thH4+Hrw4Zsb7NpuRIgfj94yht1HM/ly3X67tm1LxibG8dndC7BKyc0fLGbtsZNqS9LRAKO6xPLh+Hmcrizl9vWLKa/vnNFH3s4edPUIJaVCfUNQC4cyBE8vN26+ZxT7d6azd3uaXduePa4PowZ2460vfyQ1o8CubduSnuEhfHXfzSSGKmsVXv9hm8P2enRsx7XhMfxn7BzSyou5Y8OXnTYktadPHMcrzmCVnTMAw6EMAWDG/KGERfjz/mvrsNgx35AQgj/fOxlfLzeeens1tXWOOzkb4u3ForvmMad/L97espOHFq+kvLZzXgB0fmZsRBxvj5nN0ZJ87tr4VadcvNbTJ44qcw2ZGkl0Z28czhBcXJy4++GJnD6Zz4ZVB+zatp+3O399YBoZeSX8+/NNdm3b1rg4OfHsrEn8edo4fkw7w9z/fM6xPMft+ejYholdE3hz9CwOFuVy98avqDU77o1Pa+jjFw/AkXL7jkBoBYczBIDRk3rRo29XPn5rI7U19r2LGdIrijtmDGVF8hHW73Ds1BBCCG4b1p9FC+djtlq56b+L+WrvYYdbc6FjW6ZGJ/HqyJnsLsjm15uXY7KoU/BdDULdAglxDeBQWeecX3NIQxBCcP9jUygprrJreuxz3DdnOH3iw3nug3UU2qmIT3syoGsXlt1/K4OjI3jq2w08vmwNVXWdb7hA52euj+3Bc9dMZVPOKX67dWWnSp3dxy+ew2Vpmp1HsFjb73/TIQ0BoEefroyd0puvP9tGgZ0S353DycnIsw/PwMXZiUUb0x16PuEcAZ4evH/bDTw6bgTfHznBDf/5nEPZnXMcVUfh5sT+57OkPrFtdacpstPfL5EKc7XqeY3yaw+wu/A1dhS8TK25hLM1+/ku61d8n/1rcqp3YJW277k5rCEA3P2bSQC89+o6u7cdGujN3x+cTkFZHS987HipLZrCaDDw4JhhLFo4H4vVyi0ffsm7W3bpKS86Mff2HMpj/UexLP0Iz+xxrOy/raWvXyIAB8tSVdNQWn+K3YWvU27KoNx0mj1Fb3Ks7Au6+84l2msMe4veorje9kPWDm0IoeF+3HjXKH7ccJQDu9Lt3v6wPtFMHBjOmp9SWLax4xQ0HxgVwTcP3MakHvG8+sNP3PbREjJLytSWpaMSv+kzgru6D+LDlD28fWS72nLanSBXPyLdQ1U1hJSyrwhx78PEiFeYEvkmVQ25+LvG081nKn0C7sDbuSsFtba/5rTJEIQQLwkhjgshDgkhvhFC+DWz31QhxAkhRJoQ4sm2tHkx8++4lrAIf956aTXmBvtPfo3vH8aIfrG88tkmDpxQP3WurfB1d+OVedN5ac400gpKmP3OZyzec6hT3CHq/BIhBH8ZMpHZsb14af8WFp1wnGqCraW/fxJHyk/RYFVnOFhiwc3of/61j0sUdeafcxs5G9wxt0OJzbb2ENYDvaWUfYFU4A8X7yCEMAJvAdOAnsDNQoiebWz3PK5uzjz4u2lkpheyfPEOWx22xRiE4G8PTiM82IcnX//WIfMdNYcQguv7dmflr2+nX2QYf121kbs/XUp2qf2zMOqoi0EIXrp2OhMj4/nLznV8e9rxkj1eDQP8k6i3mjhWrs6q5QDXJM7W7KOqIY9yUwYN1mpKTWlsyv0DW88+Q7kpg3CPITZvt02GIKVcJ6U8V6R1BxDZxG5DgTQpZbqU0gQsBma1pd2LuWZ0EsNGJXL8iDp36D6ebrz021k0mC3sOtrxityH+3rz4R1z+euMCRzKyeeNTR1/2EDnUpwNRt4cPYshIZF8mXawQ/cW+/omEOoWSIVZnVookZ7DcXcKZGv+P9hd+Dpezl0YHfZ3At2642RwZVDwQ4S497H530DY6oBCiG+BL6WUn13083nAVCnlvY2vbweGSSkfbuY49wH3Nb7sDRyxicD2IwgoUltEC9B12hZdp23RddqOJCmld2s+6HSlHYQQG4CwJt76k5RyReM+fwLMwOdNHaKJnzXrQlLK94D3Go+7R0o5+Eoa1cQRNIKu09boOm2LrtN2CCH2tPazVzQEKeXEKzR+JzADmCCb7m5kA10veB0J5F6NSB0dHR2d9qetUUZTgd8DM6WUzU157wYShBCxQggX4CZgZVva1dHR0dGxPW2NMnoT8AbWCyEOCCHeBRBCdBFCrAZonHR+GFgLpABLpJRHW3j899qozx44gkbQddoaXadt0XXajlZrtNmkso6Ojo6OY+PQK5V1dHR0dGyHbgg6Ojo6OoCGDEELaTBaghBivhDiqBDCKoRoNvxMCHFGCHG4cW6l1WFgreUqdKp9PgOEEOuFECcbH/2b2U+V83ml8yMUXm98/5AQYqC9tF2FxrFCiPLGc3dACPEXe2ts1PGhEKJACNHk2iItnMtGHVfSqfr5FEJ0FUJsEkKkNP6fP9rEPld/PqWUmtiAyYBT4/N/Av9sYh8jcAqIA1yAg0BPO+vsASQBycDgy+x3BghS8XxeUadGzueLwJONz59s6u+u1vlsyfkBpgPfo6y3uQbYqUGNY4FVan0XL9AxGhgIHGnmfVXP5VXoVP18AuHAwMbn3iipg9r83dRMD0FqJA3GlZBSpkgpNV8qrYU6VT+fje190vj8E2C2ndu/HC05P7OAT6XCDsBPCBGuMY2aQEq5BSi5zC5qn0ugRTpVR0qZJ6Xc1/i8EiWCM+Ki3a76fGrGEC7ibhRnu5gIIOuC19lcehK0ggTWCSH2Nqbj0CJaOJ+hUso8UL7kQEgz+6lxPltyftQ+hy1tf7gQ4qAQ4nshRC/7SLtq1D6XV4NmzqcQIgYYAOy86K2rPp9XXKlsS+ydBqO1tERnC7hWSpkrhAhBWadxvPHOw2bYQKfq5/MqDtPu57MJWnJ+7HIOL0NL2t8HREspq4QQ04HlQEJ7C2sFap/LlqKZ8ymE8AKWAr+VUlZc/HYTH7ns+bSrIUgHSYNxJZ0tPEZu42OBEOIblK69TS9gNtCp+vkUQuQLIcKllHmN3dmCZo7R7uezCVpyftROzXLF9i+8UEgpVwsh3hZCBEkptZakTe1z2SK0cj6FEM4oZvC5lHJZE7tc9fnUzJCR6EBpMIQQnkII73PPUSbMtZi1VQvncyVwZ+PzO4FLejYqns+WnJ+VwB2NER3XAOXnhsDsxBU1CiHChBCi8flQlP/7YjtqbClqn8sWoYXz2dj+B0CKlPKVZna7+vOp5kz5RTPiaSjjXQcat3cbf94FWH3RzHkqSmTFn1TQeQOK89YD+cDai3WiRHwcbNyOalWnRs5nILARONn4GKCl89nU+QEeAB5ofC5QCkCdAg5zmcgzFTU+3HjeDqIEbIywt8ZGHV8AeUBD43fzHq2dyxbqVP18AiNRhn8OXXDNnN7W86mnrtDR0dHRATQ0ZKSjo6Ojoy66Iejo6OjoALoh6Ojo6Og0ohuCjo6Ojg6gG4KOjo6OTiO6Iejo6OjoALoh6Ojo6Og08v8BV1Ayn9vJfp4AAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
}
],
"metadata": {
"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.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}