Table of Contents
\n", "- 1 Introduction - Inversion & Optimisation
- 2 Linear algebra review
- 3 Some inversion (parameter estimation) examples
- 4 Summary
\n", "\n", "At this point we can think of $f$ in an abstract way as some mapping - it could be a mathematical operator, it could be the output of a complex numerical model (which could be very expensive to execute for each $x$), it could be a [neural network](https://en.wikipedia.org/wiki/Neural_network), ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can learn a lot about solution methodologies, pitfalls to be aware of etc, by considering $f$ to be a linear mapping, in which case we can think of $f(x)$ as the matrix-vector multiplication $A\\boldsymbol{x}$. Note that in this case I am thinking of $\\boldsymbol{x}$ as well as $\\boldsymbol{y}$ to be vectors, and will use bold font (or underlined) symbols to represent these vectors.\n", "\n", "This linear case is what my lectures will focus on.\n", "\n", "
\n", "\n", "Rhodri's lectures will focus more on practical examples where $f$ represents a complex numerical model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Forward vs inverse problems\n", "\n", "
\n", "\n", "The **forward problem** can be thought of as: given $x$ what is $y$?\n", "\n", "\n", "The **inverse problem** is the opposite: given $y$ what is $x$?\n", "\n", "
\n", "\n", "[When posing an inverse problem we will often swap the sides of the equation, i.e. $f(x)=y$, so that *known* information (e.g. given data) is on the right and all of the *unknowns* are on the left - cf. the linear case we will see later $A\\boldsymbol{x}=\\boldsymbol{b}$].\n", "\n", "
\n", "\n", "The first is in principle an easy task to achieve as we have access to $f$ (although $f$ may be very expensive and time consuming to evaluate), the second is hard as we do not have ready access to its inverse.\n", "\n", "This inversion problem has no general solution, and solving the inverse problem may be trivial, easy, difficult or impossible depending upon the form of $f$ and the value of $y$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A simple scalar linear example\n", "\n", "For example, the solution to a linear version of this problem, that is for the problem\n", "\n", "$$y = ax$$\n", "\n", "given an $x$ what is $y$, is \n", "\n", "$$ x = \\frac{1}{a}y$$\n", "\n", "where the operation of multiplying by $1/a$ is clearly the inverse of multiplying by $a$.\n", "\n", "
\n", "\n", "Here we can call \"*multiplying by a*\" our **forward model** and \"*multiplying by* $1/a$\" is then the corresponding **inverse model**. \n", "\n", "
\n", "\n", "This solution to the iversion problem is of course trivially simple at first glance. \n", "\n", "However, note that some more thinking reveals some subtle features that we will see are important in the more general, complex case.\n", "\n", "In particular, our naive solution will not work, or is not defined, if $a = 0$; in that case, there is no solution, unless $y=0$ in which case there are infinitely many $x$ solutions.\n", "\n", "
\n", "\n", "We could chose to write the inverse equation as\n", "\n", "$$ x = (a)^{-1} y$$\n", "\n", "which makes the inverse relationship explicit, that is the use of the inverse operator. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A simple scalar nonlinear example\n", "\n", "A nonlinear example is\n", "\n", "$$ y = \\cos(x) $$\n", "\n", "which of course has the solution \n", "\n", "$$x = \\cos^{-1}(y) $$\n", "\n", "but suppose we didn't know or have a way to compute $\\cos^{-1}\\equiv \\arccos$?\n", "\n", "
\n", "\n", "For this case, assuming that $x$ and $y$ are real numbers, notice that there is no solution if $|y|>1$, and there are infinitely many solutions otherwise. \n", "\n", "Note that infinitely many solutions does not imply that any real number will be a valid solution. \n", "\n", "Despite the lack of a unique solution, the solutions that we obtain are still meaningful and are likely to be useful in practical problems. \n", "\n", "In many real problems, we are likely to have additional *a priori* information that will allow us to select one or more solutions from the infinite collection. \n", "\n", "For example, we may know that possible solutions are constrained such that $0\\leq x\\leq \\pi$ , which together with our forward model ($y=\\cos(x)$), is sufficient to provide a unique solution (provided $|y|\\le 1$). \n", "\n", "\n", "
\n", "\n", "\n", "In general, the solution to our original equation ($y = f(x)$) can be written as\n", "\n", "$$ x = f^{-1}(y)$$\n", "\n", "where\n", "\n", "$$f\\left[f^{-1}(y)\\right] = y$$\n", "\n", "for all values of $y$ for which the inverse is defined.\n", "\n", "There is no general solution to this problem. It is easy to solve (in the scalar case) if $f$ is a linear function, and a much wider range of problems will have solutions we can find robustly (but maybe not cheaply) as long as we are able to begin from a reasonably good guess. \n", "\n", "
\n", "\n", "Newton's method we saw in the Numerical Methods module iteratively solves the problem $\\hat{f}(x)=0$ (I'm introducing a hat here to emphasise that this is a different function to the $f$ above, unless $y=0$ of course) from a starting guess, and can be used for multi-dimensional (i.e. vector valued) problems as well as scalar problems." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A reminder on Newton's method\n", "\n", "#### For finding roots of equations\n", "\n", "Recall that given a problem written in the form\n", "\n", "\\begin{equation}\n", "\\hat{f}(x)=0\n", "\\end{equation}\n", "\n", "The hat here is introduced as we're using $\\hat{f}$ in our definition of a problem written as a root-finding task.\n", "\n", "[What's the link between this $\\hat{f}$ and the $f$ we introduced to define our forward problem at the start?]\n", "\n", "
\n", "\n", "Newton's method is a powerful root-finding algorithm (or \"nonlinear solver\") which finds a solution $x$ (termed a root as the RHS of this equation is zero) via the iteration\n", "\n", "\\begin{equation}\n", "x_{i+1} = x_i - \\frac{\\hat{f}(x_i)}{\\hat{f}'(x_i)}\n", "\\end{equation}\n", "\n", "\n", "We wrote our own code to solve this in your Numerical Methods module from year 2, but let's remind ourselves how to use the in built SciPy function:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.567143290409784\n", "0.5671432904097835\n" ] } ], "source": [ "# An example f hat function we want to find a/the roots of\n", "def fh(x):\n", " return x - np.exp(-x)\n", "\n", "# its derivative\n", "def dfhdx(x):\n", " return 1 + np.exp(-x)\n", "\n", "x0 = -1. # initial guess\n", "print(sop.newton(fh, x0, dfhdx))\n", "\n", "\n", "# and if you don't have the derivative function it will use a quasi-Newton approach:\n", "print(sop.newton(fh, x0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### For inversion\n", "\n", "How do we turn our inversion problem into something Newton's method is appropriate to solve?\n", "\n", "Suppose we want to solve the inversion problem from above where our forward model is described by the equation\n", "\n", "$$ y = \\cos(x) $$\n", "\n", "and we have a given output of the model $y=y_d$ (the subscript $d$ is added to emphasise given data).\n", "\n", "Given these two pieces of information, what is $x$?\n", "\n", "Let's first plot our situation for a particular value of $y_d$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
\n", "\n", "In practice, we will also seldom have an analytical form for $F$.\n", "Also in general $F$ will not be a matrix; rather it is a (typically non-linear) function that turns one vector into another. You can think of it also as a large complicated numerical model/computer code that calculates $\\boldsymbol{y}$ for us given a particular value for $\\boldsymbol{x}$.\n", "\n", "
\n", "\n", "However, also note that numerical approaches to solve completely general problems such as this often approximate it as a series of linear/matrix problems (i.e. solve a nonlinear problem by considering a series of linear problems - note this is also what Newton's method does), and that is why we focus on situations where our inversion problem can be written as a linear system. \n", "\n", "Also, the concepts we learn and issues we appreciate can be more easily understood in the linear/matrix system, but these concepts are still highly applicable to/representative of the nonlinear case. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Terminology\n", "\n", "#### Parameter estimation and inversion\n", "\n", "Note that in the above abstract discussion, if we consider $\\boldsymbol{y}$ to be a vector of observations (data), and $\\boldsymbol{x}$ to be a vector of model parameters that control a model ($F$), \n", "\n", "then the problem: given $\\boldsymbol{y}$ find $\\boldsymbol{x}$ is often called a *parameter estimation problem*.\n", "\n", "
\n", "\n", "Note that the $\\boldsymbol{x}$ we want to invert for, or estimate, doesn't have to include ALL parameters that govern the model, i.e. all inputs to the model. \n", "\n", "In this case you may see notation something like\n", "\n", "$$F(\\boldsymbol{x}_1;\\boldsymbol{x}_2) = \\boldsymbol{y}$$\n", "\n", "The semi-colon helps differentiate different types of variables/parameters that are input into the model, e.g. those that we wish to invert for and those that will remain fixed during the inversion. \n", "\n", "For example, a known physical parameter may be fixed (gravity or fluid viscosity), whereas there are other parameters such as those describing the geology, the initial condition, etc, that either aren't known, or there is uncertainty over their precise values.\n", "\n", "
\n", "\n", "Sometimes when we have a relatively small number of \"parameters\" to invert for people use the term \"parameter estimation\", but when there are a large number, e.g. representing an entire unknown field, people tend to use the language \"inverse problem\", but it's really the same thing.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Model and parameters\n", "\n", "Note also that different communities use the term *model* in different contexts.\n", "\n", "Mathematicians/computational scientists would probably refer to $F$ (or $\\boldsymbol{y} = F(\\boldsymbol{x})$) as the \"model\", e.g. the underlying equations (or their discretisation in a computer code), and $\\boldsymbol{x}$ as \"parameters\".\n", "\n", "Others would refer to $F$ as the operator (and think of it as a given, and in a sense not a \"model\"; they may also refer to $F(\\boldsymbol{x})=\\boldsymbol{y}$ as the \"mathematical model\"), and would reserve the word \"model\" for the parameters $\\boldsymbol{x}$ being inverted for. \n", "\n", "
\n", "\n", "In seismic inversion problems for example this is because it represents a \"model\" for the real geological make-up of the subsurface. \n", "\n", "In fluid dynamics we would probably call this a \"parametrisation\" of some unknown, uncertain or unresolved physical process - again \"parameter\"!\n", "\n", "
\n", "\n", "In what follows I will generally try to stick to the language that $\\boldsymbol{y}=F(\\boldsymbol{x})$ is the model (e.g. a wave equation or the Navier-Stokes equations and specifically their implementation in code), with $ \\boldsymbol{y}$ the data (or variables that can be post-processed into a form that can be compared to observational data) and *some of* the $\\boldsymbol{x}$ the parameters. \n", "\n", "
\n", "\n", "
\n", "\n", "However, moving forward I will use common seismic inversion notation for the linear(ised) case we will mostly consider:\n", "\n", "$$G\\boldsymbol{m}=\\boldsymbol{d}$$\n", "\n", "where $G$ is a matrix representing the linear *model*, $\\boldsymbol{d}$ is *data* and $\\boldsymbol{m}$ are *parameters* we are inverting for." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear algebra review\n", "\n", "It turns out that to understand inversion problems linear algebra is key (even when we have nonlinear problems). So let's quickly review some notation and key results.\n", "\n", "## Notation & terminology\n", "\n", "### Vectors and matrices\n", "\n", "\n", "- We will use bold fonts to indicate vectors, e.g. $\\boldsymbol{x } $ :\n", "\n", "$$ \\boldsymbol{x} = (x_1,x_2,\\ldots, x_m)^T$$\n", "\n", "(the transpose is here since we generally assume a vector is a *column* vector (i.e. has dimension, or shape, $m\\times 1$), but it's clearly easier and more compact to display it as the transpose of a *row* vector in writing - note that when writing on paper it's common to indicate a bold symbol with an underline: $\\boldsymbol{x} \\equiv \\underline{x} $. Sometimes people will use $\\vec{x}$.)\n", "\n", "- and capital letters to indicate matrices, e.g. $A$:\n", "\n", "$$A = \n", "\\begin{pmatrix}\n", "a_{11} & a_{12} \\\\\n", "a_{21} & a_{22}\n", "\\end{pmatrix}$$\n", "\n", "\n", "- Subscript letters or numbers will be used to indicate components of arrays or matrices, e.g. $x_i$ is the $i$-th component of the vector $\\boldsymbol{x}$, and $a_{ij}$ is the entry from the $i$-th row and the $j$-th column of the matrix $A$, with of course us needing to be careful over our numbering which starts from zero in our Python implementations.\n", "\n", "\n", "- We will identify column vectors of the matrix $A$ as\n", "\n", "$$A = \\begin{pmatrix}\n", " & & & \\\\\n", " \\boldsymbol{a}_{\\,:1} & \\boldsymbol{a}_{\\,:2} & \\ldots & \\boldsymbol{a}_{\\,:n} \\\\\n", " & & & \n", "\\end{pmatrix}$$\n", "\n", "i.e. $\\boldsymbol{a}_{\\,:j}$ indicates the $j$-th column of matrix $A$.\n", "\n", "\n", "\n", "- Square or rounded (i.e. paretheses) brackets are both fine for vectors and matrices (not for intervals where they mean different things!). \n", "\n", "
\n", "\n", "\n", "See also
\n", "\n", "A homework exercise asks you to code up a version of matrix-vector multiplication based on this mathematical interpretation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll come back to this interpretation in later lectures, but for now consider briefly a $2\\times 2$ case (and so $A$ maps vectors or points in $\\mathbb{R}^2$ to $\\mathbb{R}^2$).\n", "\n", "What properties of the columns of $A$ do we require so that we know in advance that for *any* $\\boldsymbol{b}$ we choose or are given we can always find a corresponding solution $\\boldsymbol{x}$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Some inversion (parameter estimation) examples\n", "\n", "## A linear model (polynomial interpolation)\n", "\n", "
\n", "\n", "We have seen a useful example in previous modules.\n", "\n", "
\n", "\n", "Suppose that our \"model\" is given by the polynomial (of degree $N$) function\n", "\n", "$$f(x; \\boldsymbol{a} ) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \\ldots a_N x^N$$\n", "\n", "
\n", "\n", "For a given, fixed vector $\\boldsymbol{a}$, $f$ can be thought of as taking in a scalar $x$ and returning a scalar $y$.\n", "\n", "
\n", "\n", "The vector of $N+1$ parameters $\\boldsymbol{a}$ fully describes our model, and in this case our model is linear in these parameters (and hence the process of finding the parameters is sometimes termed [linear regression](https://en.wikipedia.org/wiki/Linear_regression)). \n", "\n", "We will see that we end up with a linear problem to solve, even though our model is nonlinear in the parameter(s) $x$ we do not invert for. \n", "\n", "\n", "
\n", "\n", "You can think of $x$ as locations we evaluate our model at, while $\\boldsymbol{a}$ are unknown input parameters to our model that describe in some sense the governing \"physics\".\n", "\n", "\n", "For the inversion problem we are assuming we know the $x$'s (e.g. maybe these are the locations our data is collected at) and the $y$'s (the observations of the output of our problem/model at these locations), and we want to find $\\boldsymbol{a}$ (i.e. the parameters that govern our model).\n", "\n", "\n", "We described this problem as polynomial interpolation in *Numerical Methods*, and as long as we are given $N+1$ data points $(x_i, y_i)$ (with distinct $x_i$'s) we can solve it for $\\boldsymbol{a}$ and thus fix/fully define our \"model\"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problems leading to a square linear/matrix system\n", "\n", "First of all let's consider the case where we match the number of model parameters to the number of data points we have.\n", "\n", "
\n", "\n", "### Fitting a single data point\n", "\n", "The polynomial that fits the single data point $(x_0,y_0)$ is clearly simply the *constant* function given by\n", "\n", "$$ y = f(x; \\boldsymbol{a} ) \\equiv a_0 $$\n", "\n", "where through subsitution of the data we have\n", "\n", "$$ a_0 = f(x_0; \\boldsymbol{a} ) = y_0 $$\n", "\n", "So this case is completely trivial.\n", "\n", "
\n", "\n", "### Fitting two data points\n", "\n", "The polynomial that fits the two data points $\\{(x_0,y_0),(x_1,y_1)\\}$ is clearly the *linear* function given by\n", "\n", "$$ y = f(x; \\boldsymbol{a} ) \\equiv a_0 + a_1\\,x$$\n", "\n", "where through substitution of our data into the functional form we arrive at two simultaneous equations for two unknown parameters (or a $2\\times 2$ matrix system)\n", "\n", "\\begin{align*}\n", "(1) & \\;\\;\\;\\; y_0 = a_0 + a_1\\,x_0, \\\\[5pt]\n", "(2) & \\;\\;\\;\\; y_1 = a_0 + a_1\\,x_1. \n", "\\end{align*}\n", "\n", "
\n", "\n", "[A homework exercise asks you to solve this problem via substitution and show that the solution is\n", "\n", "$$ a_0 = y_0 - \\frac{y_1-y_0}{x_1-x_0}x_0, \\;\\;\\;\\;\\;\\;\\;\\; a_1 = \\frac{y_1-y_0}{x_1-x_0}. $$\n", "\n", "]\n", "\n", "\n", "\n", "However, the substitution approach isn't easy to do with large amounts of data/large dimensional $\\boldsymbol{a}$.\n", "\n", "\n", "Instead we can note that inverting for our $\\boldsymbol{a}$ parameters in this case is equivalent to forming and solving the linear system\n", "\n", "$$\n", "\\begin{pmatrix}\n", "1 & x_0 \\\\\n", "1 & x_1 \n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "a_0\\\\\n", "a_1\n", "\\end{pmatrix}\n", "=\n", "\\begin{pmatrix}\n", "y_0\\\\\n", "y_1\n", "\\end{pmatrix}.\n", "$$\n", "\n", "
\n", "\n", "Let's implement and plot some examples to check that this is correct visually in these cases on one and two parameters." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAE5CAYAAACjyPfFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABG80lEQVR4nO3dd5gUVdbH8e8BJAoKklZQUAmKu4Oyiu4qMphWXF0JBmAMmDBhWMMaMKCAuAZWXSOioq+oiAIqKophEMyoiAsIggQJiohKjnPfP27N2NP0zPQM010dfp/nmQe6uqr63Aqnb9+6t8qcc4iIiIiISGJUCTsAEREREZFMpgq3iIiIiEgCqcItIiIiIpJAqnCLiIiIiCSQKtwiIiIiIgmkCreIiIiISAKlVYXbzPqamYv4W2NmX5lZfzOrVs515ZtZfoJCTQozaxlsh74VWDbty58IZtbNzK5M8mdWeD+m0mcEn7PQzEZWYLmkb/fyMrPcYBvmVmDZlC9fqlG+L075vvIl4rwM8sRAM0vJ+lXEedWynMvtGpSrQ4JCqxQVPdaTUb6UPCDicArwF6An8CnwX+DmUCMKx3L8dngt7EAySDcg2RUj7cdwtnt5fYHfT19UYNlupH75UpXyvac8Ufm6UfnnZS5wC+lbvyrJrvhypXSFG7g4+CuvXUlw+crVSpBCpjvn5gX/f8vMWgFXkGVJ2Dm3Cfg47DiSxcxqBGXOKNm2H9OVc2412k9hUL4n+/JEpuZ7SSzn3KywYyhJpvwC+wyoa2aNAczsODP7yMw2mNlvZjbezNqWtLCZNTWzzWZ2eYz3BprZejOrH7zON7OpZna0mX0RvPc/M+sWY9ky44hY33FmNj2Y90szO8TMqpnZ7Wa23MxWmdlIM6sTsex2lxjN7GAze9HMlgTrmhOso1YFtitmdmtQzt/MbKWZvWtmh0bNU3ipvWcQ4y9mttrMRpnZblHzOjMbYmYDImJ838wOKGG7nBhsj00Ev1rNrKOZvW1ma81snZm9Y2YdI5YdHOzPgyOm1Qm2xUclXY423xXiLKBZxGXshRHvtzWzcWb2axD3x2Z2XBzbsKmZPWVmy8xsU7A/J0Qcr7H248hg+xxoZlOC4+xbM7swxvqPDrbRRjObZ2bnBcsvjJ43xrKdg+23JtiWb5rZH8taLlj2cvNdSDaa2TQz6xRjnkZm9qiZzQ3K8L2ZPWtmzSLLSgnb3cxqmtl/zJ9ja83sBzN71cz2jSO+8hyX9czsgYh9NMfM/mlmFmN9uRHTyswHZR1XUm7K979PU75PoXxvZgPxraQAWwrXG7z3PzMbETHvLma2zcyWRK3jAzN7IeJ1mbmplHj2NrPXguP2JzO7D6gRY75ewb7+KdjOX5rZWRHvtwQWBC8fi9hefYP3jzWz14Njt/AcucrMqsYRY+G+PylYbpOZfWNmp8aYN95zLD/ideHx+o9gO64MyvmMme0aZ/n+ZmYfBp+5NtgH5f/B75xLmz+gL+CAVlHTxwBbgdrAccA2YBLwD6APMA/4CWgWsUw+kB/x+gVgVtR6qwLfAyOjllsOzARODz5vUvD5rSLmK08cPwBfA72AE4BZwI/AY8CTwN/wLTpbgDsjlm0ZbI++EdN6AjcG6+mMT1o/AM9Hla1Y+UvZ5iOAM4AuwTqfBzYDORHz5AZxfB/EexxwKbAGeC9qfYXzfYC/nHcaMAf4GWgQFd8K/ElwTvAZOcHfBuBz4OSgvJ8F09oHy1YDPgS+BXYOpo0EfgP2KqWs++Av164ADg3+Dgze2z3Yd98F+/1EYGKwj7uWsQ0nAXOBPOAI/CXyR4CWpezHkcBqYDZwAXAM8GwwX5eI+doBm4ApwfY8FX8sLQYWlnGs/B1/3L4MnBT8fQj8AuxRRpnODdZXuL/7A0uCbRx5vrQF7gv20xH4Y/wzYCFQM47tvgv+GOyFP567B9vzV6BpGTHmEsdxiW94mAKsA64Cjg1idsDtMdaXW558UFr59Kd8j/J9YXxpn++B5sE2dMBhhesN3rsf+C5i3pOA9cG8bYJpdYLtfVHwOq7cVEIs1YH5wDLgbHy+fyXYH47g+yeY94bg2DkWOBq4DX/8XRi8XwOfex1we8T2ahS8f2EQX9fg2LkmOB7uiOOYy8cfs4si4pwAFFD8u66i53puEPcCfHe0Y/HH6wbgqbLKB+yN/44dFcRwJP47+d/lzmlhJdOK/PF7Am6LP8nqBwXfBowP5pmGP/GqRSy3V3DwDItjp3SKmPaPYNqhUcttAVpHTGscxHBDxLTyxLEF2DvG574dVf6xwIKI1y2JSsBR81uwnU4PDt7dSip/nNu/arC+OcB9MbbdxKj584LpR0VMc8BKoE5UObYAg6LiKwAOiFrni/jK1q4R0+oBq4CxUev8FXgK6B18bp84yjgSWBJj+t1s/yVbNdgWX5SxzrXAZaW8v91+DOKIrlzXCLbd8Ihpz+ITTu2IaX8ANlJ2hXse8E5ULPWCz7i3lHir4JN29P4+LfiMkWUcQ3sE83Uva7uXsHxtfDL/ZxnzxnVc4isW251H+C/OTUDDqPXlxjh/y8oHcZVPf8W2f1+U7xdEvG4Z6ziNeF/5PjXy/cDg86tFTS+s0LUIXt+LrwB/C1wQTDsumGff4HVcuamEOM5n++O5Cv7HoyOiwh21XJVgvz8GfBXj+DuvjPIXHocD8I03VcqYPz9GnFWBb4ApFTzHYp3rT0V97gP470krrXz4H3oOqFee8yfWX7p2KfkGv5FXAQ/hf3mcY/7yWwdgtHNua+HMzrkF+F/YnUtaoXMuH9/ScEHE5AuAGc656H5z3zrnvo1YdgX+V/Ke4C9nlTOOuc6576LKB/BmjHI3L+1SUnD56d9mNh9/Qm4B/g9/ErQuablS1ne0mb1nZj/jE9AWoA3+SzDaC1Gvx+CT6F+ipr/unFtX+MI5txDfNzF6voXOuelR044AJjjnfo1YfjU+cXWOmLYQ/6v7THwrzNPOuWdLKmccjgA+dr/3JcU5tw14DjjAzOqVsuxnwDXmu2D8KZ5LgYH1zrn3Ij5vEz7h7Bkxz6H47bk+Yr7l+BafEplZa3wLzyjzl7KrBZde1wMfBeUtSfPgL3p/v4Q/RqI/6yLzd5dYG7y/OHirxMv+UcufamafmNmvwfLrgJ3jXT5GnNHH5RHB6+ei5nsG30oUfVxGKzUfyA5Tvi+B8n2xdaZKvi/JZPz2OTJ4fSTwbvAXOW25c67wmNiR3PQX4PvI49k5V8D2+w0za21mz5nZUvw+3wKcR/w5+g/muw4uwrfQbwEG4wciNo5jFdFxbsMfTx3NrMqOnOsRogcbf41vxGpSxnLT8eV53sxOtqArW0Wka4W7O3AwsC/+l/OZzrlV+BYQw18CjPYD0KCM9T4MnGxmu5lZC/yvzUdizLcqxrRNQM3g/+WN45eo15tLmV4N/+uvJE/iE8/9+G4IBwOXBO/VLGmhWMzfHud1fAvtufjK3cHAVyWs68fIF865zUEZmpU2X8S06Plibb8GJUz/Ab/dI72Gv3RZA/hPjGXKo7TPtRifHek0/BfEv4AZwFIzu9nKvm1U9P6H4scZ+NbsFTHmi7WNIxUmjcf5PcEW/p0A7FbCcoWfud1nBInw58hpZnYpvpL0NtAD6Ig/jiCO49HMTgRG47vW9AEOwR+DP8WzfAlxRh+XDYBVbvsBWj9EvF+asvKB7Bjl+5Ip3/8uVfJ9TMEx+xXQxcwaAn8E3gv+coPZugSvI+OoaG76AyVv+yJmtjO+m0Z74DqgE36/P0GM/t7Rgu+xV/DfG4PxPxoOBoYEs8RzHJYUZ3V8t44dPddh+/O4cJuWGl/wo+tv+Pry/wE/BA1A8VTyi0nXu5T8L/KXZ4Rf8E3/TWO815SoykAMTwND8Zcy6+P7+IyqQHw7GkeFmFlNfL+wgc65+yKm/6mCq+yJb+Xo4ZzbErG++vjLd9GK/VI0s+r47bi0tPkipkXP52LMt4qSt2v0CfUg/stqPjDczA6LLEc5lfa5LsZnFwlaxC4BLgkGeJwF3IqvND5cwXgKLSd2C0JZv9oLj8Hr8ZXhaJtjTIv8zO0+I2ghj66o98J3W7kqYr69yogtevl5zrm+EcvvRHwJtlBZx+UqoIGZVQ8qDYUK93dCzleJm/J9DMr3qZnvy/AevgGmC/64mEGQw83sMOBA4NGoOCqam5YD+8eYHr0//gK0wHevmlo40eK/1/0+wEHAGc65ZyKWPzHO5WPFVDhtM/57shYhnGOFgivN75lZDXzf/NuA18yspXNuZbzrSdcW7piCy1afA6dEjo4NWi/+ir+kU9ryq/EJ9wL8wI1ng2lJjWMH1MAnnOgk07eC66uN76tYlAjN7EhKvlQePar4FPwx9lHU9OMtavQ9vjUler5YJgN/N7O6EcvXxQ9qmRwxrQ9+8E8/fII7AH+SlGUT/uSO9bmHWsTDAoJ9exrwpXNuTRzrxjk3xzl3A/5LOq67gZThY/z2rB0R1x/wSaE0c/ADF/d3zk2L8TejlGWX4PtwR+/vnmz/I7422x+PZ8dYZ0nbvTbbd1M5g9Jb/aKVdVxODl6fEjVfHj7hR3cxqIiSyicVpHyvfB8xLZXyfWHLaaz1vodv2b8A38/YBQ0yM/GNMFXxXUwi46hobvoI2MMi7jITtEZH77fC747oH1knxVmuWMvvFMQYr+g4q+LL/KlzriBJ51hp+w3w3Tqdc+8Cd+IHuJan8ShtW7hLcxP+0tIEM3sI39fzVvyI5XviWP4hfu/XF+vyYrLiKDfn3G9m9jFwlZktxw9WOYftL93FayJ+tPxIM3sS35fvJrZvmSi0fzDf88G8Q4DJzrl3oubbgL+f7l34L41b8XfkiOcy4CD8pat3zOzf+C+Ha/En/W1Q1IL6MPC4c25MMG0AcIeZvRXZLzqGWfgWhYvwgzQ2Oue+DmLrC0wys1uCeC8Oyvn3klZmZrvgW5BH8Xtf1JPwLUFvxVHesgzGD+p408zuxm/Pm/CX4wpKWsg558zsEuDloGXqBfzx0gSfwBY754aVsGyBmd0KjIjY363wreXRFZaJwLVmdgP+oSVHBvFGK2m7TwS6mdl/8CPX/wxcRuwWt5KUdVy+AUwFHjGzRvgvv+PxfRiHlqcFoxQllU92jPK98n3K5PuIdYLfL28A25xz04Jp7+N/1BzF711/wFfE++PzbmT//h3JTU/hu4iMDfLvCnz3o+j+5x8G5XswKGsd/J1vVuLvElXoR3xLci8zm4EfS7MA391vETDEzLbhv+P+WUpcsfwIjA4+/yfgIvy2vihinkSfYyWV7xR8X/rX8Q1NDfHfdcuA/5XrE9wOjrpM5h8l3CYqxnzH4X/dbcDvjJeBtlHz5FPCqG18699nJbyXD0yNMX0hUXdnKEccU6OmtST2aNmBRIx+JvadJ1riT9I1+BPsAXyCcGx/d4WY5Y/6zEuDg24DfvDf0dHL8vso4B74Ud+/Bp//LFGjqIP5huBvQ7QEP0p4CtuPTo+5nYP3DsFXYtfiT4p3gI7Be9WCbT6H4iPjDV/BXULE6P0Y666DH6BSeJl4YcR7bYHxwb7ciG9dOK6M7VcDf4lwZhDv6mA79omYJ9Z+HEns0fPb7Td8383p+F/o3+ErEOPwLTElfkYw/S/4iuwvQZkW4r9A/xLHsXE5PtFuxH9ZHU7UeYBvLXgYn0TXBJ+1VxDLwLK2O751ZzA+ua3Ht2QcGP05JcSXS/zHZT38ubIc33I0F/+lYTHWF30elZkPSjuu9Ffi/uuL8r3yfRrl+2C5qvjuLSvwjR4u6v1PiLgTSTCt8A4mI2Osr8zcVEose+MriuvxOfg+/PeDo/htAY8Evgz2+3x8o8bAGLF3w/+g2BJ5LOKvKEwNPmcJ/sfQedGfU9o5hr9bz//w32NzgNMq41zn9+P16BLyS8vSyof/jnwZX9neFOyHMdGfG89f4e1QJGBmbfAtkec75x4PO550YGa5+F/oxzjnYvUHjpzXAUOcczcmIbSsFAyCmQe85pw7N+x4wlKe41Kyk/J9+SnfS2Uy/5Caas65w8OOJdEysUtJhZhZc/xl8Vvxv2B25JZCIkljZv/FXxZchn9gw+X4Liv3lbacSLZSvheRZEv6oEkze8LMVphZzL4v5t1v/hHVM8zfqigZzsMPVmiCv9y/IUmfK7KjagL/xl9CHY6/7Hq0K33go0ilSuHcHovyvYgkVdK7lJjZEfi+WE8757a7S4OZHY/vR3Y8vu/Wfc65Q5IapIiIlItyu4hIyZLewu2ce5/S72F5Ej5hO+efPLRrcJszERFJUcrtIiIlS8U+3M3wo0ELLQmmbfeEITPrh7/vJjVr1vzznntmx5OUCwoKqFIlo26hHlO2lBNU1kw0d+7clc65RmHHkULiyu3K65kvW8qaLeWEssu6dm01li2rRb16W2jadGMSI6t8Fc3tqVjhthjTYvZ7cc4Nx/dZpW3btm7OnDmJjCtl5Ofnk5ubG3YYCZct5QSVNROZ2aKwY0gxceV25fXMly1lzZZyQullnT0bOnaEgw+G99+HmvE87D2FVTS3p+JPryXAHhGvm+PvviAiIulLuV0ky/z2G3TrBrVrw0svpX9le0ekYoX7FeDMYET7ocBvzrntupOIiEhaUW4XySIFBXDGGfDddzBmDOyxR9nLZLKkdykxs+fwT/5paGZLgFuAnQCcc4/gn4p0PP7BHeuBs5Mdo4iIlI9yu4hEGjwYXn0V7r8fjjgi7GjCl/QKt3OudxnvO+CSJIUjIiKVQLldRApNmAC33AJnngn9+4cdTWpIxS4lIiIiIpKG5s6FvDzo0AEeeQQs1nDpLJSKdykRyRqrV69mxYoV7LLLLsyePTvscJIiU8pap04dmjdvnjW3/UpFhefPli1bwg5lh2XKeRGPevXqsXr1aurVqxd2KFLJVq/2gySrV4exY6FWrbAjSh2qcIuEZPXq1fz44480a9aMrVu3Zs2Xz5o1a6hbt27YYeyQgoICli5dysqVK2ncuHHY4WSlyPOnVq1aWJo3o2XCeREP5xw///wzP/74I0DW5L1sUFAAffv6Fu5Jk6BFi7AjSi1qmhEJyYoVK2jWrBm1a9dO+8pCtqlSpQpNmjTht99+CzuUrKXzJz2ZGTVq1KBZs2asWLEi7HCkEg0dCuPGwV13QZcuYUeTelThFgnJli1bqKXrbWlrp512YuvWrWGHkbV0/qS3WrVqZURXIPE++aQBN90EvXvDFVeEHU1qUoVbJERqmUtf2nfh0z5IX9p3mWPePBg8eD9ycmDECA2SLIkq3CIiIiJSbmvXQvfuUKWK705Su3bYEaUuVbhFREREpFycg3POgVmz4KabZrHXXmFHlNpU4RaRtPXiiy/q0rRIBen8kR1x113+ke1Dh8JBB/0SdjgpTxVuEckqLVu25O677w47DJG0pPNHAN56C66/Hk49Fa65Juxo0oMq3CJpbtSoUbRs2ZIqVarQsmVLRo0aldTP37x5c1I/T6Qy6fwRKZ8FC/zdSNq1g8cf1yDJeKnCLZLGRo0aRb9+/Vi0aBHOORYtWkS/fv0SWmnIzc3loosu4uqrr6ZRo0YcdthhAAwbNoycnBzq1KlDs2bNOO+88/j111+LlmvatCmjR48uen3YYYdRt27dolvrffvtt5gZS5cuLfGzn376aVq0aEHt2rU54YQTih6eUWj+/PmcdNJJNG3alDp16tChQwcmTJhQLPZFixZxzTXXYGZFl9N//vlnevfuTfPmzalVqxb7778/Tz755A5vK0ltqXT+PPDAAzp/JOWtX+8HSRYU+EGSO+8cdkTpQxVukTQ2YMAA1q9fX2za+vXrGTBgQEI/95lnnsE5x5QpU3j66acB/zCYe++9l5kzZ/Lss8/y6aefcumllxYt07lzZ957772iGKdNm0aNGjWYNm0aAPn5+bRq1YpmzZrF/MxPPvmEvn370q9fP6ZPn86JJ57IzTffXGyetWvX0rVrVyZNmsRXX31Fz5496dGjB9988w0AY8eOpXnz5tx8880sX76c5cuXA7Bx48aiysXMmTO5/PLLueCCC3jnnXcqd8NJSkml88fMdP5ISnMOzj8fZsyAZ5+FVq3CjijNOOcy4q9NmzYuW7z33nthh5AUmV7OWbNmFf1/9erVFVqHmTlguz8zq6wwt9O5c2f3pz/9qcz53njjDVe9enW3bds255xzDz30kGvTpo1bvXq1e+utt9y+++7rzjzzTHf77bc755zr06ePO++880pcX+/evd3RRx9dbNq5557rfBor2SGHHOIGDRpU9LpFixburrvuKjP+0047zZ177rmlzhO5D6MB01wK5MZ0/istr5e27eOVSudPdA4o6fxxzqX9+VNY1srYh6ks077Dhg1zDpwbMmT79zKtrKWpaG5XC7dIGttzzz3LNb2y/PnPf95u2rvvvssxxxxD8+bNqVu3Lj169GDz5s388MMPgL8cPXfuXJYvX05+fj5dunQhNzeX/Px8ACZPnkxubm6Jnzl79mz+8pe/FJsW/XrdunX861//ol27dtSvX5+dd96ZadOmsXjx4lLLs23bNoYMGUJOTg677bYbO++8M2PHji1zOUlvqXT+TJ48Oa7zZ9myZTp/JOnee88PjuzRww+WlPJThVskjQ0ZMoTaUU8aqF27NkOGDEno59apU6fY60WLFvH3v/+d/fbbjzFjxvD555/zxBNPAL8PCttvv/1o0qQJU6ZMKaowdOnShQ8++IBZs2axdOnSUisMvmGhdFdffTVjxoxh0KBBTJ48menTp9OxY8cyB6bdfffd3HPPPVxzzTW88847TJ8+nW7dumlAW4ZLpfPnlFNOiev8yc/P1/kjSbV4sb8bSZs2MHKkBklWVLWwAxCRisvLywN8X9TFixez5557MmTIkKLpyTJt2jQ2b97Mf/7zH6pWrQpQbLBVoc6dO/Pmm28ybdo0OnfuTOPGjWnYsCF33nlnqf1PAdq1a8fHH39cbFr066lTp3LmmWfSs2dPwPctnT9/Pm3atCmap3r16mzbtm275U488UTOOOMMwFdO5s6dy6677hr/RpC0k47nz2uvvabzR5Jmwwbfqr15M4wfD3Xrhh1R+lILt0iay8vLY+HChRQUFLBw4cKkVxYAWrduTUFBAffeey8LFizgueee4957791uvtzcXMaOHUvr1q1p3Lgx4CsRzzzzTKmtcwCXXXYZb7/9NkOHDuXbb7/lscceY9y4ccXmadOmDePGjeOLL77g66+/5vTTT2fjxo3F5mnZsiVTpkxh6dKlrFy5smi5d955h6lTp/LNN9/Qv39/FixYUPENImkj3c6f0aNH6/yRpHAOLroIPv8cRo3yLdxScapwi8gOy8nJ4b777mPYsGG0a9eOESNGxHw4RpcuXdi2bVuxykGsabEceuihPP744zz88MPk5OQwduxYBg4cWGyeYcOG0bhxYzp16kTXrl059NBD6dSpU7F5brvtNr7//nv22WcfGjVqBMCNN95Ix44d6dq1K0cccQR16tQJpeIl2SknJ4d///vfOn8kpTz4IDz1FNx6K5xwQtjRpD+Lp19XOmjbtq2bM2dO2GEkRX5+fpnJNRNkejlnz57NfvvtB8CaNWuomyXX6jKprJH7MJqZfe6cOyjJIWWU0vJ6ads+HWXSeVGWwrJm2j6Mls7fYe+/D0cdBccf7++3XaWM5tl0Lmt5VTS3q4VbRERERABYsgROOQX23huefrrsyrbER4MmRURERISNG6FnT/9Eyfx82GWXsCPKHKpwi4iIiGQ55+CSS+DTT2HsWMjg3j6h0IUCERERkSz36KPwxBMwYAB07x52NJlHFW4RERGRLPbBB3DZZX6Q5K23hh1NZlKFW0RERCRLLVsGJ58MLVr4+20Hz16SSqY+3CIiIiJZaPNmX9leswYmTQI9HDRxVOEWERERyUKXXw4ffQRjxsAf/xh2NJlNXUpEREREssyIEfDII3Dddb6VWxJLFW4REZEU0bJly5iPdRepTJ984m8BeOyxMHhw2NFkB1W4RaRCfvzxRy6//HL22WcfatSoQbNmzejatSuvv/56Uj4/Pz8fM2PlypVptW6Rvn37csIJJ8R877PPPuPiiy9OckSSTX74wT/cplkzeO45DZJMFvXhFpFyW7hwIYcddhh169Zl6NChtG/fnoKCAt555x0uvPBCFi9eHHaIImmpUaNGYYcAwObNm6levXrYYUgl27zZP7Z91Srfd7tBg7Ajyh5q4RaRcrv44otxzjFt2jROPfVU2rZty3777Uf//v356quvAFi8eDHdu3enbt261K1blx49erBkyZKidQwcOJA//vGPPP/88+yzzz7UrVuXbt26FWtV/vrrrznqqKOoV68edevWpX379rz33nssXLiQLl26AL6CYmb07dsXgIkTJ9KpUyfq169PgwYN+Nvf/sbs2bOL1rlw4ULMjJdeeoljjjmG2rVr065dOyZNmlT0fknrFkm06C4lZsbw4cM55ZRTqFOnDnvvvTfPPPNMsWWWLl1Kr169qF+/PvXr1+fvf/873377bdH78+fP56STTqJp06bUqVOHDh06MGHChO0+d+DAgZxzzjnsuuuu5OXlJbagEoorr4SpU+Hxx6F9+7CjyS5q4RZJIVdcAdOnJ/czDzgA7r03/vlXrVrFxIkTGTx4MDvvvPN279evXx/nHN26daNmzZq8++67mBn9+/enW7duvPvuu0XzLly4kNGjRzNu3DjWrVtHr169GDBgAI8++igAffr0oX379nz66adUq1aNr7/+mpo1a7LHHnvw0ksv0bNnT2bOnEmDBg2oVasWAOvWreOKK64gJyeHDRs2MHjwYE488URmzZpVrMVuwIAB3HXXXTz00EMMHjyYXr16sWjRolLXLaktHc6firjtttu44447GDp0KI8//jjnnHMOnTp1okWLFqxfv54uXbrw17/+lcmTJ1O9enXuvvtujj76aGbPnk3t2rVZu3YtXbt2ZfDgwdSqVYvRo0fTo0cPZsyYwb777lv0OcOGDePGG29k2rRpOOcSWyhJupEj4cEH4aqroHfvsKPJPqpwi0i5zJs3D+cc++23X4nzvP3223z11VfMnz+fli1bAvDss8/SqlUr8vPz+cc//gHA1q1bGTlyJLvssgsA/fr148knnyxaz6JFi7j66quLKgWtWrUqeq9BcC20cePGNGzYsGh6z549i8Xy5JNPUq9ePT799FMOP/zwoun//Oc/OfHEEwG4/fbbefrpp5k+fTqHH354iesWCcMZZ5zB6aefDsCgQYO47777mDJlCi1atOD555/HOceTTz6JmQHw6KOP0rhxYyZMmMCpp55K+/btaR/RnDlgwABeffVVXnzxRW688cai6Z07d+Zf//pXcgsnSTFtGlx4IRx5JNxxR9jRZCdVuEVSSKJbyipDPC1fs2fPZvfddy+qbAPsvffe7L777nzzzTdFFe4WLVoUVbYBdt99d1asWFH0+sorr+S8887jqaee4qijjqJnz57FWuRimT9/PjfddBOffPIJP/30EwUFBRQUFGzXrzwnJ6fY5wLFPlvSTzqcPxUReaxWq1aNRo0aFR2rn3/+OQsWLKBu3brFllm/fj3z588H/FWfW2+9lQkTJrB8+XK2bNnCxo0bi60X4KCDDkpwSSQMK1ZAjx7QtCmMHg3VVPMLhfpwi0i5tG7dGjMr1i86mnOuqLUtWuT0nXbaabv3CgoKil4PHDiQWbNm0a1bNz788ENycnJ44oknSo3vxBNP5KeffuLRRx/lk08+4csvv6RatWps3ry52HyRn10YU+Rni6SK0s6TgoICDjjgAKZPn17sb+7cuVxwwQUAXH311YwZM4ZBgwYxefJkpk+fTseOHbc7J+rUqZOcAknSbN0Kp50GP/0EY8eCLtiFRxVuESmXwoGIDzzwAGvXrt3u/V9//ZV27dqxdOlSFi5cWDT9u+++Y9myZWW2UEdr3bo1l112Ga+99hrnnnsuI0aMACjqj71t27aieX/++Wdmz57NDTfcwNFHH81+++3HmjVr2Lp1a7k+M9a6RVJRhw4dmDdvHg0bNqRVq1bF/gq7Rk2dOpUzzzyTnj17kpOTQ/PmzYtavyWzXXMN5OfD8OHQoUPY0WQ3VbhFpNweeughnHMcdNBBjBkzhjlz5vDNN9/w8MMPk5OTw9FHH0379u3Jy8vj888/Z9q0aeTl5dGhQwc6d+4c12ds2LCBSy65hPz8fBYuXMgnn3zC1KlTadeuHeC7o5gZr732Gj/99BNr166lfv36NGzYkMcee4x58+YxefJkLrzwQqqV8xpqrHWLVKbVq1cXa5GeMWNGsR+o8crLy6NJkyacdNJJTJ48mQULFvD+++9z1VVXFd2ppE2bNowbN44vvviCr7/+mtNPP52NGzdWcokk1TzzjO9mdfnlcMYZYUcjqnCLSLnttddefPHFFxxzzDFce+215OTkcOSRR/LKK6/w6KOPYmaMHz+eRo0akZubS5cuXWjatCnjx48vsatJtKpVq/LLL79w1lln0bZtW7p3785f/vIXhg0bBkCzZs249dZbGTBgAE2aNKF///5UqVKF0aNHM2PGDP74xz9yySWXMGjQIGrUqFGu8sVat0hlmjJlCgceeGDR3+GHH87VV19d7vXUrl2b999/n7333ptTTjmFfffdl7POOotffvmF+vXrA/7uI40bN6ZTp0507dqVQw89lE6dOlV2kSSFfPklnH8+dO4Md90VdjQCYGHc+sfMjgPuA6oCI5xzd0S9vwvwDLAnfmDn3c65J7dbUYS2bdu6OXPmJCji1JKfn09ubm7YYSRcppdz9uzZRXf6WLNmzXaDnjJVJpU1ch9GM7PPnXNZMwot2Xm9tG2fjjLpvChLYVkzbR9GS/Z32KhRoxgwYACLFq2jatUv2WWX+syeXYfGjRP/2Zn+fR2pork96S3cZlYVeBDoCrQDeptZu6jZLgFmOefaA7nAPWamR16JiKQg5XWRcI0aNYp+/fqxaNES4Dm2bWvIunXHMWnSqLBDk0AYXUo6AvOcc9855zYDzwMnRc3jgLrmrz3vDKwCyjfqSUREkkV5XSREAwYMYP369cBQ4GjgIjZtmsqAAQNCjkwKhXE3xmbA9xGvlwCHRM3zAPAKsAyoC5zmnNvufl1m1g/oB/4RzPn5+YmIN+WsXbs2K8qa6eXcZZddWLNmDeDvhlH4/0yXSWXduHFjRh+j5ZD0vB55/mSCTDovylJY1kw/f5L5HeafM3AqcA3+YtPIounJiCHTv68rQxgV7lgjpqI7kv8NmA4cCewDTDKzKc651cUWcm44MBx8X79s6T+ULX2lMr2cs2fPLuqzmY39NzNBzZo1OfDAA8MOIxUkPa9Hnj+ZIJPOi7IUljXTz59kfoc1bXosy5c/AUwF/lk0fc8990xKDJn+fV0ZwuhSsgTYI+J1c3yLR6SzgbHOmwcsAMp3816RNBDGoGWpHNp3xYSS17UP0pf2XeVZtQq2bRuD2W/AKcAWwN/BZsiQIaHGJr8Lo8L9GdDazPYKBsz0wl9mjLQYOArAzJoAbYHvkhqlSILttNNObNiwIewwpIK2bNlS7vt7Z7Ck53WdP+ltw4YN2z1BU8pv2zbo0wd++aUut9wygxYtamBmtGjRguHDh5OXlxd2iBJI+reFc26rmfUH3sTfPuoJ59xMM7sweP8RYBAw0sy+xl+qvNY5tzLZsYokUuPGjVm6dCnNmjVTa0+aKSgo4Mcff2SXXXYJO5SUEEZejzx/atWqFff93SVczjk2bdrEL7/8QpMmTcIOJ+3ddBO8+aZ/kuT55x/HLbcsDDskKUEozTPOudeB16OmPRLx/2XAscmOSySZ6tWrB8CyZctYs2YNNWvWDDmi5Ni4cWNGlLVOnTo0bNgw7DBSRrLzeuT5s2XLlspabWgy5byIx4YNG2jVqlXRPpSKeeklGDrUP+Dm/PPDjkbKouuhIiGqV68e9erVIz8/P6MHD0XKprJKYhWeP5kgm86L/Pz8jNlvYZk5E846Cw49FP7737CjkXjo0e4iIiIiaeLXX6F7d6hb17dy16gRdkQSD7Vwi4iIiKSBggI44wxYsADy82H33cOOSOKlCreIiIhIGrj1VpgwAR58EA47LOxopDzUpUREREQkxb38Mtx2G5x9Nlx0UdjRSHmpwi0iIiKSwr75xnclOeggeOgh0F0w048q3CIiIiIpavVq6NYNataEsWP9v5J+1IdbREREJAUVFMCZZ8K8efD227DHHmFHJBWlCreIiIhIChoyxPfdvu8+yM0NOxrZEepSIiIiIpJiXnsNbrnF992+9NKwo5EdpQq3iIiISAr59lvIy4MDDoBHH9UgyUygCreIiIhIilizxg+SrFYNxo2DWrXCjkgqg/pwi4iIiKQA56BvX38bwEmToEWLsCOSyqIKt4iIiEgKuOMOf+u/e+6BI48MOxqpTOpSIiIiIhKyiRNhwADo1Qv++c+wo5HKpgq3iIiISIjmz4feveFPf4IRIzRIMhOpwi0iIiISknXroHt3X8keNw7q1Ak7IkkE9eEWERERCYFzcO65MHMmvPEG7L132BFJoqjCLSIiIhKCe+6B0aNh6FA49tiwo5FEUpcSERERkSR7+2249lo4+WT/r2Q2VbhFREREkmjBAjjtNNhvP3jySQ2SzAaqcIuIiIgkyfr10KMHFBTA+PGw885hRyTJoD7cIiIiIkngHPTrB199BRMmQKtWYUckyaIKt4iIiEgS3H8/jBoFgwfD8ceHHY0kk7qUiIiIiCRYfj5cdZW/5/b114cdjSSbKtwiIiIiCbR4MZx6KrRuDSNHQhXVvrKOdrmIiIhIgmzY4AdJbtrkB0nWqxd2RBIG9eEWERERSQDn4OKL4fPP4eWXoW3bsCOSsKiFW0RERCQBHn7YdyG55Rb4xz/CjkbCpAq3iIiISCWbMWMXLr8cTjgBbr457GgkbKpwi4iIiFSipUth4MD92XtveOYZDZIUVbhFREREKs2mTX6Q5MaNVRg3DnbZJeyIJBWowi0iIiJSCZyDSy6BTz+F6677hnbtwo5IUoUq3CIiIiKVYPhwePxxuOEGOOKIlWGHIylEFW4RERGRHfTRR3DppXDccXDbbWFHI6lGFW4RERGRHbB8OfTsCXvsAc8+C1Wrhh2RpBo9+EZERESkgjZvhpNPht9+gzffhPr1w45IUpEq3CIiIiIVdMUV8OGHMHo0/OlPYUcjqUpdSkREREQq4PHH/dMk//UvOPXUsKORVKYKt4iIiEg5ffopXHwxHHMM3H572NFIqgulwm1mx5nZHDObZ2bXlTBPrplNN7OZZjY52TGKiEj8lNclm/z4o3+4ze67w3PPaZCklC3pfbjNrCrwIHAMsAT4zMxecc7NiphnV+Ah4Djn3GIza5zsOEVEJD7K65JNtmyBU06BVat83+3ddgs7IkkHYbRwdwTmOee+c85tBp4HToqapw8w1jm3GMA5tyLJMYqISPyU1yVrXHUVTJkCI0bAAQeEHY2kizAq3M2A7yNeLwmmRWoD1DezfDP73MzOTFp0IiJSXsrrkhWefhr++1+48kro0yfsaCSdhHFbQIsxzUW9rgb8GTgKqAV8ZGYfO+fmFluRWT+gH0CjRo3Iz8+v/GhT0Nq1a7OirNlSTlBZJe0pr++gbDov0rWsc+bszKWXduDAA3/j+ONnkJ8ffYgXl67lrIhsKmtFhVHhXgLsEfG6ObAsxjwrnXPrgHVm9j7QHiiWmJ1zw4HhAG3btnW5ubmJijml5Ofnkw1lzZZygsoqaU95fQdl03mRjmX96Sc46yz4wx/gzTfr06hR5zKXScdyVlQ2lbWiwuhS8hnQ2sz2MrPqQC/glah5XgY6mVk1M6sNHALMTnKcIiISH+V1yVhbt8Jpp8GKFTBuHDRqFHZEko6S3sLtnNtqZv2BN4GqwBPOuZlmdmHw/iPOudlmNhGYARQAI5xz/0t2rCIiUjbldclk114L773n+2936BB2NJKuQnm0u3PudeD1qGmPRL2+C7grmXGJiEjFKK9LJnr2WRg2DC69FM44I+xoJJ3pSZMiIiIiUaZPh/POg06d4J57wo5G0p0q3CIiIiIRfv4ZuneHBg1gzBjYaaewI5J0F0qXEhEREZFUtHUr9O4Ny5b5B9w0aRJ2RJIJVOEWERERCQwYAJMmweOPQ8eOYUcjmUJdSkRERESAF16AO++Eiy6Cc84JOxrJJKpwi4iISNb7+ms4+2z461/h3nvDjkYyjSrcIiIiktVWrYJu3aBePXjxRahePeyIJNOoD7eIiIhkrW3bIC8Pvv8e8vP949tFKpsq3CIiIpK1br4ZJk6ERx7x3UlEEkFdSkRERCQrjR0Lt9/uH3DTr1/Y0UgmU4VbREREss6sWXDWWf7Wfw88AGZhRySZTBVuERERySq//eYHSdap41u5a9QIOyLJdOrDLSIiIlmjoABOPx0WLIB334VmzcKOSLKBKtwiIiKSNW67DSZM8N1IOnUKOxrJFupSIiIiIlnhlVfg1lt93+2LLw47GskmqnCLiIhIxvvmG9+V5KCD/C0ANUhSkqnMCreZ9UpGICJhGzVqFC1btqRKlSq0bNmSUaNGhR2SSEIor0u2Wb0auneHmjX9IMmaNcOOSLJNPC3cT5vZu2bWLuHRiIRk1KhR9OvXj0WLFuGcY9GiRfTr10+VbslUyuuSNQoKfBeSb7+FF16APfYIOyLJRvEMmvwz8BDwpZn9FxjonFub2LDK7/vva5ObG3YUyfHrrwew665hR5F4ySznxx/vw6ZNrxWbtn49nHtuTR57LPGfny37FLKrrCksLfK6SGW4/XYYPx7+8x+ypp4gqafMFm7n3NfOuU5AP+B0YI6Z9U54ZCJJtGnTxnJNF0lnyuuSLV5/3T+6PS8PLr887Ggkm8V9W0Dn3FNmNh64Hfg/M+sH9HfOzUxUcOWxxx7ryc8PO4rkyM+fTm4W/ExPZjlbtuzLokWLtpveokUL8vMXJvzzs2WfQvaUNR0GZKV6XhfZEfPmQZ8+0L49DB+eHuekZK5y3aXEOfebc+4S4GCgIf5y5D1mVjch0YkkyZAhQ6hdu3axabVr12bIkCEhRSSSHMrrkonWrvVPkqxaFcaNg6j0LpJ0cVW4zWwnM+toZpeZ2bPAS8D++BbyS4BvzOwfCYxTJKHy8vIYPnw4LVq0wMxo0aIFw4cPJy8vL+zQRBJCeV0ylXNw9tkwezaMHg0tW4YdkUgcXUrM7EPgQKA6UAB8BbwKTAU+ANYCtwAvmtllzrlHEheuSOLk5eWpgi1ZQXldMtmdd8KLL8Jdd8HRR4cdjYgXTx/utcBQfBL+2Dm3LsY8V5nZj8ANgBKziEhqU16XjPTWW3DDDXDaaXDVVWFHI/K7Mivczrlj41zX+8AdOxaOiIgkmvK6ZKLvvoNevWD//eHxxzVIUlJLZT7a/SvgpEpcn4iIhEt5XdLCunX+SZLgB0nWqRNuPCLR4r4tYFmccxvwfQBFRCQDKK9LOnAOzjsPvv4a3ngD9tkn7IhEtldpFW4RERGRZBs2DJ5/3j9R8m9/Czsakdgqs0uJiIiISNK8/Tb861/Qsydcd13Y0YiUTBVuERERSTuLFvlBkvvuC08+qUGSktpU4RYREZG0smGDHyS5dSuMHw919VxUSXHqwy0iIiJpwzno1w+mT4dXX4XWrcOOSKRsqnCLiIhI2vjvf+GZZ+DWW+Hvfw87GpH4qEuJiIiIpIXJk+HKK+Gkk+DGG8OORiR+qnCLiIhIyvv+ezjlFGjVCp5+GqqoBiNpRIeriIiIpLSNG6FHD//v+PFQr17YEYmUj/pwi4iISMpyDi66CKZN849t33ffsCMSKT+1cIuIiEjKevhhGDkSbroJunULOxqRilGFW0RERFLS1Klw+eVw/PEwcGDY0YhUXCgVbjM7zszmmNk8MyvxYaxmdrCZbTOzk5MZn4iIlI/yulS2pUvh5JOhZUsYNUqDJCW9Jf3wNbOqwINAV6Ad0NvM2pUw37+BN5MboYiIlIfyulS2TZugZ09Yu9YPktx117AjEtkxYfxe7AjMc85955zbDDwPnBRjvkuBl4AVyQxORETKTXldKtVll8Enn8BTT8H++4cdjciOC6PC3Qz4PuL1kmBaETNrBnQHHkliXCIiUjHK61Jphg/3f9df71u5RTJBGLcFtBjTXNTre4FrnXPbzGLNHqzIrB/QD6BRo0bk5+dXUoipbe3atVlR1mwpJ6iskvaU13dQNp0XpZV15sx6XHHFARx88K8cddQM0nmTaJ9KpDAq3EuAPSJeNweWRc1zEPB8kJQbAseb2Vbn3PjImZxzw4HhAG3btnW5ubkJCjm15Ofnkw1lzZZygsoqaU95fQdl03lRUlmXL4c+fWDPPWHixAY0aLD9POlE+1QihVHh/gxobWZ7AUuBXkCfyBmcc3sV/t/MRgITopOyiIikDOV12SGbN/vHtv/2G0ycCA0ahB2RSOVKeoXbObfVzPrjR6lXBZ5wzs00swuD99W/T0QkjSivy4668kr44AN4/nnIyQk7GpHKF8qj3Z1zrwOvR02LmZCdc32TEZOIiFSc8rpU1JNPwoMPwjXXwGmnhR2NSGLoNvIiIiISis8+g4sugqOPhttvDzsakcRRhVtERESSbsUK6NED/vAH35WkWijX3EWSQ4e3iIiIJNWWLX6Q5MqV8OGHsNtuYUckkliqcIuIiEhSXX01vP8+jBoFBx4YdjQiiacuJSIiIpI0b73VhPvvh3/+0993WyQbqMItIiIiSfH553DPPW3IzYU77ww7GpHkUYVbREREEu6nn/wgyV133cILL2iQpGQXVbhFREQkobZu9ffY/vFHuO22/9GoUdgRiSSXfl+KiIhIQl17Lbz3HowcCS1arA07HJGkUwu3iIiIJMxzz8GwYdC/P5x1VtjRiIRDFW4RERFJiK++gnPPhU6dfKVbJFupwi0iIiKVbtUq6N4dGjSAMWNgp53CjkgkPOrDLSIiIpVq2zbo3RuWLvUPuGnSJOyIRMKlCreIiIhUqhtvhLfegsceg0MOCTsakfCpS4mIiIhUmjFj4I474IIL4Lzzwo5GJDWowi0iIiKV4n//g7PPhr/+Fe6/P+xoRFKHKtwiIiKyw375xQ+SrFcPXnwRqlcPOyKR1KE+3CIiIrJDtm2DvDxYtAjy8+EPfwg7IpHUogq3iIiI7JCBA+GNN+Dhh313EhEpTl1KREREpMLGjYPBg/0Dbi64IOxoRFKTKtwiIiJSIbNnw5lnQseO8MADYBZ2RCKpSRVuERERKbfffoNu3aB2bXjpJahZM+yIRFKX+nCLiIhIuRQUwBlnwHffwTvvQPPmYUckktpU4RYREZFyGTQIXn0V/vtfOOKIsKMRSX3qUiIiIiJxe+UVf1eSs86CSy4JOxqR9KAKt4iIiMRlzhzfleTPf/a3ANQgSZH4qMItIiIiZVqzxj9Jsnp1GDsWatUKOyKR9KE+3CIiIlKqggJ/+7+5c2HSJNhzz7AjEkkvqnCLiIhIqYYOhfHjYdgw6NIl7GhE0o+6lIiIiEiJ3ngDbroJ+vSBK64IOxqR9KQKt4iIiMQ0b56vaOfkwGOPaZCkSEWpwi0iIiLbWbvWD5KsUgXGjfNPlBSRilEfbhERESnGOTjnHJg1CyZOhL32CjsikfSmCreIiIgUc9ddMGYM3HknHHNM2NGIpD91KREREZEib70F118Pp54KV18ddjQimUEVbhEREQHgu++gVy/Yf3944gkNkhSpLKpwi4iICOvXQ48evv/2uHFQp07YEYlkDvXhFhERyXLOwfnnw4wZ8PrrsM8+YUckkllU4RYREcly994Lzz4LQ4bAcceFHY1I5lGXEhERkSz27rtwzTX+ntvXXx92NCKZSRVuERGRLLVoEZx2GrRpA089pUGSIokSSoXbzI4zszlmNs/Mrovxfp6ZzQj+PjSz9mHEKSIi8VFeTx+jRo2iZcuWmNWmTZuvWbduM+PHQ926YUcmkrmSXuE2s6rAg0BXoB3Q28zaRc22AOjsnMsBBgHDkxuliIjES3k9fYwaNYp+/fqxaNEi4GE2b/4TBQW9+eyzUWGHJpLRwmjh7gjMc85955zbDDwPnBQ5g3PuQ+fcL8HLj4HmSY5RRETip7yeJgYMGMD69euB/sBZwC1s2jSWAQMGhByZSGYL4y4lzYDvI14vAQ4pZf5zgTdivWFm/YB+AI0aNSI/P7+SQkxta9euzYqyZks5QWWVtKe8voOSdV4sXrwYOAL4D/Ay/mKDn56sbZ0tOSBbygnZVdaKCqPCHWtIhos5o1kXfGI+PNb7zrnhBJcl27Zt63JzcyspxNSWn59PNpQ1W8oJKqukPeX1HZSs82L33Q9h6dIXgHnAmRTupj333DNp52W25IBsKSdkV1krKowK9xJgj4jXzYFl0TOZWQ4wAujqnPs5SbGJiEj5Ka+ngY0boUaNV4HqQGdgNQC1a9dmyJAhYYYmkvHC6MP9GdDazPYys+pAL+CVyBnMbE9gLHCGc25uCDGKiEj8lNdTnHNwySXw3XcNueKKL2nRYiNmRosWLRg+fDh5eXlhhyiS0ZLewu2c22pm/YE3garAE865mWZ2YfD+I8DNwG7AQ+ZvCrrVOXdQsmMVEZGyKa+nvkcfhSeegBtvhEGDOvOf/ywMOySRrBLKo92dc68Dr0dNeyTi/+cB5yU7LhERqRjl9dT1wQdw2WVw/PEwcGDY0YhkJz1pUkREJEMtWwYnnwwtWsCoUVC1atgRiWSnUFq4RUREJLE2bfKV7TVrYNIk2HXXsCMSyV6qcIuIiGSgyy+Hjz6CF16AP/4x7GhEspu6lIiIiGSYESP8QMnrroNTTgk7GhFRhVtERCSDfPKJvwXgscfC4MFhRyMioAq3iIhIxvjhB+jZE5o1g+ee0yBJkVShPtwiIiIZYPNm333kl1983+0GDcKOSEQKqcItIiKSAa68EqZO9S3bOTlhRyMikdSlREREJM09+SQ8+CBcdRX06hV2NCISTRVuERGRNPbZZ3DRRXDkkXDHHWFHIyKxqMItIiKSplasgB49oGlTGD0aqqmjqEhK0qkpIiKShrZsgVNPhZUr4YMPoGHDsCMSkZKowi0iIpKGrrkGJk+Gp5+GDh3CjkZESqMuJSIiImnm//4P7rsPLrsMzjgj7GhEpCyqcIuIiKSRL7+Efv2gc2e4++6woxGReKjCLSIikiZWroTu3X1/7RdegJ12CjsiEYmH+nCLiIikga1b/T22f/gBpkyBxo3DjkhE4qUKt4iISBq4/np45x3/kJuDDw47GhEpD3UpERERSXHPP+/7a19yCfTtG3Y0IlJeqnCLiIiksK++gnPOgcMPh2HDwo5GRCpCFW4REZEUtWqVHyRZvz6MGQPVq4cdkYhUhPpwi4iIpKBt26B3b1i61D/gpmnTsCMSkYpShVtERCQF3XgjvPUWDB8Ohx4adjQisiPUpURERCTFvPgi3HGHf8DN+eeHHY2I7ChVuEVERFLIzJn+TiSHHgr33x92NCJSGVThFhERSRG//grdukHduvDSS1CjRtgRiUhlUB9uERGRFFBQAHl5sHAh5OfD7ruHHZGIVBZVuEVERFLAwIHw+uvw0ENw2GFhRyMilUldSkREREI2depuDBoEZ58NF14YdjQiUtlU4RYREQnRN9/A0KH7cfDBvnXbLOyIRKSyqcItIiISktWr/SDJGjUKeOklqFkz7IhEJBFU4RYREQlBQQGceSbMmwc33zyTPfYIOyIRSRQNmhQREQnBkCHw8stw332Qk/Nb2OGISAKphVtERCTJXnsNbrkFTj8dLr007GhEJNFU4RYREUmiuXOhTx844AAYPlyDJEWygSrcIiIiSbJmDXTvDjvtBOPGQa1aYUckIsmgPtwiIiJJ4Bz07etvAzhpErRoEXZEIpIsqnCLiIgkwdChMHYs3H03HHlk2NGISDKpS4mIiEiCTZwIN94IvXvDlVeGHY2IJJsq3CIiIgk0f76vaOfkwIgRGiQpko1CqXCb2XFmNsfM5pnZdTHeNzO7P3h/hpl1CCNOERGJj/J6bOvW+UGSVar4QZK1a4cdkYiEIekVbjOrCjwIdAXaAb3NrF3UbF2B1sFfP+DhpAYpIiJxU16PzTk45xyYOROeew722ivsiEQkLGG0cHcE5jnnvnPObQaeB06Kmuck4GnnfQzsamZ/SHagIiISF+X1GO6+G154AW6/HY49NuxoRCRMYdylpBnwfcTrJcAhcczTDFgeOZOZ9cO3lABsMrP/VW6oKashsDLsIJIgW8oJKmsmaht2AEmkvF6K667zf2XIlvMCsqes2VJOyK6yVii3h1HhjjVcxFVgHpxzw4HhAGY2zTl30I6Hl/qypazZUk5QWTORmU0LO4YkUl7fQSpr5smWckL2lbUiy4XRpWQJsEfE6+bAsgrMIyIiqUF5XUSkFGFUuD8DWpvZXmZWHegFvBI1zyvAmcGo9kOB35xzy6NXJCIiKUF5XUSkFEnvUuKc22pm/YE3garAE865mWZ2YfD+I8DrwPHAPGA9cHYcqx6eoJBTUbaUNVvKCSprJsqWciqvVw6VNfNkSzlBZS2TObddFzoREREREakketKkiIiIiEgCqcItIiIiIpJAaVfhzpbHB8dRzrygfDPM7EMzax9GnJWhrLJGzHewmW0zs5OTGV9liqesZpZrZtPNbKaZTU52jJUhjuN3FzN71cy+CsoZT3/elGRmT5jZipLuF50pOSmRsiWvQ/bkduX17eZJ+7wO2ZPbE5LXnXNp84cfjDMf2BuoDnwFtIua53jgDfw9Xw8FPgk77gSV869A/eD/XdOxnPGWNWK+d/EDr04OO+4E7tddgVnAnsHrxmHHnaBy3gD8O/h/I2AVUD3s2CtY3iOADsD/Sng/7XNSChwvGbENsyW3K69nXl4vR1kzIrcnIq+nWwt3tjw+uMxyOuc+dM79Erz8GH9P23QUzz4FuBR4CViRzOAqWTxl7QOMdc4tBnDOpWN54ymnA+qamQE745Py1uSGWTmcc+/j4y9JJuSkRMqWvA7Zk9uV14vLhLwOWZTbE5HX063CXdKjgcs7T6orbxnOxf/SSkdlltXMmgHdgUeSGFcixLNf2wD1zSzfzD43szOTFl3liaecDwD74R988jVwuXOuIDnhJV0m5KREypa8DtmT25XXi8uEvA7K7ZHKnZPCeLT7jqi0xwenuLjLYGZd8En58IRGlDjxlPVe4Frn3Db/ozltxVPWasCfgaOAWsBHZvaxc25uooOrRPGU82/AdOBIYB9gkplNcc6tTnBsYciEnJRI2ZLXIXtyu/J6cZmQ10G5PVK5c1K6Vbiz5fHBcZXBzHKAEUBX59zPSYqtssVT1oOA54Ok3BA43sy2OufGJyXCyhPv8bvSObcOWGdm7wPtgXRKzPGU82zgDuc7w80zswXAvsCnyQkxqTIhJyVStuR1yJ7crry+/TzpntdBuT1SuXNSunUpyZbHB5dZTjPbExgLnJGGv5IjlVlW59xezrmWzrmWwIvAxWmYlCG+4/dloJOZVTOz2sAhwOwkx7mj4innYnxrD2bWBGgLfJfUKJMnE3JSImVLXofsye3K68VlQl4H5fZI5c5JadXC7RL3+OCUEmc5bwZ2Ax4KWgi2OucOCivmioqzrBkhnrI652ab2URgBlAAjHDOxbwtUaqKc58OAkaa2df4S3PXOudWhhb0DjCz54BcoKGZLQFuAXaCzMlJiZQteR2yJ7crr2deXofsyu2JyOt6tLuIiIiISAKlW5cSEREREZG0ogq3iIiIiEgCqcItIiIiIpJAqnCLiIiIiCSQKtwiIiIiIgmkCreIiIiISAKpwi0iIiIikkCqcIuIiIiIJJAq3JL1zKyVmW0xs1ujpj9sZmvMLK2e8iYiku2U1yXVqMItWc85Nw8YAfzTzBoCmNnNwDlAd+fctDDjExGR8lFel1SjR7uLAGbWFJgPPAR8AwwHejvnXgg1MBERqRDldUkl1cIOQCQVOOd+MLN7gavw58VlSsoiIulLeV1SibqUiPzuW6AG8JFz7sGwgxERkR2mvC4pQRVuEcDMjgQeBT4CDjOz9iGHJCIiO0B5XVKJKtyS9cysAzAeP8AmF1gM3B5iSCIisgOU1yXVqMItWc3MWgFvAG8BlzrnNgO3Aseb2RGhBiciIuWmvC6pSHcpkawVjGD/EN/y8Tfn3KZgelXgf8Avzrm/hhiiiIiUg/K6pCpVuEVEREREEkhdSkREREREEkgVbhERERGRBFKFW0REREQkgVThFhERERFJIFW4RUREREQSSBVuEREREZEEUoVbRERERCSBVOEWEREREUmg/wcoDz9cYFGv4gAAAABJRU5ErkJggg==\n", "text/plain": [ "
\n", "\n", "Consider the simpler case of a *linear* (in $x$) model asked to pass through *three* data points. If the three points happended to lie along a line then we can solve the problem, i.e. we can invert for the model parameters. But if not there is no solution (at least no exact solution) that fits all the data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Least squares solution\n", "\n", "In the language of our Numerical Methods module we are now in the world of curve-fitting rather than interpolation - we can't find a model (curve) that fits all the data exactly, but can still obtain a potentially useful model if we appropriately define \"useful\".\n", "\n", "
\n", "\n", "Recall that in Numerical Methods we used the command `numpy.polyfit` to fit a polynomial to data where the polynomial wasn't of a high enough degree to go through all the data points (or some of the data points weren't at distinct $x$ values - i.e. two $y$ values for the same $x$, think a cloud of points).\n", "\n", "
\n", "\n", "Instead of requiring the curve to go through all data points exactly (as for the cases above), we used *least squares* fitting to construct a function $f(x;\\boldsymbol{a})$, e.g. a polynomial in $x$ of degree $N$, (equivalently finding the parameters $\\boldsymbol{a}$) which minimises the sum of the squares of the differences between $M+1$ data points provided and the polynomial approximation, i.e. it minimises the quantity:\n", "\n", "$$E := \\sum_{i=0}^{M} (f(x_i;\\boldsymbol{a}) - y_i)^2,$$\n", "\n", "where $f(x_i;\\boldsymbol{a})$ is the output of our model evaluated at point $x_i$, and the $y_i$ are the corresponding data values. \n", "\n", "Note that for the \"cloud of points\" example we will be in the case where $M+1 > N$.\n", "\n", "
![](\"https://upload.wikimedia.org/wikipedia/commons/b/b0/Linear_least_squares_example2.svg\")
\n",
"\n",
"*(Wikipedia: https://en.wikipedia.org/wiki/Linear_least_squares) We're computing the sum of the squares of the distances indicated in green.*\n",
"\n",
"\n", "\n", "
\n", "\n", "Note that this error $E$ can be defined equivalently as\n", "\n", "$$E=\\| \\boldsymbol{f}(\\boldsymbol{x}; \\boldsymbol{a} ) - \\boldsymbol{y}\\|_2^2\n", "= || V\\boldsymbol{a} - \\boldsymbol{y}||_2^2,$$\n", "\n", "where $\\| \\cdot \\|_2$ is the 2-norm or the Euclidean norm acting on vectors.\n", "\n", "
\n", "\n", "[Why the square of $E$ here? As the 2 norm (Euclidean norm) includes a square root. We know that $E$ can't be negative (why?) and so minimising $E$ is completely equivalent to minimising $E^2$, i.e. we are minimising the 2-norm of the difference between the observed data and the model's prediction of the data.].\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will show in a later lecture that the coefficients ($\\boldsymbol{a}$) of the polynomial that minimises $E$ are given by the solution to the linear system\n", "\n", "$$V^TV\\boldsymbol{a} = V^T\\boldsymbol{y},$$\n", "\n", "where $V$ is again the Vandermonde matrix. \n", "\n", "
\n", "\n", "This is a very important result/equation in this module we will return to.\n", "\n", "
\n", "\n", "[As we have seen, $V$ is no longer square in the case where $M+1 > N$. \n", "What is the shape of $V$, and indeed all the vectors and matrices appearing in this equation - is it \"dimensionally\" consistent?]\n", "\n", "Let's check that this is true by forming and solving the matrix system $V^TV\\boldsymbol{a} = V^T\\boldsymbol{y}$ for $\\boldsymbol{a}$ and comparing with the result we get using `numpy.polyfit`." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a = \n", " [ 0.75909819 -0.43193108 0.09889271 -0.00552147]\n", "\n", "poly_coeffs = \n", " [-0.00552147 0.09889271 -0.43193108 0.75909819]\n", "\n", "Our _a_ vector = output from np.polyfit (as long as we flip the order)? True\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhQAAAGKCAYAAABU/thLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABQ3UlEQVR4nO3dd3xT1f/H8ddp2RtkqCAtCDjw6wJBFLGCAoKi4gCtfgXBioKCE7E/FEdFvw4QBbEiOKh7i4gDKTJkqSgKDsAWGYqIILTI6vn9cdKSpukiaW/Tvp+PRx6Q25ubT05Obj456xprLSIiIiKhiPI6ABEREYl8SihEREQkZEooREREJGRKKERERCRkSihEREQkZEooREREJGQlmlAYYwYYY6wxplVJPk9xGGPGGGO6eh2HlD3GmFRjTKoHzzvCGNM3yPYxxpgyO6/bGHO0MeZzY8w/vs/5hV7HVJCino+MMbG+/QaE4TmDvreRIPDzYIw50VcnGwTZ1xpjxhzEc8T5HhsXSqxFeJ6DjW+AMeaaEggpbPzqdexBPjZsr68itlDcAyihkGBu8N1K2wgg2JfOFKBT6YZSLI8DLYHLcHHO9TacsNmEez0fhuFYIwj+3kaiE3HnzzwJBa68ppRqNKVjAFCmEwpcPe2Eq7fFNYAwvr5K4TqQeMsYY4DK1to9XsdSGkri9VprV4brWOFgrV0PrPc6jgIcA3xhrZ3ldSDhZK3dDSzyOo5IYq1VeXnEWvsn8KfXcQBgrS2xGy77sUCrQva7FvgW+BfYAjwHNAjYZxjwJbAV2Ib7wPcO2KcScD+wxu9Y84HOvr/bILcxBcTVBngH2Ow73jrgDaCS3z4nAfN8f98AjAbudUWbs0+s77kGBBw/zrc9zm9bd2AmLtvMBL4HbgWiAx6bBkzHZZc/AnuBi3x/OwF4H/gb2AUsAM4IePwpwKfAX77nWQtMKuR9qgaM88W0E/gd+AA4Op/3vQvwrm/fv4CJQPUg5XID7tfuZl8sM4DYYrzenr66sQvY7nvOo/we+5KvzsT4bTsM9yF8w29bKpAa5P25EHgGV/f+9pVBtK8M5wMZwA9AjyBl/CYuKdgF/AQ8GFAGaeStk8/7/jYGv3rk21YHeArYCOz2HfNmwASJu49v3y2+1zodqFeEz21l4AFfbHt8/z6AS+D8j5/rVsDxGvnK72ff+/sb8DLQtIjnkRa+9/B332teCzyR3/sWULbPh1AvAz+vZ+I+M9t97/m3wKAC4s73vS1KvS3guM/76lR7YKFf3ert+/stvuf+B3gPaHSQ56KccvUru8BbrN+5dYzfY8f4tv0HmON73zcB9wFRBT2vb3tf3Dk+E/fZfQNoXoSyicbV1ezzZyrQNkh8rXB16ldf+a0FngbqB7z+wNebXR4HXacpYj30O0+9iPsM7wa+A67M53ixftvScJ/3/sAqXH1dhu+7sAiv71DgBQ6cZzbhzsuNC3ptnrdQGGMewn1hTgBuB5riKsRxxpjTrLX7fbvG4prU0nCJw/nADGNML2vtR759RuJOronActzJtz0Hmug64T7Az+MqAxT8C3AGrjJfj3tDmwK98HUVGWMaAp/jTnRX4wr+dqB58Uohl5bAbOBJXJLSHvfhbATcGbDvWbhmyHtxX8ZpxpiTcQnON7hELRMYAnzmK8+vjDG1gI+BJbjKuANXvqcVEltVoDYHPrANcMnAImPM0dba3wP2nw68DkwCOgB3AzV9z+lvFO79Ggg0xn3pfmKMaWut3VvI6+2Ja/L7HOgH1MKdtOYbY0601m7wxdgJeNkY0wX3wZnuK5trC3nNAOOBt33H7wL8H64Ong08gksk/w942xgTY63d4ntcc9/reh5Xxm19ZdAS90EHuAiXQH6Le58hn18bxpgo32s92XecFUBvXDLWCLgr4CFP4OrwFcBRwP+A/bi6WpAXcN0YD+ISpk6+19fSd6yvfdveB5bikviCNMDV5VG+13Y47jO/wFdv/s3vgcaYFrh6molrbv8FOAKXeB+sotZL/zguAN7CJefX4c4HbYGYAp4n3/e2iPW2IHVwXzSP4k76icBbxpiJuB9CQ4EmuLo7Efd+huJD3Of+/4BLOXDeLKyZ/V1gKjAW6IH7wZXFgfLIwxgzBPflPg1XJrV9+881xhxvrd1RwPONwX0OHgc+wZ0/3w+y3+G+1zAC9yOhpe9xMznQzXgDrq5E495zcEkahFCn/RRYD40xNXHdiPV9sf0GXAm8ZIypYa1NLuT4Z+A+96N9sd6P+86MtdZuK+T1vYSr27f7nrcJ0A2oUeAzFuUXwsHeKKSFAvclth+4O2D76b7HXZjP46JwJ/RPgPf8ts8A3i4kJgs8UITYG/r27VPAPkm4X3DN/bbVxJ1sbMDrLNKvgoC/G9/rTMRVev/MPg13kj004DGzcRlplYCsfRXwru9+e9/zHh/i+xvtq2A7gJuDvO+TA/ZP9L3fbQLKZWXAa8t+/wcV4fUuw33J+LcatcC1YDzut62D7726F/fh3EfeVptUgrdQTA3Y72vfdv9s/3jftqsLeS+vxJ1QDwl4bdODPGZMQD06L596NAWXzDYMiPuFgP2ewp1YTLAYffscR5CWO9wXSa46gzshP3+Q9eYI3/EuKmTfF3G/4A4vYJ9c71tAuT7vd7+49XKA33uX5qtrUYW9viAxBHtvi1Rv8znm8774ugSpfz/h15qJ+2Ldm70t8LUFqetxfttylSsFnM8D6wwHWijuDNjvWdz5ol6w58UlVtvJ+5mLxX1+RxRQLvV9dSXw/R0ZrE4H7FMJ6Ozb76SAMpgf5jpd1Ho4LPA98W3/DPeDKjrgeLEB9e5vcre4ZJ/3ryjs9fnK8abi1HVrreeDMs/BJQcpxphK2TdgMS5T6pK9ozGmnTFmhjHmD9yXwV7f44/yO95SoJcxJskY09kYUyWE2P7CNYM9ZIy51hjTOsg+nYBF1tp12RustRm4boCDYow5zBjzjDEmHfcB2ov7ZVAP9+vd3yLr1ypgjKmOa5Z9A8jyK0+Dq4TZ5fkLruXlGWPMlcaYI4oR32XGmMXGmG249yEDdxI4KsjurwfcfxX3fncI2P6mtTYr+461dgHuyypwQGLg662J+7X+mrV2n9/jf8X9kjzTb9sS3C+ARFxSkWStnVfoC3Y+Crj/I5BhrZ0fsA3cSSU7vjrGmIeNMWtwX/h7cZm/AYLVp8J0wSUjrwRsnw5UIW95BQ4qXIFrZWpSyHNkHzPwOcCvTIvDGHO9MeZbY8xOXL3J/swEqzf+ugMzrLUbD+Z581HUepntKNyvtSn+9fRgFafeFiDDWvuF3/3s+veZPdCqm729Eq7p3AvByroWLnENphOu9SXwO2E97rV0yedx4LpXaubznLkYY6oYY+4yxvxojNmF+2xmnw8Kq5PZxzjYOp2tsHrYBdhgrU0N2G86rkXy2EKO/6W19m+/+yt8/xalBX0pcLsxZrgx5j++MWuF8jqhyP6CXI17Q/1vdYBDAHxfeLNxzUw34prmTwFm4fr1sz2Iaxbtg6scfxljpvm6JorFujTtHNwvibHAz8aYtcaY6/12Owz4I8jDg20rlK9J+33cL9EHcLNRTsG1hEDu1wp5mxsb4DLl0eQtz2FAfWNMlLV2O677YCOuuW2dMeZ7Y8zFhcR3PvAarrXjCqCjL74/g8QGecsh+37TQvbL3ha4X+DrrY/7cg7W7Po7eUejv8yBvsKJQR6Tn78D7u/BJWQ57IHBof7lMA3X3TQBV5dOwTVFB+5XVA2ArdYNGvT3u9/f/W0NuJ/9uIKeO/sYgWWa33MUyhhzI66efYbrG+8AnFqEWMCdA8I9MLWo9dI/BsIYR3HrbTDb/O/41b9gdRUOrr6FQ3HLOvs74TPynsP+w4H3IpjspCm/5/Q3FteKMh3XbdiBA7NxCi2rEOt0fnEFlk0D8q8j2X8vSK7Pv995oyjx9cN9F92BG7exwRhzt+87Kl9ej6H4y/dvd/J+EPz/3hOoC1xm3ch3AIwxufpzrOtvfxh42BhzKO6L+XFcs3y/4gZnrV0L/NeXnZ2A+1KeZIxJs27cxiaC/9oL3JbdnxbYYhL44TgS1yx1lbU25xei74s8aIgB97fhfsFOxDUV532A7xeWtXY5cLEv+2+P6wt83RhzgrX2+3yerz+w2lo7wC+2yuRfsZvgBiv63wc35iBwv2CPXR4YfsD9v33bDg3y+EM5UH+yk7UXcF8KtYHJlOB0PmNMNeACXDPrE37b/xPCYbcCDYwxVWzu2S3Zr/+vII85mOfIPuaaMD1Hf2C2tfbW7A2+sRFFkT12qSD/4n6ABAq1XvrHQBHiKKoi19sSUNRzUbg0wbX0+t+H/Ms6+7UPIPd7lK2g8RPZX775vb/++gMvWmsfyN7gG1tWVKHUaf+4CqqHWwne2hHOz3tQ1trNuB8/Q40xR+HGXd2L+/H4dH6P87qF4lPcF2Bza+2yILdffftlJw45A/SMMW1wfe1BWWt/t9ZOwWWQ/s1re4DqxQnSOstxo6fxO96XwKn+XQa+5szABOAP3K/DwGa+3gH3g73OykB8EePMwLXMnAB8HaxMgzxmn3VTvkbj6sMxBTxFDVzTnr+rcK0iwQQOBOuPe7+XBGy/xD/zNcacDjTDlW++fK/3K+BSY0xODMaYGFwrlv+6CKNwfaRX4GaKXGSMuY6SUxVXLnsDtg8Isu9uilYn5+Leo0sDtsfj6nU4pu5ll1n/gO3ZdfALiq8GecthYBEf+wlwnjGmoCb7dKCNfxenb/Bt7Xz2L2q9zPYzrk96cFGbfv3keW+LWW/DrajnomCyf+EW5/wZrKx34maKBbMQlzS0yuc74acCnus7XBdssOcMVNQ6md9nM5Q6na2wejgXaOY7H/q7AjeGYlUxny+YQs891tqfrLV34RLh/LqqgNJroehpjAmcAbDdWvupMeZh4ClfFjQXl0EfgWsinmKtnYNLCvYBLxpjHsM1bd2L67Py/yJ6Dzei+mvciz8J17rxDAesBHobY2b59tkYrH/WGHM8bpT8a7gumWjcl8E+3MhscNMHb8DNSBjDgVkeu/yPZa21xpjXgEHGmJ/xTfHCDUjytwp3ckwyxuzHVdibA2MrxC24k/7HxpjncFl7Q1yfbbS19k5jzHlAAm4E9q+4fsebcB/kgr7EZwEXGmPG4QbAtvM9bls++/cyxjyC+1LogOuOetFa+3PAfrWBd40xz+D6BsfixnkEbWUJMBo3VmCGMWYSrn/2XtzArscAjDEdcc2b91prv/RtmwQ8boz5wlobjg9mLtba7caYRcCtxphNuF+51xD8V+5K4Azf+/I7sMVamxZkv49wsy4mG2Ma4X7d9AIGA2PtgdklocT9gzHmFWCMr/VqIa5fezTwirX2u4M47CxgpDHmLtzJsitwSREfew/us7LQGPMg7rPYFOhprb3St8+ruPo81RjzPG5w4y24OhBMUeslkPP5HYGb6fO5MWYy7pfaMbhpdPcUEH9+722h9bYkFONcFEz2Oi1DjTEv4M5P39mC14K51vdjYSlulsdgXKvdtnzi+8cYczsw0VfHP8KVSVPc2JJUa+3L+Tx2m+/clGiM2YF7f08BBgXZfRZwtTFmBa5O9SX4LLeVwA3GmH64FrsdvqQmlDqdrbB6+DwwHDd7LBHXuhqP+268LmCszMHK8/pw9fQzIIUDU/QvwHXVfVLg0Yo7irM4N/Kfu2yB7/32uwr36yoDl72uwo1Ib+a3z2W+F/cv7kTaH1fgaX773Oo7zl8cmJs9Bt/8ed8+p+N+HfxLASN/cX15L3BgnvFWXMITuNZA9jTNfNeh8O1XDzcgb4vvWJNxH+TAkdUn4r40MnEV6D7chzDYKN48o8d9fzsGd5LdjEty1uP6w3r5/n4ULlH61Rf3n7jpUh0LeT+jcGM7Nvrim4tL2tLIf77/e773dCuFr0Pxp++4HwItAp67oNcbOJ//PXzz+XHJyhpfrP4zSarhBiktB6raAyOeU/32ifPFd3bA8z0PrA8Sh8VvBpHv9X2E+5BuxtXpYO/50b46lEnR16HYhGuV+Jn816EIjDv7fYkNVo5++2WvQ5GOO5mk47cOhd9+RZrlgfsF9LTv/d2BS0ZbUMjIe7/HH4kbiJo9F38tMC5gn+twSeguXBLUjtDr5YCA5+iKW1Nhp+/2LTCwkNiDvreF1dtCjlmk+hfwmlv5batH0c5FqQTMnsF96W3AzUbIqUuB7yUHZnkc5yuzXbgvqvsp2joUvXyP+8f32NW46afHFlI22etQ/O57XCpu8GJgfA1x58i/fbcUXPKR633HdS/MxNVby4F1Gg66The1Hvr2PczvvSr2OhT51JExBb0+XOvqM7jv2Z2+92ApfrND8rsZ30EljHytFfdYa4vbPFouGHcNhGlAa2vt6gL2i8UlNdda1z0lIuVA9jkQl4QGdpNWaEU9P0Yir8dQiIiISDmghEJERERCpi4PERERCZlaKERERCRkSihEREQkZF6vlBk2DRs2tLGxsWE7XkZGBjVr1gzb8cRRuYafyrRkqFzDT2Uafl999dUWa20jr+OAcpRQxMbGsmxZnoUgD1pqaipxcXFhO544KtfwU5mWDJVr+KlMw8+4C0mWCeryEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQlZtZHgXJyspiy5YtbNu2jf37i3bF17p167JqVdival3hqVzDr6yUaXR0NPXq1aNhw4ZERem3ikhFUyESivXr12OMITY2lsqVK2NM4RcB3bFjB7Vr1y6F6CoWlWv4lYUytdayd+9e/vjjD9avX0/z5s09jUdESl+F+BmRkZFB06ZNqVKlSpGSCREpHmMMVapUoWnTpmRkZHgdjoh4oEIkFICaYEVKgT5nIhWXPv0iIiISMiUUIiIiEjIlFCLFlJWVxXXXXcchhxyCMYbU1FQGDBjAeeedd1DH+/vvv2nSpAlr1qwJW4xFiac4MW/evJlGjRqxfv36cIQnIuVQhZjlEakGDBjAli1bmDFjRoV8/rJq5syZTJs2jdTUVFq2bEmDBg046aSTsNbm7BMXF8dxxx3HU089VejxHnzwQXr16sWRRx5ZkmHn8cQTT+SKuSCNGzfmv//9L/fccw/PPfdcCUcmIpFICYWUW3v27KFKlSphP+7q1as57LDDOO2003K2HezzZGZmMmXKFD744INwhVdkdevWLdb+AwcOpH379jzyyCM0aNCghKISkUilLo8ItnLlSnr37k3t2rVp3Lgxl19+Ob///nvO35cuXUr37t1p2LAhderUoXPnznz55Ze5jvHMM8/Qpk0bqlWrRqNGjejRowf79u1jzJgxvPDCC3z44YcYY3Ka9oNZsWIF3bp1o06dOtSuXZsTTjiBOXPm5Px91qxZHH300VSrVo0ePXrw8ssvY4whLS0NgOeff55atWrlOmZqairGGLZs2QLAX3/9xeWXX06zZs2oXr06bdu2Zdq0abkeExcXx/XXX89tt91Go0aNOP3004tUToXF72/AgAHcfPPNrFu3Lmdtk+zt2d0HAwYMYO7cuUycODGn7LJfa6CZM2cSFRWVE2u2H3/8kT59+lC3bl1q1apFp06dWLFiRZ7nyjZmzBiOO+64PMd/4IEHaNKkCbVq1WLgwIHs2rUr12vxP461lscee4zWrVtTtWpVmjVrxqhRo3L+ftxxx3H44Yfz9ttvB30tIlKxVdgWihEjYPny/P++f391oqPD+5wnngjjx4fnWJs2baJLly4MGjSIRx99lL1795KYmEifPn1YtGgRUVFR7Nixg6uuuoonnngCYwxPPfUUvXr14pdffqFhw4YsW7aMoUOH8sILL9C5c2e2bdvG559/DsBtt93GqlWr2Lp1Ky+99BJAvr9Kr7jiCk444QSWLFlCpUqVWLFiBdWqVQPgt99+48ILL+Taa69l6NChLF68mDvuuKPYr/fff//l5JNPZuTIkdSpU4fPPvuM6667jubNm9OtW7ec/aZPn05CQgLz5s3DWlukcioo/kBPPPEEMTExTJ06laVLlxIdpJI88cQT/Pzzzxx99NE8+OCDADRq1Cjo8ebNm0e7du1yrY+yceNGOnfuzOmnn86nn35KvXr1WLJkSZFXec02d+5cqlevzuzZs9mwYQPXXHMNI0eOZMKECUH3v+uuu3j66ad5/PHH6dKlC3/++SfffPNNrn06dOjA3LlzGTx4cLFiEZHyr8ImFJHu6aef5oQTTuDhhx/O2fbiiy/SoEEDli1bRocOHejatWuuxzz55JO89dZbzJo1iyuvvJJ169ZRs2ZN+vTpQ+3atYmJieGEE04AoFatWlSvXp2qVaty6KGHFhhLeno6t912G0cffTQArVq1yhVn8+bNmTBhAsYYmjZtym+//cbo0aOL9XqbNm3K7bffnnM/ISGBzz//nFdeeSVXQtGiRQsee+yxnPt33313oeVUUPyB6tatS+3atYmOjs63XOrWrUuVKlWoUaNGkcrusMMOy7Vt4sSJ1KxZkzfeeCOnK6VNmzYFHieY6Ohopk2bRq1atTjuuON4+OGHGTRoEGPHjqVmzZq59t25cyfjxo1j/PjxXHPNNYArh06dOuXa7/DDD2fp0qXFjkVEyr8Km1AU1lKwY8cuz5czLshXX33FF198kaerAGDNmjV06NCBzZs3M3r0aObMmcMff/zB/v372bVrF+vWrQPgnHPOISYmhhYtWtCjRw+6d+9O3759i/26b7nlFgYPHswLL7xAt27duPjii3O+nFetWsWpp56a6xd44JdUUezfv5+HHnqI1157jQ0bNrB792727NlDXFxcrv3atWuX635Ryqmg+Evarl27aNKkSa5t33zzDZ07dw55/Mfxxx+f63V36tSJPXv2sGbNGo4//vhc+65cuZLdu3fnSs6CqV69eq5uExGRbBpDEaGysrLo3bs3y5cvz3X75ZdfcvrFr776apYuXcq4ceNYuHAhy5cvp1mzZuzZsweA2rVr8/XXX/P666/TvHlzxo4dy9FHH83GjRuLFcuYMWNYuXIlF154IQsXLuT4449n6tSpAEWaRRAVFZVnv7179+a6/+ijj/LYY49x++23M3v2bJYvX86FF16Y81qyBf7yLko5FRR/SWvYsCF///13rm2FlVlRyqu4ijrbY+vWrfl234hIxaaEIkKdfPLJ/PDDD8TExNCqVatct+wWhvnz53PjjTfSu3dv2rZtS+3atdm0aVOu41SqVImuXbsyduxYvvvuOzIyMnKmiVapUqXI/fatW7fmpptu4sMPP2TQoEFMmTIFgGOPPZbFixfn+sJatGhRrsc2atSIzMxM/vnnn5xtywMGuMyfP5/zzz+fq666ihNPPJEjjzySn3/+OSzlVFD8B6uoZXfSSSexcuXKPDHPnz8/T7KUrVGjRnnex8DyAjfY1P+6GosWLaJKlSpBp6cee+yxVK1aldmzZxcY7/fff8/JJ59c4D4iUjEpoSjj/vnnnzy/rtPS0hg6dCjbt2+nX79+LF68mLVr1/LZZ5+RkJDAjh07ANfvPn36dFauXMnSpUvp379/rmb0GTNm8MQTT/DNN9+Qnp7Oyy+/zI4dOzjmmGMAiI2N5fvvv+enn35iy5YtQX8F79q1i6FDh5KamkpaWhqLFy9m/vz5HHvssQAMGTKEtLQ0RowYwU8//cS7777L5MmTcx2jY8eO1KxZk1GjRrF69WreeustJk2alGufNm3aMHv2bObPn8+PP/7IsGHD+PXXXwstv8LKqbD4D1ZsbCxLliwhLS2NLVu2kJWVFXS/Hj16sGrVKv7666+cbTfccAM7d+7ksssuY+nSpaxevZpXXnklJ2no2rUr33zzDVOnTmX16tWMHz+eBQsW5Dn2vn37uOaaa/jhhx/49NNPufPOO7n22mvztOKAa60aPnw4o0aNYtq0aaxZs4YlS5bw9NNP5+yTmZnJV199Rc+ePUMqGxEpp6y15eLWrl07m5+VK1fm+7f8/PPPP8V+TLhdffXVFshzu/jii6211v7888/24osvtvXq1bPVqlWzbdq0scOGDbO7d++21lq7fPly26FDB1utWjXbsmVL++KLL9q2bdvae+65x1pr7bx582xcXJxt0KCBrVatmm3btq2dOnVqzvNv3rzZnnPOObZWrVoWsHPmzMkT4+7du+3ll19umzdvbqtUqWIPO+wwe+2119rt27fn7DNjxgzbpk0bW7VqVduxY0c7ffp0C9hff/01Z593333Xtm7d2larVs12797dvvTSSxawf/75p7XW2q1bt9qLLrrI1qpVyzZq1Mjefvvt9vrrr7dnnnlmzjHOPPNMO3To0DwxFlRORYk/0COPPGJjYmLyvFe9e/fOuf/TTz/ZU0891VavXj3Paw106qmn2qeeeirXtu+//96ee+65tmbNmrZWrVq2U6dOdsWKFTl/v+eee+yhhx5q69SpYwcNGmRHjRpl27Ztmyeee++91zZq1MjWrFnT/ve//7UZGRn5xrx//347duxY26JFC1u5cmXbrFkze9ddd+X8/eWXX7ZHHXVUvq8j28F83sqiYPVdQqMyDT9gmS0D38HWWiUU+SkLCUV59M8//9ilS5cW+iVbkXz00Ue2TZs2dt++fQf1+NKqq6eccopNSUkpdD8lFKVn+vTpNiYmxhpjbExMjJ0+fbrXIRUoEso00pSlhEJdHiIe69mzJ0OHDi3T18nYvHkzl1xyCZdffrnXoYhPSkoKCQkJpKenY60lPT2dhIQEUlJSvA5NKiglFCJlwE033URMTIzXYeSrcePG3HHHHbmm/4q3EhMTyczMzLUtMzOTxMREjyKSiq7CrkMh3mnfvr3rbxORg5a9nkxRt4uUNLVQiIhEoObNmxdru0hJU0IhIhKBkpKSqFGjRq5tNWrUICkpyaOIpKJTQiEiEoHi4+NJTk4mJiYGYwwxMTEkJycTHx/vdWhSQWkMhYhIhIqPj1cCIWWGWihEREQkZEooREREJGRKKOSgDBs2LM+lw0VEpOJSQlGGDRgwAGMMxhgqVapE8+bNuf766/Nc7joSpKWlYYzh66+/9joUEREpAUooiiElJYXY2FiioqKIjY0tlSVuzz77bDZt2kRaWhpTpkzhgw8+4IYbbijx5xURESkOJRRF5NW6+VWrVuXQQw+lWbNmdO/enX79+vHJJ5/k/H3//v0MGjSIFi1aUL16dVq3bs3//ve/nMtlr1q1CmMMv//+O+CW5q1SpQrnnntuzjGeffZZWrdunW8M+/fv57bbbqN+/frUr1+fESNGsH///lz7zJo1izPOOIP69evToEGDnMtyZ2vRogUAcXFxGGNyukuWLl1K9+7dadiwIXXq1KFz5858+eWXoRWaiIiUOiUURVQW1s1fu3Yts2bNonLlyjnbsrKyaNq0Ka+//jqrVq0iKSmJBx98kGnTpgFwzDHH0KRJE1JTUwFYsGABdevWZf78+ezbtw+A1NTUAsdDPPbYYzz77LM888wzfPnll+zfvz9PIpWRkcGIESNYsmQJqamp1K1bl/PPP589e/YAsGTJEgDefvttNm3axNtvvw3Ajh07uOqqq5g3bx5LlizhxBNPpFevXmzZsiUsZSYiIqXE68udhutW0pcvN8ZYIM/NGFPsYxfV1VdfbaOjo23NmjVttWrVcp7z8ccfL/BxI0eOtN26dcu5f9lll9mEhARrrbV33XWXHTJkiI2JibELFy601lrbtGnTAi97fNhhh9kHHngg5/7+/ftt69at7ZlnnpnvY3bu3GmjoqLsvHnzrLXW/vrrrxawqampBcaelZVlDz30UPvSSy8VuJ8cUFqXLy8qXb5c8qMyDT/K0OXLtbBVETVv3pz09PSg20tSly5dSE5OZteuXTz77LOsWbOGm266Kdc+kydPZsqUKaSnp7Nr1y727t2b68qVcXFxjB8/HnCtEcOHDyczM5PU1FQaNmzIhg0b8m2h2L59O5s2baJTp04526KioujYsSO//fZbzrY1a9YwevRoFi9ezJ9//klWVhZZWVmFXqho8+bNjB49mjlz5vDHH3+wf/9+du3apQsciUipysqC7dthyxZ327oV9u4Fa93frIUePaB2ba8jLbuUUBRRUlISCQkJubo9SmPd/Bo1atCqVSsAJkyYwFlnncX999/PmDFjAHjttdcYMWIEjz76KKeddhp16tRh4sSJvPPOOznHiIuL44YbbuCXX35h2bJlxMXFkZGRwSuvvELDhg1p1aoVTZs2DSnO888/n6ZNm/LMM8/QtGlTKlWqxLHHHpvT5ZGfq6++mj/++INx48YRGxtL1apV6datW6GPExE5GFlZ8MsvsGQJrFjh/v/LL7BmDfz7b8GP/eknJRQFUUJRRNnL2yYmJrJu3TqaN29OUlJSqS97e88993DuueeSkJDA4Ycfzvz58+nYsSPDhg3L2WfNmjW5HpM9jiIpKYlWrVrRuHFjzjrrLIYNG0a9evUKHD9Rt25dDjvsMBYtWkTXrl0B1022ZMkSDjvsMAD++usvVq1axcSJEznrrLMA+Prrr3PGaABUqVIFIM9gzvnz5zNhwgR69+4NwB9//MGmTZsOsnRERHLLyoJvvoFPPoE5c2DpUti2zf2tShU48kho3Rp69oRmzaBhQzjkEGjQwP09KgqMcTe/hl8JQglFMZSFdfPj4uJo27YtDzzwAJMmTaJNmzY8//zzfPTRR7Rq1YpXX32VuXPnUr9+/VyPO/PMM5k+fTpDhgwBIDY2lkaNGvH222/zwgsvFPicw4cPZ+zYsbRp04b//Oc/TJo0iU2bNuUkFPXr16dhw4Y8++yzHHHEEWzYsIHbb7+dSpUOVK/GjRtTvXp1Zs+eTdu2balWrRp169alTZs2TJ8+nY4dO5KRkcEdd9yRk3yIiByMXbtg1ix4802XSGSP8T7+eOjXDzp0cLdjjoHoaG9jLU80yyMC3XLLLTz33HOkp6dz3XXXcdlll3HFFVdwyimnkJaWxq233prnMWeddRb79+/P1RoRFxeXZ1swt956KwMHDmTw4MF07NiRrKysXIlVVFQUr732Gt999x3HHXccQ4cO5f7776dq1ao5+1SqVIkJEybw4osvcvjhh3PBBRcAMHXqVHbu3Em7du3o378/11xzDbGxsSGVj4hUPPv2wYwZEB8PjRtD374umejZE156CTZtgm+/hcmT4Zpr4LjjlEyEm3GDRCNf+/bt7bJly4L+bdWqVRxzzDHFOt6OHTuorc6ysFO5hl9ZK9OD+byVRYVNp5biK4kyXbcOnnvO3TZscF0Vffu6loi4OKhUztvhjTFfWWvbex0HqMtDREQijLUwfz48+ih88IHb1qMHPPkknHce+C3VI6VICYWIiESE/fvhnXdcIrF4sRs8edddMHgwqKfUe0ooRESkTMvKgrffhnvugZUr3cyMiRNhwACoUcPr6CSbBmWKiEiZZK3r0jj5ZLj0Unf/1VfdehA33KBkoqypMAlFeRl8KlKW6XMm4fL113DWWdCnD2RkuJkaK1a4wZaanVE2VYiEonLlyuzatcvrMETKvV27duW6eJ1IcW3c6Loy2rd33RtPPw2rVsGVVyqRKOsqRELRuHFjNmzYQGZmpn5BiZQAay2ZmZls2LCBxo0bex2ORKB9+2DcODjqKHjlFbjtNrck9pAh5X/qZ3lRId6mOnXqALBx40b27t1bpMf8+++/VKtWrSTDqpBUruFXVsq0cuXKNGnSJOfzJlJUixfDdde5had69YIJE9zAS4ksFSKhAJdUFOdEl5qaykknnVSCEVVMKtfwU5lKpNq5E+68EyZNgsMOc0tl9+3rrpshkafCJBQiIlJ2pKa6JbDT0uDGG+GBB3Qlz0hXIcZQiIhI2ZCRAcOGuRkc0dEwdy488YSSifLAk4TCGNPTGPOTMWa1MebOIH+va4z5wBjzrTHmB2PMQC/iFBGR8Pnpp1qcfLLr4hg+3I2ZOOMMr6OScCn1hMIYEw1MBM4FjgUuN8YcG7DbUGCltfYEIA54zBija1qLiESgrCz43/9g2LCTyciA2bNh/PiyuTBVSkoKsbGxREVFERsbS0pKitchRQwvxlB0AFZba9cCGGNeBS4AVvrtY4HaxhgD1AK2AvtKO1AREQnNpk1w1VUuiejS5S/eeacRDRp4HVVwKSkpJCQkkJmZCUB6ejoJCQkAxMfHexlaRPCiy6Mp8Jvf/fW+bf6eAo4BNgIrgOHW2qzSCU9ERMIhNRVOOgkWLoQpU2DMmB/KbDIBkJiYmJNMZMvMzCQxMdGjiCKLFy0UwSYEBa421QNYDnQFjgQ+NcbMs9b+k+tAxiQACQBNmjQhNTU1bEHu3LkzrMcTR+UafirTkqFyPXhZWfDKK82ZOrUFTZvuYuLEH2jRIqPMl+m6devy3V6W4y4rvEgo1gNH+N1vhmuJ8DcQeMi6ZS1XG2N+BY4GlvjvZK1NBpIB2rdvb+Pi4sIWZGpqKuE8njgq1/BTmZYMlevB2bbNLZP94YfQvz8kJ9egdu1TgLJfps2bNyc9PT3o9rIcd1nhRZfHUqC1MaaFb6Blf+D9gH3WAd0AjDFNgKOAtaUapYiIFMuqVdChA3z8MTz1FLz8cmRNB01KSqJGwEjRGjVqkJSU5FFEkaXUEwpr7T5gGPAxsAp43Vr7gzFmiDFmiG+3+4HTjDErgNnASGvtltKOVUREiua996BjR/jnH5gzB4YOjbwVL+Pj40lOTiYmJgZjDDExMSQnJ2tAZhF5slKmtXYmMDNg22S//28Eupd2XCIiUjzWulUu777bXSH0nXegWTOvozp48fHxSiAOkpbeFhGRg/LvvzB4MKSkuHETyclQvbrXUYlXlFCIiEix/fknXHihmxKalASjRkVeF4eElxIKEREpllWroHdvt2jV66/DpZd6HZGUBUooRESkyObPhz59oEoVd2GvDh28jkjKCl1tVEREiuStt+Dss6FxY/jySyUTkpsSChERKdSTT7qujZNPhgULoEULryOSskYJhYiI5MtaSEyEm26CCy5wF/k65BCvo5KySGMoREQkqP373QJVzzwDCQkwaRJER3sdlZRVaqEQEZE89uyB+HiXTIwaBZMnK5mQgqmFQkREcsnMhIsvhlmz4JFH4LbbvI5IIoESChERybFjB5x/PsybB1OmwKBBXkckkUIJhYiIAO7S4+eeC0uXuuW0+/f3OiKJJEooRESELVuge3f4/nt48023rLZIcSihEBGp4DZvhm7dYPVqeP996NnT64gkEimhEBGpwDZvhq5dYe1a+PBD93+Rg6GEQkSkgvrjD5dA/PqrSybOOsvriCSSKaEQEamAspOJtDSYORPi4ryOSCKdEgoRkQrmzz9zJxNnnul1RFIeKKEQEalAtm6Fc85xYyaUTEg4KaEQEakgtm1zU0NXrYIPPtCYCQkvJRQiIhXAjh1u0arvvoN33nGJhUg4KaEQESnndu1yy2kvXeoWrerd2+uIpDxSQiEiUo7t2eMu9PXFFzB9ulbAlJKjhEJEpJzat89dgvyjj9xlyK+4wuuIpDyL8joAEREJv6wsSEhwXRyPPeb+L1KSlFCIiJQz1sJtt8G0aXD33XDLLV5HJBWBEgoRkXLmoYdg3DgYNgzGjPE6GqkolFCIiJQjzz4Ld90Fl18OTzwBxngdkVQUSihERMqJt96CIUPcehPPPw9ROsNLKVJ1ExEpB1JT3SyOU091AzGrVPE6IqlolFCIiES4b7+FCy6AVq3ckto1angdkVRESihERCJYWprr4qhTB2bNggYNvI5IKiotbCUiEqG2bIGePd3S2vPnwxFHeB2RVGRKKEREIlBmprs+R3o6fPYZtG3rdURS0SmhEBGJMPv3uwGYixe7mR2nn+51RCJKKEREIoq1cNNN8N578OSTcNFFXkck4mhQpohIBHnkEZg0CW6/3a2EKVJWKKEQEYkQr7wCI0dC//5ueW2RskQJhYhIBJg3DwYMgC5dtAqmlE2qkiIiZdzPP8OFF0KLFvDOO1C1qtcRieSlhEJEpAz780/o1cu1SHz4oRaukrJLszxERMqof/91LRMbNsDnn8ORR3odkUj+lFCIiJRB1sLAgbBwIbzxBnTq5HVEIgVTl4eISBk0Zgy8+iqMHQuXXOJ1NCKFU0IhIlLGTJ8O993nWihGjvQ6GpGiUUIhIlKGzJsHgwZBXBxMngzGeB2RSNEooRARKSPWrHFLacfGumt0VKnidUQiRaeEQkSkDNi+3V09NCsLZszQ9FCJPJrlISLisX37oF8/+OUX+OQTaN3a64hEik8JhYiIx269FT7+GJKT4ayzvI5G5OCoy0NExEOTJ8OECXDzzXDttV5HI3LwlFCIiHjk88/dJch79XKXJReJZEooREQ8sHq1W7DqqKPcZcmjo72OSCQ0SihERErZ9u3Qp49bY+L996FOHa8jEgmdBmWKiJSi/fvh8ssPzOjQBb+kvFBCISJSiu68Ez76CJ5+WjM6pHxRl4eISCl58UV49FEYOhSGDPE6GpHwUkIhIlIKFi+GhATXKjFunNfRiISfEgoRkRK2YYO7Rsfhh8Mbb0Dlyl5HJBJ+niQUxpiexpifjDGrjTF35rNPnDFmuTHmB2PM3NKOUUQkHHbtcsnEjh1uRschh3gdkUjJKPVBmcaYaGAicA6wHlhqjHnfWrvSb596wCSgp7V2nTGmcWnHKSISKmvhuutg6VJ45x047jivIxIpOV60UHQAVltr11pr9wCvAhcE7HMF8La1dh2AtXZzKccYskWL4LLL4N9/vY5ERLwyfjy89BLcdx9ceKHX0YiULC8SiqbAb3731/u2+WsD1DfGpBpjvjLG/LfUoguTHTtcX+kHH3gdiYh44dNP4bbboG9fSEz0OhqRkufFOhQmyDYbcL8S0A7oBlQHvjTGLLLW/pzrQMYkAAkATZo0ITU1NWxB7ty5M6TjRUXBIYd0Yvz4HTRq9H3Y4op0oZar5KUyLRmhlOuGDdW4/vp2xMTsZvDgb/jii/3hDS5Cqa6Wb14kFOuBI/zuNwM2Btlni7U2A8gwxnwBnADkSiistclAMkD79u1tXFxc2IJMTU0l1OMNGABPPFGV446Lo2HDsIQV8cJRrpKbyrRkHGy57tjhLvhVuTJ89lllWrY8I/zBRSjV1fLNiy6PpUBrY0wLY0wVoD/wfsA+7wFnGGMqGWNqAB2BVaUcZ8iuugr27YPXXvM6EhEpDVlZcPXVsGqV+9y3bOl1RCKlp9QTCmvtPmAY8DEuSXjdWvuDMWaIMWaIb59VwCzgO2AJMMVaG3H9BiecAP/5D0yf7nUkIlIaHnzQzeZ49FE4+2yvoxEpXZ5cy8NaOxOYGbBtcsD9R4BHSjOuknDllTBypLsQUOvWXkcjIiXlgw/g7rvdZ37ECK+jESl9WimzhF1xhbtEsVopRMqvH390icRJJ0FysvvMi1Q0SihKWLNm0LWrSyhs4FwWEYl427e7NSaqVnXdHdWrex2RiDeUUJSCq66CtWvhyy+9jkREwikry32+16xx6840b+51RCLeUUJRCvr2db9aXnrJ60hEJJzuu8+NnXj8cTjzTK+jEfGWEopSULu2axJ97TUtxS1SXrz3Htx7r5smOmyY19GIeE8JRSkZMAD+/ttdbVBEItuqVa6ro317mDxZgzBFQAlFqenWzQ3QnDbN60hEJBTbt7vLkVerBm+/7f4VESUUpSY62rVSfPIJbNjgdTQicjCysuC//3WDMN98E444ovDHiFQUSihK0YAB7oT04oteRyIiB+OBB1y35eOPQ5cuXkcjUrYooShFRx7pTkJTp2pNCpFIM2MG3HOPGzuhQZgieSmhKGUDB8Lq1bBggdeRiEhR/fwzxMe7lTCfeUaDMEWCUUJRyi65BGrW1OBMkUixY4eb9l25slbCFCmIEopSVqsW9OsHr78OGRleRyMiBbHWjX366Sf3mY2J8ToikbJLCYUHBg6EnTvdKHERKbseeshNDf3f/9w1eUQkf0ooPHD66e5S5lOmeB2JiORnyZL6JCZC//5wyy1eRyNS9imh8IAxcO21MH++W3FPRMqWtWvhgQeO5T//cYm/BmGKFE4JhUeuvtoN8nr2Wa8jERF/GRluJUxwgzBr1vQ2HpFIoYTCI40bu5PWCy/ogmEiZYW1rvVwxQr4v/9bRcuWXkckEjmUUHgoIQG2bnWDvkTEe+PHwyuvuBUxO3TY6nU4IhFFCYWHzjrLrZ6ZnOx1JCIyZw7cfrtrORw1yutoRCKPEgoPRUW55tW5c+HHH72ORqTiWrcOLrvMzb56/nkNwhQ5GEooPDZgAFSqpMGZIl7ZtQv69oU9e+Ddd6FOHa8jEolMSig81qSJW9ZXgzNFSp+1cP318NVX8NJLcNRRXkckErmUUJQBCQnw118anClS2iZNcsn83XdDnz5eRyMS2ZRQlAHdukGrVu7kJiKlY/58GDECzjvPXZZcREJTaEJhjOlfGoFUZFFRrtl1wQJYvtzraETKvw0b3JV/W7RwXR1R+mklErKifIxeNMZ8bow5tsSjKcdSUlKIjY0lKiqK2NhYUlJScv194EB3WeSJEz0KUKSC2L0bLr7YXaDvnXegXj2vIxIpH4qSULQDKgPfGGMeNcbUKuGYyp2UlBQSEhJIT0/HWkt6ejoJCQm5kor69SE+HlJS4O+/PQxWpJy76SZYvNiNnWjb1utoRMqPQhMKa+0Ka+0ZQAJwJfCTMebyEo+sHElMTCQzMzPXtszMTBITE3NtGzrUTWF7/vlSDE6kAklOdrdRo1wrhYiET5F7Dq21LwBHAe8CLxlj5hhjlN8Xwbp164q0/cQT4bTTXLdHVlYpBCZSgSxaBDfeCD16wP33ex2NSPlTrKFI1trt1tqhwClAQ1w3yGPGmNolEl050bx58yJvHzYM1qyBTz4p6ahEKo5Nm9ziVc2awcsvQ3S01xGJlD9FSiiMMZWNMR2MMTcZY14G3gLaApWAocCPxhjN4s5HUlISNWrUyLWtRo0aJCUl5dn34ovdYlcanCkSHnv2uBkd27e7lTAbNPA6IpHyqSjTRhcC/wBfAo8BbYAPgH5AM6Ax8CrwpjFmSMmFGrni4+NJTk4mJiYGYwwxMTEkJycTHx+fZ98qVdz1PT78ENau9SBYkXLmpptg4UKYNg3+8x+voxEpv4rSQrETGAt0B+pZa9tba4dba9+w1m601v5jrb0V+D/grpIMNpLFx8eTlpZGVlYWaWlpQZOJbEOGuCbZCRNKMUCRcujZZ+GZZ2DkSHfxLxEpOUWZ5dHdWnuftXa2tTajgF2/wLVYSIiaNoV+/eC551wzrYgU38KFbuZU9+4QpHdRRMIsnOvDfQtcEMbjVWg33+wW3nnuOa8jEYk8Gze68UhHHAGvvKJBmCKlIWwJhbV2l7X2g3Adr6Jr1w7OOMN1e+zb53U0IpFj9243o2PHDg3CFClNWsG+DLv5ZkhPdydFESmcta6bI3slTA3CFCk9SijKsD59oGVLGDfO60hEIsPTT7tuwsRErYQpUtqUUJRh0dEHprwtWeJ1NCJl2xdfwPDh0KsX3Huv19GIVDxKKMq4a66BOnXUSiFSkHXr3OJVRx6plTBFvKKEooyrXdstdPXGG5CW5nU0ImVPZiZceKEbjPnee1C3rtcRiVRMSigiwPDhYAw8/rjXkYiULdbC4MGwfLlrmTjqKK8jEqm4lFBEgCOOgCuvhClTYMsWr6MRKTseecStM/HAA9C7t9fRiFRsSigixB13wK5d8OSTXkciUjbMnAl33umW1B41yutoREQJRYQ45hg3jfSppyCjoAXQRSqAH3+Eyy+HE0+EqVNdl6CIeEsJRQS5807YutV1fYhUVNu2wQUXQLVqbtG3mjW9jkhEQAlFROnUyS3H/fjjsHev19GIlL79+13LxK+/wltvQfPmXkckItmUUESYkSPdnPtXX/U6EpHSd8cdMGsWTJwInTt7HY2I+FNCEWF69XLXJxg7FrKyvI5GpPRMnepa52680a3NIiJlixKKCGOMu07BqlWuyVekIliwAIYMgbPP1nosImWVEooIdMklbtbHffeplULKv/R0uOgiiI2F11+HSpW8jkhEglFCEYGio+H//g++/16XNpfybedON6Njzx54/32oX9/riEQkP0ooIlS/ftCmjWulsNbraETCLyvLrRC7YoUbhHz00V5HJCIFUUIRoaKj3ViKb7+FDz7wOhqR8EtMdBf7GjcOevb0OhoRKYwSigh2xRXQsqVaKaT8efFFeOghuO46N6tDRMo+JRQRrFIl9yvuq6/cdQ1EyoMFC9y00K5d3bVrtKy2SGRQQhHhrrrKtVKMHq0ZHxL51q51MzqaN4c33oDKlb2OSESKSglFhKtcGcaMgW++0boUEtm2bYPzzoN9++DDD6FBA68jEpHi8CShMMb0NMb8ZIxZbYy5s4D9TjHG7DfGXFKa8UWaK66AY4+Fu+92J2ORSLN3r7sM+S+/wNtvuxlMIhJZSj2hMMZEAxOBc4FjgcuNMcfms9/DwMelG2HkiY6GBx5wl3SePt3raIJLSUkhNjaWrl27EhsbS0pKitchSRlhLdx0E3z6KSQnQ1yc1xGJyMHwooWiA7DaWrvWWrsHeBW4IMh+NwJvAZtLM7hIdeGF0L696/7YvdvraHJLSUkhISGB9PR0rLWkp6eTkJCgpEIAGD8eJk92F74bONDraETkYHmRUDQFfvO7v963LYcxpilwETC5FOOKaMZAUpJbpnjKFK+jyS0xMZHMzMxc2zIzM0lMTPQoIikr3nkHbr0VLr4YHnzQ62hEJBRerIofbBJY4CoK44GR1tr9poA5Y8aYBCABoEmTJqSmpoYpRNi5c2dYj1caKleG448/kdGja9Cy5SKqVy8b0z7WrVuX7/ZIK+OyKBLrKsCPP9ZmxIgTOfroDK69djlffFE26mu2SC3XskxlWr4ZW8orIhljOgFjrLU9fPdHAVhrx/rt8ysHEo+GQCaQYK19N7/jtm/f3i5btixscaamphIXgZ25CxZA585w//3ueh9lQWxsLOnp6Xm2x8TEkJaWVvoBlTORWFfT0uDUU6F6dVi0CJo08TqivCKxXMs6lWn4GWO+sta29zoO8KbLYynQ2hjTwhhTBegPvO+/g7W2hbU21lobC7wJ3FBQMiEHnH469O3rVhn8/Xevo3GSkpKoUaNGrm01atQgKSnJo4jES9u2Qe/ebqzPzJllM5kQkeIr9YTCWrsPGIabvbEKeN1a+4MxZogxZkhpx1MePfSQO1mPGeN1JE58fDzJycnExMRgjCEmJobk5GTi4+O9Dk1K2e7dLuHNnh56zDFeRyQi4eLFGAqstTOBmQHbgg7AtNYOKI2YypPWreGGG+Cpp9x0vGPzTMotffHx8cTHx6vJswKzFgYNgjlz4KWX4KyzvI5IRMJJK2WWU6NHQ+3acMcdXkci4vzf/0FKipuNdOWVXkcjIuGmhKKcatjQXTjsww9h9myvo5GKLjnZTQu99loYNcrraESkJCihKMduvBFiYuC222D/fq+jkYrq/ffh+uuhVy+YNElXDxUpr5RQlGPVqsHDD8Py5fDcc15HIxXRl19C//7Qrh289hpU8mTUloiUBiUU5dxll8GZZ8Jdd8HWrV5HIxXJjz+6q4c2beq63mrV8joiESlJSijKOWNgwgT4+283UFOkNGzcCD16uBaJjz+GRo28jkhESpoSigrg+ONh6FB3Aably72ORsq7v/+Gnj1di9jMmdCypdcRiUhpUEJRQdx7LzRo4AZqlvJq61KBZGbC+ee77o6333ZjJ0SkYlBCUUHUrw9jx8L8+fDyy15HI+XR3r3Qrx8sXOjWmzjnHK8jEpHSpISiArnmGjjlFHe56L//9joaKU+ysmDwYJgxw00NvfRSryMSkdKmhKICiYqCZ56BLVtg5Eivo5Hywlq31smLL8J998EQXZFHpEJSQlHBnHQS3HwzPPsszJvndTRSHtx3H4wbB8OHu+W1RaRiUkJRAY0Z41bQTEhwV38UOVjjx7v6NHAgPP64VsEUqciUUFRANWvC00+7kfgPP+x1NBKppk51rV0XX+yu1RGls4lIhaZTQAV17rlw+eXuyo8//uh1NBJpXn3VXeirRw83o0NLaouIEooKbNw411pxzTW6eJgU3TvvuMuPd+7s1pqoWtXriESkLFBCUYE1aQJPPeUu4PTYY15HI5Fg5ky31sQpp7gpojVqeB2RiJQVSigquMsvd33go0fD9997HY2UZZ99Bn37uqXcP/oIatf2OiIRKUuUUFRwxrgBmnXrwtVXu9UORQKlpkKfPtCmjbvYV716XkckImWNEgqhUSN34bCvv3bLc4v4S02FXr3cRb4++wwOOcTriESkLFJCIYBryo6Ph/vvh2XLvI5Gyoq5c6F3b5dMfP45NG7sdUQiUlYpoZAcTz4Jhx3mxlXs2OF1NOK1uXNdy0RsrJIJESmcEgrJUb++W1Ng7VoYOtTraMRLn33m1ipRMiEiRaWEQnI54wy45x546SV3k4pn5kw47zxo1QrmzHHTi0VECqOEQvJITIQuXeD66+GXX7yORkrTu+/ChRdC27YumVDLhIgUlRIKySM62nV9VK0K/fvDv/96HZGUhldfhUsugZNPhtmzNZtDRIpHCYUE1awZTJvmppLeeKPX0UgwKSkpxMbG0rVrV2JjY0lJSTnoYyUnwxVXwOmnw6efap0JESk+JRSSrz594K67YMoUePZZr6MRfykpKSQkJJCeno61lvT0dBISEg4qqXj4YbjuOjejY9YsrYApIgdHCYUU6L77oHt3GDYMlizxOhrJlpiYSGZmZq5tmZmZJCYmFvkY1sKoUXDnnW6q8DvvQPXq4Y5URCoKJRRSoOhoePllOPxwd82PzZu9jkgA1q1bV6ztgfbtg8GD4aGHYMgQN6OncuVwRigiFY0SCinUIYfAW2/Bli1w2WWwZ4/XEUnz5s2Ltd1fRoabyTF1qrso3KRJLnEUEQmFEgopkpNPdmMp5s6FhATXXC7eSUpKokbAtcNr1KhBUlJSgY/bsgW6dXNXC5082XVpGVOSkYpIRaGEQoosPt4tevXCC7qImNfi4+NJTk4mJiYGYwwxMTEkJycTHx+f72NWr3azOL791rU4XXddKQYsIuVeJa8DkMhyzz1usavERDjySOjXz+uIKq74+Hji4+NJTU0lLi6uwH0XLIALLnD//+wzl1iIiISTWiikWIyB555zX0hXXw0LF3odkRTmlVega1do0AAWLVIyISIlQwmFFFu1am6J5iOOcJe2/u47ryOSYKx1l6O/4go49VT48kt3fQ4RkZKghEIOSsOGbkXFmjXdOhWrV3sdkfjLyHDdUXffDVddBZ98oqW0RaRkKaGQgxYb65KKffvg7LNhwwavIxKA9HTXrfHWW/DII24QbdWqXkclIuWdEgoJyTHHuOWat26Fc86BP//0OqKKLTUV2reHtDT48EO47TZNCxWR0qGEQkLWvj28/z78+ivExcHvv3sdUcVjLfzvf66l6JBDYPFi6NnT66hEpCJRQiFhERfnfhGnpbn/b9zocUAVyM6d0fTtCyNHQt++sHQpHHWU11GJSEWjhELCpmtX1/2xYQOceSb89pvXEZV/X38NQ4a0Y8YMGDcOXntNVwsVEW8ooZCwOuMMN1Dzzz+hSxf46SevIyqfrHUJxKmnwu7d0cyZAyNGaLyEiHhHCYWE3amnwuzZburiaae5VRolfDZvhvPOg1tugV69YMqUpXTu7HVUIlLRKaGQEtGunVtI6ZBD3MWo3nrL64jKhxkz4IQTXMI2cSK88w7UrbvP67BERJRQSMk58ki3NPfJJ8Oll8Ljj+sqpQdr+3YYNAjOPx8aN4YlS+CGG9TFISJlhxIKKVENG7pf0xddBLfeCklJx5CZ6XVUkeWzz+D44+H5591F2ZYudfdFRMoSJRRS4qpXhzfegAcegM8/b0ynTrBmjddRlX1btrgLsJ1zjivDhQtdGVap4nVkIiJ5KaGQUhEV5X5dP/TQCn77zS2G9d57XkdVNlnrlss++mh4+WVXbt98Ax07eh2ZiEj+lFBIqerQYSvLlkGLFnDhhXDttbBzp9dRlR3Ll7uFwQYMcItTLV/uWiWqV/c2LhGRwiihkFLXsqWbATJyJDz3HJx4ortfkW3eDAkJbgDrypXwzDMwbx60bet1ZCIiRaOEQjxRtSo89JC7mNW+fdC5s1tXYccOryMrXRkZ8OCD0Lo1TJsGw4fDL7+45CJKn04RiSA6ZYmnunSB775zXR/jx7txA6+/Xv6nl+7eDU8+6VprEhMPlMO4cVCvntfRiYgUnxIK8VydOjB5suv2aNIE+vWDHj3cQMTyJjMTnnoK2rSBm25yl39fsAA++MD9X0QkUimhkDKjY0e3YNOECbBsmRtP0K9f+bgeyNatbnBlTAzceCM0awYffwxz5rjlyUVEIp0SCilTKlVyX7hr18Lo0e6S6Mce69Zj+PZbr6Mrvq++citcNmvmXk/Hjm6w5YIF0L27VroUkfJDCYWUSfXqwX33ucRi+HB48003G6RrV9c9kJXldYT5++svN0ujQwe33sarr8KVV7qEaMYMdCEvESmXlFBImda4sbsGyPr18L//werV0KeP6zoYORJWrPA6Qmf7drcI1fnnw6GHwpAhbgbHhAmwcSMkJ2u5bBEp35RQSESoXx9uv90t2f3666614vHH3Zf08ce7mRJffAF795ZOPFlZLpl55BG3EFXDhhAf7xaiuvlm+Ppr+P57131Tt27pxCQi4qVKXjypMaYn8AQQDUyx1j4U8Pd4YKTv7k7gemttBPagS7hVruyuXHrppe5aF6+/7roUHn7YredQu7b7gu/QwV1CvV0718oRiqws+O03WLXKXZhr4UJYtAi2bXN/P+EEl+z06uUGWGr9CBGpiEo9oTDGRAMTgXOA9cBSY8z71tqVfrv9Cpxprf3bGHMukAzoSgaSS8OG7hLeN9zguhxmz3YzJ1JT3TiLbI0bQ2zsgVujRlCrlks+atZ0CcPevbBnj5vWuXkz/P67u6WlwY8/knOFVGPc6pWXXQadOsHZZ7sBlyIiFZ0XLRQdgNXW2rUAxphXgQuAnITCWrvQb/9FgE7ZUqC6daFvX3cD+Ocft47FsmUuIUhLc/fffdclDoVp0MCtidGsGQwe7NaIOOYY1xqhhadERPLyIqFoCvzmd389Bbc+DAI+KtGIpNypUwfOPNPd/FnrBkvu2OFuGRkQHe26UqpUgWrVXAuGLhEuIlI8XiQUwWbeB11o2RhzFi6hCDrRzhiTACQANGnShNTU1DCFCDt37gzr8cSJhHL95RevIyieSCjTSKRyDT+VafnmRUKxHjjC734zYGPgTsaY44EpwLnW2r+CHcham4wbX0H79u1tXFxc2IJMTU0lnMcTR+UafirTkqFyDT+VafnmxXj0pUBrY0wLY0wVoD/wvv8OxpjmwNvAVdbanz2IUURERIqh1FsorLX7jDHDgI9x00anWmt/MMYM8f19MnA3cAgwybi1ifdZa9uXdqwiIiJSNJ6sQ2GtnQnMDNg22e//g4HBpR2XiIiIHBwtwSMiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIiIiFTQiEiIiIhU0IhIiIiIVNCISIlLiUlhdjYWKKiooiNjSUlJcXrkEQkzCp5HYCIlG8pKSkkJCSQmZkJQHp6OgkJCQDEx8d7GZqIhJFaKESkRCUmJuYkE9kyMzNJTEz0KCIRKQlKKESkRK1bt65Y20UkMimhEJES1bx582JtF5HIpIRCREpUUlISNWrUyLWtRo0aJCUleRSRiJQEJRQiUqLi4+NJTk4mJiYGYwwxMTEkJydrQKZIOaNZHiJS4uLj45VAiJRzaqEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkCmhEBERkZApoRAREZGQKaEQERGRkBlrrdcxhIUx5k8gPYyHbAhsCePxxFG5hp/KtGSoXMNPZRp+MdbaRl4HAeUooQg3Y8wya217r+Mob1Su4acyLRkq1/BTmZZv6vIQERGRkCmhEBERkZApochfstcBlFMq1/BTmZYMlWv4qUzLMY2hEBERkZCphUJERERCpoQiCGNMT2PMT8aY1caYO72OJ9IZY44wxswxxqwyxvxgjBnudUzliTEm2hjzjTFmhtexlAfGmHrGmDeNMT/66mwnr2MqD4wxN/s+/98bY14xxlTzOiYJLyUUAYwx0cBE4FzgWOByY8yx3kYV8fYBt1prjwFOBYaqTMNqOLDK6yDKkSeAWdbao4ETUNmGzBjTFLgJaG+tPQ6IBvp7G5WEmxKKvDoAq621a621e4BXgQs8jimiWWs3WWu/9v1/B+4E3dTbqMoHY0wzoDcwxetYygNjTB2gC/AcgLV2j7V2m6dBlR+VgOrGmEpADWCjx/FImCmhyKsp8Jvf/fXoyy9sjDGxwEnAYo9DKS/GA3cAWR7HUV60BP4Epvm6kaYYY2p6HVSks9ZuAB4F1gGbgO3W2k+8jUrCTQlFXibINk2FCQNjTC3gLWCEtfYfr+OJdMaY84DN1tqvvI6lHKkEnAw8ba09CcgANI4qRMaY+riW3hbA4UBNY8yV3kYl4aaEIq/1wBF+95uhprmQGWMq45KJFGvt217HU06cDvQxxqThuua6GmOmextSxFsPrLfWZregvYlLMCQ0ZwO/Wmv/tNbuBd4GTvM4JgkzJRR5LQVaG2NaGGOq4AYOve9xTBHNGGNwfdKrrLWPex1PeWGtHWWtbWatjcXV08+ttfrVFwJr7e/Ab8aYo3ybugErPQypvFgHnGqMqeE7H3RDg13LnUpeB1DWWGv3GWOGAR/jRiJPtdb+4HFYke504CpghTFmuW/bXdbamd6FJJKvG4EU3w+KtcBAj+OJeNbaxcaYN4GvcbO+vkGrZpY7WilTREREQqYuDxEREQmZEgoREREJmRIKERERCZkSChEREQmZEgoREREJmRIKERERCZkSChEREQmZEgoREREJmRIKESkWY0wrY8xeY8y9AdufNsbsMMa09yo2EfGOEgoRKRZr7WpgCnCzMaYhgDHmbuAa4CJr7TIv4xMRb2jpbREpNmPMocAaYBLwI+66DJdba1/3NDAR8YwuDiYixWat/d0YMx64FXceuUnJhEjFpi4PETlYvwBVgS+ttRO9DkZEvKWEQkSKzRjTFXgG+BI43RhzgschiYjHlFCISLEYY04G3sUNzIwD1gEPehiSiJQBSihEpMiMMa2Aj4BPgButtXuAe4FexpgungYnIp7SLA8RKRLfzI6FuBaJHtba3b7t0cD3wN/W2tM8DFFEPKSEQkREREKmLg8REREJmRIKERERCZkSChEREQmZEgoREREJmRIKERERCZkSChEREQmZEgoREREJmRIKERERCZkSChEREQnZ/wPPD9np+beqGAAAAABJRU5ErkJggg==\n", "text/plain": [ "