{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hard Uncertain Convex Inequalities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robust optimization is a popular technique for dealing with uncertainty in optimization problems. Methods of reformulating uncertain problems into their robust counterparts exist for a wide range of problems. However, problems with constraints that are convex in both their optimization variables and the uncertain parameters are, in general, harder to solve. Roos et. al. (2018) (http://www.optimization-online.org/DB_HTML/2018/06/6679.html) provide a way to approximate this popular class of problems. Implementing their approximations in Fusion is going to be the goal of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adjustable Robust Optimization problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the uncertain optimization problem: \n",
    "$$\n",
    " \\mathrm{minimize} \\quad c^T x \n",
    "$$\n",
    "\n",
    "$$\n",
    " \\mbox{subject to} \\quad f(A(x) \\zeta + b(x)) \\leq 0 \\quad \\forall \\zeta \\in U\n",
    "$$\n",
    "\n",
    "where $f: \\Bbb R^p \\mapsto \\Bbb R $ is convex and closed, $A:\\Bbb R ^n \\mapsto \\Bbb R^{p\\times L}$, $b: \\Bbb R^n \\mapsto \\Bbb R^p$ are affine, and U is the uncertainty set which is a non-empty polyhedron given by:\n",
    "$$\n",
    "U = \\{\\zeta \\in \\Bbb R^L | D\\zeta = d\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for some $D \\in \\Bbb R^{q\\times L}$ and $d \\in \\Bbb R^q$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above-stated constraint will be satisfied if and only if $x$ satisfies the following adjustable robust optimization constraints:\n",
    "\n",
    "$$\n",
    "\\forall w \\in \\text{dom} f^*, \\exists \\lambda \\in \\Bbb R^q : \\begin{cases}\n",
    "d^T\\lambda + w^Tb(x) - f^*(x)\\leq 0 \\\\\n",
    "D^T\\lambda = A(x)^Tw,\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "where $f^*(w):\\Bbb R^p \\mapsto \\Bbb R$ is the convex conjugate of the original function $f$ in the constraint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note: Since we have $\\zeta \\in \\Bbb R^L$, we get an equality as the second constraint. If $ \\zeta $ has a non-negativity constraint in $U$, i.e. $\\zeta \\in \\Bbb R_+^L$, then the equality is replaced by $D^T \\lambda \\geq A(x)^Tw$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Safe Approximations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roos et. al. provide two safe approximations to the ARO given above:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Approximation 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there exists $u\\in \\Bbb R^q$ and $V\\in \\Bbb R^{q\\times p}$ for a given $x\\in \\Bbb R^n$ such that \n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "d^Tu + f\\big(b(x)+V^Td\\big) \\leq 0 \\\\\n",
    "\\delta^*\\big(A_i(x)-V^TD_i|dom\\,f^*\\big) = D^T_iu \\quad i=1,...,L.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "holds, then $x$ also satisfies the original statement of the problem. Here, $A_i$ stands for the $i$-th row of the matrix $A$. $\\delta^*$ is the support function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Approximation 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there exists $u\\in \\Bbb R^q$, $V\\in \\Bbb R^{q\\times p}$ and $r\\in \\Bbb R^q$ for a given $x\\in \\Bbb R^n$ such that\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "d^Tu + (1+d^Tr)f\\Big( \\frac{V^Td+b(x)}{1+d^Tr}\\Big) \\leq 0 \\\\\n",
    "1+d^Tr \\geq 0 \\\\\n",
    "-D_i^Tu + (-D_i^Tr)f\\Big( \\frac{A_i(x)-V^TD_i}{-D_i^Tr} \\Big) = 0 \\quad i= 1,...,L\\\\\n",
    "-D_i^Tr \\geq 0 \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "holds, then $x$ also satisfies the original problem statement. Moreover, this is a tighter approximation than the first approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roos et. al. test the above-mentioned approach by solving the following problem:\n",
    "\n",
    "$$\n",
    "\\text{Minimize}\\quad c^Tx\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{Subject to}\\quad h\\Bigg(\n",
    "\\begin{pmatrix}\n",
    "\\big( -\\mathbf{1} + B_i^{(1)}\\zeta \\big)^T x \\\\\n",
    "\\big( -\\mathbf{1} + B_i^{(2)}\\zeta \\big)^T x\n",
    "\\end{pmatrix} \\Bigg) \\leq 0  \\quad \\quad \\forall \\zeta \\in U,\\, i=1,..., m\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $c=\\mathbf{1}\\in \\Bbb R^n$ is the all ones vector, and $B_i^{(1)},B_i^{(2)} \\in \\Bbb R^{n\\times L}$ are randomly generated sparse matrices (sparsity density = 0.1), with non-zero elements that lie in the interval $[0,1]$. The uncertainty set is assumed to be a box, that is:\n",
    "\n",
    "$$\n",
    "U = \\big\\{ \\zeta \\in \\Bbb R^L \\,\\big|\\,\\, \\left\\lVert \\zeta \\right\\rVert_\\infty \\leq 1 \\big\\}\n",
    "$$\n",
    "\n",
    "which can be re-stated as:\n",
    "\n",
    "$$\n",
    "U = \\big\\{ \\zeta \\in \\Bbb R^L \\,\\big|\\,\\, D\\zeta \\leq d \\big\\}\n",
    "$$\n",
    "with\n",
    "\n",
    "$$\n",
    "D = \\begin{bmatrix}\n",
    "\\mathbf{I} \\\\\n",
    "\\mathbf{-I}\n",
    "\\end{bmatrix}  \n",
    "\\quad\n",
    "d = \\mathbf{1} \\in \\Bbb R^L\n",
    "$$\n",
    "\n",
    "where $I \\in \\Bbb R^{L\\times L}$ is the identity matrix. Therefore, we see that $q=2L$. Moreover, by making the following substitutions:\n",
    "\n",
    "$$A(x) =\n",
    "\\begin{bmatrix}\n",
    "\\big( B_i^{(1) T} x \\big) ^T \\\\\n",
    "\\big( B_i^{(2) T} x \\big) ^T\n",
    "\\end{bmatrix} \n",
    "\\quad\n",
    "b(x) = \\begin{bmatrix}\n",
    "-\\mathbf{1}^T x \\\\ \n",
    "-\\mathbf{1}^T x\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "we can express the example constraints in the form: $f\\big(A(x)\\zeta + b(x)\\big)$, where the function $f$ is now the log-sum-exp ($h$) function with $p=2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.) Exact solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve the robust optmization problem exactly, by evaluating the constraints at every vertex of the uncertainty box set. Therefore, we end up with $m\\times2^L$ log-sum-exp inequalities. Each inequality can be expressed as: \n",
    "\n",
    "$$\n",
    "\\Big(h^{(j)}_1,h^{(j)}_2,h^{(j)}_3 \\Big) \\in K_{exp} \\quad, \\quad  \\sum_{j} h_1^{(j)} \\leq 1 \\quad \\quad j=1,...,p\n",
    "$$\n",
    "\n",
    "where $h_2^{i} = 1$ and,\n",
    " \n",
    "$$\n",
    "h^{(j)}_3 = -\\mathbf{1}^T x  + \\big(B_i^{(j)T}x\\big)^T \\zeta  \\quad i=1,...,m\n",
    "$$\n",
    "\n",
    "for a given uncertainty vector $\\zeta \\in U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import random\n",
    "from mosek.fusion import *\n",
    "import sys\n",
    "import scipy\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [16, 9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function that generates a random sparse matrix B, with the specified sparsity density, and dimensions n*L.\n",
    "#Please note that we generate transposed matrices.\n",
    "def Sparse_make(n,L,s,ran_stat):\n",
    "    return scipy.sparse.random(L,n,density = s,random_state = ran_stat).todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function to set the dimensional constants, matrix D, and the vector d. \n",
    "def constants(L):\n",
    "    q = 2*L\n",
    "    d = np.ones(q)\n",
    "    D = np.concatenate((np.eye(L),-1*np.eye(L)),axis=0)\n",
    "    Dt = np.transpose(D)\n",
    "    return(q,d,D,Dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Making the vertex vectors of the uncertainty box. This will be used in the exact solution.\n",
    "def Box_vertices(L):\n",
    "    U = []\n",
    "    for i in range(2**L):\n",
    "        vertex = [2*int(b) - 1 for b in list(format(i,'b').zfill(L))]\n",
    "        U.append(vertex)\n",
    "    U_matrix = Matrix.dense(np.transpose(np.asarray(U)))\n",
    "    return(U_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Exact_sol(Bt_p,n,m,L,p):    \n",
    "    U = Box_vertices(L)\n",
    "    q,d,D,dt = constants(L)\n",
    "    \n",
    "    with Model('EXACT') as M:\n",
    "        #M.setLogHandler(sys.stdout)\n",
    "        M.setSolverParam('numThreads', 1)\n",
    "        \n",
    "        x = M.variable('x',n)\n",
    "        #Objective\n",
    "        M.objective('OBJ_exact',ObjectiveSense.Minimize,Expr.sum(x))\n",
    "        \n",
    "        #b(x),dim = [p,2^L,m].\n",
    "        bx = Expr.repeat(Expr.repeat(Expr.repeat(Expr.neg(Expr.sum(x)),p,0),2**L,1),m,2)      \n",
    "        #[A(x).(uncertainty vector)] , dim = [p,2^L,m].\n",
    "        Ax_U = Expr.stack(2,[Expr.mul(Expr.transpose(Expr.hstack([Expr.mul(b,x) \n",
    "                                                       for b in B_set])),U) for B_set in Bt_p])\n",
    "        h1 = M.variable('h1',[p,2**L,m])\n",
    "        h2 = Expr.constTerm([p,2**L,m],1.0)\n",
    "        h3 = Expr.add(Ax_U,bx)\n",
    "        \n",
    "        Cone = M.constraint('Exp1',Expr.stack(3,[h1,h2,h3]),Domain.inPExpCone())\n",
    "        Line = M.constraint('Lin1',Expr.sum(h1,0),Domain.lessThan(1.0))  \n",
    "        \n",
    "        M.solve()\n",
    "        return(M.primalObjValue(),M.getSolverDoubleInfo('optimizerTime'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.) Approximation 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to model the first approximation by using the previously mentioned expressions for $A(x)$, $b(x)$, $D$ and $d$. We see that for the first constraint in the approximation, i.e,\n",
    "\n",
    "$$\n",
    "d^Tu + h(V^Td + b(x)) \\leq 0\n",
    "$$\n",
    "\n",
    "we can write,\n",
    "\n",
    "$$\n",
    "\\sum_{k=1}^{p}e^{\\big( V_k^Td + b(x)_k + d^Tu \\big)} \\leq 1\n",
    "$$\n",
    "\n",
    "since $h$ is a log-sum-exp function. Note that k is the column number in $V$ (two dimensional variable) and the row number in $b(x)$. The above expression can be easily expressed in Fusion, using exponential cones,\n",
    "\n",
    "$$\n",
    "\\big(h_1^{(k)},h_2^{(k)},h_1^{(k)}\\big) \\in K_{exp} \\quad \\quad \\sum_{k=1}^{p} h_1^{(k)} \\leq 1 \n",
    "$$\n",
    "\n",
    "where $h_2^{(k)} = 1 $ and,\n",
    "\n",
    "$$\n",
    "h_3^{(k)} = V_k^Td + b(x)_k + d^Tu \\quad \\quad k=1,...,p\n",
    "$$\n",
    "\n",
    "The second constraint in the approximation, i.e. the equality, involves the support function of the domain of $f^*$ (i.e. the conjugate). For the case of log-sum-exp $(h)$, we find that the domain is a simplex and the support function for a simplex is the maximum function. Hence, the second constraint is easily expressed as a set of linear equations:\n",
    "\n",
    "$$\n",
    "\\max_{k} \\big\\{ A_{ki}(x)-V^T_k D_i \\big\\} = D_i^T u  \\quad \\quad i = 1,...,L\n",
    "$$\n",
    "\n",
    "where $k = 1,...,p$ and $A_{ki}$ denotes the element in $A(x)$ at the $k$-th row and $i$-th column. $V_k$ denotes the $k$-th column in $V$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_1(Bt_p,n,m,L,p):\n",
    "    q,d,D,Dt =constants(L)\n",
    "    \n",
    "    with Model('APPROX_1') as M:\n",
    "        #M.setLogHandler(sys.stdout)\n",
    "        M.setSolverParam('numThreads', 1)\n",
    "        \n",
    "        x = M.variable('x',n)\n",
    "        #Objective\n",
    "        M.objective('OBJ_approx1',ObjectiveSense.Minimize,Expr.sum(x))\n",
    "        \n",
    "        u = M.variable('u',[q,m],Domain.greaterThan(0.0))\n",
    "        Vt = M.variable('Vt',[p,q,m],Domain.greaterThan(0.0))             #Note that we define the transpose of V.\n",
    "        \n",
    "        #b(x), dim = [p,m].\n",
    "        bx = Expr.repeat(Expr.repeat(Expr.neg(Expr.sum(x)),p,0),m,1)\n",
    "        #A(x), dim = [p,L,m].\n",
    "        Ax = Expr.stack(2,[Expr.transpose(Expr.hstack([Expr.mul(bT[i],x) for i in range(p)])) for bT in Bt_p]) \n",
    "        \n",
    "        h1 = M.variable('h1',[p,m],Domain.greaterThan(0.0))\n",
    "        h2 = Expr.constTerm([p,m],1.0)\n",
    "        #Matrix mul(V^T,d), dim = [p,m].\n",
    "        Vt_d = Expr.hstack([Expr.mul(Vt.slice([0,0,i],[p,q,i+1]).reshape([p,q]),d) \n",
    "                            for i in range(m)])\n",
    "        #(d^T . u), dim = [p,m].\n",
    "        dt_u = Expr.transpose(Expr.repeat(Expr.mul(d,u),p,1))                                        \n",
    "        h3 = Expr.add([Vt_d,bx,dt_u])\n",
    "        \n",
    "        #Exponential Cone\n",
    "        Cone = M.constraint('ExpCone',Expr.stack(2,[h1,h2,h3]),Domain.inPExpCone())           \n",
    "        #Linear constraint on the vector h1.\n",
    "        LinIneq = M.constraint('LinIneq',Expr.sum(h1,0),Domain.lessThan(1.0))\n",
    "        #Matrix mult(V^T,D), dim = [p,L,m].\n",
    "        Vt_D = Expr.stack(2,[Expr.mul(Vt.slice([0,0,i],[p,q,i+1]).reshape([p,q]),\n",
    "                                      Matrix.sparse(D)) for i in range(m)])\n",
    "        #Matrix mult(D^T,u), dim = [p,L,m].\n",
    "        Dt_u = Expr.repeat(Expr.reshape(Expr.mul(Dt,u),[1,L,m]),p,0)\n",
    "        #Second constraint, i.e. the equality. \n",
    "        Eq = M.constraint('Eq',Expr.sub(Expr.sub(Ax,Vt_D),Dt_u),Domain.equalsTo(0.0)) \n",
    "        \n",
    "        M.solve()\n",
    "        return(M.primalObjValue(),M.getSolverDoubleInfo('optimizerTime'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.) Approximation 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the second approximation, we can manipulate the first constraint into the following form,\n",
    "$$\n",
    "\\sum_{k=1}^{p}e^{\\Big( \\frac{V_k^Td + b(x)_k + d^Tu}{1+d^Tr} \\Big)} \\leq 1\n",
    "$$\n",
    "\n",
    "which then gives us the following exponential cones,\n",
    "\n",
    "$$\n",
    "\\big(h_{11}^{(k)},h_{12}^{(k)}, h_{13}^{(k)}\\big) \\in K_{exp} \\quad \\quad \\sum_{k=1}^{p} h_{11}^{(k)} - (1+d^Tr) \\leq 0 \n",
    "$$\n",
    "\n",
    "where $h_{12}^{(k)} = \\big(1+d^Tr\\big)$ and $h_{13}^{(k)} = \\big(V_k^Td + b(x)_k + d^Tu\\big)$.\n",
    "Following the same approach for the third inequality gives:\n",
    "\n",
    "$$\n",
    "\\bigg(h_{21}^{(k)},h_{22}^{(k)},h_{23}^{(k)} \\bigg) \\in K_{exp} \\quad \\quad \\sum_{k=1}^{p} h_{21}^{(k)} + D_i^Tr \\leq 0  \\quad \\quad i = 1,...,L\n",
    "$$\n",
    "\n",
    "where $h_{22}^{(k)} = -D_i^Tr$ and $h_{23}^{(k)} = \\big( A_i(x) - V^TD_i -D_i^Tu \\big)$.\n",
    "\n",
    "Note that $h_{12}^{(k)}$ and $h_{22}^{(k)}$ are included in the cone domain and therefore the second and the fourth inequalities in the statement for second approximation are taken care of.\n",
    "\n",
    "##### Mosek highlight: Roos et. al. mention that interior point methods can have trouble solving the second approximation (very small denominator case). However, Mosek successfully solves every single instance of the problem without any issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_2(Bt_p,n,m,L,p):\n",
    "    q,d,D,Dt =constants(L)\n",
    "    \n",
    "    with Model('APPROX_2') as M:\n",
    "        #M.setLogHandler(sys.stdout)\n",
    "        M.setSolverParam('numThreads', 1)\n",
    "        \n",
    "        x = M.variable('x',n)\n",
    "        #Objective\n",
    "        M.objective('OBJ_approx2',ObjectiveSense.Minimize,Expr.sum(x))\n",
    "        \n",
    "        #b(x), dim = [p,m].\n",
    "        bx = Expr.repeat(Expr.repeat(Expr.neg(Expr.sum(x)),p,0),m,1) \n",
    "        #A(x), dim = [p,L,m].\n",
    "        Ax = Expr.stack(2,[Expr.transpose(Expr.hstack([Expr.mul(bT[i],x) for i in range(p)])) for bT in Bt_p])  \n",
    "        \n",
    "        u = M.variable('u',[q,m],Domain.greaterThan(0.0))\n",
    "        r = M.variable('r',[q,m],Domain.greaterThan(0.0))\n",
    "        Vt = M.variable('Vt',[p,q,m],Domain.greaterThan(0.0))\n",
    "        \n",
    "        h_11 = M.variable('h_11',[p,m],Domain.greaterThan(0.0))\n",
    "        #1 + (d^T.r), dim = [m].\n",
    "        dt_r_1 = Expr.add(Expr.constTerm([m],1.0),Expr.mul(d,r))                             \n",
    "        h_12 = Expr.transpose(Expr.repeat(dt_r_1,p,1))\n",
    "        #Matrix mult(V^T,d), dim = [p,m].\n",
    "        Vt_d = Expr.hstack([Expr.mul(Vt.slice([0,0,i],[p,q,i+1]).reshape([p,q]),d) for i in range(m)])\n",
    "        #(d^T,u), dim = [p,m].\n",
    "        dt_u = Expr.transpose(Expr.repeat(Expr.mul(d,u),p,1))\n",
    "        h_13 = Expr.add([Vt_d,dt_u,bx])\n",
    "        #First set of cones.\n",
    "        Cone1 = M.constraint('Cone1',Expr.stack(2,[h_11,h_12,h_13]),Domain.inPExpCone())\n",
    "        LinIneq1 = M.constraint('LinIneq1',Expr.sub(Expr.sum(h_11,0),dt_r_1),Domain.lessThan(0.0))\n",
    "        \n",
    "        \n",
    "        h_21 = M.variable('h_21',[p,L,m],Domain.greaterThan(0.0))\n",
    "        #Matrix mult(D_T,r), dim = [L,m].\n",
    "        Dt_r = Expr.mul(Dt,r)                                                                 \n",
    "        h_22 = Expr.repeat(Expr.reshape(Expr.neg(Dt_r),[1,L,m]),p,0)\n",
    "        #Matrix mult(V^T,D), stacked in m, dim = [p,L,m].\n",
    "        Vt_D = Expr.stack(2,[Expr.mul(Vt.slice([0,0,i],[p,q,i+1]).reshape([p,q]),\n",
    "                                      Matrix.sparse(D)) for i in range(m)])\n",
    "        h_23 = Expr.sub(Expr.sub(Ax,Vt_D),Expr.repeat(Expr.reshape(Expr.mul(Dt,u),[1,L,m]),p,0))\n",
    "        #Second set of cones. Note the stacking in 4th dimension.\n",
    "        Cone2 = M.constraint('Cone2',Expr.stack(3,[h_21,h_22,h_23]),Domain.inPExpCone())\n",
    "        LinIneq2 = M.constraint('LinIneq2',Expr.sub(Expr.sum(h_21,0),Expr.neg(Dt_r)),Domain.equalsTo(0.0))\n",
    "\n",
    "        M.solve()\n",
    "        return(M.primalObjValue(),M.getSolverDoubleInfo('optimizerTime'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximation error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate the quality of the above approximations, Roos et. al. have used the following approximation error (percentage):\n",
    "\n",
    "$$\n",
    "100\\bigg( \\frac{e^{c^T \\hat{x}}}{e^{c^T x^*}} -1\\bigg)\n",
    "$$\n",
    "\n",
    "We compare the objective values of the approximation to that of the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Making 20 instances of the problem with n=m=100 and L ranging from 1 to 20. \n",
    "#Solve some initial ones exactly and others approximately\n",
    "t, tExact, n, m = 20, 7, 100, 100\n",
    "n_list = [n]*t\n",
    "m_list = [m]*t\n",
    "L_list = np.arange(1,t+1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will generate random matrices, with sparsity density $s = 0.1$. For every inequality constraint in the problem,  there are $p = 2$ such matrices. As stated before, $m=100$ inequalities are involved. The dimension of the uncertain parameter, i.e. $L$ ($\\zeta \\in \\Bbb R^L$) ranges in $[1,20]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(n_l,m_l,L_l,p,s):\n",
    "    Ex = []\n",
    "    App1 = []\n",
    "    App2 = []\n",
    "    for j in range(len(n_l)):\n",
    "        print('\\n Iteration: {}'.format(j))\n",
    "\n",
    "        Bt = [[Matrix.sparse(Sparse_make(n_l[j],L_l[j],s,i)),\n",
    "               Matrix.sparse(Sparse_make(n_l[j],L_l[j],s,i+m_l[j]))] for i in range(m_l[j])]\n",
    "        \n",
    "        # Approximation 1\n",
    "        try:\n",
    "            print('   Solve approximation 1')\n",
    "            App1.append(approx_1(Bt,n_l[j],m_l[j],L_l[j],p))\n",
    "        except SolutionError as e:\n",
    "            print('   ' + e.toString())\n",
    "        \n",
    "        # Approximation 2\n",
    "        try:\n",
    "            print('   Solve approximation 2')\n",
    "            App2.append(approx_2(Bt,n_l[j],m_l[j],L_l[j],p))\n",
    "        except SolutionError as e:\n",
    "            print('   ' + e.toString())\n",
    "        \n",
    "        # We only solve the exact problem for reasonably small j\n",
    "        if j<tExact:\n",
    "            try:\n",
    "                print('   Solve exact version')\n",
    "                Ex.append(Exact_sol(Bt,n_l[j],m_l[j],L_l[j],p))\n",
    "            except SolutionError as e:\n",
    "                print('   ' + e.toString())\n",
    "                \n",
    "    return(np.asarray(Ex),np.asarray(App1),np.asarray(App2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Iteration: 0\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 1\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 2\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 3\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 4\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 5\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 6\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "   Solve exact version\n",
      "\n",
      " Iteration: 7\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 8\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 9\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 10\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 11\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 12\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 13\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 14\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 15\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 16\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 17\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 18\n",
      "   Solve approximation 1\n",
      "   Solve approximation 2\n",
      "\n",
      " Iteration: 19\n",
      "   Solve approximation 1\n",
      "   Solution status is Unknown but at least Optimal is expected\n",
      "   Solve approximation 2\n",
      "   Solution status is Unknown but at least Optimal is expected\n"
     ]
    }
   ],
   "source": [
    "#Sparsity density set to 0.1 and number of terms in the log-sum-exp is 2, like in the paper.\n",
    "s = 0.1\n",
    "p = 2\n",
    "Ex,App1,App2 = main(n_list,m_list,L_list,p,s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VOXZ+PHvPZnJvkDYIYGACIqIgKC4o6ivWrXaWpe6vrZYrbVqW23t4tba11q1r9b+3taq1VZF1Kq1Lq1aiVYrShCU1RWYCXuAJBOyzsz9++OchCFmnWQymcz9ua655sxZnnPPdu5znvOc54iqYowxJnV5Eh2AMcaYxLJEYIwxKc4SgTHGpDhLBMYYk+IsERhjTIqzRGCMMSnOEoHpMRE5SkQ+6qN1nS8ir/TRuo4QkU9EpEZEzuiLdSYTEXlYRH7Rg+VXicjcXgzJxMgSQRyJyNdFpMzdkGwWkZdF5MhExxVNROaKSHk3l1ERmdj8WlX/raqT4xBbibsub9S6HlPVE3t7Xe24FbhPVXNV9bmeFuZuOBvd30Pz44NeiLOjde71XbUxPV1E7hKRcjeedSLymzjE8YWkoaoHqGppb6+rN4jIJSLyVqLj6CuWCOJERL4H/C/wS2AEMBb4f8CXExmX6ZZxwKpYFoxOXq3c4SaW5sdBsYfXK24AZgGHAHnAscCyhEZk+p6q2qOXH0ABUAN8rYN5MnASxSb38b9AhjttLlAOXA9sAzYDZwCnAB8DO4EfR5V1M/A0sBAIAu8DB0VNV2Bi1OuHgV8AOUAdEHHjrQFG42wU3gEq3XXfB6S7y77plrfbnf+c5nijyt8fKHWXXwWc3mrdvwNedGN9F9innc/I766rObbDgEuAt1q9t28Dn7jl/RzYx42/GniyOXZ3/lOB5W5s/wGmtbPuz9zPpc5dd4b72Tzvfv6fAvPb+A4eddf7zTbKfBj4RTvrOwf4HMh3X58MbAGGua/vAQJu2UuBo6KWTQN+7MYcdKcXt/VdtbHeF4BrOviddvZd/sId3ut7if7dAZcBTUCjG8ff3enrgeO78X/4Pnv+D//dTrznAmWtxl0LPO8OnwKsdj+njcAP2innC+9nID8SHsBAfAAnASHA28E8twKLgeHAMHej9HN32lx3+RsBHzAf2A48jrPXdgBQD0xw57/Z/aOd5c7/A2Ad4HOnt5kIotZV3iq2g4E5gBcoAdZEbyzaKK+lDHf9n+JsmNKB49w/3eSode/ESTZe4DHgiXY+oxJ3Xd6ocXv9Qd3pzwP57ufSAPwLmICTkFcDF7vzznQ3JIfibDwvxtkYZbSz/vW4Gyr39Rs4R3WZwHT3O5nX6js4A+dIO6uN8lo+93bW95g7zxCcjeGpUdMucMd7cTaIW4BMd9p1wApgMiDAQcCQtr6rNtb5U5yE+23gQECipnXlu+w0EbT33tk7EXTl/3CrG9MpQC0wuI33k+3GuG/UuCXAue7wZtwkCgwGZrbzuXzh/Qzkh1UNxccQoEJVQx3Mcz5wq6puU9XtwC3AhVHTm4DbVLUJeAIYCtyjqkFVXYWzdzYtav6lqvq0O//dOBurObEEr6pLVXWxqoZUdT3wB+CYLi4+B8gFblfVRlV9HWev87yoeZ5R1ffcz+cxnI1qT/xKVavdz2Ul8Iqqfq6qVcDLwAx3vvnAH1T1XVUNq+ojOImj089JRIqBI4Efqmq9qi4HHmDv7+wdVX1OVSOqWtdOUT8QkcqoxyNR067E2diW4uw1v9A8QVUfVdUd7ndyF84edPN5mW8CP1XVj9Txgaru6Ow9uf4H+BXO77EM2CgiF7vTuvJd9pau/B9uVdUmVX0J58jiC+elVLUW+FtzjCKyL7Afzs5CczlTRCRfVXep6vtxeC9JxxJBfOwAhnZQTwxONcOGqNcb3HEtZahq2B1u3qhsjZpeh/MnbRZoHlDVCM6hdHR5XSYik0TkBRHZIiLVOOc5hnZx8dFAwI2h2QZgTNTrLVHDtez9PmLR+nNp73MaB3w/ekOMU4XSlc9pNLBTVYNR41q/rwCdu1NVB0U9mje6qGol8BQwFbgreiER+b6IrBGRKjfuAvZ8J8U41ULd5ibE36nqEcAg4DbgIRHZn659l72lK/+H6B2rjn43j7MnWX0deM5NEABfxTmi2CAib4jIYT2OfACwRBAf7+BU3XTU5HATzoap2Vh3XKyKmwdExAMURZVXi3PI3Gxk1HBb3c/+H7AW5/A6H6dqQLoYxyag2I2h2Vic+tju6u2ucQM4R1nRG+JsVV3QhWU3AYUikhc1rvX76lG8IjIduBRYANwbNf4o4IfA2TjVIYOAKvZ8JwGc8yI9oqp1qvo7YBcwhe59l7uJ+o2JyMhW0zv7bHrz//AKzo7YdJyE8HhLEKpLVPXLOFVQz+GcQ0p5lgjiwK2SuBH4nYicISLZIuITkZNF5A53tgXAT0VkmIgMded/tAerPVhEvuIehVyDU+Wx2J22HPi6iKSJyEnsXc2zFRgiIgVR4/JwTkrWiMh+wBWt1rUVpw6+Le/ibBSud9/zXOA0nOqt7tqOc8K2vXV11x+By0XkUHHkiMiXWm3c26SqAZx66/8RkUwRmQZ8A6dqq8dEJBPn+/8x8N/AGBH5tjs5D6eOfDvgFZEbcc6JNHsA+LmI7Ou+r2kiMsSd1tF3hYhc4zYhzhIRr1stlIfTcqg73+UHwAEiMt19Lze3mt5hHPTi/8E9cnga+DVQCLzqvtd09zqUArcKtRoIt18S4n7XLY9Y4kkGlgjiRFXvBr6HczJuO85e23dw9kLAabVTBnyIc6LvfXdcrP6G0/JkF07d6lfcHzvA1Th/4EqcutiWNvGquhbnT/i5W10yGudk89dxTrr9Eac1UrSbgUfc+c9u9b4bgdNxWr1U4JxcvchdT7e4h/O3AW+764rpnEdUeWU45wnuw/mcPsU5KdhV5+GcwN4EPAvcpKqvdjOM61tdR1Dhjv8fnBPu/6eqDTgnh3/h1nH/E+dcx8c4VSb17F0NdTfOnu0rOBu3B4Esd9rNtPNduepwqqG24HxfVwJfdc+xdPm7VNWPcU7mvobTgqt1G/wHcermK0WkrWsyevv/8DhwPPBUqyqlC4H1bpXn5Tifc3sOx/l8Wh6dVPcmLVG1G9MkOxG5Gad1Rkc/amOMaZMdERhjTIqzRGCMMSnOqoaMMSbF2RGBMcakuKQ4Az506FAtKSlJdBjGGJNUli5dWqGqwzqbLykSQUlJCWVlZYkOwxhjkoqIbOh8LqsaMsaYlGeJwBhjUpwlAmOMSXFJcY7AGJO8mpqaKC8vp76+PtGhDFiZmZkUFRXh8/liWt4SgTEmrsrLy8nLy6OkpASRrnZia7pKVdmxYwfl5eWMHz8+pjKsasgYE1f19fUMGTLEkkCciAhDhgzp0RGXJQJjTNxZEoivnn6+lghMn6oL1XH30rtZs2NNokMxxrgsEZg+VR4s508r/8SG6i5d52JMr3n22WcREdau7fatMeLmzTffZObMmXi9Xp5++umExRG3RCAiD4nINhFZ2ca0H4iIunciMinEX+0HoDi/uJM5jeldCxYs4Mgjj+SJJ2K5WV7bQqFQ5zN1YOzYsTz88MN8/etf76WIYhPPI4KHgZNajxSRYuAEwB/HdZt+yh90E0GeJQLTd2pqanj77bd58MEH90oEpaWlHH300Zx55plMmTKFyy+/nEgkAkBubi7f//73mTlzJvPmzWP79u0AzJ07lx//+Mccc8wx3HPPPWzYsIF58+Yxbdo05s2bh9/v/Ma//OUv8+c//xmAP/zhD5x//vlfiKukpIRp06bh8SS2ciZuzUdV9U0RKWlj0m+A63FurWhSjD/oZ3DGYPLT8zuf2Qw4t/x9Fas3VfdqmVNG53PTaQd0OM9zzz3HSSedxKRJkygsLOT9999n5syZALz33nusXr2acePGcdJJJ/HMM89w1llnsXv3bmbOnMldd93Frbfeyi233MJ9990HQGVlJW+88QYAp512GhdddBEXX3wxDz30EN/97nd57rnnuP/++zniiCMYP348d911F4sXL243vkTr0zQkIqcDG1X1gy7Me5mIlIlIWXMmNskvUB2waiHT5xYsWMC5554LwLnnnsuCBQtaph1yyCFMmDCBtLQ0zjvvPN56y7ndssfj4ZxzzgHgggsuaBkPtIwHeOedd1qqdi688MKW+UaMGMGtt97Ksccey1133UVhYWF832QP9NkFZSKSDfwEOLEr86vq/cD9ALNmzbK75wwQ/qCfg0ccnOgwTIJ0tuceDzt27OD1119n5cqViAjhcBgR4Y477gC+2PSyvaaY0eNzcnLaXV/0fCtWrGDIkCFs2rSpJ28h7vryiGAfYDzwgYisB4qA90VkZB/GYBKoIdzAlt1bGJs3NtGhmBTy9NNPc9FFF7FhwwbWr19PIBBg/PjxLXvu7733HuvWrSMSibBw4UKOPPJIACKRSEtLnscff7xlfGuHH354y3mHxx57rGW+9957j5dffplly5Zx5513sm7duni/1Zj1WSJQ1RWqOlxVS1S1BCgHZqrqlr6KwSTWxuBGFLWqIdOnFixYwJlnnrnXuK9+9as8/vjjABx22GH86Ec/YurUqYwfP75l3pycHFatWsXBBx/M66+/zo033thm+ffeey9/+tOfmDZtGn/5y1+45557aGhoYP78+Tz00EOMHj2au+66i0svvZTWtwZesmQJRUVFPPXUU3zrW9/igAP6/ogJ4njPYhFZAMwFhgJbgZtU9cGo6euBWapa0VlZs2bNUrsxTfIrDZRy1etX8fgpj3PgsAMTHY7pI2vWrGH//fdPdBhtKi0t5c477+SFF174wrTc3FxqamoSEFVs2vqcRWSpqs7qbNl4tho6r5PpJfFat+mfmq8hGJtvVUPG9CfW+6jpM/6gn/z0fAoyChIdijGAc03A3Llz25yWTEcDPWVdTJg+46/224liY/ohSwSmz/iDfjtRbEw/ZInA9ImmcBObd2+2IwJj+iFLBKZPbKzZSEQjdqLYmH7IEoHpE82dzdkRgUmU/tgN9d13382UKVNaOqzbsCEx3bNbIjB9IhAMANbrqEmc/tgN9YwZMygrK+PDDz/krLPO4vrrr++lyLrHEoHpE/5qPzm+HAoz+2/HW2bg6q/dUB977LFkZ2cDMGfOHMrLy+P6ObTHriMwfcIfdJqO2r1rU9zLP4ItK3q3zJEHwsm3dzhLMnRD/eCDD3LyySf3wgfSfXZEYPpEIBiwE8UmYfp7N9SPPvooZWVlXHfddb30jrvHjghM3IUiITYGN3LiuC71QG4Gsk723OOhv3dD/dprr3HbbbfxxhtvkJGR0aX31NvsiMDE3ebdmwlpyE4Um4Toz91QL1u2jG9961s8//zzDB8+vNffe1dZIjBxZ53NmUTqz91QX3fdddTU1PC1r32N6dOnc/rpp8fhE+hc3Lqh7k3WDXVyW7B2Ab9895e8/rXXGZY9LNHhmD5m3VD3jZ50Q21HBCbu/NV+srxZDM0amuhQjDFtsJPFJu4CwQDFecXWdNT0O9YNtcOOCEzcNV9DYIzpnywRmLgKR8KUB8ut+2lj+jFLBCauttZupSnSZEcExvRjcUsEIvKQiGwTkZVR434tImtF5EMReVZEBsVr/aZ/aO51dFz+uARHYoxpTzyPCB4GTmo17lVgqqpOAz4Gbojj+k0/0HwNgV1MZhIpLS2N6dOntzxuv733rnBevnw5L730Uq+VlwhxazWkqm+KSEmrca9EvVwMnBWv9Zv+IRAMkJGWwfDsxF01aUxWVhbLly+PS9nLly+nrKyMU045JS7l94VEniO4FHi5vYkicpmIlIlIWXP3ryb5bKjeQHFeMR6x01Gmf6mqqmLy5Ml89NFHAJx33nn88Y9/BOCKK65g1qxZHHDAAdx0000tyyxZsoTDDz+cgw46iEMOOYSqqipuvPFGFi5cyPTp01m4cOFe6wiHw/zgBz/gwAMPZNq0afz2t78F4NZbb2X27NlMnTqVyy67rOWK43vvvbflRjXNneTt3r2bSy+9lNmzZzNjxgz+9re/9fpnkZDrCETkJ0AIeKy9eVT1fuB+cK4s7qPQTC9rvobAGIBfvfcr1u7s3TuE7Ve4Hz885IcdzlNXV8f06dNbXt9www2cc8453HfffVxyySVcffXV7Nq1i/nz5wNw2223UVhYSDgcZt68eXz44Yfst99+nHPOOSxcuJDZs2dTXV1NdnY2t956K2VlZS1dVEe7//77WbduHcuWLcPr9bJz504AvvOd77R0WXHhhRfywgsvcNppp3H77bezbt06MjIyqKysbInluOOO46GHHqKyspJDDjmE448/vsOO77qrzxOBiFwMnArM02To38LELKIRAsEAR4w+ItGhmBTXXtXQCSecwFNPPcWVV17JBx980DL+ySef5P777ycUCrF582ZWr16NiDBq1Chmz54NQH5+fqfrfe2117j88svxep1NbXNX1IsWLeKOO+6gtraWnTt3csABB3Daaacxbdo0zj//fM444wzOOOMMAF555RWef/557rzzTgDq6+vx+/292m1HnyYCETkJ+CFwjKrW9uW6Td/bVruNhnCDdTZnWnS2597XIpEIa9asISsri507d1JUVMS6deu48847WbJkCYMHD+aSSy6hvr4eVe321fFtLVNfX8+3v/1tysrKKC4u5uabb6a+vh6AF198kTfffJPnn3+en//856xatQpV5a9//SuTJ0/utffdWjybjy4A3gEmi0i5iHwDuA/IA14VkeUi8vt4rd8knt2n2PR3v/nNb9h///1ZsGABl156KU1NTVRXV5OTk0NBQQFbt27l5ZedU5n77bcfmzZtYsmSJQAEg0FCoRB5eXkEg8E2yz/xxBP5/e9/33Jv4507d7Zs9IcOHUpNTU1LV9eRSIRAIMCxxx7LHXfcQWVlJTU1NfzXf/0Xv/3tb1vOIyxbtqzXP4d4tho6r43RD8Zrfab/se6nTX/R+hzBSSedxKWXXsoDDzzAe++9R15eHkcffTS/+MUvuOWWW5gxYwYHHHAAEyZM4IgjnKrN9PR0Fi5cyFVXXUVdXR1ZWVm89tprHHvssdx+++1Mnz695dxDs29+85t8/PHHTJs2DZ/Px/z58/nOd77D/PnzOfDAAykpKWmpagqHw1xwwQVUVVWhqlx77bUMGjSIn/3sZ1xzzTVMmzYNVaWkpKTN3lJ7wrqhNnHzm6W/4c+r/0zZ+WWkedISHY5JkP7cDfVAYt1Qm34pEAxQlFtkScCYfs4SgYkbf7XfupYwJglYIjBxoar4g347UWwAvnCLRtO7evr5WiIwcVFRV0FdqM5OFBsyMzPZsWOHJYM4UVV27NhBZmZmzGXYHcpMXDT3OmrdT5uioiLKy8uxrmLiJzMzk6KiopiXt0Rg4qKl6aglgpTn8/kYP358osMwHbCqIRMXgWAAr3gZlTsq0aEYYzphicDEhT/oZ3TuaLweO+g0pr+zRGDiwl/tt/sUG5MkLBGYXqeqBIIBOz9gTJKwRGB63a6GXdQ01VgiMCZJWCIwvc46mzMmuVgiML2uuftpOyIwJjlYIjC9zh/04xEPY3LHJDoUY0wXWCIwvW5D9QZG5YzCl+ZLdCjGmC6wRGB6XaDaWgwZk0w6vNpHRDJxbjR/FDAaqANWAi+q6qr4h2eSkT/o5+TxJyc6DGNMF7WbCETkZuA0oBR4F9gGZAKTgNvdJPF9Vf0w/mGaZFHVUEV1Y7V1P21MEunoiGCJqt7czrS7RWQ40O7xv4g8hHM0sU1Vp7rjCoGFQAmwHjhbVXd1P2zTX1lnc8Ykn3bPEajqi63HiYhHRPLd6dtUtaMbCT8MnNRq3I+Af6nqvsC/3NdmAGnpftquITAmaXR6slhEHheRfBHJAVYDH4nIdZ0tp6pvAjtbjf4y8Ig7/AhwRjfjNf2cP+hHEIryYu8b3RjTt7rSamiKqlbjbLRfwqkOujDG9Y1Q1c0A7vPw9mYUkctEpExEyuyGFskjUB1gRM4IMtIyEh2KMaaLupIIfCLiw0kEf1PVJiDu95xT1ftVdZaqzho2bFi8V2d6iT/oZ1ye3bDemGTSlUTwB5wTuznAmyIyDqiOcX1bRWQUgPu8LcZyTD8VCAas+2ljkkyniUBV71XVMap6ijp3n/YDx8a4vueBi93hi4G/xViO6YeCjUF21u+0FkPGJJl2E4GIXCAiX5iujpCI7CMiR3aw/ALgHWCyiJSLyDeA24ETROQT4AT3tRkg7Ib1xiSnjq4jGAIsE5GlwFJgO84FZROBY4AKOmj+qarntTNpXmyhmv4uUO30OmpVQ8Ykl3YTgareIyL3AccBRwDTcLqYWANcqKr+vgnRJIvmI4KiXGs6akwy6bCvIVUNA6+6D2M65K/2MzxrONm+7ESHYozpBut91PQaazFkTHKyRGB6jT/otxPFxiQhSwSmV9Q21VJRV2F9DBmThLrS19AIEXlQRF52X09xm4Ia08LuU2xM8urKEcHDwD9xbkwD8DFwTbwCMsnJeh01Jnl1JREMVdUngQiAqoaAcFyjMkmn+T4EdkMaY5JPVxLBbhEZgtvRnIjMAariGpVJOv6gnyGZQ8jx5SQ6FGNMN3V4HYHrezh9BO0jIm8Dw4Cz4hqVSTr+ar9VCxmTpDpNBKr6vogcA0wGBPjI7YramBb+oJ85o+YkOgxjTAw6TQQikgacgnOfYS9wooigqnfHOTaTJOpCdWyr3WYthoxJUl2pGvo7UA+swD1hbEy08mA5YC2GjElWXUkERao6Le6RmKRl3U8bk9y60mroZRE5Me6RmKTV3P203bDemOTUlSOCxcCz7k1qmnBOGKuq5sc1MpM0/EE/gzIGUZBRkOhQjDEx6EoiuAs4DFjh3qrSmL34g9Z01Jhk1pWqoU+AlZYETHsC1QE7P2BMEuvKEcFmoNTtdK6heaQ1HzUADeEGNu/ebInAmCTWlUSwzn2ku48eE5FrgW/idFuxAvhvVa3vjbJN39oY3IiidkMaY5JYV64svqU3VygiY4DvAlNUtU5EngTOxenl1CQZazpqTPJrNxGIyP+q6jUi8nfcDueiqerpPVxvlog0AdnAph6UZRKouddRSwTGJK+Ojgj+4j7f2ZsrVNWNInIn4AfqgFdU9ZXW84nIZcBlAGPH2kamv/IH/eSl51nTUWOSWEethq4CUNU32nrEukIRGQx8GRiPc7ObHBG5oPV8qnq/qs5S1VnDhg2LdXUmzgJBp8WQiCQ6FGNMjDpKBPHqVuJ4YJ2qbnd7MX0GODxO6zJx5q+2G9Ybk+w6qhrKFpEZOFcSf4Gqvh/jOv3AHBHJxqkamgeUxViWSaCmcBObdm/ilAmnJDoUY0wPdJQIxuBcVdxWIlDguFhWqKrvisjTwPtACFgG3B9LWSaxNu3eREQjjMsfl+hQjDE90FEi+FRVY9rYd0ZVbwJuikfZpu9YiyFjBoaudDFhTJuaryGwG9Ybk9w6SgQ/7LMoTFLyV/vJ8eVQmFmY6FCMMT3QbiJoq22/MdH8Qb81HTVmALCqIROzQDBg1ULGDAAdJgIRSRORX/dVMCZ5hCIhNgY32n0IjBkAOkwEqhoGDhY79jetbN69mZCGrMWQMQNAV7qhXgb8TUSeAnY3j1TVZ+IWlen3mu9TbFVDxiS/riSCQmAHe19ApjhdQ5gU1dL9tFUNGZP0unI/gv/ui0BMcvEH/WR5sxiWZR0CGpPsOm01JCKTRORfIrLSfT1NRH4a/9BMfxaodloM2ekjY5JfV5qP/hG4AWgCUNUPce4oZlJY8zUExpjk15VEkK2q77UaF4pHMCY5hCNh5xoCu0+xMQNCVxJBhYjsg3u7ShE5C9gc16hMv7a1ditNkSY7IjBmgOhKq6ErcbqJ3k9ENgLrgC/cUcykDrthvTEDS1daDX0OHC8iOYBHVYPxD8v0Zy3dT1vTUWMGhK60GgqLyO1AbXMSEJFY705mBoBAMEC6J53h2cMTHYoxphd05RzBKne+V0Skub9hazOYwvzVforzivGI9VlozEDQlX9ySFWvx2lG+m8RORj3xLFJTf6g31oMGTOAdCURCICqPgmcDfwJmNCTlYrIIBF5WkTWisgaETmsJ+WZvhPRCOXBcjtRbMwA0pVWQ99sHlDVVSJyJHBGD9d7D/APVT1LRNKB7B6WZ/rI9trt1Ifr7Yb1xgwg7SYCETlOVV8HxolI6399TawrFJF84GjgEgBVbQQaYy3P9C27T7ExA09HRwTHAK8Dp7UxrSe9j04AtgN/EpGDgKXA1aq6O3omEbkMuAxg7FirhugvrOmoMQNPu4lAVW9yn3u791EvMBO4SlXfFZF7gB8BP2u1/vtxLmRj1qxZdnK6n/AH/Xg9XkZmj0x0KMaYXtLpOQIRGQRcBJREz6+q341xneVAuaq+675+GicRmCQQCAYoyi0izZOW6FCMMb2kKyeLXwIWAyuASE9XqKpbRCQgIpNV9SNgHrC6p+WavuGv9lu1kDEDTFcSQaaqfq+X13sV8JjbYuhzwG5+kwRUFX/Qz+yRsxMdijGmF3UlEfxFROYDLwANzSNVdWesK1XV5cCsWJc3ifFJ5SfUheqYNHhSokMxxvSiriSCRuDXwE/Yc0Wx0sOLykzyKQ2UAnBU0VGJDcQY06u6kgi+B0xU1Yp4B2P6t9JAKdOGTmNo1tBEh2KM6UVd7XSuNt6BmP5te+12VlSsYG7x3ESHYozpZV05IggDy0VkEXufI4i1+ahJQm+UvwFgicCYAagrieA592FSWGmglDG5Y5g4aGKiQzHG9LKu3KHskb4IxPRfdaE6Fm9ezFmTzkLEbkVhzEDTUadzT6rq2SKygjbuP6Cq0+Iamek3Fm9aTEO4waqFjBmgOjoiuNp9PrUvAjH9V2l5KXm+PA4ecXCiQzHGxEG7rYZUdbM7+G1V3RD9AL7dN+GZRItohNJAKUeOORKfx5focIwxcdCV5qMntDHu5N4OxPRPKypWsLN+p1ULGTOAdXSO4AqcPf8JIvJh1KQ84O14B2b6h9JAKV7xcsSYIxIdijEmTjo6R/A48DLwP+zdTXSwJ/0MmeRSGiiVARs6AAAdYElEQVTl4BEHU5BRkOhQjDFx0tE5gipVXa+q5wFDgC8Dp+Pcl8CkgEAwwKeVn3JM8TGJDsUYE0edniMQkZ8Bj+Akg6E4t5j8abwDM4n3RsC9mrhobmIDMcbEVVeuLP46MENV6wFE5HbgfeAX8QzMJF5poJSJgyZSnG83qjdmIOtKq6H1QGbU6wzgs7hEY/qNqoYqyraWWWshY1JAV44IGoBVIvIqzhXGJwBvici9YJ3PDVRvb3ybsIYtERiTArqSCJ51H81K4xOK6U9KA6UUZhZy4NADEx2KMSbOupIIFgITcY4GPms+V2AGrqZIE29tfIvjxx2PR7pSe2iMSWbt/stFxCsidwDlOK2GHgUCInKHiPS4rwERSRORZSLyQk/LMr3r/a3vE2wKWrWQMSmio929XwOFwHhVPVhVZwD7AIOAO3th3VcDa3qhHNPLSgOlZKRlMGfUnESHYozpAx0lglOB+aoabB6hqtXAFcApPVmpiBQBXwIe6Ek5pvepKosCi5gzag7ZvuxEh2OM6QMdJQJV1bbuQxCmjfsTdNP/AtcDkfZmEJHLRKRMRMq2b9/ew9WZrvq08lM21my0aiFjUkhHiWC1iFzUeqSIXACsjXWFInIqsE1Vl3Y0n6rer6qzVHXWsGHDYl2d6abSQCkAxxRZtxLGpIqOWg1dCTwjIpcCS3GOAmYDWcCZPVjnEcDpInIKzoVq+SLyqKpe0IMyTS8pDZQydchUhmVb8jUmVXTU6dxGVT0UuBXn6mI/cKuqHqKqG2NdoareoKpFqloCnAu8bkmgf6ioq2BFxQqrFjImxXTl5vWvA6/3QSwmwd4sfxNFLREYk2K6ckFZ3KhqKXalcr+xKLCI0TmjmTR4UqJDMcb0Ibts1ABQF6pj8abFzC2ei4gkOhxjTB+yRGAAeHfzu9SH661ayJgUZInAAE5roVxfLrNGzEp0KMaYPmaJwBDRCG+Uv8ERY47Al9bjbqSMMUnGEoFhVcUqKuoqrFrImBRlicCwKLCINEnjqDFHJToUY0wCWCIwlJaXMnPETAoyChIdijEmASwRpLjyYDmf7PqEuUVzEx2KMSZBLBGkuDfK3wDg2OJjExyJMSZRLBGkuEWBRUwomEBxfnGiQzHGJIglghQWbAyydMtSay1kTIqzRJDC3t74NiENWbWQMSnOEkEKWxRYRGFmIQcOPTDRoRhjEsgSQYpqijTx743/5uiio0nzpCU6HGNMAlkiSFHLti4j2Bi08wPGGEsEqWpRYBHpnnQOG3VYokMxxiSYJYIUpKq8Uf4Gh446lGxfdqLDMcYkmCWCFPR51ecEggGrFjLGAJYIUtKiwCIAjik6JsGRGGP6gz5PBCJSLCKLRGSNiKwSkav7OoZUVxoo5YAhBzAiZ0SiQzHG9AOJOCIIAd9X1f2BOcCVIjIlAXGkpIq6Cj7c/qFVCxljWvR5IlDVzar6vjscBNYAY/o6jlT17/J/o6glAmNMi4SeIxCREmAG8G4b0y4TkTIRKdu+fXtfhzZgLQosYmTOSCYPnpzoUIwx/UTCEoGI5AJ/Ba5R1erW01X1flWdpaqzhg0b1vcBDkD1oXoWb17M3KK5iEiiwzHG9BMJSQQi4sNJAo+p6jOJiCEVvbflPepCddbJnDFmL4loNSTAg8AaVb27r9efyhYFFpHjy2HWyFmJDsUY048k4ojgCOBC4DgRWe4+TklAHCmjKdLEkx89yT/W/YPDRx9Oelp6okMyxvQj3r5eoaq+BVgFdR+IaISX173M75b/jkAwwIzhM7h6pl22YYzZW58nAhN/qsqb5W9y77J7+XjXx0waPInfzfsdR405yk4SG2O+wBLBAFO2pYx7l93Lsm3LKM4r5ldH/YqTxp+ER6w3EWNM2ywRDBBrdqzhnmX38PbGtxmWNYyfzfkZZ+57Jj6PL9GhGWP6OUsESW5D9QbuW3Yf/1j/D/LT87n24Gs5b7/zyPJmJTo0Y0ySsESQpLbs3sLvP/g9z336HOlp6cw/cD6XTL2E/PT8RIdm+thbn1Tw0srN+DxChi+N9DQPGV4P6V7nuWWcz+M+73md0TyPN40huelkp9smIRXZt54sGoKQlkFlqJYHVjzAgrULiBDh7Mlnc9m0yxiaNTTREZoEePjtddz6wmpy0r14PEJjKEJDKExEu19WhtfD3MnD+NK00czbbzg5GbZ5SBX2TfdnOz+Hj16GtS+xO7CYPw8dwSN5WdRGmjhtn9O44qArKMorSnSUJgFC4Qg/f2E1j7yzgROmjOCec6e37M2rKqGIukkh0pIcml83uK/3TIvQ0BRm1aZqXlyxmX+u2kqG18Nx+w3nS9NGcdx+w+1IYYCzb7c/iURg0zL46EVY+xJsX4MCL4zel7vGj2dHpIHjgpVcpQVMHDwbcq3T1lRU0xDiqsffZ9FH25l/1Hh+dPL+pHn2NAsWEXxpgi/NQ05G18v9GvCzU6dQtn4nL67YzEsrtvDyyi1k+jzM228EX5o2imMnDycrPa3335RJKFGN4Riyj82aNUvLysoSHUZ8hBpg3Zuw9kVn779mC4gHxh7O2pJD+GVwFct2reHAoQfyw9nXc9DOjfCvW2H7Whg9A+bdBPtY30GpYlNlHd94pIyPtwa55fQDuGDOuLitKxxR3lu3kxdXbOIfK7dQUdNIli+N4/YfzqkHjmKuJYX4UoWarZBZAL7YGn+IyFJV7bRPGUsEiVC7Ez55xdn4f/Y6NNaALwcmzoP9vkTVuMO4b+2jPPnxkxSkF3DNwddwxsQz9lwLEAnDhwth0S+hKgDjj4Hjb4IxByf2fZm4WlFexTceWUJtY5jfnT+TYyb1Xa+84Yjy7rodvPjhZv6xcgs7djeSnZ7GvP1H8KUDRzF38jAyfZYUYlK3C3Z8Djs+3fPY+Rns+MzZNlz4XMw7e5YI+ptd653qno9egg3/AQ1D7giYfDJM/hKMP5qIN51nP3mWe96/h6rGKs6ZfA5XTr+SgoyCtssMNUDZQ/Dmr6F2B+x/Ohz3Mxg2qU/fmom/V1Zt4eonllOYk85Dl8xm8si8hMUSCkd4d91OXvhwM/9ctYWduxvJSU/j+CkjOGPGGI7ZdxgeTzevYFeF+iqo3gThRkjzQVq68+xpHva6z+ng8UJ3rpIPN0FTLTTVtXq440KtXjfVgTcDMvIgPRcy8iEj1x3O2/PwdrHurbHWOefXsrH/zN3Yf+r8d5uJBwaNgyETYcg+zvOkk2BQcfc+z+biLBEkWKjB2eB/+hp8+i/YvsYZP2x/2O8UZ+M/egZ4nL38lRUruW3xbazcsZKZw2fy40N/zOTCLt48piEI7/wO/vNb54c8/XyY+yMosBPJyU5VefCtddz20hqmFQ3ijxcdzPC8zESH1SIUjrD48z3VR7tqmygZks2Fh5XwtVlF5Gf63I18JVRtdDb01dHPzcObnL3f7midIDw+N4H4IBLae2MfCcXnA/D43KSQC+l5e4Yz8sCbBdXlzka/euPey+WN2ntjX+g+Dy4Bb+91CmmJIBF2fu5s9D95Fdb/2/kBpqXD2MNg3xOdBFA4Ya9FdtXv4p737+GZT55hSNYQvnfw9zh1wqmx9Qm0uwL+fTcs+SMgcMh8OOr7kF3YO+/PUFHTwMqNVe6jmh27Gzh56ii+MnMMg7J7t1fXpnCEm55fxePv+jl56kjuPnt623XyTfXOXrQvy9kIxlsk7Gy0G4LQUOMOV9NUF2T1p5/xyScfQfUmijw72SeziiHhCjyhur3LEA/kjoSCMZA/GvLd57xR4Mt23k8k5DyHG509+nCTOz5qONxqOOIOe3zO5+HLbvUc/WhrmjvszXR25hqCUe81uPf7bqh2X9dETave85k01kLeSHeDH73Rn+Akiz5giaAvNO6G9W+5e/2vOYkAnKw+8QSYeDyUHNnmlx6OhHnq46f47bLfUttUy/n7n8/lB11Obnov/EAq/VD6K/jgcedQ9vDvwpwr+uzHN1Bsq65n5aYqVpRXs3KTs/HfXFXfMr1kSDZZ6V7WbK4m3evh5KkjOXf2WOZMKOxaIo9EoKEK6irdjUt1ywanrmYXz7+7loodFRw2xseM4WlI88aooXqveQk37ilT0vZsyJo3bs3DbT5ngy/T2XtFv7ixa71xawg6OzgdkTQac0ayMTyYNbtz2RgpJHNIMdOmTGHq/lNIKxjjVIumWaPFeLNEEA+qTmud5g3/hv84f0JvFow/2tnwT5znZP4OLN+2nNvevY21O9dy6MhDueHQG9hnUMfLxGTbWnj957D2BcgZBkdfBwd+DbIGd69+tRepKk1hJd2bwE7wQo1OvWx9JSCoJ42K3WHWbqtlzdZa1mzZzcotu9la00SYNMKSxtghuUwZU8gBYwYzdUwBU0bnU5Dl7H2vDlTw93dX8/bKT/A1VDG5IMTxJekcMtJDbiTonAyMftRXusOVQMf/v4h48WTmO1UNmfluXXVe1LM73uNz9mBDdc4RQvNzUy2E6t168OjnVvNo2FmhNzOqXtytG28Z7qC+vHk4u9DZyHucI5cdNQ08sSTAo4s3sLmqnqLBWVx02DjOnlXc60dQ5ossEfSWmu3gf2dPXX91uTN+2P7ORn/i8U7Vj6/zetuKugp+s/Q3PP/Z84zIHsF1s6/jxHEnxr9r6MAS+NctTnUVOM3RCic49ZKFE/Z+5AxtN0moKjt3N7KrtpFgfYiahhA19SGC7nNNg/PYM61p79fufKGIMrogk6ljCphWVMDUMQUcOKaAIbndaPQerakeaiucqrHaCti9Y8/r3dudjb47TXdXIA1fuEV293i87snKNEA73EOOIITTC/DmFiJZg50knDUYMgdFDRe0bOTXVsINL25gVziD/zn3CA6bPKZvkna4yXmOU9VSKBzhldVbefg/63lv3U4yfR7OnDGGiw8vYb+R1i1KvFgi6K6GoLMHvW111GONsyEBZy9owjHuXv/x3ToR2xRp4om1T/D/lv8/6sP1XHLAJcw/cD7Zvuw4vZk2qIJ/MWx636nCan5U+kEje+ZLz0MLx1OXV8J23xjW6wjWNAylrHowSyp8VNW3f9JNiDDIF2FIRoTC9AiDM8IU+CIUeEPk+8IUeEPkesNk0sSOqiBbd1VTVbObdJpIJ8yQLBiT62Fkrofh2cLQTCFDQhBucI68Qo3ucJOzZ1u309noNwbbfsseLw2+wQTTCtiheWxuyiXQkMX2SD47yScouYzMz6RkcDrFgzMpLkhnVL6XDI869eCRsFNPre5zJOI+N49z96KzBu29Yc8azPraDJ5eXcOCDyrZURtizKAszp5VzNmzixhV0Hab8JdWbObahcsZnp/Bny6ZzcThiWsZFE+rN1Xz53fW89zyjdQ3RTh0fCGXHF7CCVNG4E3r/EixvinM9mADFTUN7nPjXq937G5AFbLS08hOTyM73esM+5zXWele97l5ehpZPu+eYXeZtO62fGqDCKSJ4BHB4wGPCGkiiNAn9waxROAKhSME65290ur6JoK7awhv+wRfxRoyKz8mr+pjBu/+jMGNm1uWqSeDzz1j+ViLWdU0hg8i41nrnUx6egaZ0T8mn/OjyUpPw+dtIJK2g1DaDhrZQb1uoya8jWBoO5VNW2iM1DN18CFcNPlqJhSMJ9PnIcuXRoYvjSxfGr406fObxoTCEfzbK9m0/iMqyz+iafuneKvWU1AXoFg3Uyzb8cqeJNHoyaQmuxhPmhefNuCNNJIWbsATbkDC9Uh0XXWsMZFGo3ppxEsTXsKedMSbgdeXQUZmJpmZWXh9mU6zvexCIllD2CX5bG7MYX1dNh/vzmDFLh8f7vKyI5wNOH+6cYXZ7Dsij0kjcpk0Io9JI/KYMCyHDG982743hMK8tnobTyzx8+9PKvAIzJ08nHNnF3PcfsPxpnlQVX7/xuf86h9rmTl2EH+8aFbsR0dJpLK2kYVLAvz5nQ1srKxjdEEm588Zx34j8zrc0Acb2t4ZGZztY2huBkNzMxCB2sYwdY1hdjeGqGsMO6+bwn38LtvXUZLweJqHhXvPm87h+8TWl5glAuCup17j42VvMUkCTPaUM1kCjJfNLRu3Jk3jMx3N51KM31vCpozxVGTtQ33OGHKzMsjP8pKX6cPrEYINtexo2EJl41aqQ1upCW+jTrfTQAVNsgP17F09oOEMIk2DiTQNRpsKCe3el3DNZNq7S2eaR8j0eshKTyPD6yYYXxqZPg+ZvjQyfWmkez14PYLX4z6nifu853Va62kejzvsIc0Dmyrr+XR7DZ9urWFdxW4aw3s29KMKMpk4PJeJw3PZd3ge+w7NYN+MSgbVB2DnOvcoYh2gzobYm+U8N5949Ga6Jx4zo05Itp4vw52e4Tb9y3Cay7W0D0+jqq6JVRur+HBjFSvKq1ixsQr/zj2f79jCbCYOz2VTZR2fb9/7PRQXZjFpeB6TRjob/X2H5zFxeG6/uNgpsLOWhUsCPLU0wNbqBobnZXDWwUVsDzbw1NJyTjtoNL8+a1q/iLUvhSPKv9Zs5ZF31vP2pzv2mpaf6WVoXgbDcjNanoe1vE5nWG4mQ/PSGZKT0aXzTpGIUh8KtySJ2sYwtVGJorYpTF1jiNrG2Drua01ViagSjkBElUhEiag7HDVNVQm3M+2/jyhh0ojYjg77dSIQkZOAe4A04AFVvb2j+WNNBOV/voyizxeiCDXZRewaNJHKQROoKRxLfeEYmvIHE/Y0URvaTbAxSLAxSE1TDTWNNc7rpiA1jTVU1FWwo37vH2i6J53RuaMZkzeGotwixuSOYXTuaIpyixiZM5p0yaWuKUx9Y4TaJueHVtcUpqEp4oxvcl7XNTqdfzVPbx5f3xSmvskZXx9y5msKRwhFlFBYCUUihCPOiVfn2Xkd6uTX6xEoLsxm3+G57DN8z4Zyn2E55GX235vYVNY2snJjNR9urGTlxio+3VbD6EFZTB6R17KnP3F4blJ0jhYKR1j00XYWLvHz+tptRBS+e9xErjl+UvcvxBpg1lXspqquiaG56QzNzUi5pNjb+m0iEJE04GPgBKAcWAKcp6qr21sm1kTw+//8gufLF1ETaSTYWENIO76oxCMecnw55PnyyEvPIzc9lzxfHoMzBzMmd8xeG/0hWUP65e0fm/csQu4jHFaaInuSxJCcdPtz9SNbqurZVFXHzLGDEx2KGYC6mggSsft0CPCpqn4OICJPAF8G2k0EsRo+bApTQ9Xk+dyNenre3sPpeeT69gxne7OT/ubuIm61kG3rk8LIgkxGFvSfK4VNakpEIhgDBKJelwOHtp5JRC4DLgMYO3ZsTCv6yr5f4Sv7fiWmZY0xJlUkom6jrV3uL9RPqer9qjpLVWcNG9Z3vSwaY0yqSUQiKAeiu9IrAjYlIA5jjDEkJhEsAfYVkfEikg6cCzyfgDiMMcaQgHMEqhoSke8A/8RpPvqQqq7q6ziMMcY4EtLoWlVfAl5KxLqNMcbsrf81hDfGGNOnLBEYY0yKs0RgjDEpLik6nROR7cCGGBcfClT0MAQrY2CW0R9isDKsjHiWMU5VO78QS1UH9AMoszKsjP4ag5VhZfRFGZ09rGrIGGNSnCUCY4xJcamQCO63MqyMfhyDlWFl9EUZHUqKk8XGGGPiJxWOCIwxxnTAEoExxqS4AZsIROQhEdkmIit7UEaxiCwSkTUiskpEro6hjEwReU9EPnDLuCXGWNJEZJmIvBDj8utFZIWILBeR7t/30yljkIg8LSJr3c/ksG4uP9ldf/OjWkSuiSGOa93PcqWILBCRbt/iS0Sudpdf1dUY2vpNiUihiLwqIp+4zx3ec7KdMr7mxhERkc7vL9t2Gb92v5cPReRZERkUQxk/d5dfLiKviMjo7pYRNe0HIqIiMjSGOG4WkY1Rv5NTYolDRK4SkY/cz/aObsawMGr960VkeQzvY7qILG7+z4nIITGUcZCIvOP+d/8uIvkdlRGzeLdPTdQDOBqYCazsQRmjgJnucB7OvZandLMMAXLdYR/wLjAnhli+BzwOvBDje1kPDO3hZ/oI8E13OB0Y1IOy0oAtOBe8dGe5McA6IMt9/SRwSTfLmAqsBLJxOl58Ddg3lt8UcAfwI3f4R8CvYihjf2AyUArMijGOEwGvO/yrGOPIjxr+LvD77pbhji/G6V14Q2e/uXbiuBn4QTe+z7bKONb9XjPc18O7+z6ipt8F3BhDDK8AJ7vDpwClMZSxBDjGHb4U+Hl3futdfQzYIwJVfRPY2cMyNqvq++5wEFiDsyHqThmqqjXuS5/76NYZehEpAr4EPNCd5XqTuydyNPAggKo2qmplD4qcB3ymqrFcMe4FskTEi7Mx7+6NjfYHFqtqraqGgDeAMztbqJ3f1JdxEiTu8xndLUNV16jqR12Mvb0yXnHfC8BinBs+dbeM6qiXOXTyO+3gP/Yb4PrOlu+kjC5rp4wrgNtVtcGdZ1ssMYhzE/OzgQUxxKBA8x58AZ38TtspYzLwpjv8KvDVjsqI1YBNBL1NREqAGTh79N1dNs09tNwGvKqq3S3jf3H+WJHurjuKAq+IyFJx7gfdXROA7cCf3CqqB0QkpwfxnEsnf662qOpG4E7AD2wGqlT1lW4WsxI4WkSGiEg2zt5acSfLtGeEqm52Y9sMDI+xnN50KfByLAuKyG0iEgDOB26MYfnTgY2q+kEs64/yHbea6qHOqtvaMQk4SkTeFZE3RGR2jHEcBWxV1U9iWPYa4Nfu53kncEMMZawETneHv0bsv9MOWSLoAhHJBf4KXNNqr6lLVDWsqtNx9tIOEZGp3Vj3qcA2VV3a3fW2coSqzgROBq4UkaO7ubwX57D1/1R1BrAbpyqk28S5M93pwFMxLDsYZy98PDAayBGRC7pThqquwak+eRX4B/ABEOpwoSQhIj/BeS+PxbK8qv5EVYvd5b/TzXVnAz8hhgTSyv8B+wDTcZL9XTGU4QUGA3OA64An3b377jqPGHZYXFcA17qf57W4R9PddCnO/3UpTvV0Y4yxdMgSQSdExIeTBB5T1Wd6UpZblVIKnNSNxY4ATheR9cATwHEi8mgM697kPm8DngU6PHHVhnKgPOpo5mmcxBCLk4H3VXVrDMseD6xT1e2q2gQ8Axze3UJU9UFVnamqR+McjseyxwewVURGAbjP7VZBxJuIXAycCpyvbqVyDzxO96sh9sFJ0B+4v9ci4H0RGdmdQlR1q7vzFAH+SPd/q+D8Xp9xq2bfwzma7vDEdWtu1eNXgIUxrB/gYpzfJzg7Pd1+H6q6VlVPVNWDcRLSZzHG0iFLBB1w9yAeBNao6t0xljGsuQWHiGThbMjWdnV5Vb1BVYtUtQSnOuV1Ve3WHrCI5IhIXvMwzonFbrWmUtUtQEBEJruj5gGru1NGlJ7sZfmBOSKS7X4/83DO3XSLiAx3n8fi/Nljjed5nD887vPfYiynR0TkJOCHwOmqWhtjGftGvTydbvxOAVR1haoOV9US9/dajtPYYks34xgV9fJMuvlbdT0HHOeWNwmncUN3e/A8HlirquUxrB+ccwLHuMPHEcPORtTv1AP8FPh9jLF0LB5noPvDA+ePvRlowvlBfiOGMo7EqVv/EFjuPk7pZhnTgGVuGSvppPVBJ2XNJYZWQzj1+x+4j1XAT2Jc/3SgzH0vzwGDYygjG9gBFPTgc7gFZyO1EvgLbsuQbpbxb5xE9gEwL9bfFDAE+BfOn/xfQGEMZZzpDjcAW4F/xlDGp0Ag6nfaWYuftsr4q/uZfgj8HRjT3TJaTV9P562G2orjL8AKN47ngVExlJEOPOq+n/eB47r7PoCHgct78Ns4Eljq/sbeBQ6OoYyrcVorfgzcjtsbRG8/rIsJY4xJcVY1ZIwxKc4SgTHGpDhLBMYYk+IsERhjTIqzRGCMMSnOEoExMRKRms7nMqb/s0RgjDEpzhKBMcakOEsExhiT4iwRGGNMirNEYIwxKc4SgTHGpDhLBMbELltEyqMe30t0QMbEwnofNcaYFGdHBMYYk+IsERhjTIqzRGCMMSnOEoExxqQ4SwTGGJPiLBEYY0yKs0RgjDEp7v8Do3NIU0jRpMgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5a94849400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXd4VFXawH9nJr2RSgpJSCgBAgIBxAIoimsX66rYy9pX17Wtn1vsq2JZde19LWBFRVfWhgSkSBWQXkIKCSUTCExIMknmfH+cmTCETDJJps/5Pc88U+6dc965M3Pfe94qpJRoNBqNRgNg8LUAGo1Go/EftFLQaDQaTStaKWg0Go2mFa0UNBqNRtOKVgoajUajaUUrBY1Go9G0opWCpkcIISYIITZ4aa5LhRDfeWmucUKITUIIsxDiHG/M6Q8IIe4TQrzhpblybcfX6I35NK4hdJ6CZxBCXALcAQwG9gO/Ao9KKX/2qWAOCCEmAu9LKbO78B4JDJRSbvaYYGqePKAECJdSNntyLifz/wjMlFI+5+25gxUhxDbgD1LKH3wtS0/ozv8mkNArBQ8ghLgDeBb4J5AO5AIvAWf7Ui5Nl+gLrOnOG4UQYW6WJSDm1gQJUkp9c+MN6AWYgd93sE8kSmlU2m7PApG2bROBCuAeYBdQBZwDnA5sBGqA+xzGegD4FPgItSJZDoxw2C6BAQ7P3wEeAWKBesBqk9cMZAFjgYXAXtvcLwARtvfOtY1XZ9v/Iru8DuMPAebY3r8GmNxm7heB/9pk/QXo7+QYldnmsst2DHAV8HObz3YzsMk23sNAf5v8+4CP7bLb9j8TtWLbCywAhjuZe4vtuNTb5o60HZuZtuO/Gbiune/gfdu8f2hnzDOAFbbt5cADDtvybJ/letvvoQq4swvf8TbgL8AqoBEIc/Y9ABG2Y3Cr7bkRmA/8w2Gu99vIdbVN5j3AjcCRtrn2Ai84yNEfmA2YgGrgAyDRtu29Nsf0Hofxw2z7dHaMPwbetR2DNcAYh+1/Abbbtm0AJrXzHRwN7ACMDq+dC6yyPR4LLLV9RzuBZ5z8Pibi8JsPtpvPBQi2G3Aq0Gz/oTvZ5yFgEdAbSEOdoB62bZtoe/8/gHDgOmA3MA2IB4YCDUA/2/4PAE3ABbb97+Kg2QWcKAWHuSrayDba9ucJs/1p1wG3O2xvO17rGLb5NwP3oU4+J9r+pIMc5q6x/fnCbCeND50co0NOGLbXruJwpTATSLAdl0bgR6AfSjmvBa607TsKpWSPQp0Ir0SdTCOdzL8NOMnheTFqtRcFjLR9J5PafAfnoFbf0e2MNxE4wrZ9OOqkc06bzzodpayPsI1/kovf8TbUiT4HiHbhexiGOsEPAf6K+i0aHeZqqxResX3uk1G/vS9Qv90+tmN6vG3/AcDvUEo0DXUR8WwHx/SQ79iFY9yAujgyAo8Bi2zbBqGUVpbDuM4uNrYAv3N4/glwr+3xQuBy2+M44GgnY0xEKwV9c/mAwqXAjk722QKc7vD8FGCb7fFE1NWU/U8ab/vjHOWw/zIOnlAesP85bM8NqCvNCbbnXVIK7ch6O/C5w/OOlMIE1JWYwWH7dGxXxba533DYdjqw3sm8h5wwbK9dxeFKYVyb4/IXh+dPYzspAS9jU7wO2zdgO6G1M/82Dp6Uc4AWIN5h+2PAOw7fwdwu/k6eBf7V5rMOdtg+FXjTxe94G3CNw/YOvwfb8zuB9SjlMNDh9Qc4XCn0cdhuAi5yeP4ZDhcNbT7jOcCK9o5p2+/YxWP8g8O2QqDe9ngASjmdhE1RdnDcHwHecvhv1QF9bc/nAg8CqZ2MMZEgVgrap+B+TEBqJ7bdLKDU4Xmp7bXWMaSULbbH9bb7nQ7b61FXMnbK7Q+klFaU+clxPJcRQhQIIb4WQuwQQuxD+UVSXXx7FlBuk8FOKeqK0s4Oh8cHOPRzdIe2x8XZceoL3CmE2Gu/oU5ErhynLKBGSrnf4bW2n6ucDhBCHCWE+EkIsVsIUYsyw7Q9ro5jtP1NdPYdO77Xle/hP6iT8jdSyk0dyY6Lx1gI0VsI8aEQYrvtt/M+XfvtdHaM2/52ooQQYVIFPdyOUhy7bDI4+16nAecJISKB84DlUkr7f/FaoABYL4RYIoQ400XZgwqtFNzPQtQyt6MwxkrUScpOru217pJjfyCEMADZDuMdAGIc9s1weCzbGetl1BXkQCllAsoEIVyUoxLIsclgJxdl6+0q7cnWE8pR0V+JDrcYKeV0F95bCSQLIeIdXmv7uTqTdxrK1JUjpeyFMsm0Pa45Do/b/iY6+o7bzu/K9/AS8DVwihBifCeyu8pjNjmG2347l3HoZ+zoGLlyjJ0ipZwmpRyP+l9J4Akn+61FKZvTgEtQ34t92yYp5RSUaewJ4FMhRKwr8wcTWim4GSllLcof8KIQ4hwhRIwQIlwIcZoQYqptt+nA34QQaUKIVNv+7/dg2tFCiPNsq5PbUbb1RbZtvwKXCCGMQohTgeMd3rcTSBFC9HJ4LR7laDMLIQYDN7WZayfKZt8ev6CW4/fYPvNE4Czgw258pt0ox6SzubrK68CNtit2IYSIFUKc0eYk1C5SynKU3+cxIUSUEGI46qrygy7MH4+6Em4QQoxFnZDa8nfb72Uoyrn7kcO2jr7jtnT4PQghLkf5jq4CbgP+I4To6YoN1Gc0A3uFEH2Au9tsd/rb6ckxFkIMEkKcaLv6b0CtXlo6eMs01Oc+DuVTsI9zmRAizbbC2mt72ek4Njkdb65ePPk1Wil4ACnlM6gchb+hTm7lwB9RDjpQds2lqAiO1ahokkd6MOWXqEigPcDlwHlSyibbtj+hTgh7Uf4OuwxIKdejFNRWm0klC+XEvATlmHydQ09MoJbo/7Htf2Gbz20BJqOuwqpRV6NX2ObpElLKA8CjwHzbXEd3dYw24y1FOe1fQB2nzaiToqtMQZlbKoHPgfullN934f03Aw8JIfajLgI+bmefYptcPwJPSSkdE/U6+o4PoaPvQQiRi/JnXCGlNEspp6F+i//qwmdxxoMoh34tKsJsRpvtj6EuhvYKIe5q5/3dPcaRwOOoz7oDdaV/Xwf7T0f5BWZLKasdXj8VWCOEMAPPARdLKRucjNEHpXwcb/1dkNXv0clrAY4Q4gGU4/cyX8ui6R6dJerp71jjTfRKQaPRaDStaKWg0Wg0mla0+Uij0Wg0reiVgkaj0WhaCbjiWampqTIvL8/XYmg0Gk1AsWzZsmopZVpn+wWcUsjLy2Pp0qW+FkOj0WgCCiFEaed7afORRqPRaBzQSkGj0Wg0rWiloNFoNJpWAs6noNFogpumpiYqKipoaHBWYULTEVFRUWRnZxMeHt6t93tMKQgh3kJ1utolpRzWzvZLUd2SQBXRuklKudJT8mg0msCgoqKC+Ph48vLyCJIac15DSonJZKKiooL8/PxujeFJ89E7qAJTzihBNTgZjmqj+JoHZdFoNAFCQ0MDKSkpWiF0AyEEKSkpPVpleWylIKWcayv05Wz7Aoeni1D14TUajUYrhB7Q02PnL47ma4FZvhZCowl0DliaeW/hNkpNdb4WRROg+FwpCCFOQCmFv3Swz/VCiKVCiKW7d+/2nnAaTYDx/dqd/P3LNRz/5Bwue+MXvlldRVOLtfM3ag7j888/RwjB+vVdbgfiMebOncuoUaMICwvj008/9cgcPlUKtu5KbwBnSylNzvaTUr4mpRwjpRyTltZplrZGE7Ls3t8IwM0T+7N1t5mbP1jOMY/NZur/1lNec8DH0gUW06dPZ/z48Xz4YXcaB7ZPc/Nh7TK6RG5uLu+88w6XXNJe4z734DOlYOsANQO4XEq50VdyaDTBhKnOQrhRcPcpg5j3lxN566oxjMzpxSvFWzjuyZ+48q3FfLtmB8169dAhZrOZ+fPn8+abbx6iFObMmcNxxx3HueeeS2FhITfeeCNWqzqWcXFx3HnnnYwaNYpJkyZht2pMnDiR++67j+OPP57nnnuO0tJSJk2axPDhw5k0aRJlZWUAnH322bz77rsAvPrqq1x66aWHyZWXl8fw4cMxGDx36vZkSKq95V2qEKICuB8IB5BSvoJqSZgCvGRzjDRLKcd4Sh6NJhQwmRtJiY1ECIFRwImD0zlxcDqVe+v5aEk5Hy0p54b3lpGeEMlFY3K4aGwufRKjfS22Ux78ag1rK/e5dczCrATuP2toh/t88cUXnHrqqRQUFJCcnMzy5csZNWoUAIsXL2bt2rX07duXU089lRkzZnDBBRdQV1fHqFGjePrpp3nooYd48MEHeeGFFwDYu3cvxcXFAJx11llcccUVXHnllbz11lvcdtttfPHFF7z22muMGzeO/Px8nn76aRYtctaC27N4TN1IKadIKTOllOFSymwp5ZtSyldsCgEp5R+klElSypG2m1YIGk0PMZktpMRFHPZ6VmI0f/5dAT//5QReu3w0QzIT+PdPm5nwxGyueWcJP6zdSYtV91axM336dC6++GIALr74YqZPn966bezYsfTr1w+j0ciUKVP4+eefATAYDFx00UUAXHbZZa2vA62vAyxcuLDV/HP55Ze37peens5DDz3ECSecwNNPP01ycrJnP6QTdEazRhNEVNdZSImLdLo9zGjg5KEZnDw0g/KaA2r1sLScP7y7lKxeUVx0ZC4XHZlDRq8oL0rtnM6u6D2ByWRi9uzZ/PbbbwghaGlpQQjB1KlTgcNDPp2FgDq+Hhsb63Q+x/1Wr15NSkoKlZWVPfkIPcLn0UcajcZ9mMyNpMYevlJoj5zkGO46ZRAL7j2RVy4bRf/ecfzrh42Me2I2d368MmRXDp9++ilXXHEFpaWlbNu2jfLycvLz81uv6BcvXkxJSQlWq5WPPvqI8ePHA2C1WlsjgqZNm9b6eluOPfbYVj/FBx980Lrf4sWLmTVrFitWrOCpp56ipKTE0x+1XbRS0GiCCGfmo44INxo4dVgm7117FMV3T+Tcoj58tryC9Tvca8sPFKZPn8655557yGvnn38+06ZNA+CYY47h3nvvZdiwYeTn57fuGxsby5o1axg9ejSzZ8/mH//4R7vjP//887z99tsMHz6c9957j+eee47Gxkauu+463nrrLbKysnj66ae55ppraNsuecmSJWRnZ/PJJ59www03MHSo+1dSAdejecyYMVI32dFoDueApZnCf3zLvacN5sbj+3d7nPKaA0yY+hOPnDOMy47u60YJXWPdunUMGTLE6/O6wpw5c3jqqaf4+uuvD9sWFxeH2Wz2gVSH094xFEIsc8V3q1cKGk2QYDJbAEhx0XzkjOykaFLjIlhRttcdYmkCDO1o1miChGqzSlxL7cDR7ApCCEbmJLGifI87xAoqJk6cyMSJE9vd5i+rhJ6iVwoaTZBQU6dWCsk9XCkAFOUmsnV3HbUHmno8liaw0EpBowkSWs1HXXQ0t0dRbiIAv1ZoE1KooZWCRhMkVNcp81FKbM/MRwDDsxMRAlaUaRNSqKGVgkYTJJjMFmIjjERHGHs8VlxkGIPS47WzOQTRSkGjCRJM5sYOs5m7SlFuIr+W7z0sVj5U8MfS2c888wyFhYWtxfRKS0vdPodWChpNkGCq63riWkcU5SRRW99ESXVoNuzxx9LZRUVFLF26lFWrVnHBBRdwzz33uEmyg2iloNEECdVmi1v8CXZG2pzNoWhC8tfS2SeccAIxMTEAHH300VRUVLj9s+s8BY0mSDCZGxmR3ctt4w1IiyM+MowV5Xs4f7SPWqjPuhd2rHbvmBlHwGmPd7hLIJTOfvPNNznttNPccEAORa8UApG95VC+xNdSaPwIq1VS42bzkcEgGJGTGJIrBX8vnf3++++zdOlS7r77bjd94oPolUIg8t3fYNvPcM8WX0ui8RP2NTTRbJUku9F8BMrZ/NKcLdRbWtwS1dRlOrmi9wT+Xjr7hx9+4NFHH6W4uJjISPd+36BXCoGHlFA6Hw5UQ33oXcFp2sdky2ZOdeNKAZRSaLFKVm+vdeu4/ow/l85esWIFN9xwAzNnzqR3795u/+yglULgYdoCdcqBxR7f1FvX+B8Hi+G598pxRLbd2Rw6SWz+XDr77rvvxmw28/vf/56RI0cyefJkt39+bT4KNMoWHnxcUwJZRb6TReM3mGzF8NzpU1DjRdI3JSak/Apz5sw57LXbbrutdVtMTAwfffRRu+99+OGHefjhhzscLy8vj9mzZx/23pUrV7Y+njx5crsn/B9++KEz8XuMXikEGmULIdIWYVKz1beyaPyG6jr31T1qS1FOoq6YGkJopRBolC6AvPEQl67NR5pW7CuF5BgPKIXcJHbua6Sqtt7tYwcaEydObLfBDujS2RpfsH+HUgR9j4GkfKjZ5muJNH6CyWwhKSacMKP7/9Ijc7yfxBaqpTXcQU+PnVYKgYTdn5B7DCTn65WCphVTnXvrHjkyJDOBiDCD15zNUVFRmEwmrRi6gZQSk8lEVFRUt8fQjuZAonQhhMdA5ghI7gcrp0NTPYRH+1oyjY9RJS7cbzoCiAgzcESfXl5bKWRnZ1NRUdFaJkLTNaKiosjO7n4GulYKgUTZQsgeA8ZwZT4C2FMKvQf7Vi6NzzGZGxmckeCx8YtyEnlvUSlNLVbCPWCiciQ8PJz8/HyPzqFxjjYfBQoN+2Dnb8p0BMp8BNqEpAFU8po72nA6oyg3icZmK+ur9ntsDo1/oJVCoFC+GKT1oFKwrxRqtFIIdZpbrOw90OSRcFQ7rRVTdWhq0KOVQqBQtgCEEbKPVM9jklW+gs5VCHlqDthzFDzjaAbI6hVF7/jIkEpiC1W0UggUyhYpB3NknHouBCTnafORprXERaoHzUdCCIpyE0Oq3EWo4jGlIIR4SwixSwjxm5PtQgjxvBBisxBilRBilKdkCXiaG6Fi6UHTkZ2kfG0+0hyse+TBlQIov8I20wFqbNnTmuDEkyuFd4BTO9h+GjDQdrseeNmDsgQ2lSugpVElrTmSnA97y8Da4hu5NH6Bqc4zdY/aYk9iW1muTUjBjMeUgpRyLlDTwS5nA+9KxSIgUQiR6Sl5AhrHpDVHkvLB2gS17m/JpwkcqlvNR55dKQzP7oVBhFbF1FDElz6FPkC5w/MK22uatpQuhJSBEJt66OvJ/dS9djaHNCZzI2EGQUK0Z9OOYiLCGJyRwAq9UghqfKkU2mtX1G5euxDieiHEUiHE0pDLcrRaoXzR4aYj0LkKGkD5FFLiIpx2AHMnRbmJ/Fq2F6tVl6AIVnypFCqAHIfn2UC7PeiklK9JKcdIKcekpaV5RTi/Yfc6aKiF3GMP3xafBcZI7WwOcUx1jW5vw+mMotwk9jc2s7U6OCqCag7Hl0phJnCFLQrpaKBWSlnlQ3n8k9IF6j736MO3GQyQ1FevFEKcarPF7W04nWF3Ni/X+QpBiydDUqcDC4FBQogKIcS1QogbhRA32nb5BtgKbAZeB272lCwBTdlCiM+EpLz2tyf30yuFEKemznPF8NrSLzWWhKgwncQWxHjMMyWlnNLJdgnc4qn5gwIplZM59xiVrNYeSflQMk/t6wWbssb/MJk9Vza7LQaDYGRuko5ACmJ0RrM/s7cM9ldC33b8CXaS86GpDupCzAGvAaDe0kKdpcXjOQqOFOUksnHnfuoam702p8Z7aKXgz7TmJ7TjT7CjC+OFNPbENU/nKDgyMjcRq4RVFbVem1PjPbRS8GdKF6iid70Lne9jD0vVuQohycESF95bKYzM1hVTgxmtFPyZskWQexQYjM73ScwFYdARSCHKwRIX3lspJMVG0C81VjubgxStFPyVOhNUb+jYdAQQFgkJ2dp8FKLYS1x4K/rIzsjcRFaU7dV9lIMQrRT8lVZ/QgdOZju6hHbI4gvzEShnc7W5kYo99V6dV+N5tFLwV8oWqmzlPi5UFNcltEMWk7mR6HAjMRHebbdelJsEwK+6DlLQoZWCv1K2UCmEMBdsxcn94EC16uOsCSlMdRavrxIABmXEExVu0H6FIEQrBX/EUgdVKw8vle0MXRgvZFFKwXtOZjvhRgPD+yTqCKQgRCsFf6RiKVibO05ac0TnKoQsJnOjR9twdkRRbiJrtu+jsVk3eQomtFLwR8oWAgKyj3Rtf71SCFnsZbN9wcicRCwtVtZV7ffJ/BrPoJWCP1K6ANKHQXSia/tHxkNsmk5gCzGklJjqvFf3qC12Z7OugxRcdBiyIISIAs4EJgBZQD3wG/BfKeUaz4sXgrQ0KfNR0aVde5+OQAo59jU009QivZ6jYCejVxSZvaJYUbaXq8f5RASNB3CqFIQQDwBnAXOAX4BdQBRQADxuUxh3SilXeV7MEGLHKlXgzlUns53k/IO9FzQhgclsq3vko5UCKL+CdjYHFx2tFJZIKR9wsu0ZIURvINf9IoU4pfaktS4qhaR8WPUxNDe6FsaqCXhMdb5JXHNkZE4i36zeQbW50afKSeM+nPoUpJT/bfuaEMIghEiwbd8lpVzqSeFCkrKFqqFOQmbX3pecD0jYU+oJqTR+iH2lkOwj8xE4JLHpfIWgoVNHsxBimhAiQQgRC6wFNggh7va8aCGIlEopuFLaoi3J/dS9jkAKGex1j3x5hT4sqxdhBqFNSEGEK9FHhVLKfcA5qBaaucDlHpUqVKneBAdMnRfBaw+dqxBy2OseJcX4bqUQHWFkSGaCzmwOIlxRCuFCiHCUUvhSStkE6NKInqDM5ih2NWnNkdhUiIjTK4UQoqaukV7R4USE+TayvCg3kZXle2mx6tOCx5ASProcVn3i8alc+TW9CmwDYoG5Qoi+gC6y4wnKFkFMKqQM6Pp7hdBhqSFGtY/qHrVlZE4idZYWNu8y+1qU4GXLj7BuJlg8nyjYqVKQUj4vpewjpTxdquLpZcAJHpcsFCldoExHQnTv/cn5OoEthFAlLnwf8aOT2DyMlDDnCdU3ZeRlHp/OqVIQQlwmhDhsu1Q0CyH6CyHGe1a8EGJfJewt7Z7pyE5yvhrDqmvRhAK+LHHhSF5KDIkx4dqv4Cm2zoGKxTDhzxDm+e+7ozyFFGCFEGIZsAzYjUpeGwAcD1QD93pcwlChrJv5CY4k5UOLRSmYxBz3yKXxW0x1Fo7yA6UghKAoRyexeQQpofgJSOgDRd6J7+koT+E5YBQwHUgDJtmebwcul1KeL6Xc5BUpQ4HShRAeCxnDuz+GLowXMjS3WNlzwEKKH5iPAEbmJLFpl5n9DU2+FiW4KJmrLhjH/9lrSakd1j6SUrYA39tuGk9SthByjgRjDzpotYalboX849wjl8Yv2XOgCSkh1Q9WCqAikKSEVRW1jBuQ6mtxgofiJyA+02urBNBVUv2D+r2wc033ktYc6ZUNhnAdgRQCmOrs2cz+sVIYkaMq+mpnsxspmQel82Hc7RAe5bVptVLwB8oXAxL69sCfAGAwQlJfbT4KAeyJa/7gaAboFR3OgN5x2tnsToqfgLh0GH2lV6fVSsEfKFsAhjDoM6bnY+lchZCgurVCqn8oBcDmbN6LilzX9Iht82HbPNsqIdqrU7tS+yhdCPGmEGKW7XmhEOJaz4sWQpQuhMyREBHT87GSbUpB/zGDmhp7hVQ/MR8BjMxNpKbOQlnNAV+LEvgUPwGxvWHM1V6f2pWVwjvAt6gmOwAbgdtdGVwIcaoQYoMQYrMQ4rDwVSFErhDiJyHECiHEKiHE6a4KHjQ0NUDl8p6bjuwk91NZjwdM7hlP45eYzBaMBkGv6HBfi9JKUY6tYmq5NiH1iLJFUFIM4/7k9VUCuKYUUqWUHwNWACllM9BpdpQQwgi8CJwGFAJThBCFbXb7G/CxlLIIuBh4qQuyBweVy1VuQU/yExzRhfFCAlNdI8mxERgM3cx+9wAF6XHERBi1X6GnzHlctdcdc41PpndFKdQJIVKwFcETQhwN1LrwvrHAZinlVimlBfgQOLvNPhJIsD3uBVS6JHUwYe+W5i6loHMVQoJqs8VnbTidEWY0MDy7l45A6gnli2HrT3Dsbe4xJ3cDV5TCHcBMoL8QYj7wLnCrC+/rA5Q7PK+wvebIA8BlQogKVFluV8YNLsoWQdpgiEl2z3iJfQGhVwpBjslPO52NzElibdU+Gpp0qZVuMedxiEmBI33ntnWlIN5yVFmLY4EbgKEu9mVub13b1vs5BXhHSpkNnA681169JSHE9UKIpUKIpbt373Zh6gDB2gLlv7hvlQAqnjmhjy6MF+SY/KRCaluKchNpapGsqdSFlLtMxVJVDfXYWyEi1mdiuBJ9ZESdsCcBJwO3CiHucGHsCsCxAE82h5uHrgU+BpBSLkTVVjosHVJK+ZqUcoyUckxaWpoLUwcIO9dA4z73KgVQJiRtPgpqTGb/KXHhSJFOYus+cx6H6GQ48jqfiuGK+egr4CpUgbx4h1tnLAEGCiHyhRARKEfyzDb7lKGUDUKIISilEERLgU4oW6Tu3RV5ZCcpT5uPgpiGphbMjc1+uVLonRBFn8RoVugIpK6xfRls/h6O/SNExvlUFFcK7WRLKbtcpc1WXvuPqHBWI/CWlHKNEOIhYKmUciZwJ/C6EOLPKNPSVTKUMl/KFqga6Ym57h03OR/qdkGj2ec/MI37MbXmKPifUgBlQtIRSF2keCpEJ8HY630tiUtKYZYQ4mQp5XddHVxK+Q3Kgez42j8cHq8FxnV13KBASpW0lueBlhRJDhFIGUe4f3yNTzHZsplT/NDRDKoT29erqti1r4HeCd6r2ROwVK6Ajf+DE/8Gka4YYTyLK+ajRcDnQoh6IcQ+IcR+IYT2IvWUPSVg3uF+0xGoBDbQJqQgpXWl4IfmI3DoxKZNSK5RPBWiEmHsDb6WBHBNKTwNHAPESCkTpJTxUsqEzt6k6QS7P6GnlVHbQ+cqBDX2Ynj+0IqzPYZmJRBuFNqE5ApVK2HDN3DMLRDlH6dVV8xHm4DfQsrW7w1KF6irg7TB7h87qpeKYtArhaDkoPnIP1cKUeFGCrOCJInNvBsWPK/yiI65tWf9TtqjeKr6vx7lH6sEcE0pVAFzbAXxGu0vSimf8ZhUoUDZQsg9GgweKlSbnK9zFYIUU52FqHADMRFGX4vilKKcRD5aUk6LVWL0o1IcLtNQCwtegEUvgaUOkLAKJgoCAAAgAElEQVThf3D+6+4LDNmxGtZ/DRP/TykGP8GVM1IJ8CMQQddCUjXOMO8G02b35yc4ktxPm4+ClGpzIymxkQjhvyfbwswE6ptaKA+0iqlN9bDg3/DcCJg7FQacBH9cAue9rvKKXh4Pv33mnrmKp0Jkgl+tEsCFlYKU8kFvCBJSlC1U955UCkn56sfbbIEw/zQzaLqHyWzxqz4K7VGQoa4bN+zcT16q77JzXaalCVa8r07U+yuh/ySY9HfIKlLbUwdC9pEw4zr49BrY/COc9kT3o4V2roF1M+G4e1Qoqh/hVCkIIZ6VUt4uhPiKw8tTIKWc7FHJgpmyhRAWdfAH5wmS80FaobYcUvp7bh6N1zHVNZLmp+Godgb2VvkxG3fs55ShGd0bZN7TsPozGHgSFJwGOWNVd0F3YrXCmhnw0z+hZgtkj1UmovZCxZPz4epZSnHMe0r5Bc9/E7JHd33e4qkQEQ9H39Tzz+BmOlopvGe7f8obgoQUZQtVlzVPXsE7ltDWSiGoMJktDM7wj0gVZ8RGhpGTHM3GXebuDVCxDGY/Ar1yYOGLMP85VShu4Ckw6FTof2LPYvqlhM0/wI8PKtt+70K4eDoMOg06MssZw+HEv0K/iTDjenjrZDjhr6r3gasKa9c6WPslTLjTfYUw3UhHSuFWVIZxsbeECQkazVC1Cia4Uj6qB9jDUrWzOaiQUqq6R35uPgIYlB7Pxh37u/7G5kb48maIz4Qb56nXNv8IG2ap8M2V08AYAfnHQcGp6kTeK9v18csWwQ8PqooCiX3h3NfgiAu6tgrJGwc3/Qxf3a4Uy5bZcO6r0KttIeh2KJ6qCt4dc4vr83mRjpRCl0tbaFygYjHIFhV55Eni0iE8Rjubg4z9jc1YWqx+m6PgyMD0eOZs2I2l2UpEWBei7IqfgN3r4dLPDkblDDtP3VqaoXzRQQXxzV3qljFcKYdBp6nWtu1d7e9YDT8+DJu+Vf+P05+CUVd2f8UenQS/fwd+/QC+uQdePhYm/xsKO7Cs71oPaz6H8X/2y1UCdKwUYoQQRbRfAtteUlvTVUoXgjAo26UnEUKZkHSuQlBRY/bvbGZHBqXH02yVbDPVUZDuoqmncgX8/CyMvEz5EtpiDFP2/rzxcPIjUL0JNs5SSmLuk0qhxGdBwSkw6HS1mti3XfkMfvtUKZlJ96uIH3eUpxYCii6DnKPhs2vh48th9FVwyj/bH3/uk+pi7Zg/9nxuD9GRUuiDymZ21hfhRI9IFOyUzFUOZm9kLybnq9BXTdBgqvPvukeO2BXBhh37XVMKzRb44maI6w2nPNr5/kJAWoG6jfsT1Jlg03dqBbH6E1j2tjoBt1iUuWn8HTDuNs9E+6QOgGu/h58egfnPw7b5cMGbkDni4D67N6qIwHF/gtgU98vgJjpSCpullPrE704azbB9qWqi4Q2S8mDT9yrCwlNJchqvUm327wqpjvRLi8VoEGzc6aJfYe6TsGstXPIxRCd2fcLYFBg5Rd2aG2HbPJVwFh6lspHj07s+ZlcIi4DfPaSc4DNugNcnwUkPwNE3q//fvKcgPNp7//9u4uacbU2HlC0EazPkH++d+ZL7QUsj7K9yzQGm8Xta6x4FwEohKtxIXkoMG1xxNletVCGoI6Yo009PCYtUiWcD2jFBeZp+E+GmBTDzVvjur6qb2oS71OrlmFsg9rA+Yn5FR5ePf/GaFKFCSbFaxuYc5Z35dGG8oMNe9yg5AFYKAIMy4jtfKTRb4AvbyfKUf3pHME8TmwIXfwBnPKP8iO+cDsZIOPY2X0vWKU6VQnf6J2g6YWuxcjBHxHhnPsdcBU1QYKqzkBAV1rVoHh9SkB5Pac0B6i0tznf6+RnYuRrOfNZvI3K6hRBw5LVw/RxVDXniX5S/xM8JjF9WMHCgRoXE9fOS6QhU4o8hTOcqBBHV5saAcDLbGZQej5SwZbeTJLYdq5Uv4Yjfw+DTvSuct+g9GK6ZpcJQA4AOlYIQwiiEeNJbwgQ12+YBUoXIeQtjmFIM2nwUNJjMloBwMttprYHUnl+hpUlFG0UnwWlTvSyZxhkdKgUpZQswWvhzOcZAoWQuhMdCn27USekJyf20+SiIMNU1BkSOgp2+yTFEGA3t+xV+fhZ2rIIz/xVcZqMAx5XooxXAl0KIT4A6+4tSyhkekyoY2VoMfY9VtVO8SXK+CoPVBAUms4UxeYFzAg0zGujfO44NbZXCzrUq0WzoeTDkLN8Ip2kXV5RCMmDi0GQ1CWil4Cr7KsG0CUZf6f25k/JVw5ADNfpqLMBpsUr2HLCQGkDmI4BB6XEsLqk5+EJLM3xxk8ouPl1bp/0NV/opXO0NQYKaEltRL2/6E+wkO0QgaaUQ0Ow9YMEqAyOb2ZGCjHi++LWSfQ1NJESFw4LnoOpXVTfIz2P2Q5FOo4+EEAVCiB+FEL/Zng8XQvzN86IFESXFypmWfoT3507up+61szngMdUFTt0jRwbZSlxs2rlfFYSb8zgUng1Dz/WxZJr2cCUk9XXg/4AmACnlKuBiTwoVVEipnMx5E3xTaiIpT91rZ3PAU21LXEsJgAqpjtjrHm2s3KtKYkfEwelP+1gqjTNcOUvFSCkXt3mt2RPCBCU1W1X3M1+YjkDVWonP1CuFIOBgiYvAWin0SYwmNsJI6urXYfsy5UeIS/O1WBonuKIUqoUQ/bG15BRCXABUeVSqYKJkrrrvN9F3MiTl6wQ2D7OjtoFxj89mWWlN5zt3k0ArcWHHYBBMTNnL8ZWvw+AzYdj5vhZJ0wGuKIVbgFeBwUKI7cDtgP81FvVXSorVlXrKAN/JoHMVPM5nyyvYvreeuRurPTaHqc6CQUBiTGApBawt3NPwPAeIVLWAdNqTX9OpUpBSbpVSngSkAYOllOOllNs8LlkwYLWqyKP84337R0jOA/MOsBzwnQxBjJSSz5ZVALCuap/H5qk2W0iOjcBoCLCT6qKX6Fu/hvstV1AtulESW+NVXIk+ahFCPA4ckFLut72mu665wu51cKDad/4EO/bCeHu2+VSMYGVF+V62VtcRG2FkrQeVgsncGHBOZqo3w+xHMGVP4kvrONd7K2h8hivmozW2/b4TQtgD3V26VBFCnCqE2CCE2CyEuNfJPhcKIdYKIdYIIaa5JnaAsLVY3ftaKbTmKmi/gieYsbyCqHADVx6bR8Weemrrmzwyj6nOEljhqNYW+PIWCItEnv4vQLDRld4KGp/iilJollLegwpNnSeEGI3N6dwRQggj8CJwGlAITBFCFLbZZyAq3HWclHIoyl8RPJTMVfb8xBzfypGk+yp4isbmFr5aWcUpQzM4Ml9dM6330Gqhps4SWIlri1+D8kVw6hOkZOaSGBPOhp1OqqVq/AZXlIIAkFJ+DFwIvA30c+F9Y1EtPbdKKS3Ah8DZbfa5DnhRSrnHNscuVwX3e1qaoXS+71cJoDKZoxK1s9kD/LhuF7X1TZw/Kpuhmarvtqf8CtXmxsCpkFq9CX54EAaeDCMuRghBQboLDXc0PscVpfAH+wMp5RpgPOBK+6A+QLnD8wrba44UAAVCiPlCiEVCiFNdGDcwqPoVGvd5r/VmZyTn65WCB5ixvIL0hEjGDUglLT6SlNgIj/gVGptb2N/QHBg5Cmu+gDdOUi0xz3y2NchiUHo8G3fsR8pODQ0aH+K09pEQ4kQp5WygrxCib5vNrqwB2/M7tP01hAEDgYlANso8NUxKubeNLNcD1wPk5ua6MLUfsHWOuveHlQIoE1LlCl9LEVRUmxuZs2E3107Ib40IKsxKYF2V+6+Ga1pLXPix+ajRDP/7C6x4H7JGwflvHNIbvCAjnv2NzVTVNpCVGO1DQTUd0dFKwX6Je1Y7tzNdGLsCcDSmZwOV7ezzpZSySUpZAmxAKYlDkFK+JqUcI6Uck5YWIJmQJXMhfZj/FPxKzoe9ZaqxicYtfPlrJc1WyQWjsltfG5KZwIad+2lusbp1Lns2s9+aj7Yvg1ePgxUfwIQ74drvIKX/IbvYayAdVkZb41c4XSlIKe+33Xe3SuoSYKAQIh/YjqqXdEmbfb4ApgDvCCFSUeakwA+RaWqA8l9gzDW+luQgyf1AtqiSG8muuIQ0nfHZsgqGZ/dioO1kBzAkMx5Ls5Wt1XWtNX/cQWvdo87MR7UVsOpjOOICSPTCqtraAvOfg58ehbh0uOpryBvf7q4F6XGAKox3wiD/71UcqnRaOlsIkQhcAeQ57i+l7NCvIKVsFkL8EfgWMAJvSSnXCCEeApZKKWfatp0shFgLtAB3SylN3f0wfkPFYmhu8B/TERyMQKop0UrBDayr2sfaqn08cNYhAXUUZvYCYG3lPrcqhYMrhQ7MR5t/gM+ug/oadZIecTGMv+OwK3a3UbsdPr9BtZotPAfOelZVA3ZCYkwE6QmRbNihI5D8GVea7HwDLAJWA11aE0spv7G93/G1fzg8lsAdtlvwUDIXhBH6jvO1JAdJ1mGp7mTG8grCjYLJIw+NneiXFkuE0cC6qn2cU9Q2rqL7mOo6WClYW1QXs+Kp0LsQLnoP1s6E5f+BX6fBsAvguLsgbZDb5GHtlzDzNmWOPPtFGHmpS1n7OgLJ/3FFKURJKYPrpO1pthZDVhFEJfhakoPEZUBYlA5LdQPNLVY+X1HJCYN6H1acLtxooCAjzu0RSCazhYgwA3GRbf6yddUw4zrYMhtGTFG1hSJilAlnwp2w8N+w5E1Y/YnqYXDcXZDRg74ejWb4372w4r2DzuQurEQK0uP54JdSWqwy8Mp1hAiuhKS+J4S4TgiRKYRItt88Llmg0rhfOd36+Ukoqh2DQfVW0Eqhx8zbXE21uZHzR2e3u31IRgJrK/e5NfSy2qzacArHq/Hyxcq5u20+nPUcnPOyUgh24tPh5Efg9t9gwh2w+Ud4ZTxMn6J+o11l+3KbM/l9ZZZqx5ncGYPS42loslJeo+tw+SuuKAUL8CSwEFhmu+lO8M4oXaAcuv7kT7CT3E+bj9zAZ8sqSIoJd+osLcxKwFRnYff+RrfNWVPXeDAcVUpY9Aq8fRoYwtTJefRVzs03sSkw6R/w59Uw8T71G339RHjvPChb1Pnk1hb4+V/w5u+Ur+zKr+Ck+8EY3uXPUZChI5D8HVeUwh3AACllnpQy33bTnkpnlMwFYyTkHOVrSQ4nKV+tFHTyULeprW/iu7U7mTwii4iw9v8+Q2yZze40IbXWPWrcD59erfIBBp4MNxRD1kjXBolOgol/gdtXw6T7VYLlW6fAO2cqk2d7v4va7fDu2fDDAzD4DLjxZ8if0O3PMbC3ikDSNZD8F1cL4um1nqtsLYacsarjmb+RnA/N9bB/h68lCVi+WV2Fpdnq1HQEB5WCO5PYTGYLhcZKeO0E5eQ96UG46IMOo32cEpWgzEm3r4ZT/qlKUrw7WSmITT8cVA5rv4SXj1Vmo7NfhN//R5VM6QGxkWHkJEfrlYIf44qjuQX4VQjxE9C6Hu4sJDUkqTPBztVwwt98LUn7OBbGS8j0rSwBymfLKhjQO44j+vRyuk+v6HD6JEa7baUgpeTYuh+4fdubEJ2gzDdOcgG6REQsHHMLjLlWOY5/fhY+OF8FSaQMUM7prCI4/023hrUO0hFIfo0rK4UvgEeBBRz0KXTDSxUCbJun7v3NyWwn2SFXQdNltlXXsbR0D+ePyj7U4dsOqtyFG5RCUwPNM2/nSeOLVMcXwo3z3KMQHAmPgrHXwW0r4KznoX4PrP7U5kz+3u15DgXp8WzdXYel2b1Z3xr30OlKQUr5H28IEhSUFENEnLq68kcSc1X+hHY2d4sZK7YjBJzrQv7BkMwEfly3k4amFqLCjd2bcE8pfHwF4VW/8krzWaSNf4Tz4zO6N5YrhEXA6CtVzkHjvh6bipwxKCOeZqtkm8m9Wd8a9+B0pSCE+Nh2v1oIsartzXsiBhAlc1XCWjeiMryCMRx6ZeuVQjewWiUzllcwfkAqGb2iOt2/MDMBq4QN3XWobvxWhX/WlLD5xNd4vHkKyQkxnb/PHRjDPKYQgFZF0O1jo/EoHa0U/mS7d6X4naZ2O5g2w+juloryEsn5ugNbN1i8rYaKPfXcdbJrWcGFDhFII3K60JfY2qJKVMx7WiWZXfguW6tigGWkBlorTif0S4vFaBDar+CnOF0pSCmrbA9vllKWOt6Am70jXgBRMlfd+6s/wY7OVegWny2rIC4yjFOGuma+yU6KJj4yrGt+hQM18N45SiGMukLZ85P7YWotm+2nFVK7SGSYkbyUGL1S8FNccTT/rp3XTnO3IAFPSTFEJ0Pvob6WpGOS8pUjsX5v5/tqAKi3tPDN6ipOPyKD6AjX/AMGg2BwZjxrK11UCuZdKl+gbBGc/RJM/ndrWLPJViG1bUmNQGZQho5A8lc68incJIRYDQxq408oAbRPwREp1Uohf4IqJ+HP6MJ4XebbNTuos7Rw3ijnuQntUZiZwPod+7FaO0kWrN2uspP3lMAlH0PRpYdsNtVZiI8M677D2g8pSI+ntOYA9ZYWX4uiaUNHZ7BpqIY6Mzm0wc5oKeVlXpAtcDBtgX3b/af1Zke0ltDWfgVX+Wx5BdlJ0YzN65rzdUhmAubGZsr3dJD7WbMV3j5VrRQu/xz6n3DYLiazJWhMR3YGpccjJWzepcto+xsd+RRqpZTbpJRTgBTgbGAyqq+CxpGSYnUfEEohT93rCCSX2FHbwM+bqzlvVDaGLlb1LMyyZzY7MSHt3gBvn65KV1w5E3KPbnc3k2PdoyBB10DyXzq1dQgh/g78B6UYUoG3hRB+mrLrI0rmQkIfzzUzcSeRcapDljYfucTnK7YjJZzXjd4IBenxGATt+xWqVimTkbTCVd90mNtiMlv8tw1nN+mbHENEmEH7FfwQV8pcXAIUSSkbAIQQjwPLgUc8KVjAYLUqpVBwiktNRvyCpHyo2eZrKfweKSWfLa9gTN8k8lJju/z+qHAj/dPiWNu2BlL5ElVOIiJerRA6uZioNlsoyu1CWGsAEGY0MCAtTisFP8QVr+g2wDFbJxLY4hFpApFda1T7w0AwHdnRuQousXp7LZt3mTssftcZQzLblLsomaeqjsakwDWzOlUIVqtUZbODJEfBkUEZ8bpaqh/iilJoBNYIId4RQrwN/AaYhRDPCyGe96x4AcBWuz+h++WEvU5SPuyvhKZ6X0vi13y2rIKIMAOnH9H94oGFWQls31tP7YEm2PgdfHCBKjdy9Sx13wl765uwyuDJUXBkYHoclbUN7Gto8rUoGgdcMR99brvZmeMZUQKUkrmQ3F+VjwgUkm3tMPaUQu/BvpXFT7E0W5m5spKTC9PpFd39siX2Mto7Fn1Er59vh/RCuOxz1fjGBew5CsHmaAYVgQSwaed+RvfVzRz9BVeUwkfAAEACW+y+BQ2qaXnpfBh+oa8l6RqOuQpaKbTLTxt2sedAU49MR6ByFc4zzGXg3Ncg50i49BOIcl52uy3VZpXNnBpkjmZwrIFk1krBj+goeS1MCDEVqEBFH70PlAshpgoh/LTim5epXAEWs3+23uyIJF1CuzM+W1ZBWnwkEwak9mictPXv80zEK2yOLVJ5CF1QCKDCUSE4Vwp9EqOJjTBqZ7Of0ZFP4UkgGciXUo6WUhYB/YFE4ClvCOf32P0JeQGmFGKSITJBO5udUFNn4acNuzhnZBZhxh5kqC/4N/z3DlZEHc09EX9VTW26IQsEp0/BYBAMTI/XNZD8jI7MR2cCBVIebNwqpdwnhLgJWM/BKqqhS0kxpB/hsn3YbxBCJbHpXAWFpQ5mP6oeZw5n7o4UrC3N3TcdSQnFT8Ccx2DouXwXcydrF1TQ1GIlvItKptpsUV9XTPApBVB+hR/W7fS1GBoHOlIK0lEhOLzYIoTQnd+b6qF8sepYFYhkjYRVn8D+nRCf7mtpfEdLM3x6jepfEBYFzfWcA5wWFUHkV0dAxnDIHA4ZI5STuLPe21LCd3+DhS/AyMtg8vMMXrUDS0sZW3abGZyR0CXxTOZGkmMiMHYxmzpQKMiI56Ol5VSbG0kNQhNZINLRZctaIcQVbV8UQlyGWimENuW/QEtj4PkT7Iy7HVos8PMzPhOhYs8Brn1nCe8tKvWNAFLCN3fBxv/BGU/B/1VQetFsbrfczMbciyA8BtbMgK//DG+cCP/sAy8dAzNugIUvqpwDx2qzViv89w6lEMZeryqdGowHeyu4WjHVAZPZElTVUdtij0DSfgX/oaOVwi3ADCHENaiezBI4EogGzvWCbP5NyVzV2jL3GF9L0j1S+sPIS2DpW3DsrV4Pqf1x3U7u+HgltfVN/FZZyyVjc71/NTzvaVj2Noz/Mxz5BwCmbYvhaybwtwsnQVykUhx7S1VZih2r1P3WObDqw4PjJOWpFUVzA2z6To036f7WDPf81Fgiwwzd6tms6h4Fr1IoyIgDYOOO/Rzbv2dOfY17cKoUpJTbgaOEECcCQwEBzJJS/ugt4fyarcXQZzREdc0c4Fccfw+s/BCKp8Jk7+QhNrVYefLbDbw2dytDsxK44fh+TP3fBhZuMTF+oBdPCr9Oh9kPw/CL1AkcaLFKvlixnYmD0g6aMuz+l6Q8KJx88P3mXTZFsRKqVqrH+7bDiX+H4+46ZKowo4FBGfGs7Y5SMFsYkhXAv7FOSIuLJDEmnA07dbVUf6HTPAUp5WxgthdkCRwaaqFyOYy/w9eS9IzEXBhzNSx5E8b9yeMF/Sr31vPHactZXraXy4/uy1/PGALAyz9tYcaKCu8phS2zYeYfVWmSyS+0XtHP31zNzn2N3H+WC6umuN4w8CR1s2NtAUP7PQ8KMxP4bu1OpJSILtTIqjY3BmWOgh0hBAXpuuGOP+HRjjBCiFOFEBuEEJuFEPd2sN8FQggphBjjSXncRukCVd3S31tvusKEO8EYoaJlPMhP63dx+vPz2LjTzAuXFPHwOcOICjcSFW7k9CMy+d9vOzhgafaoDIC6ov/oCkgbDBe9B2EHT7ifLa+gV3Q4k4b07t7YThQCqMzmmjoLO/c1ujycpdnKvobmoMxRcGRQuqqB1E5ci8YHeEwpCCGMwIuo1p2FwBQhRGE7+8UDtwG/eEoWt1MyV0WqZI/1tSQ9Jz5DRVCt+hh2uT9+oKnFyuOz1nP1O0vI6hXNV7eO58zhWYfsc+6oPhywtPDdGg+HJu4tU7WHonodllm8v6GJb9fs4KwRmUSGub/DWae9FdohmHMUHCnIiGd/YzNVtbpYgj/gyZXCWGCzlHKrlNICfIhq1NOWh4GpQOD8IkrmQs5REB7V+b6BwLjbISIO5vzTrcNW1dYz5bVFvFK8hcuOzmXGzceS304J6rF5yfRJjGbGiu1unf8QDtTA+xdAUwNc9ikkHKqYZq3eQUOTtcstN11lsK2pTFf8CtX2ukdBWCHVEXsEkm644x94Uin0AcodnlfYXmtFCFEE5Egpv/agHO7FvBt2/ha4oajtEZsCx9wMa79UTlM38NOGXZz+3DzWVe3j+SlFPHLOEU57DBsMgnOKsvh502527fPAtUFTA3x4qUrWmzINeg85bJdPl1fQLzWWohzP9C2IjwonNzmmS0rBvlJIDfaVQvrBCCSN7/GkUmjPm9ZqNBRCGIB/AXd2OpAQ1wshlgohlu7evduNInaDbfPUfSD1T3CFo2+GqMSDmb3dpLnFyhP/W8/Vby8hPSGKr24dz+QRWZ2+79yibKwSZq6s7NH8h2G1wufXQ9kCOPcVyBt/2C5bd5tZXFLD+aOzu+QE7ipDMuNZ14VchWCue+RIYkwE6QmReqXgJ3hSKVQAOQ7PswHHf3w8MAyYI4TYBhwNzGzP2SylfE1KOUZKOSYtLc2DIrtASbHqmNVB+8SAJDoRxt0Gm75VmdrdYEdtA1NeX8TLc7YwZWwuX9wyjn5pcS69d0DvOIZn9+Jzd5uQvvurWgGd/AgMO7/dXd5ZsI0Io4ELx+S0u91dDMlMoMRU57JD3WQODZ8CoCOQ/AhPKoUlwEAhRL4QIgK4GJhp3yilrJVSpkop86SUecAiYLKUcqkHZeo5JXMhbxwYXak6HmAcdSPEpqn4/S5SvHE3pz8/jzWV+3ju4pE8dp5zc5Ezzi3qw5rKfe47OSx8ERa9BEfdBMf8sd1dauub+HRZBWeNyCIt3rNX5IWZCUgJ6100k1SbLYQbBfGRQfhba8Og9Hg27zLTYtURSL7GY0pBStkM/BH4FlgHfCylXCOEeEgIMbnjd/spe8tVZdFgMx3ZiYhVuRclcw9WgO2E5hYrT367nivfWkzv+Ei+unU8Z4/sepN7gLNGZGE0CGYsd8Nq4bcZ8O19MGQynPKo0/7Znywt54ClhavH5fV8zk6wN9xxNQLJZFZtOD1p0vIXCjLiaWiyUl5zwNeihDwezVOQUn4jpSyQUvaXUj5qe+0fUsqZ7ew70e9XCVtsOXzB5GRuy5hrID4LZj+iSjx0wM59DVzyxi+8+NMWLj4yh89vHkd/F81F7ZEaF8nxBWl8+et2rD25Ytz2M3x+gypBct7rTvMHWqySdxZsY2xeMsP6dK3PQXfIToomPirM5RpIpjpLSJiOwKHhjjYh+RyPKoWgY9k7kDIQ0of6WhLPER4Fx98NFYth0/dOd5NScvXbS1hdUcu/LhrB4+cPJzqi5/H95xb1oaq2gUVbTd0bYNc6+PASVZbi4mkdhg3/sG4nFXvqvbJKAJW9OyQzoWsrhSB3MtsZ2FtHIPkLWim4SsUyVdpi7PVOTRFBQ9Hl6qQ6+2EVvdMOczdVs7ZqHw+fM4xzi9wX2/+7wnTiI8O6l7Owr1LlIoRFwaWfqmZCHfD2/BL6JEbzu0LvlQ4vzExg/Y79Lq2Eqs2WoC5x4UhsZBg5ydF6peAHaKXgKotfUwleIy72tSSexxgOx9+rqoKu/6rdXd6Yt5X0hEiXwk27QlS4kdOOyGDW6irqLS2uv7FhH3zwe2jYq8OEO24AAB1dSURBVLKVk/p2uPvayn0s2lrDFcf07Vl3tS5SmJnAAUsLpZ3YzqWUQV8htS2DdASSX6CVgiuYd6u6+iMvCeyqqF1h+IWQWqDyFqyHnpzXVu5j3qZqrjo2n4gw9/+Ezi3Kps7Swndrd7j2hmYLfHQZ7F4PF74LmSM6fcvb80uIDjdy8ZG5PZS2a9jLXXTmVzhgaaGhyRoy5iNQfoWtu+uwNLe/OtV4B60UXGH5f1RDmiMDtMtadzAY4YT7oHoDrP70kE1vzNtKTISRS8Z65oR6VH4yWb2iOs9ZkBJ2rIZPr1b5I2c9DwMmdTq+ydzIlysrOX90H3rFhLtJatcY0DsOo0F06ldorXsUIuYjgEEZ8TRbJSXVdW4b8z8LtnH/l7+1lgzRdE7wB0D3lJZm1Yim30RIK/C1NN5lyNmqB/Wcx2DYeWAMp6q2npkrK7n8mL4eO6EaDIKzi/rw2tyt7N7feHj+wJ5SWP2Juu1eD4YwOOlBKLrUpfGn/VKGpdnKVcfme0D6jokKNzIgLa5TpWA/iYVSi0rHCKRBtlpRPeHbNTu4f+YaAGas2M6fTyrg8mP6drlPdqihj05nbPhGNU8Ze72vJfE+BgOc+DdVM+jXDwCV/WuVkmvGefaEel5RH1qs8mDZizoTLH4d3jwFnhuunOBRiXD6U3DnRhh/u0vjWpqtvLeolOMK0hjQu/vhsz1hSGbnDXfs2czB3IqzLf3SYjEaBJvc4FfYvMvMnR+vZER2L/5723hG5iTy0NdrOeP5eSzYUu0GaYMXvVLojMWvQa8cKDjV15L4hoJToM8YKJ6KefAFTPuljNOOyCQnOcaj0w5Mj2d0VgQ1i96H0tWw5UewNqs+CCf+HY64QEVIdZFZv1Wxa38jT1zQ9fe6i8KsBL74tZI9dRaSnJz0D9Y9Ch2lEBlmJD81lg09DEvd39DE9e8tJTLMwMuXjSYrMZp3rxnLd2t38vDXa7nk9V8444hM7jtjCH0So90kffCglUJH7FyrCuCd9ECHDVSCGiFg0t/h3bNZ/eVz7G8YyfUT+nluvpYm2PITrP6ED2u/IrylnqbKTMKPvlk5v9OH9Sgk+K352+iXGsvxA31XQ8sxs/nYAe13m6u21z0K8rLZbSlIj3M5ua89rFbJnR+vpNR0gA/+cBRZtpO+EIJThmZwfEEarxZv5aU5m/lx/U5umTiA647r1+WSLMGMNh91xJLXwRgJRVf4WhLfkn881r7jKdj4GhNyYxjh7vLSUqoifP+9C54eDNN+D5u+o3noBUxp+jvPDJsBJz8MGUf0SCEsL9vDyvK9XDUuD4PBd7kmdqXQkQnJZLYQG2F0S0JgIFGQHk9pzYGuhSM78OJPm/lu7U7+evoQju6Xctj2qHAjfzppID/ccTwTC3rz9PcbOflfc/ne1ipVo5WCc+r3qqb2R1yg+g2EMkKwoO9NpLCXv2fMd9+45t0w53F4bgS8+TtY8R7kT1CZyHdtJPq8F4gccBxf/lrVs7IXNt6ev434qDDO91AjHVdJjYukd3xkx0qhLnSymR0ZlB6PlMon0FV+2rCLZ37YyDkjszrNUs9JjuGVy0fz/rVHERFm4Lp3l3LV20vYurvr8wYbWik4Y+V0aDqgWlWGOFJKnlybyC/G0Qzc9AY01PZsQNMW+Op2eHaYimxK7gfnvAx3bYLfvwODz4AwdUI8t6gPlbUNLCrpZtkLG1W19cxaXcVFY3KI9YOqo4VZCayrcm47N5lDp+6RIwUZ3auBtK26jj9NX8GQjAQeO2+4y0UExw9MZdafJvC3M4awrHQPpzw7l8dmrcPc6IV+4X6KVgrtYbWqSJfsscHXN6EbLC6pYWVFLTVH3Y2o3wOLXu7eQOVLVJLZv0eraKbhF8ItS+CKL5wmBp5cmEFcZBif97By6nsLS7FKyZXH5vVoHHcxJDOBzbv2O03UqrZVSA01+ibHEBFm6FJmc11jMze8twyDQfDq5aO7bHILNxr4w4R+zL7reCaP6MOrxVuZ9PQcvlixPSRNSloptMfW2VCzJTTDUNvh9XklJMdGcMIJJ8OQs2DBC6rnsStYrbBhFrx1Grx5kirLPf7PcPtqmPzvTnM/oiOMnDosg1m/7ei2nbmhqYXpi8s4aUi6x6OmXKUwM4GmFunUTGKqswR9G872CDMaGJAW53IEkpSSez5bxaZd+/n3lKIefb+946N4+sIRfHbTsfSOj+L2j37lwlcXsqayhyvjAMP362h/ZPHrENsbCs8+bJOl2eqR0g7+ypbdZn5Yt5M/TRqoIjRO+Cus+xrmPwe/e9D5G5sbYdXHsODfKiu6Vw6c8hiMuhwiu5aYdF5RHz5dVsH363Z2q9bSFyu2s+dAE9eM936ymjMcnc320hd2rFbJnhAqm92WQRnxLlfJfW3uVv67qop7TxvMBDdFlI3um8QXt4zj46XlPPnt/7d33uFRlVkD/500UkgFQoYmPSaABIiKVEEsgIsLrl1Xdy1rwVX5rA/7+G37nmd11d111wJYV8VGca1rQZAI0iGAhM5AAqElUpJICnm/P+7NEOKkzE0mmTDn9zz3mXtn5j333Jn33vOW856zhZ/981sGdI5HRKx8wnbvoaoPUdWZMPY7nmPP+xYJUVae7m7toumaFE3XxCi6JUWTFBMRUDkz1CjUpHAXbP0cRj0EYaduytKKk0yfv5EPs/dx64ge3DOmN20DYGza37yUtYuIsBBuusAOMJecBgOuguUzrLzOsTUijP54xFoBvnwGFO23VkRPmQX9JluB9hwwtGc7XPGRzF+T57NRMMbw6hI3aa44zu9Rd9TU5qRH+xgiw0O8rmw+dqKcikoTlMNHYHkgzV+7l6M/lhMfVXud+XbbYZ7472YmDnDxm1FN6yYdGiJcd143JvR38a+F2zzZ8qoe3lWPcJHq+6c+O/WMP/WwLywuZcHmgz8JuRETEWoZiaRoy2jYW9ekKLokRje7u+yZ/1TzlZUvgYRA5q88bx06Xspv3ljFmj1HGNoziRcW7eD9VXk8fGkqVw7pQmgLujf6k8NFpcxbk8eVg7ucHm7hwkdh41z49hkY/4T13tE8a65h9WtQVmSFBZn8AvQc0+hQ4yEhwhUZnZmVVUvYizr4bkcBWw4c58lfNHzysTkIDRFSU+K8+uQfDqLczN5ITbFWmm8/eJwhZ3k35LmFJdz79hp6J7f1638bHx3O9InpTSqzpKyC3MIf2VNYQm5hiefVfbiYrG2HOFF++jxTx7g2lpFIjObS/ilc2i+lSfWpiRqF6pSVWG6RaT+DOKtFunHvUe749yoKS8p4/obBTBjgIjv3CH/8eBMPz13P69+5efzydM734hPd2nnju92UVlRya81hl3a9rDhDq16B3uOsGEQb51r95f5TYNi9DYpU6gtTBnfmxW928FH2Pp+GgV5Z4qZdTESTh/huCtJdcXy2MR9jzGkPtQK7JRnMPQWALfuLvBqFE+UnufPN1VRUGmbclBkQ3mS+EB0RRmpKrNf4TsYYDhWVVjMWlvHYU1jCsp0F9Ggfo0ahWdnwvuVuaU8wf7Yhn2nvZZMQHc6cO4d5UjYO7JrAnDsv4KP1+fzl0xyumbmMCQNSeGx8WsBMZDaWE+UneWPZbsalJXuPETTqYWsdx1u/gPAYK4LsBXdDgn8ip/btGEu/TnHMX7u3wUZhd0ExCzYfYOqY3gG5YjXdFcvbK/aQf/SEZ+UtWJPMELw9hc4JUcREhHr1QDLG8Ni8DWzKP8bLN2fSo31MC2joP0SE5NhIkmMjvRrE5vCGUqNQhTHWBHNyP0y3C/jngm088+VWBnVLYMZNQ0iOPT2to4gwaWAnLk7ryKysnbywaAdf5RzkthE9uPsMmG+YuyaPwuIybqstpEVCVytUddF+GHxzvVnOmoLJgzrz509y2H7wOL2T65+sfn3pbkJFuHFo3Ql3WoqqCeac/GOnG4Wi4It7VB0RoU/HWK8eSK8tdTN/7V6mXdyXsWc3X8a8QKE5hkCDx42mPvYsgwMbKBtyG/e+s45nvtzKlEGdefv2oT8xCNWJigjltxf14esHR3P5ABfPL9rBmKcW8d6q3CZZhdsSVFYaXs7axTld4uuenM24znIvbQaDADApoxMhAvMasGbh+Ily3luVy8RzXHSMq/3/a0lSU7wn3KmaU0iKDk6jAN6zsC3bWcCfP8lhXFpHpo7p3UKanfmoUahixUwq28Rz/fJufLLBcnF7+uqBDR52cMVH8cw1GXxwz3C6Jkbx8Jz1THruW1bsaqA/fwCxYPNBdh4u5raRPQNqcjY5NpKRfTrwn3X76jW4c1bnUVRawa/8HOK7MbRtE0b3dtHk7D/dKBQUl5IYHd6saUIDjb4psRQUl3k8dfKP/sjU2Ws4KymaZ64Z2KKxq850grfWVedYPmbTh7xdPoqcwxXMuimTO0f3cvRAzOiawNy7hvGPazMoKCrj6hnfcc/sNeTWk5M3kJiVtZPOCVFM6O/fCS0nTBncmb1HfmSFu3ZjW1lpeH2pm0HdEsho6uB9TUya66ceSFaIi+CcZK4i1Z5s3rr/OKUVJ7nzzTX8WHaSmb8cQlxk82bLCzbUKACbP3kWU3mSjyLGM+/u4YxLb9xYpYjlQvn1/1zI/eP6sCDnABc98w1Pfb6F4gCPqbIu9wgrdhXyq+HdA7Klekl6CjERoXWGvVi45SDugpKA7iVUke6KY3dhyWmxdgqKyoIqDac3+tpuqZv3H+fxD74nO/cIT1+d0aC5JKVxBN5d34xUVhqe/mwDSZtnsy4yk+fvvapJ0gBWERURyv3j+rLwwQuZ0D+Ffy3c7plvOBmg8w2zsnYSGxnGtX7Kv9xYrLAXLj7dkM+Jcu9hL15d4iYlLpLxAdjTqUmaKw5jYEu1IaSC4tKgSsPpjQ5t25AYHc5LWTt5d1Uu94zpxWWt4P88Ewhao1BUWsFv3lzN7qx3SJYjDJj8kN9SH7rio/j7tYOYd/cwOiVY8w2X/O0bPl5f/9h4c5JbWMJnG/K5/vxuAe09NWVwZ46XVvBVzoGffLb1wHG+3X641eTirfJA2lQtYmpBcVlQpeH0hojQt2Ms+46eYHTfDky7OLWlVQoaAv+u8QO5hSX84oWlLMg5wPQOWZjEHoT3vdjv5x3cLZH5dw/jxRsHExoiTJ29lgnPZgVMgo9XluwiRIRbAiSSaG0M7dmOlLhIr0NIry5x0yYshOsDtKdTE1d8JPFR4Z55hfKTlRwpKQ9ad9TqjOzTnr4d2/LstYPO2KgBgUjQGYWV7kKueG4Je4/8yJwroul4NBs573YrSX0zICJc1t/FZ/eN4h/XZnCi/CS3/3sVP39uCYu3Hmox43C0pJx3V+YyaWAnXPGBnbc2NES4YlAnvtl6yOPTD3CkpIz5a/OYPKhzrbmPAw0RId0V54mB9INn4VpwDx8BTB3bh8/vH0V8tE4sNydBZRTeXbmH62ctIyEqnA/uGc7gA3MgPBoybmh2XULteD5fTRvNk1eew+GiMn75ygqumbGM5Q2MENmUzF6xh5Kyk7UvVgswpgzqQkWl4aPsfZ733l6Ry4nySm6pJ+tWoJHmimPz/mOcrDSeNQrtW4lR8zeB5BIdLPjVKIjIZSKyRUS2i8ijXj6fJiKbRGS9iCwQEb8tPZ2zOo9H5m5gaM92zL97OL1iymDDHDjnGohqObfFsNAQrj63K18/OJo/XdEPd0Ex18xcxo0vLWfNnh+aRYeyikpeW7qLEb3b/ySMc6CSmhJLussKewHWsMu/v3MzrFc7zk5pHddQRXqnOE6UV+IuKKaguGo1s/YUlJbBb0ZBREKB54DxQDpwnYjUDDe4Fsg0xpwDzAGe9Jc+4/un8MhlZ/PqLeda3dG1b0DFiYBJt9kmLJSbLujO4ofH8LuJaeTkH2PK80u59bWVfk/y8VH2Pg4cK+W2kYHvwlmdKYM7k513lB2Hivj8+/3kHz3RKtxQa5LmsjzeNu07RkGQR0hVWh5/9hTOA7YbY3YaY8qAd4DTstYYYxYaY6pWdS0D/JZRPaZNGHdd2Mvyva88aYXIPmsEdOznr1M6IjI8lNtG9mTxw2N46NJUVroLmfjst9z91mq2+Zi3tiEYY5iVtZPUjrGM7ts0SUqai0kDrbAX89fs5dUlbrolRTP27OSWVstn+iTHEh4q5OQf86zgbR+kEVKVlsefRqEzkFvtOM9+rzZuBT7zoz6n2PYFHNkTML0Eb8S0CeOeMb3JemQsv72oD99sOcQlf1/MA++uw324uMnO8+32w2zef5xbR/ZodeO3yXGRjOjTgde/c7N69w/cPKx7q/RSiQgLoVeHtmzKP0ZBcRlhIUJcVOC6BCtnNv6sed7uTq+uNSJyI5AJjK7l8zuAOwC6dWsCV8MVMyG2E5w9sfGy/Ex8VDjTLu7LLcO6M2PxDl5f6ubD7H2M75/CsF7tOa9HIr06tHX8QJ+VtYsOsW24IiPw8g00hCmDOrN46yFiIkK5KtNvHU2/k+6KY8mOwyTHtqFd28BKz6gEF/40CnlA12rHXYB9Nb8kIuOA6cBoY0xpzc8BjDEzgZkAmZmZjfPZPLwNdnwNY37nOD1kS5AUE8Fj49O4dUQPnl+4g4/X7+Pj9fkAJEaHk9k9iXO7J5LZPYn+neIblEc6J/8Yi7ce4qFLU2kTFnj5BhrCJf06khAdzlVDurTqmDjpneKYt3Yv2w4WBW1yHSUw8KdRWAn0EZEewF7gWuD66l8QkUHADOAyY8xBP+pyihWzICQchtzcLKdrapJjI/n9pH7878/ScReUsNJdyMpdhaza/QNfbrJW+EaGh5DRNYFzuydxbvckBnVLINbLA/OlrF1EhYdyw/mtY6GXN6Ijwlj04IUBvQK7IaS5LI+p7NwjDO/dvoW1UYIZv91JxpgKEZkKfA6EAq8YY74XkT8Cq4wxHwJ/BdoC79vd5T3GmEn+0onS47ButpVEvm3rm5CsjojQo30MPdrHcHWm1SE7ePwEq90/sNL9AyvdhTy3cDuVBkLEeuhUGYlzuydigA+z93L9ed1IaOVx+1u7/nDKKFQagj4YntKy+LV5ZYz5FPi0xnuPV9sf58/z/4Tsd6DsuCfd5plGcmwk4we4GD/ABVjxndbtOcIKdyGr3IW8uzKX15a6AYiLDONkpfEp37HiP5JiIkiJi2T/sRO6RkFpUVp3n9sXqtJtdhoEXTJbWptmoW2bMEb0ac+IPtZwRPnJSr7fd4xV7kJWugs5OyWOs9qdWTluWzPpneJso6A9BaXlCB6jsGsxHN4CP38BgtSzIzzUmmvI6JrQasJZBBNprli+3nxQ1ygoLUrwxD6KSoABV0O/KS2tiaJ4Jd0VDxD0YbOVliV4egqugXDlrJbWQlFq5cLUDtw+sgcX9GrX0qooQUzwGAVFCXBi2oQxfWLN8GCK0rwEz/CRoiiKUi9qFBRFURQPahQURVEUD2oUFEVRFA9qFBRFURQPahQURVEUD2oUFEVRFA9qFBRFURQPYkzjctY0NyJyCNjtsHh74HAjVVAZKiPQZQSCDioj8GScZYypNxF7qzMKjUFEVhljGhUiVWWojECXEQg6qIzAlVEfOnykKIqieFCjoCiKongINqMwU2WojCCQEQg6qIzAlVEnQTWnoCiKotRNsPUUFEVRlDpQo6AoiqJ4CAqjICKviMhBEdnYCBldRWShiOSIyPcicp8DGZEiskJEsm0Zf3CoS6iIrBWRj52Ut2W4RWSDiKwTkVUOZSSIyBwR2Wz/Lhf4UDbVPnfVdkxE7negwwP2b7lRRN4WkUgHMu6zy3/fUB281SkRSRKRL0Vkm/2a6EDGVbYelSJSr+thLTL+av8n60VkvogkOJDxJ7v8OhH5QkQ6+Sqj2mcPiogRkfYO9Pi9iOytVk8mONFDRO4VkS32b/ukAz3eraaDW0TWOZCRISLLqu45ETnPgYyBIvKdfe9+JCJxdclwhDHmjN+AUcBgYGMjZLiAwfZ+LLAVSPdRhgBt7f1wYDkw1IEu04DZwMeNuB430L6Rv+vrwG32fgSQ4FBOKLAfa3GNL+U6A7uAKPv4PeAWH2X0BzYC0ViZCL8C+jipU8CTwKP2/qPAEw5kpAGpwCIg06EelwBh9v4TDvWIq7b/W+BFX2XY73cFPsdacFpnfatFj98DD/rwf3qTMcb+X9vYx8lOrqXa508DjzvQ4wtgvL0/AVjkQMZKYLS9/2vgT77U94ZsQdFTMMYsBgobKSPfGLPG3j8O5GA9lHyRYYwxRfZhuL35NNMvIl2AicBLvpRrauwWyijgZQBjTJkx5ohDcRcBO4wxTlaqhwFRIhKG9WDf52P5NGCZMabEGFMBfANMrq9QLXXqCixDif36c19lGGNyjDFbGqh7bTK+sK8FYBnQxYGMY9UOY6inntZxj/0NeLi+8vXIaDC1yLgL+IsxptT+zkGneoiIAFcDbzuQYYCqln089dTVWmSkAovt/S+BK+uS4YSgMApNjYh0BwZhtfR9LRtqdz0PAl8aY3yV8Xesm6zS13PXwABfiMhqEbnDQfmewCHgVXso6yURiXGoy7XUc5N5wxizF3gK2APkA0eNMV/4KGYjMEpE2olINFYLrquvuth0NMbk27rlA8kO5TQlvwY+c1JQRP5PRHKBG4DHHZSfBOw1xmQ7OX81ptpDWa/UNyRXC32BkSKyXES+EZFzG6HLSOCAMWabg7L3A3+1f9OngMccyNgITLL3r8J5Xa0VNQo+IiJtgbnA/TVaUw3CGHPSGJOB1Xo7T0T6+3Duy4GDxpjVvp7XC8ONMYOB8cA9IjLKx/JhWF3bF4wxg4BirCETnxCRCKxK/r6DsolYrfMeQCcgRkRu9EWGMSYHa4jlS+C/QDZQUWehVoKITMe6lreclDfGTDfGdLXLT/Xx3NHAdBwYkxq8APQCMrAM/9MOZIQBicBQ4CHgPbvF74TrcNCAsbkLeMD+TR/A7mX7yK+x7tfVWMPYZQ51qRU1Cj4gIuFYBuEtY8y8xsiyh1oWAZf5UGw4MElE3MA7wFgRedPh+ffZrweB+UCdk15eyAPyqvV05mAZCV8ZD6wxxhxwUHYcsMsYc8gYUw7MA4b5KsQY87IxZrAxZhRWd91JKxDggIi4AOzXOocp/ImI3AxcDtxg7AHoRjAb34cpemEZ62y7vnYB1ohIii9CjDEH7IZUJTAL3+spWHV1nj18uwKrl13npLc37CHKKcC7DnQAuBmrjoLVCPL5Wowxm40xlxhjhmAZpx0OdakVNQoNxG5ZvAzkGGOecSijQ5UniIhEYT3UNje0vDHmMWNMF2NMd6whl6+NMT61jO1zx4hIbNU+1sSkT55Zxpj9QK6IpNpvXQRs8lUXGtfy2gMMFZFo+/+5CGuuxydEJNl+7YZ10zvV50OsGx/79T8O5TQKEbkMeASYZIwpcSijT7XDSfhQTwGMMRuMMcnGmO52fc3DctTY76MermqHk/Gxntp8AIy15fXFcopwEml0HLDZGJPnoCxYcwij7f2xOGh8VKurIcDvgBcd6lI7TT1zHYgb1k2eD5RjVc5bHcgYgTUOvx5YZ28TfJRxDrDWlrGRejwY6pF1IQ69j7DmA7Lt7XtgukM5GcAq+3o+ABJ9LB8NFADxjfgd/oD1wNoIvIHtYeKjjCwsg5YNXOS0TgHtgAVYN/sCIMmBjMn2filwAPjcgYztQG61elqf55A3GXPt33Q98BHQ2VcZNT53U7/3kTc93gA22Hp8CLgcyIgA3rSvZw0w1sm1AK8BdzaifowAVtv1bDkwxIGM+7A8H7cCf8GOStGUm4a5UBRFUTzo8JGiKIriQY2CoiiK4kGNgqIoiuJBjYKiKIriQY2CoiiK4kGNgqI0ASJSVP+3FCXwUaOgKIqieFCjoCiKonhQo6AoiqJ4UKOgKIqieFCjoCiKonhQo6AoiqJ4UKOgKE1DtIjkVdumtbRCiuIEjZKqKIqieNCegqIoiuJBjYKiKIriQY2CoiiK4kGNgqIoiuJBjYKiKIriQY2CoiiK4kGNgqIoiuLh/wE+3FO2Gzvz+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5a93edbd30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VPXSwPHvJIQWmoROCIQqvfcuFlABKQqI2MVe77W+XrteLHhFsYCiiEpsgCKCoNJB6b0KoYXQQ4dAyrx/nA0sIZBN2M3uJvN5nvOEPXV2jTs5Z35FVBVjjDEmMyH+DsAYY0xwsIRhjDHGI5YwjDHGeMQShjHGGI9YwjDGGOMRSxjGGGM8YgnD5BkickxEqubQtT4Rkf/kxLWMySli/TCMN4jITKAhUE5VT/k5nBwlIrcDd6tqO3/HEkxEZDQQp6rP+zsW4xm7wzCXTESqAO0BBXr48Dr5fHXu3Cijzyurn6F95sadJQzjDbcCfwOjgdvcN4jIaNfjmd9F5KiIzBKRym7bVUQeEZFYEdkvIm+LSIhr2+0iMk9E/iciCcBLIhIiIs+LyDYR2SsiY0SkuGv/fq7zFHO97iYiu0WktNu1qrvF9ZGITHE9qponIuVE5D0ROSgi60WksVucz4jIZtd7WCsivVzrawOfAK1d5znkdv7X3I6/R0Q2iUiCiEwUkQrpPoP7ROQf17U/FBHJ6IN2vf+0WA6IyPciUtK1rYrrXHeJyHZgekbrXPv2EJE1InJIRGa63kfaNbaKyNMishI4nj5puP57vpNu3c8i8oTr30+LyE7XZ7VBRLpk/Gtjgo6q2mLLJS3AJuABoCmQBJR12zYaOAp0AAoAw4C5btsVmAGUBKKAjTiPdwBuB5KBh4F8QCHgTtf1qgJFgPHAV27n+8Z1zQggHrg+3bWqu8W13xVzQZwv0i04yS8UeA2Y4XbsjUAFnD+y+gHHgfJucc5N95mMBl5z/fsK17WauD6DD4DZ6eKaBJRwfQb7gK4X+Kwfw0nOka5zjQBiXNuquM41Bgh3fV4Zravpiv8qIAx4yvWZ5nedZyuwHKgEFMoghg7ADs4+0r4MOOn6fGq5tlVwi6naBd7Lmc/IluBY/B6ALcG9AO1wkkQp1+v1wONu20cD37q9LgKkAJVcr9X9yxEn8fzp+vftwPZ01/sTeMDtdS3X9fO5XpcAtgOrgBHpjk2fMD512/YwsM7tdX3g0EXe93Kgp1ucF0sYo4C30n0GSUAVt7jauW3/HnjmAtddB3Rxe10+7f27JYeqbtszWvcf4Hu31yHATqCT6/VW4M6LvHdxfcYdXK/vAaa7/l0d2AtcCYRl8rtjCSPIFnskZS7VbcA0Vd3vej2WdI+lcP7iBEBVjwEJOH+Nnrcd2HaRbbi2bUu3fz6grOv8h4AfgHrA0Exi3+P275MZvC6S9kJEbhWR5a5HOIdc5y+VyfkzjNn1GRwAKrrts9vt3yfcr51OZWCCWxzrcBJwWbd90n9m6deljyfVtb3iBfY/hzrf9t8CA1yrbsa5s0NVN+HcBb0E7BWRb90fv5ngZgnDZJuIFAJuAjq6agW7gceBhiLS0G3XSm7HFMF5/BSf0XacRzLu29I344vH+dJ03z8Z15e9iDTCeWwVA7yfjbd1HlfN5VPgISBCVUsAq3H+0s4oxvTOiVlEwnEeme3MRjg7gG6qWsJtKaiq7ufKKB73denjEZz/Bpmdw10M0Nf12bQExp05UHWsOi3GKrvO86YH78sEAUsY5lLcgPPXbR2gkWupDczBqQWkuVZE2olIfuBVYIGquv8F+6SIXCYilYBHge8ucs0Y4HERiXYlnzeA71Q1WUQKAl8DzwF3ABVF5AEvvM9wnC++fQAicgfOHUaaPUCk6/1lZCxwh4g0EpECrpgXqOrWbMTyCfB6WsMBESktIj2zeI7vgetEpIuIhAH/Ak4B8z09gaouw/k8PgOmuu7sEJFaInKF630m4typpVzkVKEiUtBtudBnaAKAJQxzKW4DvlDV7aq6O20BhgMD3VrXjAVexHkU1RQYmO48PwNLcOoCv+I887+Qz4GvgNk4RepEnPoDwH9x2vV/rE5fkFuA10SkxqW8SVVdi/N46y+c5FAfmOe2y3RgDbBbRPZncPyfOHWDccAuoBrQP5vhDAMmAtNE5ChOAbxlVk6gqhtwPpsPcIrx3YHuqno6i7HE4NQqxrqtKwAMcZ13N1AGJ4FfyDM4SSVtmZ7FGEwOso57xqckk85ZIqJADdezb2NMALM7DGOMMR6xhGGMMcYj9kjKGGOMR+wOwxhjjEdy1cBipUqV0ipVqvg7DGOMCRpLlizZr6qlPdnXZwnD1aZ+DFAOSAVGquqwdPsMBJ52vTwG3K+qK1zbtuKMQZQCJKtqs8yuWaVKFRYvXuy192CMMbmdiGzLfC+HL+8wkoF/qepSESkKLBGR311t2tNsATqq6kER6QaM5Nw25Z3dhpwwxhjjRz5LGKq6C6eTEqp6VETW4YxVs9ZtH/eepWkjcBpjjAlAOVL0FmeCncbAgovsdhcwxe214vRmXSIig30XnTHGGE/4vOjtGu9nHPCYqh65wD6dcRKG+xSXbVU1XkTKAL+LyHpVnZ3BsYOBwQBRUVHnnTspKYm4uDgSExMv/c3kMQULFiQyMpKwsDB/h2KMCQA+TRiugc3GAd+o6vgL7NMAZwCzbqp6IG29qsa7fu4VkQlAC5zxg86hqiNxah80a9bsvE4lcXFxFC1alCpVqnCBScxMBlSVAwcOEBcXR3R0tL/DMcYEAJ89knINmTwKZ1Kady+wTxTOjGmDVHWj2/pwV6E8bSjoq3GGk86yxMREIiIiLFlkkYgQERFhd2bGmDN8eYfRFhgErBKR5a51z+HMX4CqfgK8gDMvwEeuL/S05rNlcSaJSYtxrKr+lt1ALFlkj31uxhh3vmwlNZezE8xcaJ+7gbszWB8LNDz/CGOMMe7m/LOPTXuPMahVZfKF+rYdkw0NkkMmTJiAiLB+/Xp/h3LG7NmzadKkCfny5ePHH3/0dzjGmCxKSknlpYlrGPPXNlJzYFhASxg5JCYmhnbt2vHtt9967ZzJycmXdHxUVBSjR4/m5ptv9lJExpic9PXf29i87zj/d21t8ufz/de5JYwccOzYMebNm8eoUaPOSRgzZ86kQ4cO9OrVizp16nDfffeRmpoKQJEiRfjXv/5FkyZN6NKlC/v27QOgU6dOPPfcc3Ts2JFhw4axbds2unTpQoMGDejSpQvbt28HoGfPnowZMwaAESNGMHBg+knunKFUGjRoQEiI/RoYE2wOHj/Ne3/8Q/sapehSu0yOXDNXDT6YmZd/WcPa+Ay7gmRbnQrFeLF73Yvu89NPP9G1a1dq1qxJyZIlWbp0KU2aNAFg4cKFrF27lsqVK9O1a1fGjx9P3759OX78OE2aNGHo0KG88sorvPzyywwfPhyAQ4cOMWvWLAC6d+/Orbfeym233cbnn3/OI488wk8//cTIkSNp27Yt0dHRDB06lL///tur79sY41/v/bGRo4lJPH9dnRxroGJ/WuaAmJgY+vd3pnDu378/MTExZ7a1aNGCqlWrEhoayoABA5g7dy4AISEh9OvXD4BbbrnlzHrgzHqAv/7668wjpUGDBp3Zr2zZsrzyyit07tyZoUOHUrJkSd++SWNMjtm45yhfL9jOwJaVqVWuaI5dN0/dYWR2J+ALBw4cYPr06axevRoRISUlBRHhrbfeAs5vunqhvxTc14eHh1/weu77rVq1ioiICOLj4y/lLRhjAoiq8uqktYTnD+Xxq2rm6LXtDsPHfvzxR2699Va2bdvG1q1b2bFjB9HR0WfuBBYuXMiWLVtITU3lu+++o107Z3SU1NTUMy2Xxo4de2Z9em3atDlTF/nmm2/O7Ldw4UKmTJnCsmXLeOedd9iyZYuv36oxJgfM2LCXOf/s59Era1IyPH+OXtsSho/FxMTQq1evc9b16dOHsWPHAtC6dWueeeYZ6tWrR3R09Jl9w8PDWbNmDU2bNmX69Om88MILGZ7//fff54svvqBBgwZ89dVXDBs2jFOnTnHPPffw+eefU6FCBYYOHcqdd95J+ul4Fy1aRGRkJD/88AP33nsvdevm/B2YMcZzp5NTeW3SOqqWDufW1pVz/Pq5ak7vZs2aafoJlNatW0ft2rX9FNHFzZw5k3feeYdJkyadt61IkSIcO3bMD1GdK5A/P2PymlFzt/DqpLV8fnszrri8rFfOKSJLPJmgDuwOwxhjgkLC8dMM+2MjHWqWpnOtnGlGm16eKnoHmk6dOtGpU6cMtwXC3YUxJnD87/eNHD+dwn+uq+23cd7sDsMYYwLcht1H+WbBNm5pGUWNsjnXjDY9SxjGGBPA0prRFi0YxmNX5mwz2vQsYRhjTAD7c91e5m7az2NX1uCyHG5Gm54lDGOMCVCnk1N5ffI6qpUO55ZWOd+MNj1LGDkkEIc3f/fdd6lTp86ZgQu3bdvm75CMMW7G/LWVLfuP8/z1dQjz8VwXnvB/BHlEIA5v3rhxYxYvXszKlSvp27cvTz31lJciM8ZcqgPHTjHsz3/oVMt/zWjTs4SRAwJ1ePPOnTtTuHBhAFq1akVcXJxPPwdjjOfe/X0jJ06n8Px1gdNxNm/1w5jyDOxe5d1zlqsP3YZcdJdgGN581KhRdOvWzQsfiDHmUq3bdYSYhdu5tXUVqpfxXzPa9OwOIwcE+vDmX3/9NYsXL+bJJ5/00js2xmRXWjPaYoXCeOzKGv4O5xw+u8MQkUrAGKAckAqMVNVh6fYRYBhwLXACuF1Vl7q23QY879r1NVX98pKDyuROwBcCfXjzP/74g9dff51Zs2ZRoEABj96TMcZ3fl+7h/mbD/Byj7qUKOzfZrTp+fIOIxn4l6rWBloBD4pInXT7dANquJbBwMcAIlISeBFoCbQAXhSRy3wYq88E8vDmy5Yt495772XixImUKRMYRTVj8rJTySm8PnkdNcoUYWDLKH+Hcx6fJQxV3ZV2t6CqR4F1QMV0u/UExqjjb6CEiJQHrgF+V9UEVT0I/A509VWsvhTIw5s/+eSTHDt2jBtvvJFGjRrRo0cPH3wCxhhPfTl/K9sOnOD56+uQLwCa0aaXI8Obi0gVYDZQT1WPuK2fBAxR1bmu138CTwOdgIKq+ppr/X+Ak6r6zsWuY8Obe18gf37G5Cb7j52i89szaR5dks9vb55j1w2o4c1FpAgwDnjMPVmkbc7gEL3I+ozOP1hEFovI4rSmp8YYE2yGTtvIyaQU/i+AmtGm59OEISJhOMniG1Udn8EucUAlt9eRQPxF1p9HVUeqajNVbVa6dGnvBJ5DOnXqlOHdBdjw5sbkJWvjj/DdIqcZbbXSRfwdzgX5LGG4WkCNAtap6rsX2G0icKs4WgGHVXUXMBW4WkQucxW7r3aty5bcNKtgTrLPzRjfU1VembSG4oXCeLRLYDWjTc+XHffaAoOAVSKy3LXuOSAKQFU/ASbjNKndhNOs9g7XtgQReRVY5DruFVVNyE4QBQsW5MCBA0RERPht0pFgpKocOHCAggUL+jsUY3K1qWv28HdsAq/2rEvxwmH+Dueicv2c3klJScTFxZGYmOinqIJXwYIFiYyMJCwssH+JjQlWp5JTuOrd2RQMC2HyI+390jIqK0XvXD80SFhYGNHR0f4OwxhjzvPFvK1sTzjBV3e1CMhmtOkFfoTGGJML7Tt6iuHTN3Fl7TK0rxEcDXYsYRhjjB8MnbaBU8kpPHdt4DajTc8ShjHG5LDVOw/z3eId3Na6ClUDuBltepYwjDEmBznNaNdyWeH8PBzgzWjTs4RhjDE56LfVu1m4JYEnrqpJ8ULB1QLREoYxxuSQxCRnNNpaZYvSv3mlzA8IMJYwjDEmh3w+bwtxB0/yQvfAHI02M8EXsTHGBKG9RxL5cPomrqpTlrbVS/k7nGyxhGGMMTngnWkbOJ2SGlTNaNOzhGGMMT62eudhflgSxx1to4kudeEplgOdR0ODuEaMrQCcBLaqaqpPozLGmFxCVXn5lzWULJyfh66o7u9wLskFE4aIFAceBAYA+YF9QEGgrIj8DXykqjNyJEpjjAlSk1ftZtHWg7zRqz7FCgZXM9r0LnaH8SMwBmivqofcN4hIU2CQiFRV1VG+DNAYY4JVYlIKb0xex+XlitIvCJvRpnfBhKGqV11k2xJgiU8iMsaYXGLU3C3sPHSSsfe0JDQk+Ofj8Xh4cxEpDTwKFAI+VtVNPovKGGOC3J4jiXw4YxPX1C1Lm2rB2Yw2vay0khoKzAZ+A2J8E44xxuQOb0/dQHKKBnUz2vQumDBE5DcRae+2Kj+w1bUU8G1YxhgTvFbGHeLHJXHc0a4KlSOCtxltehe7w+gH9BSRsSJSDfgP8AIwBHggJ4Izxphgo6q88staShXJz0Odg7sZbXoXK3ofBv4tIlWB14GdwIOu9cYYYzIwaeUuFm87yJDe9Ska5M1o07tYP4yqwP1AEvAvoBrwvYhMwumDkXKxE4vI58D1wF5VrZfB9ieBgW5x1AZKq2qCiGwFjgIpQLKnE5QbY4w/JSalMGTKemqXL8aNzYK/GW16F3skFYNT4P4b+EpV56jqNcARYJoH5x4NdL3QRlV9W1UbqWoj4FlglqomuO3S2bXdkoUxJih8OjuWnYdO8mL3OrmiGW16F0sYBYEtrqVw2kpV/RLnzuGiVHU2kJDZfi4DsJZXxpggtvtwIh/N3Ey3euVoVTXC3+H4xMX6YTwAvA2cBu5z36CqJ70VgIgUxrkTecj9EsA0EVFghKqO9Nb1jDHGF96aup6UVOXZbrmnGW16Fyt6zwPm5UAM3YF56R5HtVXVeBEpA/wuIutddyznEZHBwGCAqKgo30drjDHpLN9xiPFLd3J/p2pERRTO/IAgdbF+GL+IyPUicl6ZX0SqisgrInKnF2LoT7rHUaoa7/q5F5gAtLjQwao6UlWbqWqz0qVLeyEcY4zxnNOMdg2lihTgwVzWjDa9i9Uw7gE6AOtFZJGITBaR6SISC4wAlqjq55dycdeIuB2Bn93WhYtI0bR/A1cDqy/lOsYY4ysTV8SzdPshnrqmFkUKeDzaUlC62COp3cBTwFMiUgUojzMfxkZVPZHZiUUkBugElBKROOBFIMx17k9cu/UCpqnqcbdDywITRCQtvrGq+luW3pUxxuSAk6edZrR1KxSjT9NIf4fjcx6lQ1XdijMkiMdUdYAH+4zGaX7rvi4WaJiVaxljjD+MnB3LrsOJvNevUa5sRpueTdFqjDHZsOvwST6ZtZlr65ejZS5tRpueJQxjjMmGt37bQIrm7ma06V00YYhIqIh8nVPBGGNMMFi6/SATlu3knvbRVCqZe5vRpnfRhOEaL6q0iOTPoXiMMSagpaY6o9GWLlqA+zvl7ma06XlS9N4KzBORicCZ1kyq+q6vgjLGmEA1cUU8y3cc4u2+DXJ9M9r0PHm38a4lBCjq23CMMSZwnTidzJAp66lfsTh9muT+ZrTpZZowVPVlAFdnOlXVYz6PyhhjAtCIWbHsPpLIBzc3JiQPNKNNL9NWUiJST0SW4fS2XiMiS0Skru9DM8aYwBF/6CQjZm/mugblaV6lpL/D8QtPmtWOBJ5Q1cqqWhlnMqVPfRuWMcYEliFT1qMKz3a73N+h+I0nCSNcVWekvVDVmUDumdXcGGMysWRbAhNXxDO4Q1UiL8s7zWjT86ToHSsi/wG+cr2+BWdSJWOMyfXSmtGWKVqA+zpW83c4fuXJHcadQGlgvGspBdzhy6CMMSZQ/LR8JyviDvN018sJz2PNaNO76LsXkVDgOVV9JIfiMcaYgHH8VDJv/raehpHF6dW4or/D8TtPeno3zaFYjDEmoIyYtZk9R07xQvc6gduMNjUVkhJz5FKe3F8tc/Xy/oFze3qP91lUxhjjZ3EHTzBidizdG1agaeUAbkY74zXYPB1umwQFivj0Up4kjJLAAeAKt3WKU88wxphcaciU9YjAM4HcjHb1OJgzFJrcCvl933jVkxrGSlX9n88jMcaYALF4awKTVu7ikS41qFiikL/DyVj8MvjpQYhqDdcOBfH9IzNPahg9fB6FMcYEiNRU5eVf1lKuWEHu61jV3+Fk7Oge+HYgFI6Am76CfDkzoLgnj6Tmi8hw4DvOrWEs9VlUxhjjJ+OX7WTVzsP8r19DCucPwGa0yafg+0FwIgHumgpFSufYpT35NNq4fr7itk45t6ZhjDFB7/ipZN76bT0NK5WgZ8MAbEarCpOegB0L4MbRUL5hjl7ek9FqO2fnxCLyOXA9sFdV62WwvRPwM2d7jY9X1Vdc27oCw4BQ4DNVHZKdGIwxJis+nrmZvUdP8cmgpoHZjHbBJ7D8a+jwJNTtleOX92S02rIiMkpEprhe1xGRuzw492igayb7zFHVRq4lLVmEAh8C3YA6wAARqePB9bJt9c7DHDpx2peXMMYEuB0JJxg5J5aejSrQJOoyf4dzvs3TYepzUOs66PScX0LwZGiQ0cBUoILr9UbgscwOUtXZQEI2YmoBbFLVWFU9DXwL9MzGeTxy6MRp+o34i6d+XImq+uoyxpgAN2TKekIEnu4agM1oD2yGH+6A0pdD7xEQ4slXt/d5ctVSqvo9kAqgqslAipeu31pEVojIFLc5NioCO9z2iXOt84kShfPz+FU1mbZ2D1//vc1XlzHGBLCFWxL4ddUu7u1QjQqB1ow28QjEDHCazfYfCwX8N/GpJwnjuIhE4BS6EZFWwGEvXHspUFlVGwIfAD+51mf04PCCf/qLyGARWSwii/ft25etQO5sG02nWqV59dd1rI0/kq1zGGOCU2qq8sqkNZQvXjDwRqNNTYFxd8OBTXDTGCgZ7ddwPEkYTwATgWoiMg8YAzx8qRdW1SNp072q6mQgTERK4dxRVHLbNRJnTvELnWekqjZT1WalS2eveVlIiPDOjQ0pXiiMh2OWcuJ0crbOY4wJPj8ujWP1ziM80+1yCuUP9Xc455r+KvwzFbq9CdEd/B1N5gnD1d+iI07z2nuBuqq68lIvLCLlRJyuiSLSwhXLAWARUENEokUkP9AfJ2H5VKkiBXivXyNi9x/n5YlrfX05Y0wAOHYqmbenbqBxVAl6NKyQ+QE5adWPMPd/0PR2aH63v6MBPOuHkVa3WJOVE4tIDNAJKCUiccCLQJjrfJ8AfYH7RSQZOAn0V6fqnCwiD+EU2kOBz1U1S9fOrrbVS3F/x2p8NHMz7WqUonug/QIZY7zqoxmb2Hf0FCMHNUVyYGgNj+1cCj8/CFFtoNvbOTLshyd81o1RVQdksn04MPwC2yYDk30RV2Yev6omf8Ue4Lnxq2hUqQSVSubd6RiNyc12JJzgs7lb6NW4Io0DqRnt0d3OsB/hZaBfzg374Qn/tM0KYGGhIbzfvzEIPByzjKSUVH+HZIzxgTcmryNUhKe61vJ3KGcln4LvboHEQzBgLISX8ndE5/AoYYhIRRFpIyId0hZfB+ZPlUoWZkjvBizfcYih0zb6OxxjjBftPpzIK7+sZcrq3dzXsRrliwdIM1pVmPQ4xC2CXp9Aufr+jug8mT6SEpE3gX7AWs72v1Bgtg/j8rvrGpRn7qYoPpm1mbbVI2hfI+cG+DLGeN/W/ccZMXsz45bsJEWVPk0iuTeQRqP9+yNY/g10fBrq+Kyv8iXxpIZxA1BLVU/5OphA88L1dViyLYHHv1vBlEfbU7poAX+HZIzJojXxh/l45mYmr9pFvtAQbmoeyb0dqgVWfXLTnzDtebj8euj4jL+juSBPEkYsTuumPJcwCuUP5YMBTegxfC5PfL+cL+9oEZgDkhljzrNoawIfzdjEjA37KFIgH4M7VOPOdlUoU7Sgv0M71/5N8OMdULo29PLfsB+e8CRhnACWi8ifuCUNVX3EZ1EFkFrlivKf6+vw/E+r+XROLPcGWk9QY8wZqsrMjfv4aMYmFm09SMnw/Pz76poMal2F4oXC/B3e+RIPQ0x/kFCnyO3jObkvlScJYyI50HEukA1sGcW8Tft5e+oGWlaNoFGlEv4OyRjjJiVVmbxqFx/P3MzaXUeoULwgL3WvQ7/mUYHXeztN2rAfB7fAoJ/gsir+jihTnsyH8aWrx3VN16oNqprk27ACi4gwpHcDVsbN4ZGYZfz6SDuKFgzAv1aMyWNOJacwYelORsyOZcv+41QtHc7bfRvQs1FF8ucL3Ec7APz5MvwzDa57F6Lb+zsaj3jSSqoT8CWwFWdgwEoicptr+PI8o3jhMIb1b0S/kX/zfxNWM6x/o8DqGWpMHnL8VDIxC7fz2Zwt7D6SSP2Kxfl4YBOurluO0GCoM674DuYNg2Z3QnNPphcKDJ48khoKXK2qGwBEpCYQAzT1ZWCBqFmVkjzWpQZDf99IuxqluKlZpcwPMsZ4zaETpxk9fyuj52/l0IkkWlUtyds3NqBd9VLB8wfcziUw8WGo3A66vunvaLLEk4QRlpYsAFR1o4jk2ecxD3SuzvzNB3jx5zU0ibqM6mUCu0hlTG6w50gin82J5ZsF2zlxOoUra5flgc7VAnNmvItJG/ajaFm46cuAGvbDE54kjMUiMgr4yvV6ILDEdyEFttAQ4b3+jeg2bA4PxyxjwgNtKBgWoEU1Y4LctgPH+WRWLOOWxJGiSvcG5bm/U3VqlfPfJELZlpToJIvEI3DXtIAb9sMTniSM+4EHgUdwahizgY98GVSgK1usIO/c2IA7Ry9myJT1vNSjbuYHGWM8tm7XET6auZlfV8aTLzSEG5s5ne2iIgKos11WqMIvj8LOxXDTV1Cunr8jyhZPWkmdAt51LcblisvLcmfbaD6ft4U21SK4um45f4dkTNBbvDWBj2ZuZvr6vYTnD+WeDlW5q200ZYoFWGe7rPprOKz8Fjo9B3V6+DuabLtgwhCR71X1JhFZRQZTpKpqA59GFgSe7laLBVsO8NS4ldSPLB44g5gZE0RUlVkb9/HRjM0s3JpAyfD8/OuqmtzaugrFC+eCcuk/f8DvL0DtHtDhSX9Hc0nEmbMogw0i5VV1l4hUzmi7qm7zaWTZ0KxZM128eHGOXjN23zGu/2Au9SoWJ+aeVsHRpM+YAJCSqkzh/QtAAAAgAElEQVRZ7XS2WxN/hPLFCzK4Q1X6Na9E4fw+m6onZ+3/Bz7tAiWi4K6pkD/c3xGdR0SWqGozT/a9YM8WVd3l+ucDqrrNfQEe8EagAePYXkjJXl/EqqWL8GrPeizcksAH0//xcmDG5D6nk1P5btF2rnx3Fg+NXcbJpBTe6tuAWU925o620bknWZw85Az7EZrPGfYjAJNFVnnyX+Yq4Ol067plsC44nUiAER2gzg3QbUi2TtGnaSRzN+3n/T//oXXVCFpWjfBykMYEvxOnk4lZuINPZ8ey+0gi9SoW46OBTbgmWDrbZUVqCoy7Cw5uhVsnOncYucDFahj349xJVBWRlW6bigLzfB1YjilcEur2csaiL98AGt2crdO8ekM9lm0/yGPfLWfyI+25LDy42lcb4yuHTpzmy/nbGD1/CwdPJNEyuiRv9W1A+xpB1Nkuq/54ETb9Ade/B1Xa+jsar7lYDaM4cBnwX8B9gPajqpqQA7FlWbZrGCnJ8HUv2L4A7pgCkdnrxL4q7jC9P55Hp1plAm9SeWNy2N4jiXw2dwvf/L2N46dTuLJ2Ge7vVJ2mlYOss11WLY+Bn+6D5nfDdUP9HU2mslLDuGDCyOCkZYAzbdtUdXsm+38OXA/sVdXzGh2LyEDOPtY6Btyvqitc27YCR3Fm+Ev29M1cUtH7+AH4tJNTyxg8y+mJmQ2fzYnltV/X8UrPutzaukr2YjEmiLl3tktOTaV7wwrc36kal5cr5u/QfC9uMXxxLVRqAYMmQGjgt/LKSsLwZPDB7jh9MCoAe4HKwDogs95qo4HhwJgLbN8CdFTVgyLSDRgJtHTb3llV92cWn9eER0D/sTDqavh+ENz2C+TL+gx7d7WLZt6m/bz26zqaVS5JnQp54H8SY3A62308czOTVsaTLySEvs0iubdDVSpHBH+x1yNH4l3DfpSDG78MimSRVZ6M//sa0ArYqKrRQBc8qGG4RrO94KMrVZ2vqgddL/8GIj2IxbfK1YeeH8KOBTD5Sad3ZhaJCO/c2JAShcJ4KGYpJ04n+yBQYwLHkm0J3DV6Ed2GzeHPdXu4p31V5j7dmTd61c87ySLppJMsTh2FATHOH6C5kCetpJJU9YCIhIhIiKrOEBFvD7F4FzDF7bUC00REgRGqOtLL17uwer1h9yqY+y6Ub5itoYcjihTgvX6NGDhqAS9NXMNbfRv6IFBj/EdVmf3Pfj6csYmFWxK4rHAYT1xVk9tyS2e7rEgb9iN+KfT7Bsrm3qGCPEkYh0SkCM4YUt+IyF7Aa382i0hnnITRzm11W1WNd9VNfheR9Reaf0NEBgODAaKivNR07YrnYc9qmPIUlKkNldtk+RRtqpfigU7V+HDGZtrVKE2PhhW8E5sxfpSSqvy2ejcfzdx0prPdf66vw4AWuaizXVbNfx9Wfgedn4fa1/s7Gp/KtOgtIuFAIs7AgwOB4sA3qnog05OLVAEmZVT0dm1vAEwAuqnqxgvs8xJwTFXfyex6Xu3pffIQfNbFmXN38EwonvUnZkkpqfQb8Rcb9xxj8iPtg3fgNGOAhOOnufnTv1m/+yhVS4VzX8dq3NA4CGa286WN02DsTVCnJ9w4GoKwZaRXenqnUdXjqpoCFAZ+Ab4mg7GlskpEooDxwCD3ZCEi4SJSNO3fwNXA6ku9XpYVKuEUwdOGJE46meVThIWGMKx/Y0Tg4W+XkZSS6oNAjfG95JRUHhq7lNj9xxnWvxG/P9GRm5pXytvJYt9Gp3NeuXpww0dBmSyyKtP/2iJyr4jsAVYCi3Hmwsj0z3gRiQH+AmqJSJyI3CUi94nIfa5dXgAigI9EZLmIpJ2zLDBXRFYAC4FfVfW3LL8zbyhdC/p8CruWO88os1EEr1SyMG/2acCKHYd4Z9qGzA8wJgD9d8p65m8+wH971adno4q5r2d2Vp086Br2Iz/0j8kVw354wpOHjv8G6ma1iauqDshk+93A3RmsjwUCp0pcqxt0/j+Y8bpTBG/9YJZPcW398tzcMooRs2JpW60UHWqW9kGgxvjGhGVxjJq7hdvbVKFPU/83ZvS7lGT48U44tN1pfl8i70zV7Mn95GbghK8DCWjt/w21u8O052HzjGyd4oXr61CzbBGe+H4F+46e8nKAxvjG6p2HeWbcKlpGl+T/rqvt73ACwx8vwubpTi/uyq39HU2O8iRhPAvMF5ERIvJ+2uLrwAJKSAjc8DGUqgU/3gEJW7J8ioJhoQy/uQlHE5N44vvlpKZechnIGJ9KOH6ae79aQkR4fj4c2ISw0Dxcr0izfKwzGVKLe6Hpbf6OJsd58hswApiO07luiduStxQo6gxRrOrqoHMsy6eoWbYoL3Svw5x/9jNyTqwPgjTGO9KK3PuOneKTQU0pVSTrox7kOjsWObXM6A5wzev+jsYvPKlhJKvqEz6PJBiUrAp9P4dv+sLPDzjd/7PYMuLmFlHM/Wc/70zdQMvokjSOyuUDsZmgNMRV5B56Y0MaRJbwdzj+d3gnfDcQilXItcN+eMKTO4wZIjJYRMqLSMm0xeeRBarqXeDKl2HtzzAn6yNRighDejegbLGCPPLtMo4kZm/iJmN85eflO/nMitxnJZ2Eb2+G08dhwLfOlAh5lCcJ42ZcdQzOPo7K2XlQA02bh6H+TTD9Ndg4NcuHFy8cxvsDGhF/KJH/m7AaT0cMNsbXVu88zFM/rrQidxpVmPgw7FoBvT91Rn7IwzzpuBedwVI1J4ILWCLQ431nwqVxdzsdeLKoaeWSPH5lDX5ZEc8Pi+N8EKQxWZNW5C5pRe6z5r0Hq35whgu6/Fp/R+N3F/yNEJErXD97Z7TkXIgBKqyQM9BYaH7ndjXxcJZPcX+n6rSpFsGLE9ewae9RHwRpjGfci9wjrMjt2PAb/PEy1O0N7f/l72gCwsX+hOjo+tk9gyV3j7DlqRKV4KYxcHALjB8MqVkb+iM0RPhfv0YUyh/KQ2OXkZiU4qNAjbm4tCL3G73qW5EbYN8G5+lB2pQHeWDYD09cMGGo6ouun3dksNyZcyEGuCptoesQ2Pib0xs8i8oWK8jQGxuyfvdR/jt5nQ8CNObi3Ivcfa3IDScSnGE/wgo6c1vkt0FD03gyltRXrvm9015XFpE/fRtWkGl+NzQeBHPegTU/ZfnwzpeX4a520Xz51zamrtntgwCNydjqnYd5etxKWliR25GS7HTOPbQD+n2drVGqczNPqlpzgQUicq2I3AP8Drzn27CCjIgzTEBkc/jpAdizJsuneKprLepVLMZTP64k/lDWR8Y1JqvSityXFc7PR1bkdvz+H4idCdf/D6Ja+TuagONJK6kROIME/gy8AnRQ1V98HVjQyVcAbvrK6REeM8C5rc2CAvlC+WBAE5JTUnns2+Uk21DoxoesyJ2BpV/B3x9By/uhySB/RxOQPHkkNQj4HLgVGA1MFpHAGU02kBQrD/2/gaO7nNvalKxNTBhdKpxXb6jHwq0JfDB9k4+CNOZskfv1G+pZkRtg+wKY9DhU7QRXv+bvaAKWJ/egfYB2qhqjqs8C9+EkDpORyGbO7WzsTGdUyyzq3SSS3o0r8sH0f/g7NtNJDY3JsrQi922tK3Njs7wzNPcFHY6D725x6hV9v4DQPDrVrAc8eSR1g6rudXu9EGjp06iCXeNboMVgZ1TLFd9l+fBXbqhH5YhwHvt2OQePn/ZBgCavWhN/tsj9/PV1/B2O/50+4fSjSjqZ54f98IQnj6QiRWSCiOwTkT0iMg4okwOxBbdr3oDK7eCXRyB+WZYOLVIgHx8MaMyB46d48seVNnSI8YqE46cZPMYpcn94sxW5UYWfH4RdK6HPZ1Dmcn9HFPA8+Y35ApgIlAcq4szr/YUvg8oVQsPgpi8hvLQzHPqxvZkf46ZexeI80602f6zbw5i/tvkoSJNXJKek8nCMa7jyW5pSuqgVuZn7LqwZD11egFpd/R1NUPAkYZRW1S9UNdm1jAZsjlFPhJdyiuAnEuD72yA5a4+X7mxbhSsuL8Prv65jTXzWhx4xJs2bv61n3ianyN2wUh4tciefgm3zYdZb8GV3+PNVqNcX2j3u78iChicJY7+I3CIioa7lFsCqsZ4q3xB6Doft8+G3Z7J0qIjwdt8GlCgcxsMxyzhxOmutrowBp8j96Zw8WOROPgVb58HMN2H09TAkCr7oBjPegJMHoe0j0OMDG/YjCzxpDnAnMBz4H6A4w5x7NDSIiHyOM+7UXlWtl8F2AYYB1+LMG367qi51bbsNeN6162uq+qUn1wxI9fs6wyPPd41w2/R2jw+NKFKA9/o3YuBnC3jx5zW8faO1aDaeO1PkrpIHitzJpyBuMWydC1vnQNwiSE4EBMrVg2Z3QpV2ENXaitvZdNGEISKhQB9V7ZHN84/GSTZjLrC9G1DDtbQEPgZauiZoehFohpOklojIRFU9mM04/O/Kl5we4L/+G0rXhijPG5q1qVaKBztVZ/iMTbSrUYqejSr6LEyTe5xT5M6NPbmTEmHnEksQOeiiCUNVU0SkJ87dRZap6mwRqXKRXXoCY9RpBvS3iJQQkfJAJ+B3VU0AEJHfga5ATHbiCAghodB3FIzsDN8PgsEznekePfTYlTX4K/YA/zdhNY0qlaByRLjPQjXBz73I/cO9rXNHkTspEXam3UHMhR0LIeUUToKoD83uchJE5dZQyKY+9gVPHknNE5HhwHfA8bSVaY+OLlFFYIfb6zjXugutD26FLnNGv/zsSqej0O2TnRExPZAvNIRh/Rtx7bA5PBKzjB/ua0P+fLnsL0bjNWlF7rf7NgjeIndmCaL53ZYgcpgnCaON6+crbusUuMIL18+o2qQXWX/+CUQGA4MBoqKivBCSj5WpDb1GOBPK//pElsbaj7ysMG/2acD93yxl6LQNPHutjS5qzpdW5L412IrcliACXqYJQ1U7+/D6cYD7b3QkEO9a3ynd+pkZnUBVRwIjAZo1axYcPdxqXw8dn4FZQ5xWVC3v9fjQbvXLM7BlFCNmx9Kmeik61rQWzuYs9yL3fwK9yJ2U6NQd0hJE3KKzCaJ8A2hxj6sG0coSRIDINGGISAROAbodzl/5c4FXVNUbTWsnAg+JyLc4Re/DqrpLRKYCb4hI2m/J1cCzXrhe4Oj4NOxeBb8969x1RHfw+ND/XF+HxVsP8q/vlzP50faUKerZYy2Tux10G648IIvcHieI1lAoSB+j5XKS2bATroLzbOBr16qBQCdVvTLTk4vE4NwplAL24CSeMABV/cTVrHY4TkH7BHCHqi52HXsn8JzrVK+raqa9y5s1a6aLFy/ObLfAkXjEqWcc3+cUwS+r7PGhG/ccpcfwuTSvUpIv72hBSIi1Jc/LklNSue2LhSzaepAf7m0dGHWLCyUICYFyDZzkYAnC70Rkiao282hfDxLGElVtmm7dYk8vkJOCLmEA7N8En14BJaLgrmlZmg5y7ILtPDdhFU93vZz7O1XzYZAm0L0xeR0jZ8fyVt8G3OSvukXSyQwSxOl0CaK96xGTJYhAkZWE4UnRe4aI9Ae+d73uC/ya3eBMOqWqQ9/P4Zu+zkBofT/3uAg+oEUl5m3az9BpG2hVtSSNo+w5b1708/KdjJwdy62tK+dsssgsQbQYbAkil/HkDuMoEA6kTQEXwtnmtaqqxXwXXtYE5R1Gmrn/gz9ecjr4ZWFsm8Mnk7h22BxEYPKj7SlWMMxXEZoAtCb+MH0+nk+DiiX45p6Wvq9bxC2Bf6aenyDKN4TKbS1BBCGv3mGoatFLD8lkqu1jzjDLf7wMZetDjUxLRAAULxTG+wMac9OIv3hu/Co+GNAYsbFx8oS0IneJQjlU5F40ymkKnpYgWt57NkEULO7ba5uA4NHUUiLSm7OtpOao6k8+jSovEnEGKdz/D/x4JwyeARGe1SWaVr6MJ66qydtTN9C+Rin6NQ+C/ijmkjg9uZex9+gpvs+JntyLPoNf/wU1roHeI+0OIo/yZAKlj3CmZV0FrAbuE5EPfR1YnpQ/3BkOPSQUYgbAqaMeH3p/x2q0rR7BixPX8P6f/7Aj4YQPAzX+9tbUDczdtJ/XbqhHI1+3iFr4qZMsanaFfl9ZssjDPKlhrAHqucZ7QkRCgFWqWjcH4suSoK5huNsyG8bc4Pof9GsI8exRw96jiTzx3QrmbtoPQMvokvRpGkm3euUoarWNXOPn5Tt59NvlDGpVmVdvOG8QaO9a+ClM/jfU7OZMCJYvF4xJZc6RlRqGJ99EGwD3ZxyVgJXZCcx4KLqDM8Xrhl9h1pseH1amaEG+vrslc5/uzL+vrsneo6d46seVNH/9Dx77dhmzN+4jJTU4OsObjK2NP8LT41bSvMplvu/JvWCEkyxqXQs3jbFkYTy6w5gFNAcWulY1B/7G1VLqEoY+97pcc4cBznzDPz0AK8ZCv2+c4USyfApl2Y5DjFsSxy8r4jmSmEzZYgW4oXFF+jSJpGZZa88QTA4eP0334XNJTlF+ebidb+sWf38Cvz0Nta6DG0dDvvy+u5bxK2933Ot4se2qOisLsflUrkoY4PSUHX0t7NsAd//hDCGSTaeSU5i+bi/jlsYxc8M+klOV+hWL07tJRXo0rEBEEfvrMZAlp6Ry+xeLWLglge/va+3busXfHzuzQ15+PfT9wpJFLufVhJHBydsCN6vqg9kJzpdyXcIAOBIPIzpCgSJwz3SvDMK2/9gpflkRz7ilcazeeYR8IUKnWmXo06QiV9QuQ4F8oV4I3HjTfyevY8TsWN7q04Cbmvuwc95fH8LU55xkceNoCLXaV27n9YQhIo2Am4GbgC3AOFUdfklR+kCuTBgA2xfA6Ouc2sbAH5xWVF6yYfdRxi+NY8Kynew9eorihcLo3rA8vZtE0rhSCevTEQAmrojnkZhlvi9yzx8O0/4PavdwRhywZJEneCVhiEhNoD8wADiAM4HSv1XV8xHycliuTRgAS0bDL49C20fhqlcy3T2rUlKVuZv2M35pHFPX7CYxKZWqpcLp3aQivZpEUrFEIa9f02RubfwRen88j/oVi/PN3a18N2nW/A9g2vNQpyf0GWXJIg/xVsJIBeYAd6nqJte6WFWt6rVIvSxXJwyASU/A4lHO/9D1+/rsMkcTk5iyajfjlsaxYEsCAK2rRtC7SUW61S9PkQIe9fc0l8i9yD3x4ba+G8Z+3jD4/QWocwP0+cySRR7jrYTRC+cOow3wG/At8JmqRnsrUG/L9Qkj+TSM6QHxy+Guqc7wDD62I+EEE5btZPzSOLYeOEGhsFC61itHnyaRtK4WQagNq+4T7kXu7+5t5buBJee+B3+8CHV7Qe/PINT+GMhrvN1KKhy4AefR1BXAl8AEVZ12qYF6W65PGADH9sLITs54PoNnQnipHLmsqrJ0+0HGLd3JJFcT3XLFCnJD44r0bVqR6mWsia435UiRe8678OfLULc39P7UkkUe5bNWUiJSErgR6Keq3pjT26vyRMIA2LkUvugGkc1h0IQcf4SQmJTCn+v2Mn5pHDNdnQEbRBanT5NIujesQMlwa4Z5KdKK3Le0iuK1G+r75iJzhsKfr0C9PtBrpCWLPMynzWoDWZ5JGAArvoMJg6HFvXDtW34LY9/RU0xcEc+4JXGs3eU00e18eRn6NInkisvL+K5Im0vlSJF79jsw/VWofyPc8IklizzO2xMomUDUsB/sXgl/DXfmQ258i1/CKF20AHe1i+audtGs23WE8Uvj+Gl5PL+v3UOJwmH0aFiB3k0iaRhZ3JroZuLg8dPc+/XiM8OV+yRZzHobZrwG9W+CGz62ZGGyxO4wgllKMnzTB7bNhzumQGRgzJqbnJLK3E37Gbd0J9PW7OZUcirVSofTu0kkvRpXpII10T1Pckoqd4xexIJYHxa5Z70FM16HBv2cZOHF/jwmeNkjqbzkRIJTBE857RTBi5bzc0DnOpKYxJRVuxi3ZCcLtyYgAm2qRdC7cSRd65Uj3JroAmeL3G/2qe+b+UxmDoGZ/4UG/eGGjyxZmDMCJmGISFdgGBCK0yR3SLrt/wM6u14WBsqoagnXthScOTgAtnsyyGGeTBgAe9bAZ1dC2Xpw+6SAHVV0+4ETjF8Wx/ilO9mecILC+d2a6FaNICSPNtH9ZUU8D/uyyD3jvzBrCDS82Zmky5KFcRMQCUNEQoGNwFVAHLAIGKCqay+w/8NAY1W90/X6mKoWyco182zCAFjzE/xwGzS5Fbq/78zgF6BUlSXbDjJuaRyTVuzi6KlkKhR3muj2bhJJ9TJZ+s8e1NbGH6HPx/OpW6EYY+/xcpFb1bmrmPUmNBoIPT6wZGHOEygJozXwkqpe43r9LICq/vcC+88HXlTV312vLWFk1Z+vwpx3nDuN5ndDg5ucWfwCWGJSCr+v3cP4pXHM/mc/KalKw0ol6NOkItfVL5+rR9E9ePw0PT6cy+nkVH55uJ13e3Krwow3YPZb0OgW6PG+JQuToUBJGH2Brqp6t+v1IKClqj6Uwb6VcebYiFTVFNe6ZGA5kAwM8WQe8TyfMFJTYfk3sHAE7F4FBYpDo5ud5FGqur+jy9Teo4lMXB7Pj0viWL/bmZ42Ijw/VUuHU7VUEaJLh1O1VDhVSxchqmThoG6ym5Kq3P7FQt8UuVWd4vbst53Wc90/8HjWRpP3BEqz2oyeiVwoO/UHfkxLFi5RqhovIlWB6SKySlU3n3cRkcHAYICoKB8UC4NJSAg0GeR8SexY4EyvuegzWPAxVLsCmt8DNa8J2L80yxQtyN3tq3J3+6qsiT/M3H/2E7vvOFv2H+fP9XvYv/j0mX1DQ4RKlxWiaukiRJcKP5NUqpYOp0zRAgHfhPetqeuZ889+3uxT3/vJYvqrTse8JrfC9cMsWRiv8WXCiMOZzjVNJBB/gX37A+fMr6Gq8a6fsSIyE2gMnJcwVHUkMBKcO4xLjjo3EIGoVs5y9A1YOgYWfw7fDoDiUdDsDufLJIeGFcmOuhWKU7dC8XPWHT6ZxJb9x4ndd8z18zib9x1j/ub9JCalntmvSIF8RJcKP5NIokuFU82VWAKhVdYvK+IZMSuWW1pFebdFlKrTe3vuu9DkNrj+PUsWxqt8+UgqH07RuwuwE6fofbOqrkm3Xy1gKhCtrmBE5DLghKqeEpFSwF9AzwsVzNPk+UdSF5OSDBsmw8KRsHUOhBaAer2du47Ipv6O7pKkpiq7jiSel0i27D/OzkMncf8VL1uswJk7kbREUrV0OBVLFCJfqO+/XH1W5FaFP16Cee9B09vhuv9ZsjAeCYgahiuQa4H3cJrVfq6qr4vIK8BiVZ3o2ucloKCqPuN2XBtgBJAKhADvqeqozK5nCcNDe9c7j6pWxMDpY1ChsZM46vWGsNzVqS4xKYVtB04Qu+8Ysa5kErv/GLH7jnP4ZNKZ/cJChcoRTo0kunQ41dxqJiXD83vlEdehE85w5V4vcqs6I87OGwbN7oRrh1qyMB4LmISR0yxhZFHiEVj5nVPr2L8BCpV0aiDN7oLLAnaeLK9QVQ6eSDo3kbjuSrYdOMHplLOPuIoXCjvzeKuaW82kSkQ4BcM8qwe5F7m/vbcVTbxVt1B15rKY/77z3+3adyxZmCyxhGGyRtV5TLXwU1j/K2iqUxxvfo9TLM9jX0DJKansPHTyvEQSu+84u48kntlPBCoUL3ReIokuFU6F4oXO6Yj43ynrGDErliG969O/hZfqFqrOLHl/DXdawl37TkD3vzGByRKGyb7DO2HJF7DkSzi+F0pWdb6MGt0MhXw0iU8QOX4q2Uke+4+zxe3xVuy+Yxw/fbaRX8GwEKpEOAmkeKH8xCzczsCWUbzey0s9uc9JFvfAtW9bsjDZYgnDXLrk07BuonPXseNvyFcIGtzofDmVb+Dv6AKOqrLv6Ck2u5oBpz3q2rL/ONsTTtCmWgSjbmvunSK3Kkx9Dv7+yBnevtublixMtlnCMN61ayUs+hRW/gDJJ6FSK2hxD9TuAflssqTMJKWkki9EvNM3RBV+e9bpW9PyPug6xJKFuSSWMIxvnDwIy8c6dx0Ht0B4GacJZ7M7oFgFf0eX+6nCb8/Agk+g5f3Q9b+WLMwls4RhfCs1FTZPd+46Nk515he//DrnrqNKe/sS8wVVmPK0M+xLqwfhmtftczZeEShDg5jcKiQEalzpLAe3Or3Il45xah6lL3eK5A37Q4Gi/o40d1CFyU86Cbr1Q3D1a5YsjF/YHYbxjqSTsHq886UWvwzyF3WSRot7oHQtf0cXvFRh8r+djpZtHoarXrVkYbzKHkkZ/4pb4gxBsma8MxNgdAendVWta20O6axITXWSxeJR0OYRuOoVSxbG6yxhmMBwfP/ZgQ8P74BiFaHpHdD0NihSxt/RBbbUVPj1CadPTNtH4cqXLVkYn7CEYQJLagps/M1pXRU7A0LCoE5PaDEYKrWwL8L0UlPh18dhyWho9zh0edE+I+MzVvQ2gSUk1GlFdfl1sH+T8zx++VhY/SOUq+88rqp/I+Qv7O9I/S81FSY9Bku/hHZPQJcXLFmYgGF3GMY/Th2DVd/Dws9g7xooWBwaD3JGW42o5u/o/CM1FSY96jzGa/9vuOJ5SxbG5+yRlAkeqrD9L+dx1bqJkJoM1a907jpqXBWwswN6XWoq/PIwLPsaOjwJnf/PkoXJEfZIygQPEajcxlmO7nae2y/+AmL6QZGyUL4hlK0LZes5PyOqQ2iYv6P2rtRUmPgwLP8aOj4NnZ61ZGECkt1hmMCTkgTrJ8H6ybB3LezbAKmuyY5C80OpWq4kUvdsMilSJji/ZFNTXMniG+j4DHR+1t8RmTzG7jBMcAsNg7q9nAWckXMP/AN71pxdtsyGld+ePaZwKShb5+ydSNm6Tq/zQJ5BMDUFfn4IVox17io6PZP5Mcb4kSUME/jy5T+bBNydSHBLIqudu5HFXzgj6oIzxlVEdSiTLpGUiPL/3UhqCvz0gJP0Oj0HndGe1dMAAAbhSURBVJ72bzzGeMAShglehUtCdHtnSZOa4oxvtWf12WSyazms/ensPgWKuZJInbOPtMrUgYLFcibu1BT46X5netzO/wcdn8qZ6xpziSxhmNwlJNRplhtRzekcmObUUdi7/txEsmqc0ws9TYmos3ciaXclJat6dziT1BSYcJ/TpPiK550WUcYECZ8mDBHpCgwDQoHPVHVIuu23A28DO12rhqvqZ65ttwHPu9a/pqpf+jJWk8sVKAqVmjtLGlU4svPsI620RLJxKqhrutV8BZ1ayJlHWq5EEl4q6zGkJMNP98GqH+CK/0CHf3vnvRmTQ3zWSkpEQoGNwFVAHLAIGKCqa932uR1opqoPpTu2JLAYaAYosARoqqoHL3ZNayVlvCIpEfZvOLfIvmeNM8d5miJlz22lVaaOMypvvgIZnzMlGSYMhtXjnKE+2j+RM+/FmEwESiupFsAmVY11BfUt0BNYe9GjHNcAv6tqguvY34GuQIyPYjXmrLCCTv+P8g3PXX9s77kJZO8aWDASUk4520PyQUSNcxNJ2TpQpNzZZHHlS874UMYEIV8mjIrADrfXcUDLDPbrIyIdcO5GHlfVHRc4tqKvAjXGI0XKOEu1zmfXpSRDwuZzH2ntWOiMk5UmXyGn5daVL0O7x3I+bmO8xJcJI6N2i+mff/0CxKjqKRG5D/gSuMLDY52LiAwGBgNERUVlP1pjsiM0n/MoqnQtqNfn7PqTh2DvurPNfSNbQKMB/ovz/9u7m1Cp6jiM49/Hl0gtUVBCuqIF4aZFihghhGRJkkhLg4La1KIXw0VUm8hWbaJdYGoYmVGaIBG9QEm16EXNMNNNYXSzupaE3TZhPS3mSINEHc/M+L9z7vOB4c5c5tzz+3G595nzO/+ZE9EHgwyMUWBh1+MR4GT3E2z/0vXweeDprm1Xnbft/n/bie0twBbonMPopeCIvpkxBxbd0LlFtMSUAf7sz4BrJF0l6RJgA7Cv+wmSFnQ9XA8cq+6/DayRNFfSXGBN9b2IiChkYEcYts9KeoDOP/qpwHbbRyVtBg7Y3gc8JGk9cBY4DdxdbXta0lN0Qgdg87kT4BERUUY+fDAiYhK7kGW1gxxJRUREiyQwIiKilgRGRETUksCIiIhaEhgREVFLq1ZJSToFfNtw83nAz30sp6S29NKWPiC9TERt6QN662WR7fl1ntiqwOiFpAN1l5ZNdG3ppS19QHqZiNrSB1y8XjKSioiIWhIYERFRSwLjH1tKF9BHbemlLX1AepmI2tIHXKRecg4jIiJqyRFGRETUksCIiIhaJn1gSNouaUzSl6Vr6YWkhZLel3RM0lFJG0vX1JSkSyV9KumLqpcnS9fUC0lTJX0u6Y3StfRC0glJRyQdljTUHwstaY6k3ZKOV38zQ3mlK0lLqt/HudsZSQO7DvCkP4dRXU98HHjR9rWl62mquhjVAtuHJF0OHARut/1V4dIumCQBs2yPS5oOfARstP1x4dIakbQJWA7Mtr2udD1NSToBLLc99G92k7QD+ND21uoCbzNt/1q6rl5Imgp8D1xvu+kbmP/TpD/CsP0BnYs3DTXbP9g+VN3/jc7VC68sW1Uz7hivHk6vbkP5ykbSCHAbsLV0LdEhaTZwI7ANwPYfwx4WldXA14MKC0hgtJKkxcBS4JOylTRXjXEOA2PAu7aHtZdngUeAv0oX0gcG3pF0UNK9pYvpwdXAKeCFalS4VdKs0kX1wQZg1yB3kMBoGUmXAXuAh22fKV1PU7b/tH0dMAKskDR040JJ64Ax2wdL19InK20vA9YC91fj3GE0DVgGPGd7KfA78GjZknpTjdXWA68Ncj8JjBap5v17gJ22Xy9dTz9Uo4L9wK2FS2liJbC+mv2/Atwk6aWyJTVn+2T1dQzYC6woW1Fjo8Bo11HrbjoBMszWAods/zTInSQwWqI6UbwNOGb7mdL19ELSfElzqvszgJuB42WrunC2H7M9YnsxnXHBe7bvLFxWI5JmVYspqMY3a4ChXFlo+0fgO0lLqm+tBoZucch57mDA4yjoHJpNapJ2AauAeZJGgSdsbytbVSMrgbuAI9XsH+Bx228WrKmpBcCOatXHFOBV20O9JLUFrgD2dl6XMA142fZbZUvqyYPAzmqU8w1wT+F6GpM0E7gFuG/g+5rsy2ojIqKejKQiIqKWBEZERNSSwIiIiFoSGBERUUsCIyIiaklgRAyQpPH/f1bEcEhgRERELQmMiIioJYERERG1JDAiIqKWBEZERNSSwIiIiFoSGBGDNVPSaNdtU+mCIprKp9VGREQtOcKIiIhaEhgREVFLAiMiImpJYERERC0JjIiIqCWBERERtSQwIiKilr8BAWormo7XesAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5a8c3a6da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot for the computation time of the exact case.\n",
    "L_exact_comp = np.arange(1,len(Ex)+1,1)\n",
    "L_comp = np.arange(1,20,1)\n",
    "plt.plot(L_comp,App1[0:19,1],label='Approx 1')\n",
    "plt.plot(L_comp,App2[0:19,1],label='Approx 2')\n",
    "plt.plot(L_exact_comp,Ex[:,1],label='Exact case')\n",
    "plt.legend()\n",
    "plt.xlabel('L')\n",
    "plt.xticks(L_comp)\n",
    "plt.ylabel('Optimizer Time (s)')\n",
    "plt.title('Computation time for Exact Solution vs L')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "#Plot for the computation times of the approximations.\n",
    "\n",
    "plt.plot(L_comp,App1[0:19,1],label='Approx 1')\n",
    "plt.plot(L_comp,App2[0:19,1],label='Approx 2')\n",
    "plt.legend()\n",
    "plt.xlabel('L')\n",
    "plt.xticks(L_comp)\n",
    "plt.ylabel('Optimizer Time (s)')\n",
    "plt.title('Computation time for approximations vs L')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "#Error in the approximation vs L\n",
    "L_err = np.arange(1,tExact+1,1)\n",
    "error1 = 100*((np.exp(App1[0:tExact,0])/np.exp(Ex[0:tExact,0]))-1)\n",
    "error2 = 100*((np.exp(App2[0:tExact,0])/np.exp(Ex[0:tExact,0]))-1)\n",
    "plt.plot(L_err,error1,label='Approx 1')\n",
    "plt.plot(L_err,error2,label='Approx 2')\n",
    "plt.legend()\n",
    "plt.xlabel('L')\n",
    "plt.xticks(L_err)\n",
    "plt.ylabel('Approximation error (%)')\n",
    "plt.title('Approximation error vs L')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More results\n",
    "\n",
    "Results from an extended run of the above experiment are shown below. They were performed on a desktop with 15.6 GB of RAM and an Intel$^{\\circledR}$ Core$^{\\text{TM}}$ i7-6770HQ CPU @ 2.6 GHz $\\times$ 8.\n",
    "\n",
    "![](fullrun1.png)\n",
    "![](fullrun2.png)\n",
    "![](fullrun3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}