![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open"){kind=icon}
![\"Download\"](\"https://img.shields.io/badge/Github-Download-blue.svg\" "\\"Download"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"id": "f3e9a52a",
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.9 Algorithms Homework 4: Interior Point Methods](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9-Algorithms-Homework-4:-Interior-Point-Methods)",
"section": "4.9 Algorithms Homework 4: Interior Point Methods"
}
},
"source": [
"# 4.9 Algorithms Homework 4: Interior Point Methods"
]
},
{
"cell_type": "markdown",
"id": "ba74c034",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.9.1 Tips and Tricks](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1-Tips-and-Tricks)",
"section": "4.9.1 Tips and Tricks"
}
},
"source": [
"## 4.9.1 Tips and Tricks"
]
},
{
"cell_type": "markdown",
"id": "a8715f98",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.1 Background](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.1-Background)",
"section": "4.9.1.1 Background"
}
},
"source": [
"### 4.9.1.1 Background"
]
},
{
"cell_type": "markdown",
"id": "4844e9f9",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.1 Background](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.1-Background)",
"section": "4.9.1.1 Background"
}
},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "11a91e88",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.2 Problem Formulation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.2-Problem-Formulation)",
"section": "4.9.1.2 Problem Formulation"
}
},
"source": [
"### 4.9.1.2 Problem Formulation"
]
},
{
"cell_type": "markdown",
"id": "164bb509",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.2 Problem Formulation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.2-Problem-Formulation)",
"section": "4.9.1.2 Problem Formulation"
}
},
"source": [
"Consider the following nonlinear program:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\t\\min_{x} \\quad & f(x) \\\\\n",
"\t\\mathrm{s.t.} \\quad & c(x) = 0 \\\\\n",
"\t& x_i \\geq 0, \\quad i \\in \\mathcal{I}\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $x \\in \\mathbb{R}^{n}$, $f(x): \\mathbb{R}^{n} \\rightarrow \\mathbb{R}$, $c(x): \\mathbb{R}^{n} \\rightarrow \\mathbb{R}^{m}$ and $|\\mathcal{I} |$ = $r$ (i.e., there are $r$ variables with a lower bound and $r \\leq n$). This is an extension of (6.48) in Biegler (2010).\n",
"\n",
"This has the corresponding log-barrier approximation:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\t\\min_{x} \\quad & \\phi_{\\mu_l}(x) := f(x) - \\mu_l \\sum_{i \\in \\mathcal{I}} \\log(x_i) \\\\\n",
"\t\\mathrm{s.t.} \\quad & c(x) = 0\n",
"\\end{align}\n",
"$$\n",
"\n",
"which is an extension of (6.49) in Biegler (2010).\n",
"\n",
"Let the $n \\times r$ matrix $G$ encode which variables are bounded. If variable $i$ corresponds to the $j$th bound, then $G_{i,j} = 1$ and otherwise $G_{i,j} = 0$. $G$ is assembled as follows:\n",
"\n",
"\n",
"\n",
"Notice that $G$ is the gradient of $x \\geq 0$. $G^T G = I$ but the converse does not hold unless $r = n$."
]
},
{
"cell_type": "markdown",
"id": "c2fed8b8",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.3 Reformulation Example](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.3-Reformulation-Example)",
"section": "4.9.1.3 Reformulation Example"
}
},
"source": [
"### 4.9.1.3 Reformulation Example"
]
},
{
"cell_type": "markdown",
"id": "fae61884",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.3 Reformulation Example](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.3-Reformulation-Example)",
"section": "4.9.1.3 Reformulation Example"
}
},
"source": [
"Start with:\n",
"\n",
"$$\n",
"\t\\begin{align}\n",
"\t\t\\min_{x} \\quad & x_1^2 + x_2 \\\\\n",
"\t\t\\mathrm{s.t.} \\quad & x_1 + x_2 = 1 \\\\\n",
"\t\t& x_1 + 1 \\geq 0\n",
"\t\\end{align}\n",
"$$\n",
"\n",
"Add slack variable $x_3$ and convert the inequality constraint to an equality constraint and bound:\n",
"\n",
"$$\n",
"\t\\begin{align}\n",
"\t\t\\min_{x} \\quad & x_1^2 + x_2 \\\\\n",
"\t\t\\mathrm{s.t.} \\quad & x_1 + x_2 = 1 \\\\\n",
"\t\t& x_1 + 1 - x_3 = 0 \\\\\n",
"\t\t& x_3 \\geq 0\n",
"\t\\end{align}\n",
"$$\n",
"\n",
"Now assemble $G$:\n",
"\n",
"$$\n",
"\tG = \\begin{bmatrix}\n",
"\t\t0 \\\\\n",
"\t\t0 \\\\\n",
"\t\t1\n",
"\t\\end{bmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "00e568fa",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.4 Primal Dual Optimality Conditions](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.4-Primal-Dual-Optimality-Conditions)",
"section": "4.9.1.4 Primal Dual Optimality Conditions"
}
},
"source": [
"### 4.9.1.4 Primal Dual Optimality Conditions"
]
},
{
"cell_type": "markdown",
"id": "b923aa13",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.1.4 Primal Dual Optimality Conditions](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.1.4-Primal-Dual-Optimality-Conditions)",
"section": "4.9.1.4 Primal Dual Optimality Conditions"
}
},
"source": [
"$$\n",
"\\newcommand{\\dims}[3]{\\underbrace{#1}_{#2 \\times #3}}\n",
"\\newcommand{\\I}[0]{\\mathcal{I}}\n",
"\\newcommand{\\dfx}[0]{\\dims{\\nabla f(x^k)}{n}{1}}\n",
"\\newcommand{\\cx}[0]{\\dims{c(x^k)}{m}{1}}\n",
"\\newcommand{\\dcx}[0]{\\dims{\\nabla c(x^k)}{n}{m}}\n",
"\\newcommand{\\dcxT}[0]{\\dims{\\nabla_x c(x^k)^T}{m}{n}}\n",
"\\newcommand{\\W}[0]{\\dims{W^k}{n}{n}}\n",
"\\newcommand{\\G}[0]{\\dims{G}{n}{r}}\n",
"\\newcommand{\\GT}[0]{\\dims{G^T}{r}{n}}\n",
"\\newcommand{\\vk}[0]{\\dims{v^k}{m}{1}}\n",
"\\newcommand{\\Xinv}[0]{\\dims{(\\hat{X}^k)^{-1}}{r}{r}}\n",
"\\newcommand{\\X}[0]{\\dims{\\hat{X}^k}{r}{r}}\n",
"\\newcommand{\\e}[0]{\\dims{e}{r}{1}}\n",
"\\newcommand{\\uk}[0]{\\dims{u^k}{r}{1}}\n",
"\\newcommand{\\dx}[0]{\\dims{d_x^k}{n}{1}}\n",
"\\newcommand{\\dv}[0]{\\dims{d_v^k}{m}{1}}\n",
"\\newcommand{\\du}[0]{\\dims{d_u^k}{r}{1}}\n",
"\\newcommand{\\Uk}[0]{\\dims{U^k}{r}{r}}\n",
"\\newcommand{\\Sk}[0]{\\dims{\\Sigma^k}{n}{n}}\n",
"\\newcommand{\\dphi}[0]{\\dims{\\nabla \\phi_{\\mu_l}}{n}{1}}\n",
"$$\n",
"\n",
"Next we extend (6.51) in Biegler (2010):\n",
" \n",
"$$\n",
"\\begin{gather}\n",
"\t\\dfx ~+~ \\dcx \\vk ~-~ \\G \\uk = \\dims{0}{n}{1} \\\\\n",
"\t\\X \\uk = \\mu \\e \\\\\n",
"\t\\cx = \\dims{0}{m}{1}\n",
"\\end{gather}\n",
"$$\n",
"\n",
"We now extend (6.56) in Biegler (2010):\n",
"\n",
"$$\n",
"\t\\begin{bmatrix}\n",
"\t\t\\W & \\dcx & \\G \\\\\n",
"\t\t\\dcxT & \\dims{0}{m}{m} & \\dims{0}{m}{r} \\\\\n",
"\t\t\\Uk \\GT & \\dims{0}{r}{m} & \\X\n",
"\t\\end{bmatrix} \n",
"\t\\begin{bmatrix}\n",
"\t\t\\dx \\\\\n",
"\t\t\\dv \\\\\n",
"\t\t\\du\n",
"\t\\end{bmatrix} = -\n",
"\t\\begin{bmatrix}\n",
"\t\\dfx ~+~ \\dcx \\vk ~-~ \\G \\uk \\\\\n",
"\t\\cx \\\\\n",
"\t\\X \\uk ~-~ \\mu_l \\e\n",
"\t\\end{bmatrix}\n",
"$$\n",
"\n",
"where $U^k = \\mathrm{diag}\\{u^k\\}$ and $W^k = \\nabla_{xx} L(x^k,v^k)$. Notice that $W^k$ does NOT include a contribution from the barrier term:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\t\\W = \\dims{\\nabla^2 f(x^k)}{n}{n} + \\sum_{j=1}^{m} \\left( \\dims{\\nabla^2 c_{j}(x^k)}{n}{n} \\underbrace{v^k_{j}}_{\\mathrm{scalar}} \\right)\n",
"\\end{equation}\n",
"$$\n",
"\n",
"We can verify that $W^k$ does not include $\\hat{X}$ by showing that the KKT system above is a Newton step to solve the nonlinear system for the primal dual conditions.\n",
"\n",
"\n",
"The Newton step can be simplified, similar to (6.57) and (6.58) in Biegler (2010):\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\t\\begin{bmatrix}\n",
"\t\t\\W + \\Sk & \\dcx \\\\\n",
"\t\t\\dcxT & 0\n",
"\t\\end{bmatrix} \n",
"\t\\begin{bmatrix}\n",
"\t\t\\dx \\\\\n",
"\t\t\\dv\n",
"\t\\end{bmatrix} = -\n",
"\t\\begin{bmatrix}\n",
"\t\\dphi ~+~ \\dcx \\vk \\\\\n",
"\t\\cx\n",
"\t\\end{bmatrix}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\begin{equation} \\label{eq:du}\n",
"\t\\du = \\mu_l \\Xinv \\e ~-~ \\uk ~-~ \\GT \\Sk \\dx\n",
"\\end{equation}\n",
"$$\n",
"\n",
"where\n",
"$$\n",
"\\begin{equation} \\label{eq:Sk}\n",
"\t\\Sk = \\G \\Xinv \\Uk \\GT\n",
"\\end{equation}\n",
"$$\n",
"\n",
"and \n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\t\\dphi = \\dfx - \\mu_l \\G \\Xinv \\dims{e}{r}{1}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Notice that the equation for $\\du$ simplifies by substituting $G^T G = I$:\n",
"\n",
"$$\n",
"\t\\du = \\mu_l \\Xinv \\e ~-~ \\uk ~-~ \\Xinv \\Uk \\GT \\dx\n",
"$$\n",
"\n",
"Finally, inertia correction can be applied to simplified KKT step similar to (6.59) in Biegler (2010). "
]
},
{
"cell_type": "markdown",
"id": "1ca1ed21",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.9.2 Basic Interior Point Method for Inequality and Equality Constraint NLPs](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2-Basic-Interior-Point-Method-for-Inequality-and-Equality-Constraint-NLPs)",
"section": "4.9.2 Basic Interior Point Method for Inequality and Equality Constraint NLPs"
}
},
"source": [
"## 4.9.2 Basic Interior Point Method for Inequality and Equality Constraint NLPs"
]
},
{
"cell_type": "markdown",
"id": "3fd841f9",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.9.2 Basic Interior Point Method for Inequality and Equality Constraint NLPs](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2-Basic-Interior-Point-Method-for-Inequality-and-Equality-Constraint-NLPs)",
"section": "4.9.2 Basic Interior Point Method for Inequality and Equality Constraint NLPs"
}
},
"source": [
"Implement a basic interior point method for inequality and equality constrained nonlinear programs. See pg. 154 – 155 in Biegler (2010). You may skip the line search, i.e., always take a full step."
]
},
{
"cell_type": "markdown",
"id": "a7a244c5",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.1 Pseudocode](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.1-Pseudocode)",
"section": "4.9.2.1 Pseudocode"
}
},
"source": [
"### 4.9.2.1 Pseudocode"
]
},
{
"cell_type": "markdown",
"id": "23e86ba7",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.1 Pseudocode](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.1-Pseudocode)",
"section": "4.9.2.1 Pseudocode"
}
},
"source": [
"Write detailed pseudocode on paper or a whiteboard. Scan/take a photo and turn in."
]
},
{
"cell_type": "markdown",
"id": "d1ec363f",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"source": [
"### 4.9.2.2 Python Implementation"
]
},
{
"cell_type": "markdown",
"id": "2ad37a24",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"source": [
"Implement in Python. Hints: Reuse code from Algorithm 5.2 example."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ec182a38",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"outputs": [],
"source": [
"### Load Python libraries\n",
"import numpy as np\n",
"from scipy import linalg\n",
"\n",
"### Define helper functions\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",
"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",
"## Calculate gradient with central finite difference\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 Hessian using central finite difference\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\")\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f51b5c25",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"outputs": [],
"source": [
"## Assemble KKT matrix (equality constrained only)\n",
"def assemble_check_KKT(W,Sk,A,deltaA,deltaW,verbose):\n",
" \n",
" # YOUR SOLUTION HERE\n",
" \n",
" return KKT,inertia_correct,pos_ev,neg_ev,zero_ev"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "ae9be074",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"outputs": [],
"source": [
"def barrier_subproblem(x0,v0,u0,calc_f,calc_c,var_bounds,mu,eps,max_iter=100,verbose=False):\n",
" '''\n",
" Basic Full Space Newton Method for Solving Barrier Subproblem Constrained NLP\n",
" \n",
" Input:\n",
" x0 - starting point (vector)\n",
" calc_f - function to calculate objective (returns scalar)\n",
" calc_c - function to calculate constraints (returns vector)\n",
" var_bounds - list of indicies for variables with lower bound\n",
" mu - barrier penalty\n",
" eps - tolerance for termination\n",
" \n",
" Histories (stored for debugging):\n",
" x - history of steps (primal variables)\n",
" v - history of steps (duals for constraints)\n",
" u - history of steps (duals for bounds)\n",
" f - history of objective evaluations\n",
" c - 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",
" S - history of sigma matrix\n",
" \n",
" Outputs:\n",
" x - final value for primal variable\n",
" v - final value for constraint duals\n",
" u - final value for bound duals\n",
" \n",
" Notes:\n",
" 1. For simplicity, central finite difference is used \n",
" for all gradient calculations.\n",
" ''' \n",
" \n",
" ### Specifics for Algorithm 5.2\n",
" # Tuning parameters\n",
" delta_bar_W_min = 1E-20\n",
" delta_bar_W_0 = 1E-4\n",
" delta_bar_W_max = 1E40\n",
" delta_bar_A = 1E-8\n",
" kappa_u = 8\n",
" kappa_l = 1/3\n",
" \n",
" # Declare iteration histories as empty lists\n",
" x = []\n",
" v = []\n",
" u = []\n",
" f = []\n",
" L = []\n",
" c = []\n",
" df = []\n",
" dL = []\n",
" A = []\n",
" W = []\n",
"\n",
" # YOUR SOLUTION HERE"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6549fabb",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.2.2 Python Implementation](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.2.2-Python-Implementation)",
"section": "4.9.2.2 Python Implementation"
}
},
"outputs": [],
"source": [
"def interior_point(x0,calc_f,calc_c,var_bounds,max_iter=20):\n",
" \n",
" \n",
" # YOUR SOLUTION HERE\n",
" \n",
" return x, v, u, mu, E"
]
},
{
"cell_type": "markdown",
"id": "4fb21995",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.9.3 Test Problems](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3-Test-Problems)",
"section": "4.9.3 Test Problems"
}
},
"source": [
"## 4.9.3 Test Problems"
]
},
{
"cell_type": "markdown",
"id": "2e7ab8ef",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.3.1 Problem 1: Convex](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3.1-Problem-1:-Convex)",
"section": "4.9.3.1 Problem 1: Convex"
}
},
"source": [
"### 4.9.3.1 Problem 1: Convex"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "4c290f53",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.3.1 Problem 1: Convex](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3.1-Problem-1:-Convex)",
"section": "4.9.3.1 Problem 1: Convex"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 3.0000e+00 \t1.0000e+00 \t 2.00e+00 \t 9.06e-01 \t 2.00e+00 \t 3.56e-04 \t 0.00e+00 \t 0.00e+00\n",
"1 \t 1.0996e+00 \t1.3978e-10 \t 3.20e-04 \t 5.03e-04 \t 1.78e-04 \t 3.56e-04 \t 0.00e+00 \t 0.00e+00\n",
"2 \t 1.1000e+00 \t3.9524e-14 \t 1.26e-07 \t 1.70e-07 \t 3.94e-08 \t 5.91e-08 \t 0.00e+00 \t 0.00e+00\n",
"3 \t 1.1000e+00 \t0.0000e+00 \t 6.87e-11 \t 1.48e-11 \t 1.11e-10 \t 1.04e-10 \t 0.00e+00 \t 0.00e+00\n"
]
}
],
"source": [
"f = lambda x: x[0] + 2*x[1]\n",
"c = lambda x: (x[0] + x[1] - 1)*np.ones(1)\n",
"\n",
"# Indices of variables with lower bound of zero\n",
"vb = [1]\n",
"\n",
"x0 = np.ones(2)\n",
"u0 = np.ones(1)\n",
"v0 = np.ones(1)\n",
"\n",
"x_, v_, u_, E_ = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,verbose=False)"
]
},
{
"cell_type": "markdown",
"id": "1f44a1a9",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.3.2 Problem 2: Convex](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3.2-Problem-2:-Convex)",
"section": "4.9.3.2 Problem 2: Convex"
}
},
"source": [
"### 4.9.3.2 Problem 2: Convex"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "08e17352",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.3.2 Problem 2: Convex](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3.2-Problem-2:-Convex)",
"section": "4.9.3.2 Problem 2: Convex"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"*** Barrier Subproblem 0 ***\n",
"mu = 10.0\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 3.0000e+00 \t1.0000e+00 \t 9.00e+00 \t 1.35e+01 \t 2.00e+00 \t 8.45e-03 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 1 ***\n",
"mu = 2.0\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.1008e+01 \t1.3978e-10 \t 7.92e+00 \t 1.14e+01 \t 4.22e-03 \t 8.45e-03 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 2 ***\n",
"mu = 0.4\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 2.9318e+00 \t0.0000e+00 \t 1.53e+00 \t 2.17e+00 \t 1.70e-04 \t 3.40e-04 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 3 ***\n",
"mu = 0.08000000000000002\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.3993e+00 \t2.5520e-10 \t 3.19e-01 \t 4.51e-01 \t 1.70e-04 \t 3.40e-04 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 4 ***\n",
"mu = 0.016000000000000004\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.0801e+00 \t5.3163e-11 \t 6.41e-02 \t 9.07e-02 \t 1.07e-05 \t 1.78e-05 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 5 ***\n",
"mu = 0.0020238577025077633\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.0160e+00 \t0.0000e+00 \t 1.40e-02 \t 1.98e-02 \t 1.15e-05 \t 1.86e-05 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 6 ***\n",
"mu = 9.104790579399288e-05\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.0020e+00 \t7.7571e-13 \t 1.93e-03 \t 2.73e-03 \t 7.76e-07 \t 8.83e-07 \t 0.00e+00 \t 0.00e+00\n",
"1 \t 1.0001e+00 \t0.0000e+00 \t 1.07e-07 \t 2.43e-09 \t 0.00e+00 \t 1.07e-07 \t 0.00e+00 \t 0.00e+00\n",
"\n",
"*** Barrier Subproblem 7 ***\n",
"mu = 8.687702517211205e-07\n",
"Iter. \tf(x) \t\t||c(x)|| \tE \t\t||dx|| \t\t||dv|| \t\t||du|| \t\tdelta_A \tdelta_W\n",
"0 \t 1.0001e+00 \t0.0000e+00 \t 9.02e-05 \t 1.28e-04 \t 1.11e-10 \t 5.12e-09 \t 0.00e+00 \t 0.00e+00\n",
"1 \t 1.0000e+00 \t9.9920e-15 \t 4.95e-09 \t 6.52e-13 \t 3.88e-17 \t 4.95e-09 \t 0.00e+00 \t 0.00e+00\n"
]
}
],
"source": [
"f = lambda x: x[0] + 2*x[1]\n",
"c = lambda x: (x[0] + x[1] - 1)*np.ones(1)\n",
"\n",
"# Indices of variables with lower bound of zero\n",
"vb = [1]\n",
"\n",
"x0 = np.ones(2)\n",
"#u0 = np.ones(1)\n",
"#v0 = np.ones(1)\n",
"\n",
"# x_, v_, u_, E_ = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,verbose=False)\n",
"\n",
"x, v, u, mu, E = interior_point(x0,f,c,vb)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "2efe01ca",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.9.3.2 Problem 2: Convex](https://ndcbe.github.io/CBE60499/04.09-Algorithms4.html#4.9.3.2-Problem-2:-Convex)",
"section": "4.9.3.2 Problem 2: Convex"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy70lEQVR4nO3deXxU9b3/8dcnO1nYsgBJWGXfEwMIiHVDQREktnXHrVLvFdvb3ntb79L2ent/9/bXX3t7b6tdEFGw1KUKFisq1qqgsoWwb4JAICFA2AmQ/fP740xkDJPJkMzMmSSf5+Mxj5kz58w5bxTmM+ec7yKqijHGGNOYKLcDGGOMiWxWKIwxxvhlhcIYY4xfViiMMcb4ZYXCGGOMXzFuBwiFtLQ07dOnj9sxjDGm1Vi/fv0xVU33ta5NFoo+ffpQUFDgdgxjjGk1RKSosXV26ckYY4xfViiMMcb4ZYXCGGOMX23yHoUxxoRSdXU1xcXFVFRUuB3lsiUkJJCdnU1sbGzAn7FCYYwxl6m4uJiUlBT69OmDiLgdJ2CqyvHjxykuLqZv374Bfy7kl55EZL6IHBWRrV7vdRWR90Rkt+e5SyOfnSIiu0Rkj4g8GeqsxhgTiIqKClJTU1tVkQAQEVJTUy/7TCgc9yheAKY0eO9J4H1VHQC871n+EhGJBp4BpgJDgbtFZGhooxpjTGBaW5Go15zcIS8UqroCONHg7RnAAs/rBcDtPj46FtijqntVtQp42fO59mfnMijb5XYKY0w75Varp26qWgrgec7wsU0WcNBrudjznk8iMltECkSkoKysLKhhXXWmFF65D974W7C5Q4wxfixYsIABAwYwYMAAFixY0PQHAhTJzWN9nR81+k2pqnNVNU9V89LTffZCb502/B60FkoK4OAat9MYYyLUiRMneOqpp1izZg1r167lqaee4uTJk0HZt1uF4oiI9ADwPB/1sU0x0NNrORs4FIZskaOuFta/AL3GQ4cu8Omv3E5kjIkA69atY+TIkVRUVHDu3DmGDRvGM888w+TJk+natStdunRh8uTJvPPOO0E5nlvNY5cCDwA/8Tz/ycc264ABItIXKAHuAu4JW8JIsPs9OFMMU/4TDm+BFT+DY3sgrb/byYwxHk+9uY3th84EdZ9DMzvyo9uGNbp+zJgxTJ8+nX/913/lwoUL3HfffcTGxtKz58Xf1tnZ2ZSUlAQlTziax74ErAIGiUixiDyCUyAmi8huYLJnGRHJFJFlAKpaA8wB3gV2AK+q6rZQ540oBfMhuRsMugXGzoboWFj9jNupjDER4Ic//CHvvfceBQUFfO9730N93MMMVsuskJ9RqOrdjay6wce2h4BbvJaXActCFC2ynToAu5fDNf/gFIjkDBh1F2z8A1z3L5CU5nZCYwz4/eUfSidOnKC8vJzq6moqKirIzs7mww8//GJ9cXEx1157bVCOFck3s9u39QtABHIfuPje+DlQUwHr5rmXyxgTEWbPns2Pf/xj7r33Xr7//e9z8803s3z5ck6ePMnJkydZvnw5N998c1COZUN4RKLaaihcCANuhs5e9/PTB8HAKbB2Lkz8NsR2cC+jMcY1CxcuJCYmhnvuuYfa2lomTJjAxo0b+cEPfsCYMWMA59JU165dg3I8KxSRaOdbcO4o5D186boJT8ALt8Kml3yvN8a0ebNmzWLWrFkAREdHs2bNxabzDz8c/O8Fu/QUiQrmQ6de0P+S2zjQeyJk5sCnT0NdXfizGWPaHSsUkebYHtj3EVz5AERFX7pexDmrOPE5fPZ2+PMZY9odKxSRZv3zEBUDOfc3vs2QGc4Zh3XAM8aEgRWKSFJdARsXweBpkNKt8e2iY2D838KBVXBwXfjyGWPaJSsUkWT7n+DCycBuUufcBwmdYJWdVRhjQssKRSQpmA+p/aHvNU1vG5/iFJQdb8KJvaHPZoxpt6xQRIoj2+DgarjyIeeGdSDGfhMkGlb/JrTZjDGtwpQpU+jcuTPTpk0L6n6tUESKguchOh5GX8a4hx17wMivO0ORn284N5Qxpr35x3/8R1588cWg79cKRSSoLIdNL8OwmZB4mT0px8+B6vNQ8FxoshljIo6vYca3bt3KDTfcQEpKStCPZz2zI8HW16HqbPN6WncbCv1vhDVzYfwTEJsQ/HzGmMa9/aQzDUAwdR8BU3/S6Gpfw4wPHz48uBm82BlFJCiYDxnDoOfY5n1+whPOkB9bXg1uLmNMxGo4zHgo2RmF20oKoXQj3PKzwG9iN9T3K84vkE+fhtH3QZTVf2PCxs8v/1BqOMx4UlJSyI5l3yhuK5gPsUkw8s7m70MEJnwLju2CPe8FL5sxJmI1HGY8lKxQuOnCKef+xIivQkLHlu1r2EzomGXDehjTDngPM/7kk0+ybt06/vrXvzJp0iS+9rWv8f7775Odnc27774blOO5dulJRAYBr3i91Q/4oar+j9c21+LMp73P89ZiVf33MEUMvc2vOi2W8h5q+b6iY+Gqv4Hl/+pczsrKbfk+jTERqbFhxq+//vqQHM+1MwpV3aWqo1V1NHAlcB5Y4mPTlfXbtakioepcdsrMdYYND4bcByC+I6x6Ojj7M8YYIufS0w3A56pa5HaQsDmwGsp2BHfyoYSOzvDk296Ak+3nP6UxJrQipVDcBbzUyLrxIrJJRN4WkUZnMReR2SJSICIFZWVloUkZTAXzIb4TDM8P7n7HPebc3F7z2+Du1xjzJarqdoRmaU5u1wuFiMQB04E/+lhdCPRW1VHAr4A3GtuPqs5V1TxVzUtPTw9J1qA5dxy2vwGj7oK4IDdp65QNw++A9QuckWiNMUGXkJDA8ePHW12xUFWOHz9OQsLldcyNhH4UU4FCVT3ScIWqnvF6vUxEfi0iaap6LKwJg23jIqitCs5NbF/Gz4HNr8D6F+Dq74TmGMa0Y9nZ2RQXF9Mqrl40kJCQQHZ29mV9JhIKxd00ctlJRLoDR1RVRWQszhnQ8XCGC7q6OmcWu14TIGNIaI7RYyT0uxZW/xauehxi4kJzHGPaqdjYWPr27et2jLBx9dKTiCQCk4HFXu89JiKPeRa/CmwVkU3AL4G7tLWd6zW07yNn/ohg3sT2ZcITUH4Ytr4W2uMYY9o8V88oVPU8kNrgvd96vX4aaFttPQvmQ2IqDJ0e2uNccQNkDHU64I26u/nDgxhj2j3Xb2a3K2dKYedbMPpeiIkP7bFEnLOKo9vh8/dDeyxjTJtmhSKcNvwetBaufDA8xxv+VUjpYcN6GGNaxApFuNTVOq2Q+l0HqVeE55gxcTDum7D3QyjdHJ5jGmPaHCsU4bL7PThTHPqb2A1d+RDEJdtZhTGm2axQhEvBfEjuDoOmhve4HTpD7ixnlNrTxeE9tjGmTbBCEQ6nDsDu5c4XdnRs+I8/ztPaePVvwn9sY0yrZ4UiHNYvcFoh5c5y5/hdesOw250cFafdyWCMabWsUIRabTUULoQBN0Pnnu7lGD8Hqs46xcIYYy6DFYpQ2/kWnDsa/pvYDWXlQp9JzuWnmip3sxhjWhUrFKFWMB869YL+N7idxOmAd/YQbPM1P5QxxvhmhSKUju1xxna68gGIinY7DfSfDGmDnKayrXzILGNM+FihCKX1z0NUDOTc73YSR1QUTJgDR7Y4nfCMMSYAVihCpbrCmXdi8DRI6eZ2motGfB2SMqwDnjEmYFYoQmX7n5wZ5ty+id1QbAKMm+0MFHh4q9tpjDGtgBWKUCmYD6n9oe81bie5VN4jEJsIq55xO4kxphWwQhEKR7bBwdXOOEuROA9EYlfIuQ+2/BHOHHI7jTEmwlmhCIWC5yE6Hkbf43aSxl31N86Q52t+53YSY0yEc3sq1P0iskVENopIgY/1IiK/FJE9IrJZRHLdyHlZKsth08swbKbzyz1Sde0HQ25zilrlWbfTGGMiWCScUVynqqNVNc/HuqnAAM9jNhD5o9ptfd0ZKiPSbmL7MuFbUHkaCl90O4kxJoJFQqHwZwawUB2rgc4i0sPtUH4VzIeMYdBzrNtJmpadB73Gw+pfQ22N22mMMRHK7UKhwHIRWS8is32szwIOei0Xe967hIjMFpECESkoKysLQdQAlBRC6UbIi9Cb2L5MeAJOH4Ttb7idxBgTodwuFBNVNRfnEtPjItKwLamvb1ufY0+o6lxVzVPVvPT09GDnDEzBfIhNgpF3unP85hg4FbpeAZ/+0ob1MMb45GqhUNVDnuejwBKg4fWaYsB7bO5sIDLbc1445dyfGPFVSOjodprA1Q/rUboJ9n/sdhpjTARyrVCISJKIpNS/Bm4CGnYVXgrM8rR+ugo4raqlYY4amM2vQvV557JTazPqbkhMtWE9jDE+uXlG0Q34WEQ2AWuBt1T1HRF5TEQ8c3eyDNgL7AGeBf7WnahNUHUuO2XmQmaO22kuX2wHGDsbdr8LR3e6ncYYE2Fi3Dqwqu4FRvl4/7derxV4PJy5muXAaijbAdOfdjtJ8435Bnz8C1j1NMxoxX8OY0zQuX0zu20omA/xnWB4vttJmi8pzelJvvkVOHvE7TTGmAhihaKlzh13mpaOugviktxO0zJXPe7M8b12rttJjDERxApFS21cBLVVrfMmdkNp/WHwrbBuHlSdczuNMSZCNFkoROSnItJRRGJF5H0ROSYi94UjXMSrq3Nmses1ATKGuJ0mOCY8ARWnYMMit5MYYyJEIGcUN6nqGWAaTr+GgcA/hjRVa7HvIzixt3WM6xSonuMge4xzU7uu1u00xpgIEEihiPU83wK8pKonQpindSmY7/Q/GDrd7STBI+KcVZwqgh1vup3GGBMBAikUb4rITiAPeF9E0oGK0MZqBc6Uws63YPS9EBPvdprgGjwNuvSxYT2MMUAAhUJVnwTGA3mqWg2cxxnVtX3b8Htn4p8rH3Q7SfBFRcP4OVCy3ukjYoxp1wK5mZ2I0+mtfi6ITJyzi/arrhbWvwD9roPUK9xOExqj74EOXWxYD2NMQJeengeqgAme5WLgP0KWqDXY/R6cKW5bN7EbiktyemvvWgbHdrudxhjjokAKxRWq+lOgGkBVL+B7+O/2o2A+JHeHQVPdThJaY2dDdCysesbtJMYYFwVSKKpEpAOeeSBE5AqgMqSpItmpA7B7OeTOcr5E27LkDKfH+aaXoNylyaCMMa4LpFD8CHgH6Ckii4D3ge+FNFUkW7/AaUKaO8vtJOExfg7UVDi9tY0x7VIgrZ7eA/KBB4GXcFo/fRjaWBGqthoKF8KAm6Fzz6a3bwvSB8HAKbDuWag673YaY4wLGi0UIpJb/wB6A6U4s8v18rzX/ux8C84dbds3sX2Z8AScP+5cgjLGtDv+5qP4uZ91Clwf5CyRr2A+dOoF/W9wO0l49Z7oTMi06hmn30hUtNuJjDFh1GihUNXrwhkk4h3b44ztdP0P2t8XZf2wHq89DLvehiHT3E5kjAmjQDrcJYjId0VksYi8LiJ/JyIJLT2wiPQUkQ9EZIeIbBORb/vY5loROS0iGz2PH7b0uM22/nmIioGc+12L4KohM6BzL+uAZ0w7FEirp4XAMOBXwNPAUODFIBy7Bvh7VR0CXAU8LiJDfWy3UlVHex7/HoTjXr7qCmfeicHTIKWbKxFcFx3jTGx0cDUcXOt2GmNMGAVSKAap6iOq+oHnMRtnqPEWUdVSVS30vD4L7ACyWrrfkNj+J7hwsv3dxG4o5z5I6GRnFca0M4EUig0iclX9goiMAz4JZggR6QPkAGt8rB4vIptE5G0RGeZnH7NFpEBECsrKgtw5rGA+pPaHvtcEd7+tTXwy5D3iDD9+Yq/baYwxYeKveewWEdkMjAM+FZH9IrIPWAUE7RtTRJKB14G/80yQ5K0Q6K2qo3Aufb3R2H5Uda6q5qlqXnp6erDiwZFtzuWWKx9ybuq2d+O+6dyrWfVrt5MYY8LEX/PYkDdtEZFYnCKxSFUXN1zvXThUdZmI/FpE0lT1WKizfaHgeYiOd0ZTNZDSHUbe6Qyzft0/Q2JXtxMZY0Ks0TMKVS3yfgAXcPpP1D9aREQEeA7Yoar/3cg23T3bISJjPXmPt/TYAassh00vw7CZ9oXobcIcqLkA655zO4kxJgwCaR47XUR2A/uAj4D9wNtBOPZE4H7geq/mr7eIyGMi8phnm68CW0VkE/BL4C7VME65tvV1qDprN7EbyhgC/SfD2t85LcKMMW2av0tP9X6M03z1L6qaIyLXAXe39MCq+jFNDFeuqk/jNMkNP1UoeA4yhkHPsa5EiGgTnoCF02Hlz2HS30Nsi7vWGGMiVCCtnqpV9TgQJSJRqvoBMDq0sSLAoUIo3QR5dhPbp77XODP8rfgp/GIo/OXfnCHYjTFtTiBnFKc8LZNWAItE5ChOZ7m2rWA+xCY5N27NpUTgvsWw70NYOw8++V/nMXAqjP2GU0SswBrTJgRSKGYAFcB3gHuBToA7PaTD5cIp2PI6jPw6JHR0O03kioqCK653HqcOOMW1cCHsegtSBzhTqY6+2+mkZ4xptSSc94bDJS8vTwsKCpq/gzW/g7e/B7M/dEZNNYGrroDtb8DaZ6GkwHNW9nUY+yh0a7S/pDHGZSKyXlXzfK1r9IxCRD5W1atF5Cxfbg4rgKpq2/yprer8Ms7MtSLRHLEJzvSpo+6CQxucy1KbXnIGVew1wSkYQ25r+9PIGtOG+OtHcbXnOUVVO3o9UtpskQA4sArKdlqT2GDIzIHbn4Hv7oDJ/w5nSuC1h+AXw+GD/4IzpW4nNMYEwG+rJxGJEpGt4QoTEQrmQ3wnGJ7vdpK2I7ErTPw2fGsD3PMqdB8BH/0E/mc4/PFB2P+JcyZnjIlIfm9mq2qdZ0C+Xqra9ts+njvmjBR75UMQl+R2mrYnKhoG3uw8jn/uFOUNL8K2JU5/lTGPOK3M4pPdTmqM8RJIP4oewDYReV9EltY/Qh3MFRsXQW2V03fChFbqFXDz/4Hv7oTpv3KKyFvfhf8eAm9/H47tdjuhMcYjkOaxT4U8RQQ4V1HF6fd/g3TKoUfGELfjtB9xiZA7y5k5sHid01pq3XOw5rfQ71oY8ygMnOJMnNTK1dYpFdW11NQqNXV11NSp86h1XtfWKdW1dZ5nZ7l+XU1dnedzX/5MTa1SW3fxdU3Dz3yxTcN911Fdp9R6PlOniqo6A7mpZ0A3z+VAZ1mdZ+/XAA2WG+6DLy1776OR/XvtA8XvoHJNtdj0/1n//6+05cPZuaJrYhx/mnN10Pfb5L8+Vf0o6EeNQElSxYdxeRTWjOAHbodpj0ScoVJ6jnXONAoXOCP3vnIvdOrpnOXlPgBJaW4n/ZILVbUcK6/kWHklx8urnOdzVZ73qjju9f6J81VhvRUTHSXE1D+iozzPQkxUlLMuun59FFFRIAginnF1RBDnyfPsvey8KYBEgRDlvC8X98Eln/nyMl7H8v7cJfv3fKYxTfXp9P9Z/x9ujd1FUxJC84OqyX4UnkmLfgUMAeKAaOBcJLd8am4/ioWr9vPDP21j2bcmMTQzYv947UdtDXz2tnOWse8jiI5zRvId8yhk54Wk53ddnXLyfFWjX/bHyqs4fu5iYThfVetzP8nxMaQlx5GaHE9qUhxpKfGkJcWRFB/T4Evb+aL2/gKPjRbPs/fypZ+5dJuLyzFR0uQXoTHemtWPwsvTwF3AH4E8YBYwIHjxIse0kZn8+M/bWbKhmKGZvqbvNmEVHeP0uRhyG5TtgnXzYONLsPkV6DHa6ZMx/A6I7eB3NxXVtZd86Zd5no+f+3IROHGukjofv52io4SuSXHOl35yPL17JTpFINlZTkuOIzUpnrQUpzAkxEaH5r+JMS4I5IyiQFXzRGSzqo70vPepqk4IS8JmaEnP7NkLC9hw8BSrnryemOhA7vWbsKo86xSKtfOgbAd06OLM5Z33CHTty/ZDZ1i0pogdpWecs4KzlZxr5Fd/Ulw0qfVf8p7nNM8ZQH0RSE+OJzU5ns4dYomKsl/opu1q6RnFeRGJAzaKyE+BUqDNth3Nz81i+fYjfPL5cb4yMIhTqprgiE9xxpDKewSKPoG1z6Krfg2fPk1h3Bh+VX4ta6JHk9MrlVHZnZ0v/uS4i0Wg/lJQcjwd4uxXvzGB8DeER56qFuBMLhQFzMEZGLAncEd44oXfdYMz6NQhlsWFxVYoIpkIh7vk8YfO6SyPupUplW9zv3zAC3Frqe3cl+jMGyApA5LTIcn70QHiO9jItsZcBn9nFM96hhd/CXhZVbfTDprKxsdEM21kD14vLKa8sobk+NbfLLMtUVVW7T3Oi6uKWL79CHWqXD+oL6PH/1+69OsEu/5MdMHzsOU1qDjleycxCZ6ikeZ5zrj4Ojnjy+8npraJprnGtESj/wI8s9kNwrmR/ZqIVHGxaBQF4+AiMgX4X5yWVPNU9ScN1otn/S3AeeBBVS0MxrH9yc/NZtGaA7y9pZSv5fUM9eFMAMora1hSWMzCVUXsPlpOl8RYvjGpL/eN603ProkXNxx+h/MAqKmC88fgXJnzKC+7+PqL947AkW1QfhTqqn0fvENX30XEV3GJS7azFdPmNDWExy6cs4inRGQUTtH4q4gcVtWJLTmwiEQDzwCTgWJgnYgs9Zy51JuK08JqADAO+I3nOaRye3WmT2oiSzaUWKFw2e4jZ3lxdRGLC0sor6xhZHYnfva1UUwb2aPplkUxcdAx03k0RRUqTjvDuJw72nhxObzFea443cgxO1w8W7mksKQ7w5NExTo90aNjISrGx7Ln4W/ZipEJo4DOqUUkCsgAuuHcyC4LwrHHAntUda/nGC/jTJLkXShmAAvVaZq1WkQ6i0gPVQ3psKMiwsycbP7n/c8oOXWBrM7+m1+a4KqpreO97UdYuKqIVXuPExcTxbSRPZg1vg+je3YOzUFFoENn55HWP4CQlZ6i0vAM5ejF98+UONPpniuDuiBPCinRXsUj2lNsfC3HeBUj7+UGBUqicHq5SSPPNLE+kGea//lG/zs0VTD9fbaJjza1QSQW69gkuOqxoO/Wb6EQkUnA3cDtwFbgZeA7qtrIz6nLkgUc9Fou5tKzBV/bZOG0vGqYdTYwG6BXr14tDjczJ4tf/OUz3thQwuPXBfDFYVqs7GwlL689wB/WHqD0dAVZnTvw/SmDuXNMT7omxbkd78ti4qFTlvNoiipcOOkUjKpzTtGoq4Ha6ouvv/RerXMZzOdy/bYNl+vfq/Xab8Nlz/6qzl963LpanDEzlItjZ3gvX+5zCz9vmicpI7yFQkQOAgdwisNTqnokyMf2VY4b/g0JZBvnTdW5wFxw+lG0LBr0Sk1kbJ+uLNlQwt9ee4X1cg0RVWV90UkWriri7a2lVNcqkwak8e8zhnP94Ayi20LfBRFnqPXErm4naV3qB4ZqfIOmPx+Szwbw+TbG3xnF1cG6ad2IYpymtvWygUPN2CZkZuZm8U+Lt7Cl5DQjszuH67DtwvmqGv608RALVzmd41ISYrj/qj7cd1Uv+qXbMOMGvhj8ybjOX6unUBYJgHXAABHpC5Tg3Ci/p8E2S4E5nvsX44DTob4/4e2WET340dJtLC4ssUIRJPuOneP3q4v4Y8FBzlTUMLh7Cv85cwS352SSGGfNUI2JRK79y1TVGhGZA7yL0zx2vqpuE5HHPOt/CyzDaRq7B6d5bFgniujUIZbJQ7qxdNMh/uXWIcTakB7NUlunfLDzKAtXF7HiszJiooSpI3owa3xv8np3sct6xkQ4V3/CqeoynGLg/d5vvV4r8Hi4c3nLz83irS2lfLSrjBuHdnMzSqtz4lwVrxYc5Periyg+eYFuHeP5zo0DuXtsTzI6JrgdzxgTIH83s3+Fnzs2qvqtkCSKMNcMTCc1KY4lG0qsUARo08FTLFxVxJubD1FVU8dV/bryz7cMYfLQbnZWZkwr5O+MonnDr7YxsdFR3DYqkz+sPcDp89V0Sox1O1JEqqiu5c+bS3lx1X42FZ8mKS6aO/N6cv/43gzsluJ2PGNMC/i7mb0gnEEi2R252bzw6X7e2lLKPeNa3kejLTl44jyL1hzglXUHOHm+mivSk3hq+jDyc7NISbCiakxb0OQ9ChFJB74PDAW+uLCsqteHMFdEGZ7Vkf4ZySzZUGyFwsuP/rSVhauLEGDy0G48ML4P469ItZvTxrQxgdzMXgS8AtwKPAY8QHCG8Gg1RIT83Cx++s4uio6fo3dqm52OI2BbS06zYFUR+TlZ/MPNg8i0YU6MabMCubOYqqrPAdWq+pGqPgxcFeJcEef20VmIwJINJW5HiQjzVu4lKS6aH00fZkXCmDYukEJRP/ZyqYjcKiI5OD2k25XMzh0Y3y+VJRtKaGr62Lbu0KkLvLm5lLvG9qJTB7sPYUxbF0ih+A8R6QT8PfAPwDycme7anfzcbIqOn6fwwEm3o7jq+U/2AfDQxD7uBjHGhEWThUJV/6yqp1V1q6pep6pXqurScISLNFOGd6dDbDSLC9vv5aczFdW8tPYgt47oQXaXxKY/YIxp9ZosFCLSV0T+W0QWi8jS+kc4wkWa5PgYbh7WjTc3HaKyptbtOK54ee0ByitreHRSP7ejGGPCJJBWT28AzwFvAnUhTdMK5Odm88bGQ/x1x1GmjujhdpywqqqpY/7H+xnfL5UR2Z3cjmOMCZNACkWFqv4y5ElaiYn908hIiWfxhpJ2Vyje2nKIw2cq+K/8EW5HMcaEUSA3s/9XRH4kIuNFJLf+EfJkESo6Srg9J4sPdh7lxLkqt+OEjaoyd8U+BmQk85WB6W7HMcaEUSCFYgTwKPAT4Oeex89CGSrS5edmUVOnvLkpbHMoue6TPcfZUXqGRyf1I6otzDpnjAlYIJeeZgL9VLX9/HxuwuDuHRnSoyOLN5TwwIQ+bscJi7kr95KWHM+MnEy3oxhjwiyQM4pNQOcQ52h17sjNYtPBU3xeVu52lJDbUXqGFZ+V8dDEPsTHRLsdxxgTZoEUim7AThF5t703j/U2fVQmUQJL2kGfinkr99EhNpp7bUBEY9qlQC49/SjYBxWR/wfcBlQBnwMPqeopH9vtB84CtUCNquYFO0tzZXRMYNKAdJZsKOG7kwe22ev2h09XsHRTCfeO603nxDi34xhjXBBIz+yPfD1aeNz3gOGqOhL4DPgnP9tep6qjI6lI1MvPzaLk1AXW7DvhdpSQeeHT/dTWKQ9P7Ot2FGOMSxotFCLysef5rIic8XqcFZEzLTmoqi5X1RrP4mpa6SCDNw3tTnJ8DEs2FLsdJSTKK2tYtKaIqcN70CvVhuswpr1qtFCo6tWe5xRV7ej1SFHVjkHM8DDwdmMxgOUisl5EZvvbiYjMFpECESkoKwvPdBkd4qKZOrw7y7Yc5kJV2xvS45V1BzlbUcM3JtnZhDHtmd9LTyISJSJbm7NjEfmLiGz18Zjhtc2/ADU4kyP5MlFVc4GpwOMick1jx1PVuaqap6p56enh6xCWn5tNeWUNy7cfDtsxw6Gmto75H+9jbJ+u5PTq4nYcY4yL/N7MVtU6EdkkIr1U9cDl7FhVb/S3XkQeAKYBN2gjEzyo6iHP81ERWQKMBVZcTo5QG9e3K1mdO7BkQwkzRme5HSdolm09TMmpC/zb9GFuRzHGuCyQ5rE9gG0i8n6wmseKyBScebinq+r5RrZJEpGU+tfATUCzzm5CKSpKuD0nkxWflXH0bIXbcYLCGa7jc/qlJXHD4Ay34xhjXBZI89inQnDcp4F44D0RAVitqo+JSCYwT1Vvwem/scSzPgb4g6q+E4IsLTYzJ5tnPvicpRsP8Y02MPz26r0n2Fpyhv+cOaLNNvs1xgSu0UIhIgnAY0B/YAvwnFdLpRZR1f6NvH8IuMXzei8wKhjHC7X+GcmMyu7E4sKSNlEonl25l9SkOPJz286lNGNM8/m79LQAyMMpElNxBgM0jcjPzWZ76Rl2Hm5Ry2HX7T5ylr/uPMqs8X1IiLXhOowx/gvFUFW9T1V/B3wVmBSmTK3SbaMyiYmSVj+kx7yV+4iPieL+8b3djmKMiRD+CkV1/YtgXXJqy7omxXHtoAyWbCihts5nI66Id/RsBUs2lPC1vGy6JtlwHcYYh79CMcq7NzYwMlg9s9uq/Nwsjp6t5JM9x9yO0iwLPy2iuq6OR65u/fdZjDHB469ndnSD3tgxIeqZ3WZcPziDjgkxLNnQ+i4/na+q4cXVRdw0tBt905LcjmOMiSCB9KMwAUqIjWbaqEze2XqY8srWdbXujwXFnL5Qzexr7GzCGPNlViiCLD8niwvVtbyztfUM6VFbp8z7eC+5vTpzZe+ubscxxkQYKxRBdmXvLvTqmtiqRpR9d9thDp64YGcTxhifrFAEmYgwMyeLTz8/TunpC27HaZKq8rsVe+mdmsjkod3djmOMiUBWKEIgPzcLVXhjwyG3ozSpoOgkmw6e4htX9yXahuswxvhghSIEeqcmkde7C4sLi2lkYNyIMXfFXrokxvLVK3u6HcUYE6GsUITIzNwsdh8tZ2tJ5HY5+bysnL/sOML9V/WmQ5wN12GM8c0KRYhMG5FJXHQUiyP4pvZzH+8jNjqK+8f3cTuKMSaCWaEIkU6Jsdw4NIOlGw9RXVvndpxLHCuv5PX1xdyRm0V6SrzbcYwxEcwKRQjNzMnm+LkqVnwWnjm8L8eLq4qorLHhOowxTbNCEUJfGZhO16Q4FkfYkB4Xqmp5cXURNw7JoH9GsttxjDERzgpFCMXFRHHbyB68t/0Ipy9UN/2BMHm9sJgT56p4tA1MsmSMCT1XCoWI/JuIlIjIRs/jlka2myIiu0Rkj4g8Ge6cwZCfm01VTR1vbyl1OwrgDNfx3Mf7GJXdibF9bbgOY0zT3Dyj+IWqjvY8ljVcKSLRwDM4s+sNBe4WkaHhDtlSI7M7cUV6EosjZEKjv+w4wr5j53j0mn545iM3xhi/IvnS01hgj6ruVdUq4GVghsuZLpuIkJ+bzdr9Jzh44rzbcXh2xV6yu3RgyjAbrsMYExg3C8UcEdksIvNFpIuP9VnAQa/lYs97PonIbBEpEJGCsrLIamV0e44T2+15KtYXnaSg6CSPXN2XmOhI/o1gjIkkIfu2EJG/iMhWH48ZwG+AK4DRQCnwc1+78PFeo+NhqOpcVc1T1bz09PRg/BGCJqtzB67q19X1IT3mrdxLx4QYvp5nw3UYYwIXE6odq+qNgWwnIs8Cf/axqhjw/kbLBiJ/lL1G5Odm873XNlN44BRX9vZ1AhVaRcfP8c62w/zNV64gKT5k/9uNMW2QW62eengtzgS2+thsHTBARPqKSBxwF7A0HPlCYerw7iTERrk2T8VzH+8jJkp4cEIfV45vjGm93LpQ/VMR2SIim4HrgO8AiEimiCwDUNUaYA7wLrADeFVVt7mUt8VSEmK5aWh33txUSmVNbViPffJcFa8WHOT20VlkdEwI67GNMa2fK9cgVPX+Rt4/BNzitbwMuKTpbGuVn5vF0k2H+GBnGVOGh6/V0e9XF1FRXcejNoOdMaYZrOlLGF3dP430lHgWF4bv8lNFdS0LVu3n2kHpDOyWErbjGmPaDisUYRQTHcWMUZl8sOsoJ89VheWYb2wo4Vh5FbNtuA5jTDNZoQiz/NxsqmuVP28OfQOuujrl2ZV7GZbZkfFXpIb8eMaYtskKRZgNzezI4O4pvB6GIT0+2HWUz8vOMduG6zDGtIAVChfk52ax8eApPi8rD+lx5q7YS2anBG4Z0aPpjY0xphFWKFwwY3QWUeLcPwiVTQdPsWbfCR6+ui+xNlyHMaYF7BvEBd06JjCxfxqLC0uoqwvNkB7PrtxLSnwMd46x4TqMMS1jhcIld+RmU3LqAuv2nwj6vg+eOM+yLaXcM64XKQmxQd+/MaZ9sULhkpuGdSMxLjok81TM/2QfUSI8OLFP0PdtjGl/rFC4JDEuhqnDe7BsSykV1cEb0uP0+WpeWXeQ6aMy6dGpQ9D2a4xpv6xQuOiO3CzOVtbw3vYjQdvnorVFnK+q5RvWwc4YEyRWKFx0Vb9UenRKCNqQHpU1tbzwyX4mDUhjaGbHoOzTGGOsULgoKkq4PSeLFbuPUXa2ssX7W7rxEEfPVvKonU0YY4LICoXL8nOyqK1Tlm5q2ZAeqs5wHYO7pzBpQFqQ0hljjBUK1w3olsKIrE4tvvz00WdlfHaknEcn2XAdxpjgskIRAfJzs9h26Ay7Dp9t9j6eXbmXbh3juW1UZhCTGWOMFYqIcNuoTKKjhMXNnCZ1a8lpPtlznIcm9iUuxv6XGmOCy605s18RkY2ex34R2djIdvs9U6ZuFJGCMMcMm7TkeK4dmM4bG0qobcaQHvNW7iUpLpq7x/YKQTpjTHvn1lSod9a/FpGfA6f9bH6dqh4LfSp35edm8/7Oo6z6/DhXX8bN6EOnLvDm5lIenNCHTh1suA5jTPC5ep1CnLuuXwdecjNHJLhhSAYpCTGXfVP7+U/2AfCQDddhjAkRty9oTwKOqOruRtYrsFxE1ovIbH87EpHZIlIgIgVlZWVBDxpqCbHRTBvZg7e3HuZcZU1AnzlTUc1Law9y64geZHdJDHFCY0x7FbJCISJ/EZGtPh4zvDa7G/9nExNVNReYCjwuItc0tqGqzlXVPFXNS09PD9KfIrxm5mRzobqWd7cdDmj7l9ceoLyyxjrYGWNCKmT3KFT1Rn/rRSQGyAeu9LOPQ57noyKyBBgLrAhmzkiS17sLPbt2YHFhCfm52X63raqpY/7H+xnfL5UR2Z3ClNAY0x65eenpRmCnqvq8KC8iSSKSUv8auAnYGsZ8YRcVJczMyeaTz49x+HSF323f2nKIw2cqmH2NnU0YY0LLzUJxFw0uO4lIpogs8yx2Az4WkU3AWuAtVX0nzBnDbmZOFqrwxsbG56lQVeau2MeAjGS+MrB1XmYzxrQerjSPBVDVB328dwi4xfN6LzAqzLFc1zctidxenVlcWMw3r/E9HMcne46zo/QMP71jJFFRNlyHMSa03G71ZHzIz83msyPlbDt0xuf6uSv3kpYcz4wcG67DGBN6Vigi0LSRPYiLjvI5TeqO0jOs+KyMhyb2IT4m2oV0xpj2xgpFBOqcGMf1gzNYuqmEmtq6L62bt3IfHWKjuXecDddhjAkPKxQRamZuFsfKq1i5++LoJYdPV7B0Uwl3julJ58Q4F9MZY9oTKxQR6rpBGXROjOV1ryE9Xvh0P7V1ysMT+7qYzBjT3lihiFBxMVFMH5XJe9uPcKaimvLKGhatKWLq8B70SrXhOowx4WOFIoLNzMmisqaOt7eU8sq6g5ytqOEbk+xswhgTXq71ozBNG92zM/3Skni1oJjDpysY26crOb26uB3LGNPO2BlFBBMRZuZksb7oJCWnLvCoDddhjHGBFYoId3tOFgD90pK4YXCGy2mMMe2RXXqKcD27JvIvtwxhWFZHG67DGOMKKxStgF1yMsa4yS49GWOM8csKhTHGGL+sUBhjjPHLCoUxxhi/rFAYY4zxywqFMcYYv6xQGGOM8csKhTHGGL9EVd3OEHQiUgYUNfPjacCxJrcKP8t1eSzX5bFcl6ct5uqtqum+VrTJQtESIlKgqnlu52jIcl0ey3V5LNflaW+57NKTMcYYv6xQGGOM8csKxaXmuh2gEZbr8liuy2O5Lk+7ymX3KIwxxvhlZxTGGGP8skJhjDHGLysUHiIyRUR2icgeEXnS7Tz1RGS+iBwVka1uZ6knIj1F5AMR2SEi20Tk225nAhCRBBFZKyKbPLmecjuTNxGJFpENIvJnt7N4E5H9IrJFRDaKSIHbeeqJSGcReU1Ednr+ro2PgEyDPP+d6h9nROTv3M4FICLf8fy93yoiL4lIQtD2bfconH/AwGfAZKAYWAfcrarbXQ0GiMg1QDmwUFWHu50HQER6AD1UtVBEUoD1wO1u//cSEQGSVLVcRGKBj4Fvq+pqN3PVE5HvAnlAR1Wd5naeeiKyH8hT1YjqQCYiC4CVqjpPROKARFU95XKsL3i+N0qAcara3A6+wcqShfP3faiqXhCRV4FlqvpCMPZvZxSOscAeVd2rqlXAy8AMlzMBoKorgBNu5/CmqqWqWuh5fRbYAWS5mwrUUe5ZjPU8IuKXkIhkA7cC89zO0hqISEfgGuA5AFWtiqQi4XED8LnbRcJLDNBBRGKAROBQsHZshcKRBRz0Wi4mAr74WgMR6QPkAGtcjgJ8cXlnI3AUeE9VIyIX8D/A94A6l3P4osByEVkvIrPdDuPRDygDnvdcrpsnIkluh2rgLuAlt0MAqGoJ8DPgAFAKnFbV5cHavxUKh/h4LyJ+iUYyEUkGXgf+TlXPuJ0HQFVrVXU0kA2MFRHXL9eJyDTgqKqudztLIyaqai4wFXjcc7nTbTFALvAbVc0BzgGRdO8wDpgO/NHtLAAi0gXnKkhfIBNIEpH7grV/KxSOYqCn13I2QTxta4s89wBeBxap6mK38zTkuUzxITDF3SQATASme+4FvAxcLyK/dzfSRap6yPN8FFiCcynWbcVAsdcZ4Ws4hSNSTAUKVfWI20E8bgT2qWqZqlYDi4EJwdq5FQrHOmCAiPT1/FK4C1jqcqaI5blp/BywQ1X/2+089UQkXUQ6e153wPnHs9PVUICq/pOqZqtqH5y/W39V1aD92msJEUnyNEjAc2nnJsD1Fnaqehg4KCKDPG/dALjeuMTL3UTIZSePA8BVIpLo+fd5A869w6CICdaOWjNVrRGROcC7QDQwX1W3uRwLABF5CbgWSBORYuBHqvqcu6mYCNwPbPHcDwD4Z1Vd5l4kAHoACzytUaKAV1U1opqiRqBuwBLnu4UY4A+q+o67kb7wBLDI8+NtL/CQy3kAEJFEnBaS33Q7Sz1VXSMirwGFQA2wgSAO52HNY40xxvhll56MMcb4ZYXCGGOMX1YojDHG+GWFwhhjjF9WKIwxxvhlhcIYP0Sk3PPcR0TuCfK+/7nB8qfB3L8xwWKFwpjA9AEuq1B4+nP486VCoapB60lrTDBZoTAmMD8BJnnmIPiOZ/DB/yci60Rks4h8E0BErvXM1fEHYIvnvTc8A+5tqx90T0R+gjPS50YRWeR5r/7sRTz73uqZJ+JOr31/6DVHwyJPL1xjQsp6ZhsTmCeBf6ifR8LzhX9aVceISDzwiYjUj9Y5Fhiuqvs8yw+r6gnPsCLrROR1VX1SROZ4BjBsKB8YDYwC0jyfWeFZlwMMwxmL7BOcXvIfB/sPa4w3O6MwpnluAmZ5hjBZA6QCAzzr1noVCYBvicgmYDXO4JMD8O9q4CXPSLhHgI+AMV77LlbVOmAjziUxY0LKziiMaR4BnlDVd7/0psi1OENiey/fCIxX1fMi8iHQ1BSV/i4nVXq9rsX+DZswsDMKYwJzFkjxWn4X+BvPcOuIyMBGJtbpBJz0FInBwFVe66rrP9/ACuBOz32QdJyZ3tYG5U9hTDPYrxFjArMZqPFcQnoB+F+cyz6FnhvKZcDtPj73DvCYiGwGduFcfqo3F9gsIoWqeq/X+0uA8cAmnAm0vqeqhz2Fxpiws9FjjTHG+GWXnowxxvhlhcIYY4xfViiMMcb4ZYXCGGOMX1YojDHG+GWFwhhjjF9WKIwxxvj1/wHGkZAa0gZrKgAAAABJRU5ErkJggg==\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": 5
}