{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Sequential Subspace Optimization Method for Large-Scale Unconstrained Problems\n",
"## (ECE 602 Final Project)\n",
"\n",
"## Table of Contents\n",
"* [Introduction](#first-bullet)\n",
"* [SESOP (Sequential Subspace Optimization)](#second-bullet)\n",
"* [Fast SESOP](#third-bullet)\n",
"* [Newton method](#forth-bullet)\n",
"* [Modified Newton method](#fifth-bullet)\n",
"* [Implementation details](#sixth-bullet)\n",
"* [References](#seventh-bullet)\n",
"\n",
"## Introduction \n",
"\n",
"Cosider the function $f$ is the objective we want to minimize over all $x$, belonging to $R^n$. When the number of variables is very large, let say $n$ is greater than $10^4$ and more, we need \n",
" for a optimization algorithm, for which storage requirement and computational cost per iteration grow not more than linearly in $n$. Conjugate gradient (CG) method is one the first methods which was proposed to tackle computational costs. It is known that CG worst case\n",
"convergence rate for quadratic problems. However, CG to nonlinear\n",
"functions is not worstcase optimal. In the paper, a method is presented ,which is equivalent to CG in the quadratic case, and often outperforms CG and Trancated Newton (TN) in non-quadratic case, while preserving worst-case optimality.\n",
"\n",
"This method is used for large scale unconstraint optimization problems, in which at each iteration the minimum of the objective function is searched over a subspace spanned by the current gradient and by directions of few previous steps.\n",
"\n",
"\n",
"## SESOP (Sequential Subspace Optimization) \n",
"\n",
"In this section we describe the SESOP algorithm, which is a general method for\n",
"smooth unconstrained optimization. Generaly, in this algorithms , we minimize the objective function $f$ in a subsapace of $R^n$ at each iteration($k$). This subspace is constructed using the points $x_{k-1}, x_{k-2},...,x_{k-M}$ from previous iterations. Also, the gradient of $f$ (g^k) at point $x_k$ beside two Nemrowski directios are used to biuld the subspace. \n",
"\n",
"\n",
"To construct the subspace, following steps need to be done at each iteration($k$):\n",
"\n",
"\n",
"\n",
"\n",
"$min\\:\\:f(x) \\:\\:\\: where f\\:\\:is\\:\\:smooth$\n",
"\n",
"\n",
"\n",
"$1.\\:\\: g^k = \\nabla f(x^k)$\n",
"\n",
"\n",
"$2.\\:\\: d^{k-1} = x^k - x^{k-1}$\n",
"\n",
"\n",
"\n",
"$3.\\:\\: P^k = (g^k|d^{k-1}...d^{k-M})$\n",
"\n",
"\n",
"\n",
"$4.\\:\\:\\alpha ^k = argmin\\:\\:f(x^k+P \\alpha)$\n",
"\n",
"\n",
"\n",
"$5. \\:\\:x^{k+1}=x^k+P\\alpha ^k$\n",
"\n",
"\n",
"$6. \\:\\: update \\:\\:matrix\\:\\: P\\:\\: with \\:\\:new \\:\\:directios$\n",
"\n",
"\n",
"In step six, optimization is done on the subspace which is spanned by the column of matrix $P^k$. Note that the number of the columns of matrix $P^k$ is $M$ which is much smaller than original space dimenstion, $n$. To do the optimization on this subsapce, we are looking for a $\\alpha$ which minimize $f(x^k+ P\\alpha)$. This optimization can be done using Newton ,Truncated Newton(TN), CG or even exact line search method. Newton method will be explained in the sections below.\n",
"\n",
"\n",
"\n",
"## Fast SESOP \n",
"\n",
"In some applications:\n",
"\n",
"\\begin{equation}\n",
"f(x)=\\Phi (Ax)\n",
"\\end{equation}\n",
"\n",
"where, $\\Phi$ is cheap to compute, but $Ax$ is expensive to compute, but\n",
"\n",
"\\begin{equation}\n",
"f(x^k+P\\alpha)+\\Phi(Ax^k+AP\\alpha), \\:\\:\\: AP = R\n",
"\\end{equation}\n",
"\n",
"Hence, in each iteration we only need to compute $Ax$ once. Furthermore, for each $x^k$, we change only two columns of P. Considering $Ax^{k+1}=Ax+Rx^k$, we do not need to compute $Ax^{k+1}$ for each k, just update it. We can compute it directly in each let say 10 iteration. Therefore, if\n",
"\n",
"\\begin{equation}\n",
"argmin\\:\\:||Ac-y||^2 + \\mu _0 * smooth\\_function(c)\n",
"\\end{equation}\n",
"\n",
"Then,\n",
"\\begin{equation}\n",
"\\nabla _\\alpha \\Phi(Ax^k+R\\alpha) = R^T\\nabla _u\\Phi(u) \\Longrightarrow Gradient\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"\\nabla ^2 _{\\alpha \\alpha} \\Phi(Ax^k+R\\alpha) = R^T\\nabla ^2 _{uu}\\Phi(u)R \\Longrightarrow Hessian\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"x^{k+1}=x^k+\\alpha ^kd^k,\\:\\:f'\\:^k_d(x^k)\\:<0\\:\\:or\\:\\: <\\nabla f(x^k),d^k>\\:\\:<0\n",
"\\end{equation}\n",
"\n",
"\n",
"## Newton method \n",
"In this method, to find the best direction $d$ in order to optimize function $f$ around point $x^k$, Teilor exanstion is applied.\n",
"\n",
"\\begin{equation}\n",
"f(x_k,d) \\cong f(x_k) + d^Tg(x_k) + \\frac{1}{2}d^TH(x_k)d\n",
"\\end{equation}\n",
"\n",
"Define:\n",
"\n",
"\\begin{equation}\n",
"S(d) = d^Tg(x_k) + \\frac{1}{2}d^TH(x_k)d\n",
"\\end{equation}\n",
"\n",
"Then:\n",
"\n",
"\\begin{equation}\n",
"\\nabla S(d) = 0\n",
"\\end{equation}\n",
"\n",
"So we will have:\n",
"\n",
"\\begin{equation}\n",
"g^k + H^kd = 0\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"H^kd = -g^k\n",
"\\end{equation}\n",
"\n",
"## Modified Newton method \n",
"\n",
"In Newton method we have:\n",
"\n",
"\\begin{equation}\n",
"d^k = -H^{k^{-1}}g^k\n",
"\\end{equation}\n",
"\n",
"To find the minumum value of $f$ in direction $d$ and around point $x^k$, line search method can be applied. In 'Line search' method, we would like our direction to be direction of descent:\n",
"\n",
"\\begin{equation}\n",
"\\Longrightarrow f'_{d^k}=\\:\\:<0\\:\\:\\Longrightarrow \\:\\:g^{k^T}H^{k^{-1}}g^k\\:\\:>0\n",
"\\end{equation}\n",
"\n",
"Therefore $H^{k^{-1}}$ or $H^k$ should be positive semidefinite:\n",
"\n",
"\\begin{equation}\n",
"H^k \\succeq 0\n",
"\\end{equation}\n",
"\n",
"But sometimes $H^k \\nsucceq 0$, so we need to modify it in order to have no eigenvalues less than zero. As a result we will have modified Newton method:\n",
"\n",
"\\begin{equation}\n",
"H^k + E = LDL^T \\Longrightarrow Cholesky\\:\\:modified\\:\\:factorization\n",
"\\end{equation}\n",
"\n",
"Regarding Newton direction, we will have:\n",
"\n",
"\\begin{equation}\n",
"LDL^Td^k = -g^k\n",
"\\end{equation}\n",
"\n",
"In modified Newton method, Hession matrix of $f$ at point $x^k$ is modified in way that all eighen values are positive.\n",
"\n",
"After we find the best direction($d^k$), to minimize $f$ around the point $x^k$ we can do exact line search on direction $d^k$ or can apply Back tracking method.\n",
"\n",
"Exact Line search\n",
"\\begin{equation}\n",
"\\Longrightarrow\\alpha ^k = argmin\\:\\:f(x^k,\\alpha d^k)\n",
"\\end{equation}\n",
"\n",
"Back tracking in exact Line search:\n",
"\\begin{equation}\n",
"\\Phi (\\alpha) = f(x^k+\\alpha d^k) - f(x^k)\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"\\Phi '(0) = f'_{d^k}(x^k) < 0 \\equiv c\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"0<\\sigma <1,\\:\\:\\sigma = 10^{-4}:10^{-1} \n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"0<\\beta <1,\\:\\:\\beta = 0.2\n",
"\\end{equation}\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Implementation details \n",
"In this project, we implemented Fast Sesop with modified Newton method usig Back tracking search strategy. We implemented a Basis Pursuit(BP) to check the efficincy of SESOP method.\n",
"Basis Pursuit is a way for decomposing a signal into a sparse superposition of\n",
"dictionary elements, using $l_1-norm$ penalization. When the signal is sparse enough,\n",
"it is equivalent to $l_0$ minimization, where $l_0(\\alpha)$ stands for the number of non-zero\n",
"elements in $\\alpha$. BP in highly overcomplete dictionaries leads to large-scale\n",
"optimization problems. The noisy signal $y$ is described as:\n",
"\n",
"\\begin{equation}\n",
"y=x+n\n",
"\\end{equation}\n",
"\n",
"where $x$ is the original signal, and $n$ is Gaussian noise with\n",
"variance $σ_n^2$\n",
". We assume signal $x$ can be represented as linear combination of some basis(dictionary), $x=Ac$( $c$ is the cofficents and $A$ is the basis matrix).\n",
"\n",
"The MAP estimation for the coefficients is achieved by maximizing the log-likelihood:\n",
"\\begin{equation}\n",
"min_c \\:\\:\\: 0.5\\times||Ac-y||_2^2 +\\lambda \\times||c||_1\n",
"\\end{equation}\n",
"\n",
"Now, the aim is to find coefficents $c$ and reconstruct our original signal $x$ from noisy signal $y$.\n",
"\n",
"\n",
"The main parts of code are modified Newton optimization method and SESOP subspace constructing and updating. To use this code for other use( not Basis Pursuit), one needs to implement the cost function which returns the functions derivitev at input poit $x$ and function value at this point.\n",
"\n",
"## References \n",
"\n",
"1. Guy Narkiss and Michael Zibulevsky (2005). \"Sequential Subspace Optimization Method for Large-Scale Unconstrained Problems\", Tech. Report CCIT No 559, EE Dept., Technion.\n",
"2. Stephen G. Nash, “A Survey of Truncated-Newton Methods”, Journal of Computational and Applied Mathematics, 124 (2000), pp. 45-59.\n",
"3. Michael Elad, Boaz Matalon, and Michael Zibulevsky, \"Coordinate and Subspace Optimization Methods for Linear Least Squares with Non-Quadratic Regularization\", Applied and Computational Harmonic Analysis, Vol. 23, pp. 346-367, November 2007.\n",
"4. M. ELAD AND M. ZIBULEVSKY, A probabilistic study of the average perfor- mance of the basis pursuit, submitted to IEEE Trans. On Information Theory, (2004)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"#global constants\n",
"global flag_is_New\n",
"flag_is_New = True\n",
"\n",
"Nem = True #two Nemirovski directions\n",
"\n",
"newton_iters = 7\n",
"isHess = True\n",
"weights_pen = 0.01\n",
"eps_smooth_abs = 0.01\n",
"all_r = []\n",
"sesop_iters = 200\n",
"last_iters = 5\n",
"sub_grad_norm = []\n",
"glob_grad_norm = []\n",
"\n",
"#tolerance values\n",
"diff_x_thresh = 1.0e-8\n",
"func_thresh = 1.0e-16\n",
"glob_thresh = 1.0e-8"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def mult_norm(x1, x2):\n",
" return np.matmul(x1, x2)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def mult_adj(x1, x2):\n",
" return np.matmul(x1.T, x2).T"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def l2norm(r, y, Z=0, isNew=True, nouts=1): # same as ls_fgh\n",
" if isNew:\n",
" l2norm.tmp = r - y\n",
" o = 0.5 * np.matmul((l2norm.tmp.flatten()).T, l2norm.tmp.flatten())\n",
" l2norm.f_prev = o\n",
" else:\n",
" o = l2norm.f_prev\n",
" \n",
" if nouts > 1:\n",
" \n",
" grad = l2norm.tmp\n",
" \n",
" if nouts > 2:\n",
" grad = l2norm.tmp\n",
" Hes = Z\n",
" \n",
" return o, grad, Hes\n",
" \n",
" return o, grad\n",
"\n",
" return o"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def fast_smooth_approx(x, Hes=np.array([]), isNew=True, nouts=1): #same as sum_fast_abs_smooth_new\n",
" if isNew:\n",
" fast_smooth_approx.out_f = True\n",
" fast_smooth_approx.g_f = True\n",
" fast_smooth_approx.h_f = True\n",
"\n",
" m = weights_pen\n",
" epsilon = eps_smooth_abs\n",
" \n",
" if fast_smooth_approx.out_f:\n",
" fast_smooth_approx.t = (1.0 / epsilon) * x\n",
" fast_smooth_approx.abs_t = np.absolute(fast_smooth_approx.t)\n",
" fast_smooth_approx.backup = 1.0 / (fast_smooth_approx.abs_t + 1)\n",
" f = epsilon * (fast_smooth_approx.abs_t + fast_smooth_approx.backup - 1)\n",
" out = m*(f.sum(axis=0))\n",
" fast_smooth_approx.out_prev = out\n",
" fast_smooth_approx.out_f = False\n",
" else:\n",
" out = fast_smooth_approx.out_prev\n",
" \n",
" if nouts > 1: #Gradient\n",
" if fast_smooth_approx.g_f:\n",
" fast_smooth_approx.tmp = fast_smooth_approx.backup * fast_smooth_approx.backup\n",
" g = m * (fast_smooth_approx.t * (fast_smooth_approx.abs_t + 2) * fast_smooth_approx.tmp)\n",
" fast_smooth_approx.g_old = g\n",
" fast_smooth_approx.g_f = False\n",
" else:\n",
" g = fast_smooth_approx.g_old\n",
" \n",
" if nouts > 2: # Hessian-matrix product\n",
" if fast_smooth_approx.h_f:\n",
" fast_smooth_approx.to_f = m * (2.0/epsilon) * fast_smooth_approx.tmp * fast_smooth_approx.backup\n",
" fast_smooth_approx.h_f = False\n",
" \n",
" if Hes.size > 0:\n",
" HZ = np.expand_dims(fast_smooth_approx.to_f.flatten(), -1) * Hes\n",
" else:\n",
" HZ = []\n",
" \n",
" if nouts > 3: #diagonal of Hessian\n",
" d_comp = fast_smooth_approx.to_f\n",
" \n",
" return out, g, HZ, d_comp\n",
"\n",
" return out, g, HZ\n",
" \n",
" return out, g\n",
"\n",
" return out"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def ux_body(x, Bx, y, R, P, diff_alpha, alpha, a, costf1, costf2):\n",
" global flag_is_New\n",
" new_aplha = alpha + a * diff_alpha\n",
" r = Bx + np.matmul(R, new_aplha)\n",
" out = x + np.matmul(P, new_aplha)\n",
" \n",
" flag_is_New = True\n",
" \n",
" r_func, g2 = costf1(r, y, isNew = flag_is_New, nouts = 2)\n",
" r2_func, g2_func = costf2(out, isNew = flag_is_New, nouts = 2)\n",
" \n",
" flag_is_New = False\n",
"\n",
" g_alpha = np.matmul(g2.T, R).T + np.matmul(g2_func.T, P).T\n",
"\n",
" return r_func + r2_func, g_alpha, new_aplha, r, out\n",
"\n",
"#Newton iterations\n",
"\n",
"def newton_body(r, R, P, x, x0, Bx, y, alpha, i, costf1, costf2):\n",
" global flag_is_New\n",
"\n",
" first_func, first_x, new_R_mult = costf1(r, y, R, isNew = flag_is_New, nouts = 3)\n",
"\n",
" if isHess:\n",
" second_func, second_x, new_P_mult = costf2(x, P, isNew = flag_is_New, nouts = 3)\n",
" \n",
" flag_is_New = False\n",
" out = first_func + second_func\n",
"\n",
" a_g = np.matmul(first_x.T, R).T + np.matmul(second_x.T, P).T # compute a_g= R'*first_x + P'*second_x in efficient way\n",
" \n",
" if i == 1:\n",
" sub_grad_norm.append(np.linalg.norm(a_g, 2))\n",
" \n",
" n_a = np.matmul(R.T, new_R_mult) + np.matmul(P.T, new_P_mult)\n",
" \n",
" Lambda, eig_vecs = np.linalg.eig(n_a)\n",
"\n",
" l = np.abs(Lambda)\n",
" l = np.maximum(l, 1e-12 * np.amax(l))\n",
" \n",
" a_d = np.matmul(-eig_vecs, np.matmul(np.diag(1./l), np.matmul(eig_vecs.T, a_g)))\n",
"\n",
" f0 = out\n",
" step = 1\n",
" a_d_abs = np.absolute(np.matmul(a_g.T, a_d))\n",
"\n",
" #r, x calculations loop\n",
" for j in range(0, 100):\n",
"\n",
" out, a_g, upd_alpha, r, x = ux_body(x0, Bx, y, R, P, a_d, alpha, step, costf1, costf2)\n",
" if (out < f0) or ((out < (f0 + 1e-9 * np.maximum(np.abs(f0), 1)) ) and (np.matmul(a_g.T, a_d) <= 0.1 * a_d_abs)):\n",
" break\n",
" \n",
" step = 0.25 * step\n",
"\n",
" alpha = upd_alpha\n",
" all_r.append(step)\n",
" sub_grad_norm.append(np.linalg.norm(a_g, 2) + 1e-20)\n",
" \n",
" return r, x, a_g"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def NewtonOOS(x0,y,Bx,P,R,costf1,costf2):\n",
" alpha = np.zeros([P.shape[1], 1])\n",
" \n",
" x = x0\n",
" r = Bx\n",
" \n",
" for i in range(0, newton_iters):\n",
" u_n, x_n, a_g = newton_body(r, R, P, x, x0, Bx, y, alpha, i, costf1, costf2)\n",
" \n",
" if np.linalg.norm(a_g, 2) < 1e-12:\n",
" break\n",
"\n",
" return u_n, x_n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def sesop(basis, y, coeffs, costf1, costf2, mult1, mult2):\n",
" p = 0 #Previous step\n",
" x0 = coeffs\n",
"\n",
" AC = mult1(basis, coeffs)\n",
" \n",
" n = np.max(coeffs.shape)\n",
" m = np.max(AC.shape)\n",
" P1 = np.zeros([n, 1])\n",
" indP = 0 # Columns of P span the subspace of optimization\n",
" R1 = np.zeros([m, 1])\n",
" indR = 0 # R=AP\n",
"\n",
" x = coeffs\n",
"\n",
" for i in range(0, sesop_iters):\n",
" flag_is_New = True\n",
" \n",
" first_norm_f, first_norm_g = costf1(AC, y, isNew=flag_is_New, nouts=2)\n",
" second_norm_f, second_norm_g = costf2(x, isNew=flag_is_New, nouts=2)\n",
"\n",
" flag_is_New = False\n",
"\n",
" f = first_norm_f + second_norm_f\n",
" first_grad_x = mult2(first_norm_g, basis)\n",
"\n",
" x_grad = first_grad_x + second_norm_g\n",
"\n",
" #Check stopping criteria\n",
" if i > 0 and (np.linalg.norm(x - x_prev, 2) < diff_x_thresh or np.abs(f - f_old) < func_thresh or np.linalg.norm(x_grad, 2) < glob_thresh):\n",
" break\n",
"\n",
" f_old = f\n",
" d_new = -x_grad #New direction: gradient\n",
"\n",
" if i > 0:\n",
" p = (1.0/dx_norm) * dx\n",
"\n",
" #r = (1.0/dx_norm) * (AC-AC_prev)\n",
"\n",
" r = mult1(basis, p)\n",
"\n",
" if last_iters > 0:\n",
" indP = indP + 1\n",
" if indP > last_iters:\n",
" indP = 1\n",
"\n",
" P1[:, indP-1] = np.squeeze(p)\n",
" R1[:, indP-1] = np.squeeze(r)\n",
"\n",
" d_new_norm = (1.0 / np.linalg.norm(d_new, 2)) * d_new\n",
" Ad_new_norm = mult1(basis, d_new_norm)\n",
"\n",
" i2 = min(i, last_iters + 1) - 1\n",
" i2 = i2 + 1\n",
"\n",
" if i2 == last_iters + 1:\n",
" P1[:, i2-1] = np.squeeze(d_new_norm)\n",
" R1[:, i2-1] = np.squeeze(Ad_new_norm)\n",
" else:\n",
" P1_temp = np.zeros([n, i2+1])\n",
" P1_temp[:, 0:i2] = P1\n",
" P1_temp[:, i2] = np.squeeze(d_new_norm) #Put normalized [preconditioned] gradient at the last free column\n",
" P1 = P1_temp\n",
"\n",
" R1_temp = np.zeros([m, i2+1])\n",
" R1_temp[:, 0:i2] = R1\n",
" R1_temp[:, i2] = np.squeeze(Ad_new_norm) #Put normalized [preconditioned] gradient at the last free column\n",
" R1 = R1_temp\n",
"\n",
" if Nem:\n",
" P2 = np.zeros([n, 2])\n",
" R2 = np.zeros([m, 2])\n",
"\n",
" if i == 0:\n",
" NemWeight = 1\n",
" d2Nem = -x_grad\n",
" else:\n",
" NemWeight=0.5 + np.sqrt(0.25 + NemWeight ** 2)\n",
" d2Nem = d2Nem - NemWeight * x_grad\n",
"\n",
" d1Nem = x - x0\n",
" Ad1Nem = mult1(basis, d1Nem)\n",
"\n",
" tmp = 1.0/(np.linalg.norm(d1Nem, 2) + 1e-20)\n",
"\n",
" P2[:, 0] = np.squeeze(tmp * d1Nem)\n",
" R2[:, 0] = np.squeeze(tmp * Ad1Nem)\n",
"\n",
" Ad2Nem = mult1(basis, d2Nem)\n",
" tmp = 1.0/(np.linalg.norm(d2Nem, 2) + 1e-20)\n",
"\n",
" P2[:, 1] = np.squeeze(tmp * d2Nem)\n",
" R2[:, 1] = np.squeeze(tmp * Ad2Nem)\n",
"\n",
" if i == 0:\n",
" P2 = tmp * d2Nem\n",
" R2 = tmp * Ad2Nem\n",
"\n",
" #p=(1/dx_norm)*dx\n",
"\n",
" x_prev = x\n",
" AC_prev = AC\n",
"\n",
" P = np.concatenate([P1, P2], axis=1)\n",
" R = np.concatenate([R1, R2], axis=1)\n",
"\n",
" AC, x = NewtonOOS(x_prev, y, AC_prev, P, R, costf1, costf2)\n",
"\n",
" flag_is_New = True\n",
"\n",
" AC = mult1(basis, x)\n",
"\n",
" dx = x - x_prev #If TN, (x_prev+d_tn) - current approx. minimizer of quadratic model\n",
" dx_norm = np.linalg.norm(dx, 2) + 1e-30\n",
"\n",
" glob_grad_norm.append(np.linalg.norm(x_grad, 2))\n",
" \n",
" return x, AC, f"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"def plotting(x, y, A, n1, c00):\n",
" plt.figure(1)\n",
" plt.plot(y[0:50], label='original')\n",
" plt.plot(np.matmul(A, x)[0:50]-20, label='reconstructed')\n",
" plt.title('Orignal and reconstructed signal')\n",
" plt.legend(loc='best')\n",
"\n",
" plt.figure(2)\n",
" plt.plot(np.arange(0, n1), c00[0:n1], label='original')\n",
" plt.plot(np.arange(0, n1), x[0:n1]-3, label='reconstructed')\n",
" plt.title('Coefficients')\n",
" plt.legend(loc='best')\n",
"\n",
" plt.figure(3)\n",
" plt.title('Norm of gradients')\n",
" plt.semilogy(glob_grad_norm)\n",
"\n",
" plt.figure(4)\n",
" plt.title('Norm of subspace gradient with Newton iterations')\n",
" plt.semilogy(sub_grad_norm)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnWd0XNW5sJ896r3bltVlG8tYtiVbtlww4ADBQMDBhE7AaVw+ElJuGiT3JiGBtEuAhCSAA8SEQIJpSeiJ6Ta4W7LlbqsXF5UZ9Tazvx9nzmgkzUgjTZW0n7W0JM0pe5+ZM+9591uFlBKFQqFQTH4M/p6AQqFQKHyDEvgKhUIxRVACX6FQKKYISuArFArFFEEJfIVCoZgiKIGvUCgUUwQl8CcxQogfCCGe8ME4G4QQW709jnWs94UQX/bFWFMNIUS2EEIKIYLHeNyku88mK0rgTxCsN/sBIUSnEOKUEOJRIUT8SMdIKX8upVTCMUDwtsASQvxECPFXb53fGeo+mzgogT8BEEJ8G/gV8F0gDlgOZAH/EUKEOjlmTFraZGQivgdCiCB/z0ExeVECP8ARQsQC9wJ3SSnfklL2SSkrgeuAbOAW634/EUK8KIT4qxCiFdgwVOMTQtwqhKgSQjQJIf5XCFEphLjY7vjNQoi/CCHahBAHhRBFdsfeLYQ4ad12SAhx9Riu4QXrqsQkhPhQCDHfbtsmIcQfhBCvW8+9Qwgxy277JUKII9Zjfw+IEcZx9B4Y7ObeZL3GRLtjzhNCfCyEMAohaoQQG6yvx1nfi7PW9+x/hBAG67YNQoitQogHhBAtQogKIcRldufcIIQot15PhRDiZiHEPOAxYIUQol0IYbS7/keFEG8IITqANUPNVkNXBkKI+UKI/wghmoUQp60mlbXAD4DrrecvtbuOJ4UQDUKIOiHEffpDRQgRZL2GRiFEOXDFKJ/j963naBNCHBVCXGT3vvv9PlOMjhL4gc9KIBx42f5FKWU78AZwid3L64AXgXjgWfv9hRDnAn8EbgZS0VYKaUPGugr4u/X4fwG/t9t2ElhtPe5e4K9CiFQXr+FNYA4wDdg7dG7ADdZzJgAngPutc05Gu+7/AZKtc1g1ylhD34O7gM8CFwAzgRbgD9bzZ1nn9giQAhQAJdbzPGK91lzrsbcCX7Abpxg4ap3Xr4EnhUYU8DvgMillDNrnVyKlPAzcAXwipYyWUtqb426yXnMMMKLJRwgRA2wB3rJez2zgHSnlW8DPgeet519kPWQT0G/drxD4NKA/TL4CfMb6ehHwuRHGnQt8DVhqva5LgUoH+/nzPlOMghL4gU8y0Cil7HewrcG6XecTKeU/pJQWKWXXkH0/B7wqpdwqpewFfgQMLaS0VUr5hpTSDDwD6EIDKeULUsp667mfB44Dy1y5ACnlU1LKNillD/ATYJEQIs5ul1eklDut1/gsmuAFuBw4KKV8UUrZBzwMnBpluKHvwR3AD6WUtXbjf05o5p6bgC1Syr9ZV05NUsoSqwZ8A3CPdd6VwG+Az9uNUyWl/JP1vXoaTbhNt26zAPlCiAgpZYOU8uAoc/6nlHKbdc7do+z7GeCUlPI3Uspu6/x2ONpRCDEd7T38ppSyQ0p5BnjIem2grRIfllLWSCmbgV+MMK4ZCAPOFUKESCkrpZQnHeznt/tMMTpK4Ac+jUCycGyPTrVu16kZ4Twz7bdLKTuBpiH72AvTTiBcH9e6TC+xmj6MQD6DHzYOsZoNfmldprcyoBXaHzt03Ggnc5aMfI042J4FvGI378Nowms6kIGmUQ4lGQgBquxeq2Kwpmqbs/W9BIiWUnYA16M9aBqspqq8Mc55JJzN2RFZaNfRYHf9j6OttGDI+8vg6x2ElPIE8E20B+YZIcTfhRAzHezql/tM4RpK4Ac+nwA9wHr7F4UQ0cBlwDt2L49U+rQBSLc7PgJIcmUCVtPHn9CW9ElWc0QZI9jT7bgJzcxyMdoyPVs/rQvHNqAJOH0ewv5/Jwx9D2rQzCvxdj/hUso667ZZw09BI9CHJjB1MoE6F+aMlPJtKeUlaA/kI2jvnaO5OZtzBxBp9/+MIdeT6+J5atDunWS7a4+VUuo+lEHvL9o1OkVK+ZyU8jy090WiBRIMxV/3mcIFlMAPcKSUJjRb5iNCiLVCiBAhRDawGahFWxK7wovAlUKIlUKL7PkJrn+RotC+4GcBhBBfQNO8XCEGTeg0oQmxn7t4HMDrwHwhxHqrBvh1Bgs/V3gMuN8qTBBCpAgh1lm3PQtcLIS4TggRLIRIEkIUWE0Nm63HxViP/W9g1JBHIcR0IcQ6qy2/B2hHM/EAnAbShZPIKjtKgPVCiEghxGzgS3bbXgNShRDfFEKEWedXbHf+bGF1LkspG4B/A78RQsQKzYE9SwhxgXX/zcDXhRDpQogE4O4RrmuuEOJTQogwoBvosrsue/x1nylcQAn8CYCU8tdoERgPAK3ADjTt7SKrXdqVcxxEc2D+HU0LawfOoAml0Y49hGbD/gRNqCwAtrk4/b+gmQrqgEPAdhePQ0rZCFwL/BLtgTFnDOPq/BbNMfhvIUSbdfxi6/mr0Wzc3waa0QStbk++C03TLkdzpD4HPOXCeAa0h0O99ZwXAP/Puu1d4CBwSgjR6PhwQLOz96K9109j5+SWUrahOeqvRDONHAfWWDe/YP3dJITYa/37ViAU7b1vQRPIuhP0T8DbQCmaM31QYMAQwtA+h0bruNOAe4bu5Mf7TOECQjVAmZpYTUJGYI6UssLf81FMTtR9FlgoDX8KIYS40momiEJbLRzAQWidQuEO6j4LXJTAn1qsQzM11KOZR26Qaomn8DzqPgtQlElHoVAopghKw1coFIopQkAVl0pOTpbZ2dn+noZCoVBMKPbs2dMopUwZbb+AEvjZ2dns3r3b39NQKBSKCYUQwmmWtD3KpKNQKBRTBCXwFQqFYoqgBL5CoVBMEQLKhq9QKPxLX18ftbW1dHePVqVZ4Q/Cw8NJT08nJCRkXMd7ROALIZ5Cq9N9RkqZb30tEXgerTpiJXCdlLLFE+MpFArvUFtbS0xMDNnZ2WjFSRWBgpSSpqYmamtrycnJGdc5PGXS2QSsHfLa3WideOaglfB1WolPoVAEBt3d3SQlJSlhH4AIIUhKSnJr9eURgS+l/BCtMqA969Aq/WH9/VlPjKVQKLyLEvaBi7ufjTedttOt9bhBK6c63dFOQojbhRC7hRC7z54968XpKACOnGplR/nQBkQKhWIq4JMoHWvhJIdFe6SUG6WURVLKopSUURPFFG7yf28d5Z6XD/h7GgqFW1x++eUYjcYR9/nRj37Eli1bxnX+999/n8985jPjOjaQ8WaUzmkhRKqUssHadf6MF8dSuEhNSydn213qmaJQBBxSSqSUvPHGG6Pu+9Of/tQHM5pYeFPD/xdwm/Xv24B/enEshQtIKalr6aKtu5+efrO/p6NQOOTBBx8kPz+f/Px8Hn74YSorK5k7dy633nor+fn51NTUkJ2dTWOj1jTsZz/7GXPnzuW8887jxhtv5IEHHgBgw4YNvPjii4BWtuXHP/4xixcvZsGCBRw5cgSAnTt3smLFCgoLC1m5ciVHjx71z0X7CE+FZf4NuBBIFkLUAj9Ga4e2WQjxJbQWd9d5YizF+Gnt6qejVxP0LR19zIgL8vOMFIHMva8e5FB9q0fPee7MWH585Xyn2/fs2cOf//xnduzYgZSS4uJiLrjgAo4fP87TTz/N8uXLB+2/a9cuXnrpJUpLS+nr62Px4sUsWbLE4bmTk5PZu3cvf/zjH3nggQd44oknyMvL46OPPiI4OJgtW7bwgx/8gJdeesmj1xxIeETgSylvdLLpIk+cX+EZao2dtr8b23uYERfux9koFMPZunUrV199NVFRUQCsX7+ejz76iKysrGHCHmDbtm2sW7eO8PBwwsPDufLKK52ee/369QAsWbKEl1/W2veaTCZuu+02jh8/jhCCvr4+L1xV4KAybacQdS1dtr+bO3r9OBPFRGAkTdzX6A8AdwgLCwMgKCiI/v5+AP73f/+XNWvW8Morr1BZWcmFF17o9jiBjKqlM4WoNw4I/KYO5bhVBB6rV6/mH//4B52dnXR0dPDKK6+wevVqp/uvWrWKV199le7ubtrb23nttdfGNJ7JZCItLQ2ATZs2uTP1CYES+FOIOmMXBmveRlO70vAVgcfixYvZsGEDy5Yto7i4mC9/+cskJCQ43X/p0qVcddVVLFy4kMsuu4wFCxYQFxfn8njf+973uOeeeygsLLRp/ZOZgOppW1RUJFUDFO/x1Wf3crihlermTr5yfi7fX5vn7ykpAozDhw8zb948f09jTLS3txMdHU1nZyfnn38+GzduZPHixf6eltdw9BkJIfZIKYtGO1bZ8KcQtcYu0hIiaO/pp1lp+IpJwu23386hQ4fo7u7mtttum9TC3l2UwJ9C1LV0MW/eNBrbe5UNXzFpeO655/w9hQmDsuFPEbr7zDS29zAzPoKkqFCaVJSOQjHlUAJ/itBg0kqqpsVHkBQdqpy2CsUURAn8KYIeg5+WEEFiVChNqp6OQjHlUAJ/iqDH4KfFR5AcHUZHr5nuPlVPR6GYSiiBP0Wotcbgz4gLJykqFEDZ8RUKB2zatIn6+nqPne/hhx+ms7Nz9B3t8FZ5ZiXwpwh1LV1Mjw0nJMhAolXgq9BMRSAjpcRisfh83JEEvtk89lXxeAS+t1ACf4pQb+xiZnwEAEnRWk2RRhWaqQgwhpZCfuaZZ1ixYgWLFy/m2muvpb29HdCqZK5cuZJFixaxbNky2tra6O7u5gtf+AILFiygsLCQ9957D9AE+Pr161m7di1z5szhe9/7HqAJ7w0bNpCfn8+CBQt46KGHePHFF9m9ezc333wzBQUFdHV1kZ2dzfe//30WL17MCy+8wIUXXoieINrY2Eh2drbtfN/5znfIz89n4cKFPPLII/zud7+jvr6eNWvWsGbNGgD+/e9/O7ymt956i7y8PBYvXmwr7uZpVBz+FKHO2EVBRjyAzaSjNPzJw9bjjdSburiuKMNzJ33zbjjl4e5oMxbAZb8ccRe9FPLs2bNZv349W7ZsISoqil/96lc8+OCD3H333Vx//fU8//zzLF26lNbWViIiIvjtb3+LEIIDBw5w5MgRPv3pT3Ps2DEASkpK2LdvH2FhYcydO5e77rqLM2fOUFdXR1lZGQBGo5H4+Hh+//vf88ADD1BUNJC4mpSUxN69ewF47LHHHM5748aNVFZWUlJSQnBwMM3NzSQmJvLggw/y3nvvkZycTGNjI/fdd9+wa/re977HV77yFd59911mz57N9ddf74l3exhK4E8BLBZJg6mLKxamApAUrdvwlYY/WXjwP0c5eqqN9YVpBAdN7IW7Xgr5tdde49ChQ6xatQqA3t5eVqxYwdGjR0lNTWXp0qUAxMbGAlpp5bvuuguAvLw8srKybAL/oosustXYOffcc6mqqmL+/PmUl5dz1113ccUVV/DpT3/a6ZxcEcBbtmzhjjvuIDhYE6uJiYnD9tm+fbvDazpy5Ag5OTnMmTMHgFtuuYWNGzeO/maNESXwpwBn23voM0ubSSc6LJjQIINy2k4SOnv72V9rot8iOdTQysL0eM+ceBRN3FvopZCllFxyySX87W9/G7T9wIGxrzr00sgwUB45ISGB0tJS3n77bR577DE2b97MU089NeKcAIKDg22+he7u7jHNw9k1lZSUjOk848XrqoAQolIIcUAIUSKEUJXR/ECtNQY/3SrwhRAq+coL/LOkjmd3VPl83L1VRvotWhHEnRXNPh/fWyxfvpxt27Zx4sQJADo6Ojh27Bhz586loaGBXbt2AdDW1kZ/fz+rV6/m2WefBeDYsWNUV1czd+5cp+dvbGzEYrFwzTXXcN9999lMNjExMbS1tTk9Ljs7mz179gDYWigCXHLJJTz++OO2qpvNzc3DzufsmvLy8qisrOTkyZMAwx4InsJXa781UsoCV6q5KTyPLQY/IcL2WmJUqGqC4mEeff8kj39Q7vNxd1Q0EWQQzIgNn1QCPyUlhU2bNnHjjTeycOFCm+kjNDSU559/nrvuuotFixZxySWX0N3dzZ133onFYmHBggVcf/31bNq0aZBmP5S6ujouvPBCCgoKuOWWW/jFL34BaL1w77jjDpvTdijf+c53ePTRRyksLLT11QX48pe/TGZmJgsXLmTBwoW2h8/tt9/O2rVrWbNmjdNrCg8PZ+PGjVxxxRUsXryYadOmefjd1PB6eWQhRCVQJKVsHG1fVR7ZOzz2wUl++eYRyu69lOgwzYp361M7MXX28s+vnefn2U0OuvvMzP/x2wAc/dlan9rRr3vsE3r6zcyZHsM7h0+z538uwaA3PhgjE7E8cqDR3t1HeWMHs1KiiQrzvNXcnfLIvrgrJfBvIcQeIcTtQzcKIW4XQuwWQuw+e/asD6Yz9ahr6SIuIsQm7AFVQM3DHDnVhtkiMVsk9cax2XXdobvPTEmNkeLcJJblJNLS2cfJs+0+G18xnJ5+q30/ADPZfSHwz5NSLgYuA74qhDjffqOUcqOUskhKWZSSkuKD6UwepJRs2lbB2baRo23sY/B1kqKUDd+THKw32f6ubvZdks2+aiO9ZgvFOYksy9aiQnZMIrPORKTXrAl8XfAHEl4X+FLKOuvvM8ArwDJvj2nP6/sbBvVynUwcOdXGT149xDPbR3YU1hm7SBsq8KPD6Ooz09k7+du6+YKyulaCrWYUXwr8HRVNCAFF2YlkJUUyLSaMXZXuCfxA6oI3Eent957Ad/ez8arAF0JECSFi9L+BTwNl3hzTnsb2Hr763F4eff+kr4b0KfuqjQDsGkWjq2vpIj1huIYPqretpzhYb6IoO4HQIANVzR0+G3dHeTPnpsYSFxGCEIKlOYnsrGget2AIDw+nqalJCX03GNDwPWvSkVLS1NREeHj4uM/h7Tj86cArQgh9rOeklG95eUwbuiDcXdXiqyF9yr5q7br21bTQ228hNHj487u1u4+2nn5mxg++SQaSr3rJSIz0/mQnMX1mC0ca2tiwKpszrT3U+EjD7+k3s7e6hZuLs2yvFeck8vr+Bmpbusb1uaanp1NbW4vyp42fBmMXZgkCsDSHY5V/HiE8PJz09PRxH+9VgS+lLAcWeXOMkdBtmUdPtdLW3UdMeIhb55NSevTDc5d9NUYiQ4Po7DVzoM7EkqyEYfvY6uDHD/7y2wqoqWxbtzl+up1es4X5M2M5drrNZyad/bUmevotLMsZyOhcarXj76psHpfADwkJIScnx2NznGp09PRz2Y/fJm9GDEdOtfH2N89n7owYf0/LxsTOwR6FnRXNxIQFY5ED5o/xsu1EI4U/+w/HTjtPyPAlpq4+Tpxpt9VOcWa3dRSDD5CsF1BTJh23KbM6bPPT4shMjKSqqdMnJpEd5U0AgwT+3OkxxIYHT6p4/IlEnfX7dv45WgBKeYBFTE1agW/q6uPwqVZuKs7EINw363x47CzGzj7+5x9lAWHf3F+rPcAunjed3JQop3Z8/QYcatIZ0PCVwHeXg3UmokKDyEmKIjMxkrbufkxdfV4fd0dFM3Onx9g+SwCDQbA0O5GdbjpuFeOjtkVb3a2ekwxAeaPv/DmuMGkF/p6qZqSENXnTyJsRy54q974AZfUmQoMM7Kxo5qW9dR6a5fjZV21ECFiYEcey7ER2VTZjsQx/ENW1dBEabCA5anDGYWRoEOEhhoBodfhqaT3rfr+Vjp6JGTF0sL6Vc2fGYjAIMq1mlKom75p1+swW9lS1UJw7vEDX0pxEys92jBquq/A8Nc2agjV3RgzTY8MCLidi0gr8HRXNhAYZKMiIpyg7gX3VRvrN4wuTklJSVtfK+sVpLM6M5+dvHMbY6V/NeF91C3OmRRMbHsLS7ERau/s56sDcpIdkDs28FEKQFBXm9+Sr7eVNfHtzKaW1pgmplZqtBcvmz9QqMWYmaQLf23b8A3UmOnvNFOckDdumm3h2T8D3c6JT29JJWLCBlOgwcpOjKT+rNHyfsLOimUUZcYSHBLEkK4HOXjNHTo3P/l7T3IWpq48F6XHcf/UCTF19/Oqtox6esetIKdlXY6QwQ3PS6l9wR3b8OmPXMHOOjr8LqJ08285/PbOH9MQIQoIEO8onnoCqaOygs9dMfppV4Cf6RuDrNnp7+71O/sw4wkMMKgHLD9RaQ6CFEOSmRFF+tj0gTMA6k1Lgd/b2c6DWZPsyFGW7p/EcqNOccgvS4piXGssXV2Xzt53V7PFTuGdlUyfGzj4KM7UyuOkJEaTGhTv8gte1DE+60vFnAbWm9h6+8OddBBsEmzYsY2F6PDsqmvwyF3c4aHPYajXZI0ODSY4Oo9rLJp0d5U3MSokiJWZ4cbDQYAOLMxPcTsBSjB1N4GsP/dyUaFq7+/2+irZnUgr8fdVaudhl1uVuWrwmEMfruD1QZyIkSNjCq7558TmkxoXzP/8oG7eZyB1KarTrKLAKfCE0R92uIQk3Pf1mzrT1DAvJ1EmKCvOLDb+7z8yX/7Kb063dPHFbEZlJkRTnJHKg1jTh7PhldSZCgw3MSom2vZaZGOFVDd9skeyubKE4d7g5R2dpdiKHG1pp7fa+81gxQE1Lpy3JMTdFq6EfSGadSSnwd1Q0YxAMiktfkpUwbo28rM7EOdNjCAsOAiAqLJgfXzmfww2tbPq40hNTHhP7qo1EhQYxZ9pAfO+ynETOtPUMEjSnTFoRL2cmneRorYCaL5ecFovkW8+XUFJj5Lc3FFCYqX1GxblJ9Fuk31ZN46WsrpV5M2IIsauOmZUU5VWBf6i+lbaefoodmHN0inMSsUgm3Ps5kWnr7sPY2WfT8Gcla0pAIIVmTkqBv7Oiify0uEHVIYuyEmgwddvCFF1FSklZvYkFVhutzqXzp/OpvGk89J9jNJh8W6tnX7WRRRnxBNk5YnXzlb1Zp85JDL5OYlQoPf0WOnp9V9XvV28d4c2yU/zw8nmszU+1vb4kK4Egg5hQZh0pJQfrTcwfcm9kJEZSb+qy1VTxNPp7tHwEDb8wM4Fggxi17IbCc+jft4xE7fuWlhBBaLAhoEIzJ53A7+k3s6/aaKscqDNeO35tSxfGzr5hX2ohBPdeNR+zlPz01UPuTXoMdPWaOdzQarPf68xOiSY+MmTQF3wgy9axwE+yJl/5qpn5i3tqefzDcj6/PIsvnTc4mzM6LJgFaXETynFb29JFa3c/+TMH3xtZiZFIyZiVC1fZXt5MdlIk02Od11SJCA1iQXpcQCZg9Zst/OWTyklnbqq1hmTqGn6QQZCdFKk0fG9ywEG6OUDejBgiQ4PGvMQts3PYDiUjMZK7PjWHN8tOsfX4qP1dPEJZvda7VI/Q0dETbuwddXXGLoSA1DgnAt+asNPog/IKZovk4S3HKMiI58dXnuuwREVxbiKltUa6fLjicAf93tAdtjp6aGZVk+c1O4tFsquy2WE45lCWZSeyv9YUcHXZ/1Vaz4/+eZBN2yr9PRWPUmNNurIvVBhooZmTTuDrJo2lQzT84CADhZnx7K4cm8A/UGci2CDIc1IP4yurc0mNC+eRd4+Pb8JjpMRaIqIgc3ij6mXZiVQ2dXKmVbPd1xu7mBYT5rCoGgwUUPOFhv/+0TPUtnTx5dU5TrtBLc9Jos8s2Vs9MezOZfXavXHO9MH3hh6a6Y0iakdOtWHq6nOYcDWUZTmJ9JotlNS4V1bEk0gpeeKjCgBe2FPjMFlwolLb0kVESJBNkQLNcVvd3EmfH4I7HDHpBP5Oa7p5gt2brrMkM4Ejp1ppH0MkyIE6E3OmxxAeEuRwe2iwgS+vzmVHRbNPHGT7alrISIyw1cKxZ6l1VaMnMNU5aHxij56S3+QDDf+Z7VVMiwnj0vkznO5TlJ2AQQzUiPEFLR29/LOkblz29rK6Vof3xrSYMMKCDV5x3Or2e0fx90MpykpEiNHLZ/uS7eXNHGpo5bzZydQ0d02qXIFaa4SO/eo1NyWafov0aY+EkZhUAr/fmm7u7MuwJDvRWkjNNcGsZdiaWDBkyT6UG5dlkBAZwqPvnxjznMfKvmrjMHOOzvyZsUSEBNm+4CPF4IMWlgl4PU64qqmDD46d5cZlmYOiWYYSEx7C/JlxbPeRENh6vJG1v/2Qb/y9hPWPbuPEGddtrfq9kT9z+L0hhLAVUfM0O8qbSYuPsNmJRyIuMoS502MCKoP5ya3lJEWF8vubCokJC+aF3TX+npLHqHXQdyLQQjMnlcA/3NBGe0+/U4FfmBmPELhs1qk3ddPS2efQfm9PZGgwG1bmsOXwGY6OM5vXFRpMXTSYuoc5bHVCggwsyUpgR4VWV6fe1O00Qgc0x15kaJDXs23/ur2KICG4qThz1H2X5yZSUmMct925qqmDsjrTiKGmPf1m7n/9ELc8uYPosGDuvWo+dS1dfOaRj3huR7VLYaqnW3to6uhlvgOBD5CVFOkVrW5/rZGibMcPfEcsy0lkT1VLQOQ3lJ9t550jZ7h5eRbxkaFcWTCTN8oaJo3z1j7pSkcPzaxoDAzH7aQS+KMtd2PDNY3HVdPLgdqBsrejcdvKLCJDg3jsA+9119Lt93rsuiOWZidy9HQb5Y0d9PZbRtTwQS+v4D2TTlevmc27a7l0/owRo0p0inOS6O0fn925paOXax79hM88spVPP/Qhj31wktOtgxuKHzvdxmf/8DF/+qiCzy/P4rW7VnPbymze+ub5FGUl8oNXDnDHX/fQMsqq52D9yPdGRqIm8Ed6eJTVmSi6b4utwuJo9JktnGrtJmsMde7XFcykq8/Mva8edPkYb/HnbZWEGAx8frnWsOW6ogy6+yy8Vtrg55m5T2t3H6auvmEaflxkCElRoVNHwxdCrBVCHBVCnBBC3O3NsXZUjB6uphVSa3EpQ7aszkSQQTAvdWSTDkB8ZCg3LcvkX6X1Xut4VFJjJDTYwLkjzGdpTgJSwr9KtIqeown8RC8XUHu1tB5TVx+3LM8afWc0P4QQWlG1sXLvqwcxdvby3UvnEhMezC/fPMKKX7zDrU/t5J8ldWzaVsGVj2zlTGs3T95WxM8+m09EqGZ/nx4bzl++uIwfXj5rFUquAAAgAElEQVSPd4+cYe1vP2TbCeeRV2V1rQiB03sjMzGSzl7ziO/t6wcaaGzvoayu1aXrO93ajUUyol9mKEuyErnzwlls3l3La/vrXT7O0xg7e3lxTy3rCmbaykEsSo/jnOnRbJ4EZh09JNNR0xmtps4UEPhCiCDgD8BlwLnAjUKIc70xlh6uNpozqygrkQ4XC6kdqDMxZ1q0U4ftUL68OheDgI0flru0/1jZV21k/sxYp1E3AIUZCYQECf5Ron25RzLpACRHea+AmpSSv2yv5Jzp0Sx3IaoEIC4ihHkzYsccj/+fQ6f5R0k9X10zm6+umc3Ld67i3W9fwJ0XzubE6Ta+8fcSfvLqIVbOSuKtb57PRfOmDzuHwSD4yvm5vHLnKqLDgrnlyR38eVuFw/HK6k3kJkcRFea4aVyWC1UzPzymtRF0NV6/3qhnTrsu8EErBVKYGc89Lx/wWfvFoTy3s5quPjNfWj2QfyGE4LqiDEpqjBz3YWOhs209Hs8ur3UQkqmTmxxN+RQx6SwDTkgpy6WUvcDfgXXeGOj4mXaMnX22+jnO0MstjGbWsTnlXDDn6MyIC+eaxels3l3j8VrkfWYL++ucO2x1IkKDWJAWZxM0owkHbxZQK6kxUlbXyueXZ42pNWRxbiJ7q1tcbgJt7OzlB68cIG9GDF9dM9v2em5KNN+5dC5bv/8pnvtyMX+6tYinNix1WHDMnvy0OF67azWXzJvOva8e4omPhj/AD45yb9iqZjpx3Da293CwXtPs9QS50ai3NbMZm8APCTLwuxsKQcI3/r7P5/WfevstPP1xJefNTiZvxuAV0WcL0wg2CF7YU+uTufzpw3KW3r+Ff5R4tqdFbcvgpCt7clOiaGzv9UlTnNHwtsBPA+zXa7XW12wIIW4XQuwWQux2p3HyTqv9fqT6IqA9gafHho1aSK3B1E1TR++oDtuh3H5+Lr1mC0850QzHy9FTbXT3WZw6bO3RwzNjwoOJHaWPb1J0GE0dntd4AJ75pIrosGCuXjy2psvLc5Po6bew3+pDGY2fvnaI5o5eHrh2kcPVj8EgWDk7mUvOne7ygyciNIg/3LyYy/JncN/rh9n44YBvpqm9h3pT97AMW3v0L74zDV9P1AsNNlBndE3rtpXKGKPAB83UcN/V+eytNvK7d3yTM6LzxoEGTrf2DMuuBq3V5kXzpvHy3lqvxqpbLJKfvXaI+984DMB7RzzbpL2mpZPI0CASIod/33JTAqemjt+dtlLKjVLKIillUUpKyrjPs6OimdS4cIdLKnuEEBRlJbJnlFC1gSzKsQn83JRoLs9P5a+fVHk0+kAPJXVF4OsPPVcEQ3J0KH1mSZuHozia2nt4bX8D6xenDapp5Ap6WYztJ0e347975DQv763jzgtnjfmzGo2QIAO/u7GQKxak8vM3jvDo+5rQ1zXz+SOE64aHBDEjNtypwP/w+FkSIkMozkl02aRTZ+wiMSrU5ncYK+sK0vjcknQeee/EuHwk40FKyRNby5mVEsUF5zj+fl+7JIPG9l7eO3LGK3Po6TfzjedLeHJrBV9Ylc0VC1LZUdHkUSWntqWLjIRIhwpFTnLghGZ6W+DXARl2/6dbX/MoUkp2Vmj2e1c0uCVZCdSbum1LZEeU1ZkwCEZ0kDrj/104i7aefv66vWrMxzpjX7WRlJgwl4T4EmvCjSv72pKvPGzHf353Db1mi8vOWnsSokLJmxEzalKOqauPe14+wDnTo/nap2aPuO94CQky8NsbCrhy0Ux+9dYR/vDeCVvT8vkjaPigmXUcmXSklHx0vJFVs5PJSIwck0nHWeVTV7n3qvlkJ0XxredLRo1E8gQ7K5opq2vli+flDOu6pnPh3BRSYsLYvNvzZp227j6+8OddvFpaz92X5fGjz5zL8llJnG7t8WiehKMYfJ3MxEiCDCIg7PjeFvi7gDlCiBwhRChwA/AvTw9S1dTJmbYel7IPAVsc80hmHc1hGzMubSo/LY7Vc5J5amuFx+qYlNQYKciId+mBFhcRwg1LM7k033lWq46tgJoHs23NFsmz26tZnps4rOyAqxRb48dHWubf//ohGts1U45eutobBAcZeOi6RawrmMn/vX2UJz+qICMxgriIkc1lmU5i8Y+cauNsWw/nn5NCWnwELZ19dPaOvsKqN3Yx00ldJFeJCgvmkRsLaWzv4fsv7fd6aewnt1aQEBnC+kLnZr3gIAPrF6fx3tEznGnrdrrfWDnT2s31j29nZ0Uzv7l2EXdcMAshBMutcsKTq5xauzr4QwkNNpCZGDn5NXwpZT/wNeBt4DCwWUrp8YBgs5RcXZjGihHKxdozL1XLSHVWOVNKyYG6VrdMBHdeOJvG9l6e3+V+yFlLRy/ljR0umXN0frF+AdcVZYy6n62A2hg1fIu1GNqv3zrCszuq+ODYWU6ebae7z8x7R85QZ+zi1hXZYzqnPctzk+jqMzu1479/9Aybd9fyX+fnsjDd9fdlvAQHGXjwugLWF6bR1NE7ov1eJzMxklOt3cMe+np0zvlzUmxCYjQtX0pJXcvIpTJcJT8tju+vzePfh07zh/e8lx1e2djBfw6f5ubirFEVp2uXZGC2SP6xzzMGgNOt3ax/9GMqmzp44rYirlky8MCZPS2apKhQj5V1MHX20dbdP2L2c25yYIRmjs24Og6klG8Ab3hzjFkp0Tx0fYHL+4cEGbjgnBSe31XDdUUZwwT76dYeGtt7hlVBHAvLcxMpzknk4S3HWFcwk/jI4bV9XOU/h08DjBqhMx5sBdTGuLx/dX89D285jkHA0PpXoUEGpseGccm5w0MfXWWgvn/ToEY2vf0WXthTw2/+fYzZ06L5+kVzxj3GWAkyCP7v2kXMS411qXiZHqlT29LJbLtmNR8db+Sc6dHMiAu3md1qjV3MGWE11NrdT0eveVQflat86bwcDta38sC/jzEtNtwl5WCsvH3wFFLCzctHz7CePS2aJVkJbN5dy1dW544pqssRT22roMHUzUv/byUFGYMVAiEExbmJ7CjX7PjujqVXydTr4DsiNyWKj040YrbIQX0sdLp6zYSHGNyey2j43WnrL+67Op+kqFD+65k9wzJND4xQEtlVhBD85Kr5mLr6+M2/j437PKdbu7n/9cMUZsa7bLIaCwM2fNdNOt19Zn791lHmz4zl6H2X8ck9n+KFO1bw0PWL+PYl57B+cRo/W5c/Yt2c0UiKDmPOtGi2W+PxzRbJS3tquejB9/nhK2VkJ0Xy6M2LXc6R8BRB1lh9V1YVmQ5i8bt6zeysbOb8OZoDM81FDV/f7gkNH7T781fXLGT1nGTuefkA7x457ZHz2lNaayQjMcJpee6hXFeUzokz7exzs7pnb7+FF3fXclHetGHCXqc4J4l6U7ctnNIdRgrJ1MlNiaa33+LQb9hvtnDbn3fyw3+UuT2X0ZiyAj85OozHPr+Es+09fO25wbHJNoetkzoprjIvNZZbV2Tz7I4qWyr+WJBScvdL++npN/Obaxc51AzcJSw4iJiw4DFl2/55WyV1xi5+ePk8QoIMpMZFsDQ7kasL07nrojn88pqFfHqEqpiuUpyrRVP9q7SeTz/0Ad9+oZTY8BD+vGEpL/2/lSNqxIGAruHbOwe3VzTR229htTViZVpMOMEGMWqkznhj8EciNNjAo7cs4dzUWO58dq/LRQVdpbTGRMEYVqVXLJxJREgQL+91z3m75fBpmjp6uXGZ85WFvkLzhB1/pKQrnVxrpM5JB6GZv3vnODsrminK8vwKfihTVuADLEyP5xdXL+CT8iZ+/sYR2+tldSZmpUQTGeq+xetbF59DfGQoP/nXwTE7yDbvruG9o2e5e22eLZbXGyRGu55t29Tewx/fO8FFedNYOTvZa3MCTQvr6DXz9b/twyAEj92ymNfuOo81edO8vvT1BElRoUSFBg3S8D861khosMEWOhtkEKTGh4+q4debdIHvXpTOUKLDgnlqw1KmxYTzxU27PBYrfqZNaye6KN31VXJ0WDBr8lJ4++BpzG7Uyf/bzmpmxoVzvpMwUIBzpsUQHxniETt+bUsX0WHBIzrxB2LxB9vxt51o5JH3TvC5JemsH2O+yniY0gIf4Jol6WxYmc1T2ypsmsWBuuE9bMdLXGQI3187l12VLfyzxPVaJjXNnfz01UOsyE1yy/npCkljyLb97TvH6ewzc8/l87w6J4A1edNYX5jGw9cX8NY3z2dtfuqEEPQ6QggyEiMHlTP48PhZinMSB5mi0uIjRtXw64xdhAYZSI4aOUt4PKTEhPGXLy7DIAS3PrXT1kDHHUprtBWtM5OKM9bmp3K2rWfcvSVqmjv56Hgj1y3NGHFFbDAIlmUneqSHsqM6+ENJjg4lJjx4UGjm2bYevvH3EmalRPPTdfPdnocrTHmBD/DDK+ZRnJNos2WeaevxaBLPtUsyWJQex8/fOOxS8xWLRfLdF0sRQvDrzy10Gr/sKRKjwmh0wYZ/4kw7z+6o5qZlmcye5r0Vh050WDAPXl/AZwvTvGLO8gX2dfHrjV2cONNus9/rpMWPHotfb+xmZny41+6F7OQontqwlOaOXjb8eRdtbiYNltYYCTKIMX+PPpU3jdBgA2+Wja+C5vO7ajAIXHJCF+cmUdPcNerDtrTGOGIhPUdlkYcihCA3ZaDdodki+dbzJbT39PGHmxZ7xJrgCkrgo0Xt/OHmxSRFhXLHM3sBWDCGpehoGAyCe9flc6ath0dcSGt/+pNKtpc387+fmeew+p6nSY4OdcmG/8s3DxMREsQ3LvZdZMxER6+LryVbWcMxh5ga0hIiON3WPWLXrbqWTo/a7x2xKCOeP968mKOn2/j9u+6Fa5bUGMmb4bxTnDOiw4I5f04Kb5edGrMJtN9sYfPuGi6cO82l90o3q43UYa3fbOHOZ/dy+192O1wFSylHTLqyZ5ZdaOaj759g64lG7r1qPnOdtE/1BkrgW0mODuPxzxchBIhxZtiOREFGPNcXZfDk1ooROyudPNvOL988wpq5KV4JlXNEUnQoLR29I/YX/fhkI1sOn+HONbMctldUOCYzMZKefgtn2nr48Hgj02PDOGf64NVRenwEUsIpk3NTiqbhe1fgA1w4dxoX5U3jJTdq21gsktJa45jNOTqX5c+g3tRNqYu1lHTePXKGM2093LDUte/NvNRYYsODR6zM+vqBBuqMXXT0mh0W0TN19dHe0++SYpabEsWp1m7eO3KGB/+jhWv76juuowS+HQvS4/j9TYv52prZTsveusN3184lIjSIe1917MDtN1v4zgulhIcE8ctrFvrMXp0YFUa/RTqt/WOxSO5//TBp8RF8cdXwAlgK5+iCoKKxg63HG1k9J2XY56qHZtY6KaLWZ7Zwus03Ah/g+qXu1bYpb+ygrbufReMU+BfPm06wQfDmgbGZdf6+q4ZpMWF8Km+aS/sHGQTLcpzb8aWUPP6BVgfoigWpPP1x5bByFDXNekjm6J+N7ri989m9ZCVFcf/VC3zuk/KN4WgCccm5091KGBqJ5Ogwvn3JOfzk1UP87p0TRIUFccrUTUNrN6dNWlRDg6mb395Q4FJ3KM/NS29m3uswQezlfXUcrG/ltzcU+DzufaKTlaSF471xoAFTVx+r5wyPbNKTr5zZ8U+ZupES0jwcoeOMC84ZqG0znvDaUmsc/Xg1/LjIEFbOTubNslPcfVmeS0Kx3tjF+0fPcOeFswkeQ/5HcU4SWw6f4XRr97Dv3LYTTRxqaOXX1yykIDOeN8oaeGJrOd+9NM+2jyshmTp6f1uzlPz+psIxFxX0BErD9zG3LM8ib0YMD205xn2vH+avO6o4VN9KcJBgeW4S9302n6sWzfTpnEYqoNbVa+aBt4+yMD2OKxf6dl6TgbT4CAwCXt5bhxCwes7wUMFUqyB35jyst5VF9r4/B7QyEtcsTh93bZvSWiPRYcHMciOU+PL8GVQ3d3KowbVuYJt312CR2upkLIwUj//4hyeZFhPGusKZnDM9hssXpPL0x1WDtHxXkq50cpKjyJsRw33r8kctvOctlIbvY4KDDDx/+woaWrtIjY0gNiLY76GGSVHOC6j98s3DnGrt5nc3Fno9WmgyEhqsJabVGbtYkBZne7jaExYcxLSYMKcafp3ROzH4I3FtUTqPfXCSl/fWcccFs8Z0bEmNkQVpcW5FVl1y7nR+8MoB3io7NapwNFskm3fVsHpO8piDHM5NjSU6LJgdFc2sKxho1VFWZ+Kj4418f22erTDf1z81h9f3N/Dk1gq+c+lcQNPwY8NHjsHXCQsO4q1vnj+m+XkapeH7gbjIEPJmxBIXGeJ3YQ8D9XSGFlB77+gZnv6kii+uyvFKWYepgp5xe/45zhPV0hKcx+J7I8t2NGalRFOUlcDm3TVjipbp7jNzuKF13PZ7naToMIpzknjDBTv+h8fOUm/qHjGz1hnBQQaKshOGReps/LCc6LBgbioeOOfcGTFcsSCVTR9XYuzUvis1LoRkBhJK4CtIiBxeQK2pvYfvvbifudNj+N7auf6a2qTAJvAdmHN0Rkq+qjN2kxQV6nP/yXVLMyg/28HeMZRcONzQSp9Zjtt+b89lC2Zw8mzHqP1u/7azmqSoUC520KfYFYpzkjh5tsPWlrSmuZPXDzRwU3HmMM39rotm097Tz5NbtY52I5VFDkSUwFcQGmwgNjzYVkBNSsndLx/A1NnHw8pR6zbLZyUye1o0hZnOa6WkJUTQYOx2GBqrNT7xvVC5YkEqkaFBbN7lem2bEjcdtvZcOn8GQsCbZaec7nOmtZt3jpzhc0vSHba3dAXdjq9H6zy5tQKDgC+syh62b96MWC5fMINN2zQt35Wkq0BCCXwFoEUQ6clXz++q4T+HTvO9tXOZ5+F8hKnI1YXpbPnvC0YUSOnxEfSaLZx1kPFcb+waVx9bd4kKC+YzC1N5bX89HS62wCytMTI9NowZce77G6bHhrMkM8GpwLdYJL/59zHMFjlmZ609C9LiiAwNYkd5My0dWg+LdQVpTqt8fv2iObT19PPAv4/S2WsesSxyoKEEvgLQInWa2nupaOzg3lcPsWp2koq59yG2WPwhjlspJXV+0vBBK1HQ0WvmdRdj4vXObJ5ibf4MDje0Utk4uOhYv9nCd14s5fndNdx+fq5bxQVDggwsyUpgR0UTz2yvoqvPzO3n5zrdP29GLJflz+DZHdWAaxE6gYLXBL4Q4idCiDohRIn153JvjaVwn6ToUE63dfOt50sIDTbwwLWLVFSOD9FDLofa8U1dfXT2mn0aoWPPkqwEclOieGH36J3bjJ29VDZ1uu2wtWettU2nvZbf3Wfmjr/u5eW9dfz3Jedwz2V5zg53meW5SRw73c6TWyu4KG/aqK05v37RHHRftrLhD/CQlLLA+uPVrlcK90iKDqP8bAclNUZ+fvUCl5tWKDyDs0YodbYYfP98HkIIrivKYFdly6ilk/VSCJ7U8NMTIlmUHsdb1mJqelPyLYdPc+9V8/n6RXM8Eumm19UxdfXxXy6Eoc5LjWWt1cegBL5iwqH3tl2/OI0rFqb6eTZTD72eet2Q8gr1Ri3xyV8mHYD11mqlL+wZ2XlbUm1ECPc6xTlibX4qpbUmyupM3PzEDnZWNvPw9QXctjLbY2MsTI8nIiSIgox4lma71ojkvqvzefyWJcSEjx6DHyh4W+B/TQixXwjxlBDC4bsohLhdCLFbCLH77NmzXp6OwhmrZifzqbxp3HuVb+pyK4aTFh8xTMO3Zdn6UYucFhvOmrkpvLSndlBnuKGU1hqZnRLtcQF4mdWsc82jH3P0VBuP37KEzxamjXLU2AgNNrDx1iU8eN0il1cMydFhHuns5kvcEvhCiC1CiDIHP+uAR4FZQAHQAPzG0TmklBullEVSyqKUFOdxygrvsjw3iac2LJ1Q2spkw1HyVZ2xi9Bgg20F5i+uLcrgTFsPHxxzrJRJKSn1sMNWJzs5ivy0WEKCDDz9xWVc7KVaV6vnpHi1s1wg4FZpBSnlxa7sJ4T4E/CaO2MpFJOdtPgIPj7RiJTSpmXWWUMy/Z2R/am8aSRHh/K7d0+wclYyEaGDczNqW7po6uj1qMPWnidvW4pFSuVbchNvRunYG4KvBrzfkl2hmMCkJ0TQ0WvG1DVQplpLuvJPhI49IUEGfrYun/21Rr723N5hph1PJlw5YnpsuBL2HsCbNvxfCyEOCCH2A2uAb3lxLIViwqNH4tjH4vsr6coRly1I5afr8nnnyBnufvnAoBo7pTVGwoINPu3epBg7XquWKaX8vLfOrVBMRmyhmcYu8tPi6LV2yvJnhM5QPr88i8a2Hn77znGSokO55zKtmX1JjZH8tDhCxlCLXuF7VHlkhSJAGNoIRW98EkgCH+CbF8+hqaOHxz8oJyU6jNtWZlNWb+KmZVn+nppiFJTAVygChMSoUMJDDLZIHX8nXTlDCMG9V+XT3NHLfa8f5nRrN919FgoyvWO/V3gOJfAVigBBCDEoFt8fdfBdJcggeOj6Aoydu/jTR1qp4IJ0JfADHWVwUygCiLSESJtmrwv8VA9UnvQGYcFBPP75JeSnxTI9NmxCVY2cqigNX6EIINLiIyir02rS1Ju6SI4OC+h+BDHhIbzwXysxdvX6PVdAMTpK4CsUAUR6QgTNHb109vZT29JFWgDE4I9GRGgQEaFKu58IKJOOQhFA6A7aemOX3zpdKSYvSuArFAGEfSOUemO3EvgKj6IEvkIRQOga/sH6Vrr6zAEXkqmY2CiBr1AEENNjwwk2CHZWNAOBGZKpmLgoga9QBBBBBsGMuHD2VLUAgZd0pZjYKIGvUAQYafERtPf0AwREpUzF5EEJfIUiwNAdt+EhBhL93PhEMblQAl+hCDDSrWacmQHQ+EQxuVACX6EIMHQNX9nvFZ5GCXyFIsBIi48EYKbq8KTwMErgKxQBhq7hq5BMhadxS+ALIa4VQhwUQliEEEVDtt0jhDghhDgqhLjUvWkqFFOHzMRIblyWyaX50/09FcUkw93iaWXAeuBx+xeFEOcCNwDzgZnAFiHEOVJKs5vjKRSTniCD4BfrF/h7GopJiFsavpTysJTyqINN64C/Syl7pJQVwAlgmTtjKRQKhcI9vGXDTwNq7P6vtb42DCHE7UKI3UKI3WfPnvXSdBQ2qnfAkTf8PQuFQuEHRjXpCCG2ADMcbPqhlPKf7k5ASrkR2AhQVFQk3T2fYhQ++BU0l0Pe5f6eiUKh8DGjCnwp5cXjOG8dkGH3f7r1NYW/MdVAR6O/Z6FQKPyAt0w6/wJuEEKECSFygDnATi+NpXAVKcFUC71t0Nft79koFAof425Y5tVCiFpgBfC6EOJtACnlQWAzcAh4C/iqitAJADqboa/T+rfS8hWKqYZbYZlSyleAV5xsux+4353zKzyMqXrg746zEJfuv7koFAqfozJtpxKm2oG/O5r8Nw+FQuEXlMCfShjtImU7VAisQjHVUAJ/KmGqAYPViqds+ArFlEMJ/KmEqQYSZ0FQqNLwFYopiLu1dBQTCWMNxGdAT5uy4SsUUxCl4U8lTLVaZE5UstLwFYopiBL4U4XeTs1uH5ehCXxlw1cophxK4E8V9JDM+EyISlEavkIxBVECf6pgsoZkxqVDZLKy4SsUUxAl8KcKNoFvNen0dUBvh3/npFAofIoS+FMFYw2IIIhJ1QQ+qKqZCsUUQwn8qYKpFmJnQlCwZsMH5bhVKKYYSuBPFUw1mjkHBgS+0vAViimFEvhTBT3pCiAySfutBP7k4eA/YOvD/p6FIsBRAn8qYDFDa91AOWSbhq9CMycNWx+C938B/T3+nokigFECfyrQ1gDSPGDSCY2C4HBlw58sdJvg1H7o74b6ff6ejSKAcbfj1bVCiINCCIsQosju9WwhRJcQosT685j7U1WMG1vSlVXgC2FNvlIC36Ps/BN8+IDvx63eAdKi/V31se/HV0wY3NXwy4D1wIcOtp2UUhZYf+5wcxyFOxjtYvB1IpOUwPc0O/8Ee572/bhVW8EQAgnZUP2J78dXDMZihkP/AovF3zMZhlsCX0p5WEp51FOTUXgJvbWhfUtDVV7Bs/R2QOMxaK2F/l7fjl25DdKWQO6FUL1dEzgK/3HyXdj8eah0pAf7F2/a8HOEEPuEEB8IIVY720kIcbsQYrcQYvfZs0oAeQVTrabRh0YNvBaVDJ2qvILHOFUGSM20YqoZdXeP0dOu2e2zV0HWKuhphdMHfTe+YjjN5drvs4GnC48q8IUQW4QQZQ5+1o1wWAOQKaUsBP4beE4IEetoRynlRillkZSyKCUlZXxXMVWxmOGdn0FL1cj7GWuGNyzXSyRL6b35TSUaSgf+bqnw3bg1OzSHfNYqyFyhvabs+P7FaF1RNx737zwcMKrAl1JeLKXMd/DzzxGO6ZFSNln/3gOcBM7x3LTHwI6NcPaYX4b2Oqf2w0cPwJ4/j7yffdKVTmSyFtWh6ul4hoZSCArT/m6p9N24Vdu0khkZxZpTPi4TqpXA9ytGqwLWdMK/83CAV0w6QogUIUSQ9e9cYA5Q7o2xRqS1Ht78LnzyiM+H9gm1u7XfVSM46qS0Nj4ZIvBVLL5naSiFrJVauGuzDzX8ym0wsxDCorX/s1ZoGr5aufmPlkkq8IUQVwshaoEVwOtCiLetm84H9gshSoAXgTuklM3uTXUc6Evb6h0+H9on1O7Sftfvhb4ux/t0tUBv+0BIpo6tno6y47tNfw+cPawJ3vgs32n4vZ1Qt0ez3+tkrdQe4k0nfTMHxXCM1YDQVta9nf6ezSDcjdJ5RUqZLqUMk1JOl1Jean39JSnlfGtI5mIp5aueme4Yqdqm/W48Cp1uPm+k1BxkgUTtLgiPA3Ov9sV3hMlBSCZAlF5eQWn4bnPmEFj6IXURJOb4TuDX7gJLn2a/18lcqf3W732Fb+k2QbcRZhZo/zcH1oN3cmfaVoj2W8AAACAASURBVH08UDemZqd75zr2Fvw6x7lg9TWdzVo0wJIvaP87M+voSVfDnLbKpOMxdIdt6iItFr6l0jcmlaptIAyQuXzgteQ5mn9GxeP7Bz3nZdZF2u8Ac9xOXoHf0Qhnj0DRl8AQDDXb3TtfxUeaJv3atwIjzlm338+5BKad69xRp9+A8ZmDX49UNfE9RkMphMVpwj4hRzOh+cJUVrkNZizQVnk6Qljt+ErD9wu6w3bWGu13gNnxJ6/A1zWc2RfDjIXu2/EbSiA0Wvty73rC/fm5S+0uTbubWajZbWt2grl/+H6mGgiOGFjp6IRGQkhUYAj8nX+Ch/Khy+jvmYyPhlJIXagJ24Rs7TVvO277e7R7IOu84dsyV2p2ZFOdd+egGI4ekpmSp5lRlYbvI6o+1iImZhZqS976vePPgLRYoGE/LLoRctfAu/dB2ynPznes1O6C6fO1ZKrMFZpWeWr/8P1M1hh8IYZvi0ryfwG1I2/AG9/V5jkRtVJzn5Z0lbpI+18X+N6249ftAXPPYIetTpbVjq/MOr7HWA0hkZqClTQbmpTA9w1V2yB9KQSHajHK/d2Dk2PGQvNJ6G3THDFX/EbTrt7+gWfnOxYsFu0Ln75U+3+kL7h9Hfyh+Lu8Qt1eeOlL2vsaFKaZKCYajcc0wWsT+Fnab28L/MptgBhItrJnxgIIjZmYD9CJjrFai9QSQvOnNJ4IqBDZySnwu01w6sBA9ILu1BqvHb++RPudugiSZsHq/4ayl7SaGf6g8ZiWQq8L/NiZmmbpKMPSUQy+TmSy/0w6xmp47not4/emzZCxTCsCNtGwd9gChERofYO9nW1btVVb4UUmDt9mCILM4pHzMxTewVg14C9LmqMpiu2n/TsnOyanwK/ZqdU00TXfmBnaU7d6nAK/oUTTQFPytP9XfRMSZ8Hr34a+bs/MeSzUWR22usAHzW5b/clgbaKvGzrOOBf4/iqR3GWEZ6/TVko3vQDR07SH86kDE8+O37BfW8InzR54LcHLoZnmPu0ez3JgztHJXKHlBrgbjqwYGy3VAwI/2XpPBJAdf3IK/KptWmTOIIG43Fp3ZBzLq4ZSmJEPQSHa/yHhmmmnuRy2+aGtnB5/nzhr4LWsFVpkSKNdGYmhdfCHotvwfbnk7O+Fzbdqts3rn4Fp1odo9irtIT3eh7K/aCjVTCiGoIHXErK967St3wd9nY7t9zr6w0DZ8X1HlxF6TIM1fAgoO/4kFfgfw8zFWiSKTkaxZq9uHmOFB4vFGoVRMPj1WWsg/xr46EHfZzXW7oa0IjDYfXz6F9zerOMs6UonKkULNe1p9c48hyKlFtZa8QFc9QjkXjCwLX0pBIVOLLOOxaI5ynVzjk5CNrTVe2/1V2l9j0bS8NMWa6tSVUjNd+gROrofJzZNi5BrDJzQzMkn8Hs7NWegbs7Rsdnxxxie2VKhCcSZBcO3XfpzCA6DN77jOy25p03L7LRfvQAk5kLUNCcCf0jSlY4t+cpHZp0dj0PJX+GC70PBTYO3hURoNd0nkuO2uVyLjhoq8BNztN+6APA0Vdsgea7m/3BGcJj2fgaiwO/rhvd+Du2TLOlP/7x1Dd9gCLhInckn8Ot2D083B0iZpyXHjNVkoPcIHarhg+Yb+NT/as7bI6+Pb75jpX6fZvoYKvD1hBv7JbyxRovVj53p+Fy+TL4y92mNtnPOhwvvcbxP9nnaaqqnzfvz8QQNds58e2yhmV4w65j7tXt4JHOOTtZK6/sZYCVBSp6FD34FOx7190w8i550FZ818FrSLGXD9ypVH6OFqxUPft1ggIylYy+x0FCimRqmzXO8femXNO36owd8o+XrBdPSFg/flrlS0+p1TcNUCzEzB3wPQ9E1RF/E4h95DdpPwfKvOs4JAO0hLc0Tp9hdQ6l2b+jOfB1vxuKfKtVWFSOZc3SyVmjvZ62bZUU8icUC262CvuS5wMha9xTGai05MyJh4LXkOdqDoL/Hf/OyYxIKfAfp5joZy7XIha4W18/XUKqFvzkTmoYgWPUNTfMuf39cUx4TtXu0ZaKjcDzdjKWH45kcND6xRxf4vojF3/mEttSdc4nzfTKWac72yo+8Px8dYw188gfoHocfo6FUK2sx9N6IStGymL3huNVNXtkOMmyHklGsrfACyaxz4j+aiWPeVdDW4L/QZm9gH4OvkzRHW5H7smT2CEwugd/fCzW7nGs/utZfs8u180np2GE7lEU3QvQM2Pqg63MdD1JqGv5Qc47O9PkQFjtQV8dY7TxCB3xn0jl9SHPGFn1pcDTLUEKjNGe7rxKGDrwIj67SkugeO29sqz8pHTtsYaDEgjc0/KptWnRWzIzR9w2L0cqKBFI8/id/0FadVz+mZaPue8bfM/IcxurhNav00MwAseNPLoHfUAL9XcMdtjppS7TuQK4mYLVUaElcjhy29gSHwYqvQsWHmgbuLYzVWlx9epHj7QZr56OqT7Slcmv9yBp+SLiWkeltgb/rCS1ipPDzo++bvUpbLY23E1f9Pji+xXFdIZ1uE7x8u5blmzIXrn0akPDUWnj/VyMfq2Oq0VaKjgQ+uF4meaxL/YZSbSXkKlmrNL9WIMTjnzqgRWgV36493Bder5XW6JgEPRmkdCzw9dDMALHjTy6Br2uGzgR+aJRW5MpVG7Etw3YUgQ9Q9AUIj/eulq/b751p+KBde+NRrZG1pc95SKaO3tvWW3S3wv7nIX/9QA3+kcg+T6stP9ZoKtAecE+vg2evgYfOhbd/OLyhd9Un8Oh5mnZ/4T3whTdh/mfhjq1amO37P4dNV4zeJ9iWYevk3nClTHL1dvh5muthvf29Wg0ne6fgaBTcqD38/3WX/1P8tz+qJakt2aD9X3iLdo8e2OzXaXmErhYtmm+owA+PhejpAVM1092OV/8nhDgihNgvhHhFCBFvt+0eIcQJIcRRIcSl7k/VBao+1hxoI4WrZSy3Fp7qG/18NoftuaPvGxYDy27XnJPe6lZfu1uL65023/k++sPuwAva76E34FCikr3rtN3/vOZkXPoV1/bPKNZWYWMNz5QSXv2mlldw1SPaQ3HHY/DoSnhstSZs3r0PNl2uOfC/+BZceDcEBWvHh8fBNX+CqzdqD4nHrA8FZzSUavOc7uTeSMjWVpsjpdUf+pcm8E4dcO0aW+sAObKZbigzFsDFP9buy91PuX6cp2k7rd2TBTcPODWnz9eKG+77q/8fRu4yNAbfnqQ5k0bD/w+QL6VcCBwD7gEQQpwL3ADMB9YCf9R73HoNi1nTmJxp9zqZxdoXscFBZcmh1Jdowj441LU5FN+hCeStXsq+rdttTagJdr7PzELNfKILq5FMOuDd8gpSaqWPZxZC+hLXjgmL0UxolWNMwNr/PBx/Gy76ESy+lf/f3pmHR1Vfffx7shBCwASSgAqBsIQCZYsigiyyiAKylYpFBNFaEatW61K12vpYy1Ot79MWitoHW4u2KCirIloWsVKUCmGXxQIFISChYlhkDZz3j3NvcjPMzL137p3MMHM+z5MnM3e2353cnHvuWb4Ho2cAD28HBj4vMfUPHgc+fkHyLRP/FTos0ukHwMQV4jjMuRNY+kxwY3RggzwnPTP4+9Q3avHDhXXMhKXZL2FHqGE2dnS7F2jZT3IVZVvdvdYvVv9ZnKxu91TfXjwWOLg5cmHDSPhqs1QL+UlgDb6VvPipxfc64nAxM5sBz1UAzCNxOICZzHyamf8LYAcAF4HHCDi4WS6p7MrVChwKqVUmbEPEaIORlQtcOV4uUcsd/hM7peK0rKexjeFMyxDv9th+uW8X0qmTGz2Dv3uFhJecevcmzXrIVZjTeaDHvgLef0yuDq6+u2p7Vh7QbSJw98fAPZ8Cdy4FRrwkJ5VwNGguoZ4rxkuIbunTFxp9u2PDThf/SKlUjAHOj5VKg+/CwwfkimbEn2S/Z98Zev5xtDh7EljzF+A7g6Qu3Ur7m0TGfN3fo7+O8+clzPenHsDa6f6+dziDn1skIZ84yFX4GcP/IYD3jduNAViP4n3GtgsgoglEtIaI1hw65CGWbIYAgsnFWrnkMvmj2DVgle+pPpvSKd3vk9+f/NHd6+z4apOEK8LF702aGd9BZn0go27452blR09P57NXZA3tR7p7XWFPCXXsc1BNxQwsfEjkr4e/GLoKqFE76cNwSmoaMOQPUlm0cjKw5BdV39GxryRUc1nH0K/PaQqAQnv4u5bL7/Qs9x5+qEa6cNRrBIx4GSj7HFj8C/ev98LGWaLz1P3eCx/LzAHaDhUnKZpChBVngLl3AZ9OlVDczuX+vn/5HqmQq51z4WN58aOpY2vwiWgpEW0O8jPc8pwnAVQAmOF2Acw8jZm7MHOX/Px8ty+vYs9K8aqyg55XqlPgQEjNTcLWSk6BVB+sfd1fz9lJwtbEPOk5ufTPypck6SmfVSqPlEr3cfHY0GGPUDTtZtSPO4jjb54DbH8P6Ptk1T+WX6SkiEjeVXfJCfwfTxpXfkY4MJyHn1ZLvv9Q3bY7lkkyr7CHCw9/r1Hj7/L7NCkaIOGd1a9IdUxNwAx8+pKUh4a6+i4eK5VT2xZGZw2njgJvjAI2zwb6Pw10GCX5Pj+dnGA1+Ca58aOaaWvwmfk6Zm4f5GcBABDR7QCGALiVufIbLAVgve5sYmyLDszyB3TSfQhIHP/4wfDx1QPrgZR0SSy5pceD4nGu8rF1fN9q4JImcoViR0FX8WKybRK2gKX5yuewTsl0aTjpcqf719bOFgNhF8c/XiY6Rk2uCu49+gERMPgFyc+sehH44IkqSYVLO4R/baha/PPnxMNv2U+uBI441Nw5std9OCeQ656W73bBj6WqKdrsWCZhve73he6wLuwtx2o0wjrHDkrV1X9XAMNfklkWhT3kqtbP4opgJZkmOc3EllwMHn44iGgggJ8BGMbM1oDrOwBGE1EGETUHUAQgev3dh7YDJw/bJ2xNChwIqe1fL3IKaRnu15PfGmg7RDypSDo4g7Fvdej6+0Ay6gE9fwp0Gm3/3GgY/IozYvCLBlQJibmlsKdUJYW6zGcG3ntI4vzhQjl+QAQMfE6843+/LEn53Fb2uYBQBv/AeonptuwvBvzUEWfHyZF97hO2gaRlADe9KjmhuROiL23w6VRpSvzu90I/JyVFxPR2feSv4NzXO4G/DJCSyDGzgOJbZXulsqxPyqyhavBNUtNEfiUOVDO9xvCnAqgHYAkRrSeiPwEAM38O4C0AWwB8AOBeZo7ekZWeKR5E897Ont+wrcTbQnmQzPJP6TZ+b6XnT+UfefUrkb+HyfEyOaCchHNM+v8CaDfM/nl1IpRXOFchUsfzfwx89Jzoouz+l6xzy3xpEHObrLXSrIeMDjSHvQTy+Txg67tA3yekeSraEAE3TJLj7Oy34iXbUb9QriQDm8jM6pwWfapKLO3i+Mzhp5e5Ia9Irlp2r5DKnWiVRB7cIlcyV0+wr3TrPAYAA+vf9OezD+8SY3/mODB+YXVJjwYtZCqZ20qwUJw4LJ8TrgQ6ryguPPww9X32MHOrMI9NAjDJy/s7pn4z+Wd0Skoq0OZGMVIdfwA071X98fIvjS5KDwa/8ZVA60Gil9/5Vmet8KHYPFd+u+mwdIopkey2Fn/d61LXndXQOFkEGI36hUCr6yJfV7PuAEiS8VbdmFNHJVT2yRSRYeh+f+Sf4RYi4PpfS+z+8mL751eKqO2pXq+/40N5j7r5VWG38r3hw4cnv5GhJ149fJPOt4pBXvWiGL+eD/rzvla2G3mCK8bbP7d+M6D5taKk2fvR6rMeImHVy6K6es8nF+Z2iOSY2vVPOdmFCjU5xVTJDFaDb5LbCvjiH+IoBSur/vZ/QGYD7/ttQ2J12rph0G+lROzt8Rcmzey6KJ1ywySprFnyy8jf4/AuYNkzEu914+E7pY7R/eompHPqqOiZN70GeOQL4KmDwP1rgXHzgaFT5B92+EveDt7M+jJlzLzsPnNCQimTO0o3bIs+wM2vh+9JiAZEQMebnSWIGwSpxT91VNQrW/aT+049fPNxN01X4TBPXu2/L2WnG2b6875WSteKoQvXCGmleJwYT6+hlrMnpTKo7bDQf6dmPeQq1I8O2HAlmSZ5RVJ5Vh6kg/vsKeBvI4C5P/K+FhuS1+DXvgQY/YY0g8y6tXpt8oH1otoYScLWSm5L4JqfyMEXiYDV+XMSMklJB4ZN9e6JBCOtliRJ3Rj8lZPFq7/h17KmtAzZ15Z9pQ+h31PO9NrtaGYImq16GZjcSQxT4yuBu5ZLY5Vfxi9aVDZfWSp1dq+QqqiW/eV+VkPp5raLXZtOiV8ePmDU578MFPYCFtwL7Fjq33szG42CDvNOgOS9Mi4BNszy9tlbFkg41ZRwCIZ51ehHWMf824ULt4XT1Fn8lJRddxjlfS02JK/BB+SsO3KaePTvPlgVy9y/XgampNf2/hm9HpLqmkWPuk+QrXpJBpoM/q2zctNIqeNCT+fIPknEdRhl3wTmlcIeUu30weNAXmvgjg+AsXOCzwKIRzLry9Adq4e/80OpvS8wlFtTUsSI23r4ETZd2ZGWISfP/LbArNvEK/eDo6WSv3BzjKRnSnPW9vecSZ+EomS6KIqGk5DObWVMiPNBmbX8S3GaMoPU4JuEqsXfskDyfN3vk32PMslt8AH5kvv8HNg4U7RXKhO2Ljpsw1ErSzzhg5vcaZmUbQOWPQu0GSJ5hmiSle/c4C97Vr6j/h7CVE5pdR1wzf0SKrp9YVVD2cUCkcR1rd22O5ZJzsiaxMwusK/FP7JXOlLrOBCgc0vtbGDsbHnvN252P/c5GKWGaqxTSQ2TtsMkXxHpTIRD28VJunJ8+CtiInEodq/0nrQu32MvaFengcTorR7+N7uBBffLSbH/097W4BA1+IDEnNsMkcaaDTOlK9Br/N5KuxFSQfThr521V587C8y7W7pkh/w+OqEcK1l5ss927F8nJ8Zu99iLsvlBeqbEmVv2jf53EC2spZmHd0l4x4zfm+QUOPPws5tE73uodykwbq5chf5tpPd5s6UlEqpq1N7d61r1lyugLe9E9rklr0kItNMY++c26yESJHajKHcul6qwUIQrybSSV1SVM6g4A7x9h9y+6VXnel0eUYMPVMUyc1sC8w1xJydVGE4hkiTx6WPAh7+yf/6K38lVxo2/A+o29G8doXAikcwsLfl1ciVMpTijQXPxAM+fryrHNOP3JtlNJfwRTlrAr5LMcOQVAWPeknDMx7/19l77SqQxzW0fS3om0Pp66bp1GwI9ewrY8CbQZrBUQNlRGccPE9apOC3O19y7gaMHLny8sgbfgWS1VTVz2TPA/rXA8KlV1Vw1gBp8EzOJm1HPkL31mLANpGFb6dYsea1qMHow9q+Xf7YOo0SnvSbIyhcPP5yC4Pb35TK7zxPBx0cqwalfKJVax/ZLOWZO0wsFxMzk89Ewzeh24yr9ouAqCats9KBtc/6cHONuErZW2g4TB8RO7yqQbQulATNcstZKfhtxYMLF8TfOkpPxudPAyiAquCe+lnLZcCWZJnmtpDJowyzJg3Wd4KxXxkfU4FvJK5KOvIHPRa5XEo4+j4lxXfRocONacRqYN1GSqIM8elhuqJMnMgihZv2eOyviYblFzv+ZFMGs1Pnff2QiWst+F4ZlTM89VKVOxWkxOtH28E2Kx4q2UqTaNoe2SXNapEn9ouslX7FlgbvXlUwXT7t5H2fPJ5KwTqhKnfPngZVTpMmu81hgzV8v9PLNQTlOQjpmpc6CH8t7DnjW2Tp9RA1+IM2ukc7AaFA7GxjwjMgkvP+o1OfP+RHw1xuByZ2B55qKZO6wPwYfUh4tzDrpUM1XJdMl9nj9s6GHuSvBMS/XN80Gzhy7MJwD2Nfim55/TXj4gDRAedG2MRO2kRr8jLqSsN/6rnPd+q93yhXoFePc9X8U9pTvPdiEsy/el6qaHg8AvR+RctqVk6s/p9yFwTcrddJqA6Om+1MF6BI1+DVNx9GiZrn6z1JfvvczgM9JqWHXu4Af/F1imDVJVhh5hZPlwEe/kVrt1gNrdl2JQHYTCRFuniO/g8l/XNJY1EFDVepEOvgkUlJSRHcmUm2b0hJxbgJDV25oO0zCYObJw461r8n323msu8+p1NUJEtZZOVkMebsRkovpdAtQ8leRxzZx0nRl0qAF0KKvDHD38t14oIbbFBWkpAC3LZCOyzq5UW+ldoQprxCs+WrRo2L0b5h08VbKxJLUdPHgv9kttffBarVT00XeIJSHbxr8mmw06zymSiOpz+PuXruvRLx7L8dL6xuk2mbrAvs5BhVnZJ2tBzpTk7XSsJ30S+xeaej5GHy5SsQVB71Q1c3d+2FJCq+cDAz8jWwr/1LKLe2E9AD5O9823936fCYOrE0SkpYhVQTxYOyB0AJqm2bLYIprH3M3+UupjhnWCSzHtBKuFt/cfkkUm+8CyWkq8hXrZrgbB3jmW6Bsi/emvMwc+fwtC+zr5LcvkmP3SgeaPYGkpIhESKCcw8opYshNhU1APPROo6WfxvTyy/fUTImyT8SJxVFiitnMY63FL98rk6SaXAX0ejg260oUKg1+kPi9SU5BaF38I3tlWEokUt1eKB4ra9r9sfPXHNhohCgjrNCx0m64eNB2827XviYnw0jF+gp7yBXYESNXcugL6fbtOkEaJ630eliKGMxYvtMa/DhBDb4il6yZ9as8/PPnpB+Bz4n0RE0LlCUarQaILlC43o7sAhlIEqz23A8d/EhoM0Ri8W6St6actR/yF21ulLj81jBNWN/slsao4nGRz0QIjON/MgVIy5ScWiC5LaXz3fTy1eArFyVZ+VUx/E+nSsXDoOflMlbxRtshwB3vhT9x5hRIFcixIM09sTL46bWBDjdL12uokt1ASkukwsePhsE6DaSKJlRYp+KMCI8BcjUSKZd2EM2j3f+SssuNxrCUUCqfvR8RL3/xU6L1VIONU15Rg68IdfLE4B/YKHo5bYeKZrpSM1h18a34OfgkEq4YJ01Hm2Y7e/6+Evf6OeFoN0xKgsu2Vt9+5gQwc4yUbg74lbeEdkqq6DTtWSl6Wucrwo/MzG0pEtmb3pb7yeLhE9ELRLSNiDYS0TwiyjG2FxLRSWMKVuUkLCWOycoDju4D5t4lMf2hU7QqpyYJVYt/4mug4mTsDP5lncQDdhLWOV4mMX8/VVTbDAVA1cM6J8uBv48UOeehk4EeP/H+Oc16yInls1ckd2B3Zdv7USmlBZLH4ANYAqA9M3cE8AWAJyyP7WTmzsbPRI+fo0SbrHyJhx7aBox4qWYbv5SqkE1g3bt5AohFSMekeJxoO321KfzzTGllPw1+vUbSt2KKqR0vA6YPkXnHN73qX+e3Ob/h7Lcyw8IOM5afkha7k3EEeDL4zLyYmSuMu6sAxPCoVDxhxiuvvkcUC5WapVaWXFkFevg13XQVjA6jRPly3YzwzytdI0lWv0t42w0Dyj4X8blXbwAO7wTGzATaj/TvMy7tJHH8wl7OE86DXwBuXySdwRcJfsbwfwjgfcv95kS0joj+SUS9Qr2IiCYQ0RoiWnPokEdJViVyvjNIPLnrakaXWwlCsFr8yqarGIYN6jSQip2NM0XXJxSlJdLIFFjK6JW2Q+X330ZKiGvcfG/zkoORmgaMXwCMfMX5azLqAU2v9ncdUcbW4BPRUiLaHORnuOU5TwKoAGC6AAcANGXmYgAPAXiDiC4J9v7MPI2ZuzBzl/x8B5KmSnS4vFikWqMhGqc4I5gufvleIL2OlM3GkuKxUqljDiYPhFkMvp8JW5PsJlLWmpUvHnW0jOzlxe47dS8ybAusmTnsqZSIbgcwBEB/ZqmdYubTAE4bt0uIaCeA1gDWeF2woiQs2U2B/ywV42kmzE1Z5Fgn0Fv0kSuQj56T24EnoK93yhzZaI29vOUNSZI6kTBQQuK1SmcggJ8BGMbMJyzb84ko1bjdAkARAB/mpilKApNTIBU51o7nWNXgB5KSKleAh3cBb94CnD1Z/XGvCpl21M5WY+8DXmP4UwHUA7AkoPyyN4CNRLQewGwAE5n5sMfPUpTEJpgufrwYfEA8+5HTRFjs7TuAcxVVj5WukdGE+W1itTrFAZ565pm5VYjtcwDM8fLeipJ0WGvxG18hE6e+LatqyooHvvs9adBb9Ajw7gPi9ROJh395ceTyBkqNoCIpihIvVHr4RuK2pgefOKXrXaK79M/npZy378+lRv9qbbeJd9TgK0q8kFkfqFW3qlInHpquQtHnCTH6K/8gImLnzgBNfFDIVKKKGnxFiReIqtfix0PTVSiIgMH/JwnmjTNlW7QStopvqMFXlHjCqot/ZB8AqtnBJ25ISZVGpZPlkmiO13UqlajBV5R4IrtA5hwDEtKpdymQViu2awpHWgYwbh5w9kTsewUUW1QeWVHiiZwC4FQ5cPqYhHbiMZwTSEqq1shfJKjBV5R4wlqpE081+EpCoAZfUeIJUySt/Es1+IrvqMFXlHjC9PD3r5NJUxeR1roS/6jBV5R4om4j0Z7/8hO5rwZf8RE1+IoST6SkSHnjPkNYVkM6io+owVeUeCOnQMocATX4iq+owVeUeMMUS0vPiv3gEyWhUIOvKPGGqZqZU6DNTIqvqMFXlHjDTNRqOEfxGTX4ihJv5KjBV6KDZ4NPRM8S0UZj4tViIrrc2E5ENIWIdhiPX+F9uYqSBKiHr0QJPzz8F5i5IzN3BrAQwC+N7YMgs2yLAEwA8LIPn6UoiU/9QuDax4D2N8V6JUqC4Vktk5mPWu5mAWDj9nAArzMzA1hFRDlEdBkzH/D6mYqS0BDJFClF8Rlf5JGJaBKA2wAcAdDX2NwYwF7L0/YZ29TgK4qixABHIR0iWkpEm4P8DAcAZn6SmQsAzABwn5sFENEEIlpDRGsOHTrkfg8URVEURzjy8Jn5OofvNwPAIgBPAygFYBUCaWJsnxbMUAAABCVJREFUC3zvaQCmAUCXLl048HFFURTFH/yo0imy3B0OYJtx+x0AtxnVOt0AHNH4vaIoSuzwI4b/HBF9B8B5AHsATDS2LwIwGMAOACcA3OHDZymKoigR4keVzvdDbGcA93p9f0VRFMUftNNWURQlSVCDryiKkiSQRF7iAyI6BMkDREoegP/5tJyLCd3v5EL3O7lwst/NmDnf7o3iyuB7hYjWMHOXWK+jptH9Ti50v5MLP/dbQzqKoihJghp8RVGUJCHRDP60WC8gRuh+Jxe638mFb/udUDF8RVEUJTSJ5uEriqIoIVCDryiKkiQkhMEnooFEtN0Yp/h4rNcTLYjoVSIqI6LNlm0NiGgJEf3H+F0/lmuMBkRUQETLiWgLEX1ORA8Y2xN634moNhF9RkQbjP1+xtjenIj+bRzvs4ioVqzXGg2IKJWI1hHRQuN+suz3biLaZIyNXWNs8+VYv+gNPhGlAngRMlKxHYBbiKhdbFcVNaYDGBiw7XEAy5i5CMAy436iUQHgYWZuB6AbgHuNv3Gi7/tpAP2YuROAzgAGGsqzzwP4PTO3AvANgDtjuMZo8gCArZb7ybLfANCXmTtb6u99OdYveoMPoCuAHcy8i5nPAJgJkWlOOJj5YwCHAzYPB/Cacfs1ACNqdFE1ADMfYOa1xu1jECPQGAm+7ywcN+6mGz8MoB+A2cb2hNtvACCiJgBuBPBn4z4hCfY7DL4c64lg8EONUkwWGlnmDHwFoFEsFxNtiKgQQDGAfyMJ9t0Ia6wHUAZgCYCdAMqZucJ4SqIe738A8DOI7DoA5CI59huQk/piIiohognGNl+OdV9m2irxATMzESVsnS0R1QUwB8CDzHxUnD4hUfedmc8B6ExEOQDmAWgT4yVFHSIaAqCMmUuIqE+s1xMDejJzKRE1BLCEiLZZH/RyrCeCh+9olGICc5CILgMA43dZjNcTFYgoHWLsZzDzXGNzUuw7ADBzOYDlALoDyCEi01lLxOO9B4BhRLQbEqLtB2AyEn+/AQDMXGr8LoOc5LvCp2M9EQz+agBFRga/FoDRkPGKycI7AMYbt8cDWBDDtUQFI377FwBbmfl3locSet+JKN/w7EFEmQAGQPIXywHcZDwt4fabmZ9g5ibMXAj5f/6QmW9Fgu83ABBRFhHVM28DuB7AZvh0rCdEpy0RDYbE/FIBvMrMk2K8pKhARG8C6AORSz0IGRY/H8BbAJpCpKVvZubAxO5FDRH1BLACwCZUxXR/DonjJ+y+E1FHSIIuFeKcvcXMvyKiFhDPtwGAdQDGMvPp2K00ehghnUeYeUgy7Lexj/OMu2kA3mDmSUSUCx+O9YQw+IqiKIo9iRDSURRFURygBl9RFCVJUIOvKIqSJKjBVxRFSRLU4CuKoiQJavAVRVGSBDX4iqIoScL/A9IL5lDWGBonAAAAAElFTkSuQmCC\n",
"text/plain": [
"