{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Introduction to Computational Physics - Lecture 6: Least-squares line fitting ##\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"TOC"
]
},
"source": [
"\n",
"## Table of contents\n",
" \n",
"Introduction to Computational Physics - Lecture 6: Least-squares line fitting: [>>](#Introduction-to-Computational-Physics---Lecture-6:-Least-squares-line-fitting) \n",
"-Recap of Lecture 5: [>>](#Recap-of-Lecture-5) \n",
"-Introduction to Lecture 6: [>>](#Introduction-to-Lecture-6) \n",
"-Least squares fitting: [>>](#Least-squares-fitting) \n",
"--An imaginary experiment: [>>](#An-imaginary-experiment) \n",
"--Loading the data: [>>](#Loading-the-data) \n",
"--Lecture 6 formative exercise 1: [>>](#Lecture-6-formative-exercise-1) \n",
"--Fitting the data: [>>](#Fitting-the-data) \n",
"--Visualising the result: [>>](#Visualising-the-result) \n",
"--Lecture 6 formative exercise 2: [>>](#Lecture-6-formative-exercise-2) \n",
"--Lecture 6 formative exercise 3: [>>](#Lecture-6-formative-exercise-3) \n",
"--Lecture 6 summative exercise: [>>](#Lecture-6-summative-exercise) \n",
"-Summary: [>>](#Summary) \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
" ## Recap of Lecture 5\n",
" \n",
" Last week we looked at how to input and output information:\n",
" * String formatting using f-strings or the `format` function\n",
" * Taking input from the keyboard using the `input` function and from files using the `open` function\n",
" * Writing output to files using the `write` function\n",
" \n",
" \n",
" An example of this is the following code, which reads in comma-separated (csv) data and prints it out as a formatted table, as you might want to for data from a lab experiment for example.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| Measurement | $x$ value | $x$ error | $y$ value | $y$ error |\n",
"--------------------------------------------------------------------\n",
"| 1 | 1.50 | 0.21 | 14.3 | 2.1 |\n",
"| 2 | 2.31 | 0.11 | 20.2 | 1.7 |\n",
"| 3 | 2.78 | 0.43 | 30.1 | 3.3 |\n",
"| 4 | 3.58 | 0.13 | 36.5 | 1.1 |\n",
"| 5 | 4.08 | 0.17 | 42.7 | 0.9 |\n",
"| 6 | 4.76 | 0.18 | 47.1 | 1.1 |\n",
"| 7 | 5.62 | 0.15 | 52.9 | 1.5 |\n",
"| 8 | 6.02 | 0.19 | 78.8 | 0.9 |\n",
"| 9 | 8.45 | 0.17 | 85.2 | 1.2 |\n",
"| 10 | 9.65 | 0.11 | 99.4 | 2.9 |\n"
]
}
],
"source": [
"# \n",
"\n",
"# Print out the header\n",
"print (\"| {:12s} | {:10s} | {:10s} | {:10s} | {:10s} |\".format(\"Measurement\", \"$x$ value\", \"$x$ error\", \"$y$ value\", \"$y$ error\"))\n",
"print (\"-\"*68)\n",
"\n",
"# Open the file\n",
"f = open(\"fitdata.csv\")\n",
"\n",
"# Print out the data, tabulated\n",
"imeas = 0\n",
"for line in f:\n",
" imeas = imeas +1 \n",
" data = line.strip().split(\",\") # remove \\n from end of line (strip) and split by comma to get a list\n",
" print (f\"| {imeas:<12d} | {float(data[0]):^10.2f} | {float(data[1]):^10.2f} | {float(data[2]):^10.1f} | {float(data[3]):^10.1f} |\")\n",
" \n",
"# Don't forget to close the file\n",
"f.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
" ## Introduction to Lecture 6\n",
" \n",
"We now have all the basic python we need to look at real physics use cases. We will start by looking this week at how to fit data i.e. find the parameters of a function that best describe it. You will use this frequently in labs (for fitting various data), replacing the black-box GUI you have used up to now to give you more control over the fit and the results."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Least squares fitting\n",
" \n",
"It is very common in physics, and scientific fields in general, to have a set of data points ($y_{i}$, $x_{i}$) for $i = 0 ... N$ (where $x$ is the independent variable and $y$ is the dependent one) that we need to determine the relationship between i.e. $y = f(\\beta, x)$ where $\\beta$ are the parameters of the function. I am sure you've already come across such cases at A-level and in your PHYS106 labs. Often we know the functional form of $f(\\beta, x)$ (linear, polynomial, power, exponential, etc) but we need to determine the actual $\\beta$ parameters that best describe the data. \n",
"\n",
"To do this we must minimise the (absolute) distance (often called the residual $r_i$) between the function prediction $f(\\beta, x_{i})$ and the data points $y_i$ over the full set of points. The most common method to do this is the [\"least-squared\"](https://en.wikipedia.org/wiki/Least_squares) technique (that was briefly introduced in your PHYS106 lectures - if you do labs), so-called because it minimises the sum of the squares of the residuals\n",
"\n",
"$$ S = \\sum_{i=0}^{N} r_{i}^{2} = \\sum_{i=0}^{N} \\left [y_{i} - f(\\beta, x_i) \\right ]^{2} $$\n",
"\n",
"Of course, generally we don't measure the points ($y_{i}$, $x_{i}$) exactly, but with some uncertainty on each, which also needs to be taken into account. This is done by dividing the residual by the (total) uncertainty on $y_i$, $\\sigma_{y_{i}}$, to get the residual at each point in terms of standard deviations i.e\n",
"\n",
"$$ S = \\sum_{i=0}^{N} r_{i}^{2} = \\sum_{i=0}^{N} \\left [ \\frac{ y_{i} - f(\\beta, x_i)}{\\sigma_{y_{i}}} \\right ]^{2} $$\n",
"\n",
"\n",
"For linear equations, this can be solved analytically (by setting the differential $dS/d\\beta_{i} = 0$ to find the minimum). However, in the general, non-linear case it cannot. Instead we need to use a numerical method to find the solution iteratively, starting from an initial guess of the parameters. Luckily, python contains a powerful scientific python module called [SciPy](https://www.scipy.org), which includes a [least squares](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) fitting routine as part of its [optimization](https://docs.scipy.org/doc/scipy/reference/optimize.html) package. \n",
"\n",
"We will now proceed to use the SciPy `least_squares` function, along with the NumPy `array` data structures and matplotlib `pyplot` methods we have already seen, to:\n",
"\n",
"* Load experimental data\n",
"* Fit the data with a given functional form to find the best fit parameters\n",
"* Visualise the result.\n",
"\n",
"Note: SciPy is a powerful module, that does much more than just fitting, and you will use other features of the module later in the course and in PHYS205 next year. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### An imaginary experiment\n",
"\n",
"Let's start with an imaginary experiment that performed 10 measurements of two quantities: $x$ and $y$. These are thought to be linearly related i.e. to lie on a straight line, $y = mx + c$. The $x$ and $y$ values (with their measurement errors) are \n",
"\n",
"|Measurement | $x$ value | Error in $x$ | $y$ value | Error in $y$ |\n",
"|------------|-----------|--------------|-----------|--------------|\n",
"| 1 | 1.50 | 0.21 | 14.3 | 2.1 |\n",
"| 2 | 2.31 | 0.11 | 20.2 | 1.7 |\n",
"| 3 | 2.78 | 0.43 | 30.1 | 3.3 |\n",
"| 4 | 3.58 | 0.13 | 36.5 | 1.1 |\n",
"| 5 | 4.08 | 0.17 | 42.7 | 0.9 |\n",
"| 6 | 4.76 | 0.18 | 47.1 | 1.1 |\n",
"| 7 | 5.62 | 0.15 | 52.9 | 1.5 |\n",
"| 8 | 7.02 | 0.19 | 68.8 | 0.9 |\n",
"| 9 | 8.45 | 0.17 | 85.2 | 1.2 |\n",
"| 10 | 9.65 | 0.11 | 99.4 | 2.9 |\n",
"\n",
"**Table 1** *A number of data points which are expected to lie on a straight line.*\n",
"\n",
"and are stored in the fitdata.csv file we printed above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Loading the data\n",
"\n",
"The first thing we need to do is to load the data into NumPy arrays. We could do this manually by creating a NumPy array and copying the data from the file by hand e.g.\n",
"\n",
"```python\n",
"nPoints = 10\n",
"xData = np.zeros(nPoints)\n",
"xData[0] = 1.50\n",
"xData[1] = 2.31\n",
"```\n",
"etc, but this is error prone. Better is to read the data directly from the csv file, using the [np.loadtxt](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.loadtxt.html) function (which we covered in week 3) along with the optional `unpack` argument to transpose the data so that it can be read into separate arrays:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.5 2.31 2.78 3.58 4.08 4.76 5.62 6.02 8.45 9.65]\n",
"[0.21 0.11 0.43 0.13 0.17 0.18 0.15 0.19 0.17 0.11]\n",
"[14.3 20.2 30.1 36.5 42.7 47.1 52.9 78.8 85.2 99.4]\n",
"[2.1 1.7 3.3 1.1 0.9 1.1 1.5 0.9 1.2 2.9]\n"
]
}
],
"source": [
"# \n",
"\n",
"import numpy as np\n",
"xdata, xerror, ydata, yerror = np.loadtxt(\"fitdata.csv\", delimiter = ',', unpack = True)\n",
"print(xdata)\n",
"print(xerror)\n",
"print(ydata) \n",
"print(yerror)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"Once we have loaded the data, the first thing we should do is check it and the best way to do this is to plot it. This allows us to:\n",
"\n",
"* Check that it is correct\n",
"* Confirm that it lies on the purported straight line\n",
"\n",
"### Lecture 6 formative exercise 1\n",
"\n",
"Plot the data, including its errors, and check it looks sensible. Use `matplotlib.pyplot` (which we introduced in weeks 3 and 4) with the [errorbar](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.errorbar.html) function. Does the data look like it lies on a straight line? If not, double check the values in the file compared to Table 1 above (which we assume is the definitive source) and fix any errors. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Fitting the data\n",
"\n",
"The data should now look compatible with a straight line, as expected, so we are ready to fit it. As mentioned above, we will do this using SciPy's [optimise.least_squares](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html). Looking at the documentation the function definition is \n",
"\n",
"```python\n",
"scipy.optimize.least_squares(fun, x0, jac='2-point', bounds=(- inf, inf), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options={}, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs={})[source]\n",
"```\n",
"\n",
"This looks scary, but for now there are only a couple of arguments that are important to us: the first two, which have no default value, and the `args` tuple. \n",
"\n",
"The first,`fun`, is a function that calculates the residuals $r_{i}$ (with error) defined above, which will be minimised to find the best fit. Note that this function (which you will sometimes see referred to as a _cost_ or _merit_ function) defines the residuals, not the sum of the squares of them, since the `least_squares` method takes care of the latter automatically. The function must have the form ```fun(params, arg1, arg2, ...)``` where `params` are the parameters to be minimised and `args` are the \"data\", which will be passed in indirectly via the `args` parameter of `least_squares`. In our case these are `args = (xdata, ydata, xerror, yerror)` so the minimisation function will have the form ```fun(params, xdata, ydata, xerror, yerror)```. It should return an array of residuals, one for each point.\n",
"\n",
"The second, `x0`, is simply the initial guess of the parameters that we will fit. It is a tuple/array of the same length as the number of parameters we will fit, 2 in our case: $m$ and $c$. \n",
"\n",
"Let's first define the minimisation function to calculate the residuals:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# \n",
"\n",
"def straight_line(params, xdata):\n",
" \"\"\"\n",
" Function for a straight line\n",
" - params array of function parameters\n",
" - xdata array of x data points\n",
"\n",
" \"\"\"\n",
" \n",
" f = params[0] + params[1]*xdata\n",
" return f\n",
"\n",
"def straight_line_diff(params, xdata):\n",
" \"\"\"\n",
" Differential of function for a straight line\n",
" - params array of function parameters\n",
" - xdata array of x data points\n",
"\n",
" \"\"\"\n",
"\n",
" df = params[1]\n",
" return df\n",
"\n",
"def minimise(params, xdata, ydata, xerr, yerr):\n",
" \"\"\"\n",
" Calculation to minimise to find the best fit using chi2: difference between each y point and its \n",
" prediction by the function, divided by the sum in quadrature of the error on y, both from the \n",
" y error and from the related error in x.\n",
" - params array of function parameters\n",
" - xdata array of x data points\n",
" - ydata array of y data points\n",
" - xerr array of x data errors\n",
" - yerr array of y data errors\n",
" \"\"\"\n",
"\n",
" residuals = (ydata - straight_line(params, xdata)) / (np.sqrt(yerr**2 + straight_line_diff(params, xdata)**2 * xerr**2))\n",
"\n",
" return residuals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"The minimisation function simply calculates the residuals defined in the second equation above, using NumPy's array-at-a-time functions we covered previously. \n",
"\n",
"In order to do this, it depends on the mathematical function we are fitting $f(\\beta, x) = mx + c$, which we have coded in the python function `straight_line`. The `straight_line` function takes the parameters (${\\rm params} = \\beta$) as a array, in this case a two-component array whose first component `params[0]` is the intercept $c$ and whose second component `params[1]` is the gradient $m$, along with the independent variable $x$\n",
"\n",
"It also depends on the derivative of the straight line function, `straight_line_diff`. You should be able to work out why this is the case from the error analysis you do in the PHYS106 lectures or by thinking about the effect of an error on the $x$ variable on the result of the function.\n",
"\n",
"Note that we are following good practice and writing _doc strings_ for our functions. In this case we have not only explained what the function does but also what the inputs are. Since you will be using this code often in your labs, this will help you remember what everything does.\n",
"\n",
"Now, that we have the function we need an initial guess of the intercept and gradient, which from the plot above look like roughly 0 and 10, respectively."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# \n",
"init_params = [0.0, 10.0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now we have the functions and initial parameters defined we can call the `least_squares` method, passing in the data as `args`:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" message: `ftol` termination condition is satisfied.\n",
" success: True\n",
" status: 2\n",
" fun: [ 1.567e-01 -9.440e-01 5.743e-01 7.912e-01 1.247e+00\n",
" -5.584e-02 -1.457e+00 -7.314e-01 8.735e-02 6.729e-01]\n",
" x: [-1.535e+00 1.024e+01]\n",
" cost: 3.2745972129027687\n",
" jac: [[-3.327e-01 -5.068e-01]\n",
" [-4.903e-01 -1.105e+00]\n",
" ...\n",
" [-4.729e-01 -4.002e+00]\n",
" [-3.214e-01 -3.110e+00]]\n",
" grad: [ 2.254e-07 -6.923e-07]\n",
" optimality: 6.923129177509927e-07\n",
" active_mask: [ 0.000e+00 0.000e+00]\n",
" nfev: 5\n",
" njev: 5\n"
]
}
],
"source": [
"# \n",
"from scipy.optimize import least_squares\n",
"result = least_squares(minimise, init_params, args=(xdata, ydata, xerror, yerror))\n",
"print(result)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"OK, we have run the fit, now what? We have a results object, which contains several pieces of information documented [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.OptimizeResult.html#scipy.optimize.OptimizeResult)\n",
"\n",
"The first thing is to check if the fit was successful, using the `success` option of the result, and the status at termination, using the `status`. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fit succeeded\n"
]
}
],
"source": [
"# \n",
"\n",
"if not result.success or result.status < 1:\n",
" print (\"ERROR: Fit failed with message {}\".format(result.message))\n",
" print (\"Please check the data and initial parameter estimates\")\n",
"else:\n",
" print (\"Fit succeeded\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"OK, the fit succeeded, so we can then get the fitted parameters via `result.x`, which is an array containing the resulting intercept (index 0) and slope (index 1). \n",
"\n",
"The first thing we should do with these is to check if they are a good fit to the data by calculating the [chi-squared](https://en.wikipedia.org/wiki/Chi-squared_test) ($\\chi^2$) test statistic, that was introduced in your PHYS106 lectures and in this week's lecture. To do this, we first evaluate our residuals function `minimise` for the fitted parameters. We can do this manually by passing the `result.x` values into `minimise`\n",
"\n",
"```python\n",
"final_params = result.x\n",
"residual_array = minimise(final_params, xdata, ydata, xerror, yerror)\n",
"```\n",
"\n",
"but the `least_squared` result already contains this as `result.fun`. We then square this to get the chi-squared per data point and sum to get the total chi-squared. \n",
"\n",
"More useful, is the reduced chi-squared ($\\chi^2$/NDF) by dividing by the number of degrees of freedom (i.e. number of points - number of fit parameters). The closer the $\\chi^2$/NDF is to one the better the fit. Generally, for small numbers of degrees of freedom, values in the range $0.25 < \\chi^2$/NDF $ < 4$ represent good agreement. \n",
"\n",
"A somewhat more robust fit statistic is the chi-squared probability, which is the probability that we get the resulting chi-squared given the number of degrees of freedom and should be greater than 5% (ie. within 95% or $\\approx 2\\sigma$). The `scipy.stats.chi2` [module](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html) has a function (called the \"survival function\" or \"sf\") to calculate this."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"=== Fit quality ===\n",
"chisq per point = \n",
" [0.025 0.891 0.33 0.626 1.556 0.003 2.123 0.535 0.008 0.453]\n",
"chisq = 6.5492, ndf = 8, chisq/NDF = 0.81865, chisq prob = 0.58596\n",
"\n"
]
}
],
"source": [
"# \n",
"\n",
"from scipy.stats import chi2 as stats_chi2\n",
"\n",
"# Get fitted parameters \n",
"final_params = result.x\n",
"c = final_params[0]\n",
"m = final_params[1]\n",
"nparams = len(final_params)\n",
"\n",
"# Calculate chi2\n",
"chi2_array = result.fun ** 2\n",
"chi2 = sum(chi2_array)\n",
"npoints = len(xdata)\n",
"reduced_chi2 = chi2 / (npoints - nparams) \n",
"chi2_prob = stats_chi2.sf(chi2, (npoints - nparams))\n",
"\n",
"# Print chi2\n",
"np.set_printoptions(precision = 3)\n",
"print(\"\\n=== Fit quality ===\")\n",
"print(\"chisq per point = \\n\",chi2_array)\n",
"print(\"chisq = {:7.5g}, ndf = {}, chisq/NDF = {:7.5g}, chisq prob = {:7.5g}\\n\".format(chi2, npoints-nparams, reduced_chi2, chi2_prob))\n",
"\n",
"if reduced_chi2 < 0.25 or reduced_chi2 > 4:\n",
" print(\"WARNING: chi2/ndf suspiciously small or large. Please check the data and initial parameter estimates\")\n",
"\n",
"if chi2_prob < 0.05:\n",
" print(\"WARNING: chi2 probability for given degrees of freedom less than 0.05 . Please check the data and initial parameter estimates\") "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"Now that we have checked that it is a good fit, we can print out the parameters and their errors. We have the parameters already but we need to calculate the errors from the [jacobian matrix](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant), which is the matrix of the first order partial derivatives. We square this and take the inverse to get the [covariance matrix](https://en.wikipedia.org/wiki/Covariance_matrix), whose diagonal elements are the variance (i.e. error squared) of the parameters. In doing so, we check that the squared jacobian is not singular (i.e. has a finite determinant). You will learn the details of this in PHYS108 in semester 2 so for now you will just have to trust me. Note: NumPy has all the functions we need to deal with the matrix calculations."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"c = -1.5352 +- 1.7438\n",
"m = 10.243 +- 0.3193\n"
]
}
],
"source": [
"# \n",
"\n",
"# Calculate errors\n",
"jacobian = result.jac\n",
"jacobian2 = np.dot(jacobian.T, jacobian)\n",
"determinant = np.linalg.det(jacobian2)\n",
" \n",
"if determinant < 1E-32:\n",
" print(f\"Matrix singular (determinant = {determinant}, error calculation failed.\")\n",
" param_errors = np.zeros(nparams)\n",
"else:\n",
" covariance = np.linalg.inv(jacobian2)\n",
" param_errors = np.sqrt(covariance.diagonal())\n",
"\n",
"print (\"c = {:7.5g} +- {:7.5g}\".format(final_params[0], param_errors[0]))\n",
"print (\"m = {:7.5g} +- {:7.5g}\".format(final_params[1], param_errors[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Visualising the result\n",
"\n",
"Finally, we plot the data with error bars (using the routine `plt.errorbar` from `matplotlib.pyplot`) and draw the fitted line, which we obtain by evaluating our `straight_line` function with the `final_params` on the plot (using `plt.plot`). We also use the [plt.text](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.text.html) function to print the fitted parameter values on the plot so everything is in one place."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAIqCAYAAAAq8PwLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACOn0lEQVR4nOzdeZyN5f/H8deZlVmMGIaxb5UII1IoWmilvV+hqG9aqEh9JYpRslUoorQho/KtFBXxrUxJStb2sqZkjRkzmLPdvz/u5nydMwdnOGfuM+e8n4/HPHJfZ7uuOe/OfOae67pum2EYBiIiIiIiUSLG6g6IiIiIiJQlFcAiIiIiElVUAIuIiIhIVFEBLCIiIiJRRQWwiIiIiEQVFcAiIiIiElVUAIuIiIhIVFEBLCIiIiJRRQWwiIiIiEQVFcAiImVsxowZ2Gw2v18PPfQQW7ZswWazMWPGDM9jli9fTnZ2Nvv377es3ycjOzsbm83m1TZ16lSvMYqIlJU4qzsgIhKtXnvtNU4//XSvtszMTDIyMvjqq69o1KiRp3358uWMHDmSPn36ULly5TLuaWhMnTqV9PR0+vTpY3VXRCTKqAAWEbFI8+bNadOmjd/bzjnnnKC/3sGDB0lKSgr684qIlDeaAiEiEmZ8p0BkZ2fz73//G4AGDRp4pkssXbr0qM/Rp08fUlJS+O677+jatSupqalcdNFFANjtdkaNGsXpp59OYmIi1apV47bbbmP37t1ez/Hpp5/SuXNnqlatSsWKFalbty7XXXcdBw8eBGDp0qV+++FvCoev+vXr88MPP5Cbm+sZT/369Uv1fRIROVE6AywiYhGXy4XT6fRqi4sr+bF8xx138PfffzN58mTeffddatasCcAZZ5xxzOe32+10796du+66iyFDhuB0OnG73Vx11VV88cUXDB48mPbt27N161ZGjBhB586d+fbbb6lYsSJbtmzhiiuu4LzzzuPVV1+lcuXK/PnnnyxatAi73X7SZ5LnzZvH9ddfT1paGlOnTgUgMTHxpJ5TRCRQKoBFRCzib5qDw+Eo0Va7dm3q1q0LQFZWVsBnSh0OB8OHD+e2227ztL355pssWrSId955h2uvvdbT3rJlS9q2bcuMGTO45557WLVqFYcPH+app56iZcuWnvv16NEj0OEdU1ZWFhUrVqRSpUohme4hInIsKoBFRCwya9YsmjZt6tXm7wzwybjuuuu8jj/44AMqV65Mt27dvM4+t2rViho1arB06VLuueceWrVqRUJCAnfeeSf9+vXjvPPOo2HDhkHtm4iIVVQAi4hYpGnTpkddBBcMSUlJVKpUyatt586d7N+/n4SEBL+P2bNnDwCNGjXiv//9L+PHj6d///4UFhbSsGFD7r//fgYMGBCyPouIlAUVwCIiEcp3312A9PR0qlatyqJFi/w+JjU11fPv8847j/POOw+Xy8W3337L5MmTGThwIBkZGdx0001UqFABgKKiIq/nKC6iRUTClXaBEBEpB4oXiB06dOiknufKK69k7969uFwu2rRpU+LrtNNOK/GY2NhY2rVrx/PPPw/A6tWrATxzkdevX+91//nz5wfUl8TExJMej4jIidAZYBGRcuDMM88E4Nlnn6V3797Ex8dz2mmneZ2xDcRNN91ETk4Ol19+OQMGDODss88mPj6eP/74g88++4yrrrqKa665hhdeeIFPP/2UK664grp163L48GFeffVVAC6++GIAatSowcUXX8yYMWM45ZRTqFevHp988gnvvvtuwGN68803eeutt2jYsCEVKlTwjFNEJJRUAIuIlAOdO3fmkUceYebMmbz00ku43W4+++wzOnfuXKrniY2NZf78+Tz77LO8/vrrjBkzhri4OGrXrk2nTp08BWirVq1YvHgxI0aMYMeOHaSkpNC8eXPmz59P165dPc/3+uuvc9999/Hwww/jcrno1q0bb7zxRkBzm0eOHMlff/1F3759OXDgAPXq1WPLli2lGo+IyImwGYZhWN0JEREREZGyojnAIiIiIhJVVACLiIiISFRRASwiIiIiUUUFsIiIiIhEFRXAIiIiIhJVVACLiIiISFTRPsABcrvdbN++ndTUVL+XFxURERERaxmGwYEDB8jMzCQm5ujneVUAB2j79u3UqVPH6m6IiIiIyHFs27aN2rVrH/V2FcABKr7c6LZt26hUqZLFvQlvRUVFTJgwgUGDBpGYmGh1dyQMKBPij3IhvpQJ8VXaTOTn51OnTp3jXiZeBXCAiqc9VKpUSQXwcRQVFVGhQgUqVaqkDzABlAnxT7kQX8qE+DrRTBxvuqouhRyg/Px80tLSyMvLUwF8HG63mz179pCenn7M+TcSPZQJ8Ue5EF/KhPgqbSYCrdeULgk6m81GWlqaFguKhzIh/igX4kuZEF+hyoQKYAk6u93O2LFjsdvtVndFwoQyIf4oF+JLmRBfocqECmARERERiSoqgEVEREQkqqgAFhEREZGool0gAqRdIAJnGAZ2u52EhAQtZBBAmRD/lAvxpUyIr9JmItB6TfsAlwGHw4HL5bK6G2XG7Xazb98+TjnllKjfxiY2Npb4+Hiru2E5wzDIy8sjPT1dP9TEQ7kQX8qE+ApVJlQAh1B+fj579uyhqKjI6q6UqeKw5uXl6QMMSExMJD09Par/cuBwOJg2bRpDhgzR5vbioVyIL2VCfIUqEyqAQyQ/P58///yTlJQU0tPTiY+Pj5piUBuZmwzDwOFwkJeXx59//gkQ1UWwiIhIuFABHCJ79uwhJSWF2rVrR03hW8ztdhMXF0eFChWiugAGqFixIqmpqfzxxx/s2bNHBbCIiEgYiO7qJEQcDgdFRUVRfTWbaB23P8VXsSkqKsLhcFjdHcskJCRY3QUJQ8qF+FImxFcoMqFdIAJUml0gDh8+zObNm6lfvz4VK1Ysox5KODt06BBbtmyhQYMGVKhQweruiIiIRKRA6zWdAQ6haD0LahgGhw8fRr9b/U+0ZqGY2+1mw4YNuN1uq7siYUS5EF/KhPgKVSZUAEvQGYbB33//rQJYPBwOBzk5OVE9BURKUi7ElzIhvkKVCRXAIiIiIhJVVACLiIiISFRRASwhERcXvB32OnfujM1mY+nSpUF7TilbNpuNatWqRf1caPGmXIgvZUJ8hSoTKoAl6GJiYqhevfpx9wCuX78+NpvtmF+TJk066uOXLl1Kdna2CuNyICEhgX79+ml7I/GiXIgvZaJ8KSws9Py8LiwsDMlrhCoTuhCGBJ1hGBw8eJCkpKSAfmNr0qQJ1atX93tbrVq1qFu3LqeddhpJSUlety1dupSRI0cC5lliCV8ul4t169bRsmVLYmNjre6OhAnlQnwpE+IrVJlQASxBZxgGeXl5VKxYMaACeOjQofTp0+eot99www1B7J1Ywel0smDBApo1a6YfauKhXIgvZUJ8hSoTmgIhIiIiIlFFBbCEPX+L4Gw2m2f6w8iRI73mDR/rbLKIiIiIpkBISCQmJob0+Tt06MDvv//Otm3bqFOnDnXr1vXcduqpp4b0taX0bDYbjRo10spu8aJciC9lQnyFKhMqgCXoYmJiqFq1akhfY9myZWRnZzNy5Ehuv/12srOzQ/p6cnISEhLo1auX1d2QMKNciC9lQnyFKhNhPwXi888/p1u3bmRmZmKz2Xjvvfe8bjcMg+zsbDIzM6lYsSKdO3fmhx9+8LpPUVER9913H+np6SQnJ9O9e3f++OOPMhyFD8OAwsLy+RXA5Y0Nw+DAgQMBXwr5tttu87sFmnZ2iBxOp5OlS5fidDqt7oqEEeVCfCkT4itUmQj7M8CFhYW0bNmS2267jeuuu67E7ePHj2fChAnMmDGDU089lVGjRtGlSxd++eUXUlNTARg4cCALFizgzTffpGrVqjz44INceeWVrFq1yppVpgcPQkpK2b9uMBQUQHLyMe9SXAAnJyef1DZoZ5555gl3U8KLy+UiNzeXc889N6gXSZHyTbkQX8qE+ApVJsI+XZdddhmXXXaZ39sMw2DSpEkMGzaMa6+9FoCZM2eSkZHBnDlzuOuuu8jLy+OVV17h9ddf5+KLLwZg9uzZ1KlTh//+979ccsklfp+7qKiIoqIiz3F+fn6J9piYGOLj43E4HLjdbs99i39LcbvdXu3FZzbdbnf4n3o/CrfbDUcZU7Ejz/we2V58f9/7DBkyhNtuu61Eu7/XLn6+mJgYz30NwyjR7vs8MTExfvvi2/cTbffX9yPb3W43hmFgt9upUKECbrcbh8Phdd+EhARcLpfXb7nF7U6nE5fL5TUef9mLjY0lLi4Ou93u1Z+4uDhiY2NLtMfHxxMTE+OV9eJ2m82G3W73ak9ISMAwDK++gznn+1hjKn7+oqKiiBlTJL5PZT2mYkeOq7yPKRLfp7Ic05GvEyljOrI90sZ0pKKiIuLi4kI2puL/Hm9Mvt+Down7AvhYNm/ezI4dO+jataunLTExkU6dOrF8+XLuuusuVq1ahcPh8LpPZmYmzZs3Z/ny5UctgMeMGePZZeBIEyZMoEKFCgBkZWXRvXt3Fi5cyJo1azz3ueCCC6hWrRp5eXns37/f056WlkZycjJ7Dh7E9dtvnvZTTjmFChUqsGPHDq83Pz09ndjYWHbu3OnVh4yMDFwuF3v27PG02Ww2atSoweHDh9m3b5+nPTY2lurVq1NYWOgp4sEMStWqVTlw4AAFBQWe9ooVK1K5cmX279/PoUOHPO0pKSmkpqay79Ahig4cKDmmPXv8/nli586dXmOqVq0asbGx7Nixw/MhkJeXh2EYuFwudu/e7TWmmjVrev4H+vvvv9mxYwdxcXFUr17d8z9oQUEBO3bsIDExkapVq1JQUMCBI/qYlJRE5cqVyc/P5+DBg5721NRUc0z79nn9D3O0MVWpUoUKFSocc0xHqlGjhmdMTqeTvLw8lixZwoABA9i0aRM5OTlez9GvXz/WrVvHggULPO2NGjWiV69eLFu2jNzcXE/70bLXqVMnOnfuzNy5c9m4caOnvVu3brRu3ZqXX37Z63vcs2dPGjduzIQJE7w+TO655x7S0tIYO3as15iGDBlCXl4e06ZN87QlJCTwyCOPBDSmiRMnRtyYIPLep7IaU506dQAzF5Eypkh8n8p6TMUiaUyR+D41atSIa665xnP8zDPPkJCQEPQxPf/888D/PieON6YjP0+OxWYEOlEzDNhsNubNm8fVV18NwPLly+nQoQN//vknmZmZnvvdeeedbN26lY8//pg5c+Zw2223lfiNoGvXrjRo0IAXX3zR72v5OwNcp04ddu3aRaVKlYCj/+bmcDjYtm0b9erV8xTLxf0P9pnFIx3r7Gdp20/mrGjxFIji79PR+t6wYUO2bt3KK6+8ctQzwDExMXTu3Jnc3Fw++eQTz7zgmJgYRowYweOPP87w4cMZMWJESMd0vHZ/fT+y/fDhw2zZsoXatWtTqVKlqDtrcOjQIZYsWUKXLl1ISEiIiDFF4vtU1mNyuVx88MEHdOnSxXNGuLyPKRLfp7Ick8PhYMmSJVx55ZXYbLaIGNOR7eX5fcrPzy8xJofDQUZGBgC///47ycnJfseUnJx8wmMqKCjw/PyIj48/7pj27t1L9erVycvL81uHeF73qLeUI76n4Q3DOO7c0+PdJzEx0e9WXv7aj/xTXvFzgxnsmJiSkx38tR2r/Wj99NdeXKSdbHtp++jbXrlyZb/3O/J1j3xs8fGx3hPf72fxv202m1d7qMZ0vPZjvU/FYyy+lnlMTIzffMXGxvqdlx4XF+d37pNv9ood7ZrpR2s/2rZ1/tptNpvf9mONKSUlxetMAZT/MUXi+1TWY4qJiSmRi2LldUwQee8TlN2YEhMTvTIRCWM6Unl+n9LS0vw+d7EjtyP1dWRhW9ox+fv5AaUbkz/ldSoqYP6JGSjxp+ddu3Z5fiOpUaMGdrvda1qA730kuNxuN/v37y9xpjTYKlasCOA1VUPCk8PhYP78+SXONkh0Uy7ElzIhvkKViXJdADdo0IAaNWqwZMkST5vdbic3N5f27dsDcNZZZxEfH+91n7/++ovvv//ecx8JviPn2oZKw4YNAXMqjLbMCW9ut5s1a9aE/JciKV+UC/GlTISvgoKCEl9HrlHauXOn3/scuc7oRIQqE2E/BaKgoIANGzZ4jjdv3szatWupUqUKdevWZeDAgYwePZomTZrQpEkTRo8eTVJSEj169ADMU/b/+te/ePDBB6latSpVqlThoYce4swzz/TsCiHlU9euXTnllFNYtmwZdevWpWHDhsTFxXHppZcyZMgQq7snIiISMZKPswVqcnLyce8TTsK+AP7222+54IILPMeDBg0CoHfv3syYMYPBgwdz6NAh+vXrx759+2jXrh2LFy/27AEM5srBuLg4brzxRg4dOsRFF13EjBkzrNkDWIKmUqVKLF68mOHDh/P111/z1Vdf4Xa7qV+/vtVdExERkTBWrnaBsFJ+fj5paWnHXVUIcPjwYTZv3kyDBg28doGIFoZhUFBQQEpKiq7n/o9oz4TT6WTZsmV07NhRm9uLh3IhvpSJ8qWwsJCUfy7sVVBQEJIzwKXNRKD1mtIlQWez2bzOwIvExcV5trATKaZciC9lQnyFKhPlehGchCe3283evXu1iEE87HY7s2fPLrFvo0Q35UJ8KRPiK1SZUAEsIRHopQglOhiGwcaNG495qWuJPsqF+FImxFeoMqECWEREREROWDurO3ACVACLiIiISOnl5zMLWAHEzptndW9KRQWwBJ3NZiMtLU07QIhHXFwc3bp106pu8aJciC9lohz56isqtm/PLYALsG3bFpKXCVUmVABL0NlsNpKTk1UAi0dsbCytW7fW3tviRbkQX8pEOeBywRNPwHnnEbNlC5uB8wDn/feH5OVClQkVwBJ0brebXbt2aRcI8bDb7UydOlUru8WLciG+lIkw9/vvcMEFMHw4uFw4b7yRVsBXIXzJUGVCBbCEhNPptLoLEkYMw2D37t1a2S1elAvxpUyEsblzoWVL+OILSE2F11+n6NVXyQ/xy4YqE5pkIyIiIiL+FRTA/ffDa6+Zx+3aQU4ONGpEMpTbX1Z0BlhERERESvr2W2jd2ix+bTYYNsw8A9yokdU9O2k6AyxBZ7PZqFKlihbBiUd8fDw9e/YkPj7e6q5IGFEuxJcyESbcbnjqKXj0UXA6oXZtmD0bOnUq866EKhMqgCXobDYbFSpUsLobEkZiYmJo3Lix1d2QMKNciC9lIgz8+Sfceit8+ql5fP31MH06nHKKJd0JVSY0BUKCzu1289dff2kXCPEoKipizJgxukS2eFEuxJcyYbH33oMWLcziNykJXnnFXPxmUfELocuECuAIUlhYiM1mw2azUVhYaGlfyuukeAkdbWsk/igX4kuZsMDBg3D33XDNNfD333DWWbBmDdx+uzn312KhyIQKYBEREZFotXatWfC++KJ5PHgwLF8Op55qabdCTQWwWKZ+/fqeM9Y2m42YmBgqVapEnTp16NKlC48++ig//vhjUF9z//79ZGdnM2nSpKA+r4iISLnidsPEiea2Zj//DDVrwpIlMG4cJCRY3buQsxn6W3VA8vPzSUtLIy8vj0qVKh3zvocPH2bz5s00aNCgTBeDFRYWkpKSAkBBQQHJycll9tpHMgwDp9NJXFzcMXeCqF+/Plu3bqVJkyZUr14dML93e/bsYevWrZ77XXfddbz44otUrVr1pPu2ZcsWGjRoQL169diyZctJP1+grMpEuHC73ezZs4f09HRiYvR7t5iUC/GlTJSRHTugTx/4+GPz+Kqr4OWXIT3d0m75U9pMBFqvKV0SEqW5ZvfQoUNZtmwZy5Yt49tvv2XLli3s3r2bSZMmkZ6ezjvvvEPHjh3Jy8sLYY8llGw2G2lpadoaT7woF+JLmSgDH35oLnT7+GOoWBGmTYN588Ky+IXQZUIFsASdYRjs2LHjpBbCpaenM2DAAL799ltq1qzJzz//zMCBA4PXSSlTdrudsWPHanGLeFEuxJcyEUKHD5tXdLvySti92yyCv/3WXPwWxr9whCoTKoAlrNWrV4+pU6cCMHv2bLZt2+a5bdOmTYwbN47OnTtTp04dEhMTqVatGpdeeikffvhhiefq06cPDRo0AGDr1q1e84+P/M3y0KFDvPHGG9x0002cdtpppKSkkJKSQqtWrRg1apTlO2yIiIiUyvffQ9u2MHmyeTxwIHz9NZxxhqXdspIKYAl73bt3JzMzE6fTyeLFiz3to0ePZsiQIaxatYqkpCRatGhBfHw8H3/8MVdeeSXjxo3zep5TTz2VNm3aAJCYmEiHDh28voqtWrWKHj168M4773Dw4EGaNm1KZmYmP/zwA4899hjnn38+hw4dKpvBi4iInCjDgClToE0bswiuXh0WLjQXv0XhepQjqQCWsBcTE8O5554LwMqVKz3t1113HStWrCA/P59ffvmFlStXsn37dj7//HNq1qzJsGHD2Lhxo+f+Q4cO5T//+Q8ANWrU8Mw7Lv4qVqdOHebOncu+ffvYtm0bK1eu5Ndff2Xbtm1cf/31rF69mvHjx5fR6EVERE7A7t3QvTvcdx8UFcFll8H69XDppVb3LCyoAC4nCgsLA/o60fsHk81mo0aNGkGdsF6nTh0Adu3a5Wm77LLLaNeuXYnXOe+883jiiSdwuVy89dZbpX6tevXqccMNN3h21ChWo0YNZs2aRUJCAjk5OScwiuiVkJDAkCFDSIiCrXUkcMqF+FImgmTxYnOO7wcfmFuaPfusufgtI8PqnpVaqDIRF9Rnk5DxLcaOJyPAkIdqFzyXy0VcXPDiVbyl24EDB7zad+/ezZw5c/j666/ZtWsXhw8fBvDsGLFu3boTej23282CBQtYvHgxmzZtoqCgwPO9stls/Pbbbxw8eJCkpKQTHVJUMQyDvLw80tPTtbpbPJQL8aVMnKSiIhg6FCZMMI/POAPeeMMshsupUGVCBbAEnWEY7N69O6hngQsKCgC89vRbvHgxN9544zG3R/v7779L/Vr79+/n8ssv56uvvjrm/fbt26cCOEAOh4Np06YxZMgQEhMTre6OhAnlQnwpEyfh55/h5pvNK7sB9OsHTz9tbnVWjoUqE5oCUU4UFBQc92vnzp2e++/cuTOgx5QXv//+O4Dnghn79+/npptuIi8vj1tvvZUVK1awb98+XC4XhmGwZMkSwPwfp7QGDRrEV199xWmnncY777zDn3/+SVFREYZhYBgGtWrVOuHnFhERCSrDgOnToXVrs/itWhXefx+ef77cF7+hpDPA5URpr+qWnJxs2ZXggs3tdnvOxp599tkALFy4kH379nHuuecyY8aMEmeaj9wurTScTidz584F4P333+e0004rcfuOHTtO6LlFRESCau9e6NvXvJAFwMUXw8yZkJlpbb/KAZ0BlpAI5jyd9957jx07dhAfH0/Xrl0BPJcxPvfcc/2+1tHm/h6vX7t376awsJAqVaqUKH4Bvv/+e1wuVylHIIAWtYhfyoX4UiYC9Nln0LKlWfzGx5vTHT7+OCKL31BkQmeAJehiYmKoWbNmUJ5r69at3HvvvQDceuutnukHFf/5s86R0z6K7d27l1deecXv8xU/7mj7+Bbfnp+fz6FDhzzHxbT92YlJTEzkkUcesbobEmaUC/GlTATA4YDhw2HcOHP6w6mnmgvdWre2umchEapM6AywBJ1hGBw+fPikdpjYs2cPzz33HG3atOGvv/7ijDPOYELxqlbMrc4A5s6dy3//+19P+19//cV1112H0+n0+7zVqlUjNTWVXbt28dNPP5W4vXLlyjRr1gyn08kDDzzgufSiy+Vi3LhxvPXWWzo7cQLcbjcbNmzA7XZb3RUJI8qF+FImjmPDBujQAcaONYvfO+6A1asjtviF0GVCBbAEnWEY/P333wEXwKNHj6Zjx4507NiRtm3b0qBBA6pVq8aAAQPYs2cPN9xwA1988YXXDhBnnXUW119/PQ6Hgy5dutCkSROysrKoW7cuq1evZuzYsX5fy2azccMNNwDQunVr2rZtS+fOnencubPnPmPGjMFms/Hiiy9Ss2ZN2rZtS40aNRgyZAjDhg0L2tntaOJwOMjJydHCQfGiXIgvZeIoDANmzIBWrWDlSjjlFHj7bXjpJYiQ9T5HE6pMaAqEWO63337jt99+A8z9jitXrszFF19Mu3bt6NmzJ02bNvX7uJycHJo2bcrrr7/O1q1bqVq1Ktdffz3Z2dn89ddfR329Z599ltTUVN5//33WrVtX4n+qbt26sXDhQh5//HHWrFnDL7/8QrNmzZg0aRI9e/Zk1qxZwRu8iIjIsezfD3ffDcUXdurcGWbNgn8uECUnRgWwWKZ4IduJSkhI4PHHH+fxxx8vcdtpp5121DPQKSkpTJo0iUmTJh31uS+55BIuueQSv7edbL9FREQCsmwZ9OwJv/8OsbHwxBMweLD5bzkpKoAlJIJ5FTgp/2w2G9WqVdOVncSLciG+lIl/OJ1msTtqFLjd0LAhzJkD7dpZ3bMyF6pM2IxQXQs3wuTn55OWlkZeXp7XXFR/Dh8+zObNm2nQoAEVKlQoox5CYWGh55LJBQUFEbMPcCSwKhMiIlLObN5snvUtvhpp794weTKkplrbr3Ii0HpNi+AiSHJysudqZVYWv4ZhUFhYeFK7QEhkcblcrF69WnsoixflQnxFfSbmzDEXun31FVSqZG5vNmNGVBe/ocqECmAJOsMwyMvLUwEsHk6nkwULFhx1ezqJTsqF+IraTOTnw623mmd+8/OhfXtYtw5uusnqnlkuVJlQASwiIiJila+/hqwseP11iImB7GzIzYX69a3uWUTTSiURERGRsuZymRe0GDHC/He9epCTY17oQkJOBbCERGJiotVdkDBis9lo1KiRVnaLF+VCfEVNJn7/HW65BT7/3Dy+6SaYNg0qV7a0W+EoVJnQLhABKg+7QEj4UiZERASA//wH7rzTvMBFSgo8/7xZDEd60V9GtAtEGIjW3y0Mw+DAgQNRO35/ov174XQ6Wbp0afQtbJFjUi7EV0RnoqAA/vUvuPFGs/g9+2xYu9Zc/Kbi96hClQkVwCEQE2N+W6N1GxcVwCUVZ6E4G9HG5XKRm5sbtf9PiH/KhfiK2EysWgWtW8Orr5rF7tCh5lXeGjWyumdhL1SZiM6fxiEWHx9PbGwshw4dsrorEiYOHTpEbGws8fHxVndFRETKitsN48fDuefCb79B7drw6afw5JOgnweWUgEcAjabjaSkJPLy8iLvt1gpNZfLRV5eHklJSZG/sENEREx//gldu8LDD4PDAdddZ+7t27mz1T0TtAtEyFSvXp0tW7awdetWqlSpQmJiYtQUP263m7i4OA4fPhy1f/IHcypIUVERf//9N263m+rVq1vdJcvExMSQlZUV1XmQkpQL8RUxmXj/fXO+7969kJQEzz0Ht9+uub4nIFSZ0C4QASrNLhDFDh48yJ49eygsLAxx7yScJScnk56eTlJSktVdERGRUDp4EB58EF54wTzOyjIvZ3zaadb2K4oEWq/pDHAIJSUlUbduXZxOZ2SuaD0Kh8PBF198wXnnnRf1c17j4uKIi9P/Zg6Hg4ULF3LZZZdFfSbkf5QL8VWuM7FuHdx8M/z0k3n80EMwahRoX/yTEqpM6CdzGYi2Ishms7Fy5UouuugiXRBDAHNazJo1a7jkkkus7oqEEeVCfJXLTLjd5hSHhx8Gux1q1oSZM6FLF6t7FhFClYnoqcpEREREgmnnTujTBxYtMo+7dYNXXoFq1SztlhxfOZ9lLiIiImKBjz6CM880i98KFWDqVHPxm4rfckFngCXoYmNj6dSpE7GxsVZ3RcKEMiH+KBfiq1xk4vBhc7rDc8+Zx2eeaS50a9bM2n5FqFBlQrtABOhEdoEQERGRCPLDD+ZCt+++M4/vvx/GjTPPAEtYCLRe0xQICTq73c7s2bOx2+1Wd0XChDIh/igX4itsM2EY5hSHNm3M4rd6dfjwQ3j2WRW/IRaqTGgKhASdYRhs3LgR/XFBiikT4o9yIb7CMhO7d5sXtViwwDy+9FKYMQMyMiztVrQIVSZ0BlhERESiVmFhITabDZvNVvLCVUuWQIsWZvGbkACTJplnflX8lns6AywiIiJypKIiGDYMnnnGPG7a1Fzo1rKltf2SoFEBLEEXFxdHt27douriH3JsyoT4o1yIr7DIxC+/mAvd1qwxj++5B55+GnQ5e0uEKhPaBSJA2gVCRESk/CgxneEY98v4Z0pD3jPPkPrYY9gOHsSoUoWiqVNxXXllicckJycHta8SPNoFQixjt9uZOnVq+K3iFcsoE+KPciG+gpmJlJSUgL4yMjI4BXgbqPTgg9gOHuS/QK2//6biTTf5fYyUnVB9TqgAlqAzDIPdu3eH1ypesZQyIf4oF+LLikx0AtYD1wEO4N9AV+CvMuuBHEuoMqGJVyIiIhJxCgoKjn0Hh4P4J58k/plnsBkGvwLpixeT3b492WXRQbGUCmARERGJOMecp7thA/ToAStXAvAyMBDY2b695vdGCU2BkKCLj4+nZ8+exMfHW90VCRPKhPijXIivkGfCMGDmTMjKMovfypU5/Prr9AUCWzInZS1UmdAuEAHSLhAiIiLl2P795pZmb75pHp9/PsyeTWGVKp6FbQUFBToDXM5pFwixTFFREWPGjKGoqMjqrkiYUCbEH+VCfIUsE8uWQatWZvEbGwujRsGnn0KdOsF9HQm6UGVCc4AlJLStkfhSJsQf5UJ8BTUTTqdZ7D7xBLjd0LAh5OTAOed47pKcnKydSMJcKD4nVACLiIhI5NmyBXr2hOXLzeNbboEpU0DTGAVNgRAREZFI88Yb0LKlWfxWqmSe9Z01S8WveGgRXIC0CC5wbrebPXv2kJ6eTkyMfscSZUL8Uy7E10ln4sABuPdes9gFaN8eZs+GBg2C21EpM6XNhBbBiWVsNhtpaWnYbDaruyJhQpkQf5QL8XVSmfj6a3Oh26xZEBMDI0ZAbq6K33IuVJ8TKoAl6Ox2O2PHjtXiFvFQJsQf5UJ8nVAmXC548kno0AE2bYK6dc3CNzsb4rTUqbwL1eeEkiEiIiLl07Zt0KsXfP65efx//wcvvACVK1vaLQl/OgMsIiIi5c/bb5sL3T7/HJKTYcYMc/Gbil8JgM4Ai4iISPlRWAgDBsArr5jHbdvCnDnQuLG1/ZJyRbtABEi7QATOMAzsdjsJCQla3CKAMiH+KRfi67iZWLUKevSAX38Fmw2GDIGRIyE+vuw7K2WitJ8T2gVCLGMYBnl5ebqyjngoE+KPciG+jpoJtxueegrOPdcsfmvVMi9lPHq0it8IF6rPCRXAEnQOh4Np06bhcDis7oqECWVC/FEuxJffTGzfDl27wuDB4HDAtdfCunXQubNl/ZSyE6rPCRXAIiIiEp7mz4cWLeCTTyApCaZPNxe/Va1qdc+knNMiOBEREQkvBw/CAw/AtGnmcVaWudDt9NOt7ZdEDBXAEhIJCQlWd0HCjDIh/igX4qvW3r3Ed+gAP/1kNjz4oHmhi8REazsmlgnF50S53wXC6XSSnZ1NTk4OO3bsoGbNmvTp04dHH33Uc81owzAYOXIk06dPZ9++fbRr147nn3+eZs2aBfw62gVCREQkhAwDnnvOnOtrt0ONGuZljbt0sbpnUo5EzS4Q48aN44UXXmDKlCn89NNPjB8/nqeeeorJkyd77jN+/HgmTJjAlClTWLlyJTVq1KBLly4cOHDAwp5HLrfbzYYNG3C73VZ3RcKEMiH+KBfisXMnXHEFDBwIdjvGFVfA+vUqfiVknxPlvgD+6quvuOqqq7jiiiuoX78+119/PV27duXbb78FzLO/kyZNYtiwYVx77bU0b96cmTNncvDgQebMmWNx7yOTw+EgJydHK7vFQ5kQf5QLAWDhQnOh28KFGBUq8NHll2N/+22oVs3qnkkYCNXnRLmfA9yxY0deeOEFfv31V0499VTWrVvHsmXLmDRpEgCbN29mx44ddO3a1fOYxMREOnXqxPLly7nrrrv8Pm9RURFFRUWe4/z8/BLtMTExxMfH43A4vH4ziY2NJS4uDrvd7rVvXVxcHLGxsSXa4+PjiYmJ8Xq94nabzYbdbvdqT0hIwDCMEmFITEzE7XZ7tdtsNhISEnC5XDidzhLtTqcTl8vlaQ/GmI4cR6SMqbjvkfQ+leWYip+/qKgoYsYUie9TWY+p2JHjKu9jisT3KWRjKioiYfhwePZZANzNmnHw5ZdZ+fHHXGSzlc8xReL7FCZjKv7v8cbk+z04mnJfAD/88MPk5eVx+umnExsbi8vl4sknn+Tmm28GYMeOHQBkZGR4PS4jI4OtW7ce9XnHjBnDyJEjS7RPmDCBChUqAJCVlUX37t1ZuHAha9as8dynU6dOdO7cmblz57Jx40ZPe7du3WjdujUvv/wyu3fv9rT37NmTxo0bM2HCBK839J577iEtLY2xY8d69WHIkCHk5eUxrXh1LOYb/8gjj7Bp0yZycnI87dWqVaNfv36sW7eOBQsWeNobNWpEr169WLZsGbm5uZ72YI3pyO9XpIwpEt+nsh7TxIkTI25MEHnvU1mNqU6dOoCZi0gZUyS+T6EYU7Vdu7jx/fdJ//NPAL5u144lF1+M6+OPPY8rb2OCyHufwmFMzz//PPC/z4njjenIz5NjKfeL4N58803+/e9/89RTT9GsWTPWrl3LwIEDmTBhAr1792b58uV06NCB7du3U7NmTc/j+vbty7Zt21i0aJHf5/V3BrhOnTrs2rXLM6lav7n5H5PdbmfWrFn07du3xJVbyuuYivseSe9TWY7p4MGDzJw5k969e5OYmBgRY4rE96msx+R0Opk+fTq9e/f2rPIu72OKxPcpqGOy24l56SXiBg/GdvgwVKuG6+WXcV5yCeD98yMmJqZ8jCkS36cwGtOBAwc8Pz8SEhKOO6a9e/dSvXr14y6CK/cFcJ06dRgyZAj9+/f3tI0aNYrZs2fz888/s2nTJho1asTq1avJysry3Oeqq66icuXKzJw5M6DX0S4QIiIiJ2HPHvjXv8yLWwBccgnMmGHu9iASJFGzC8TBgwc9250Vi42N9fzW0aBBA2rUqMGSJUs8t9vtdnJzc2nfvn2Z9jVauFwuVq9e7fUboUQ3ZUL8US6iyH//ay50mz8fEhJg4kT46KMSxa8yIb5ClYlyXwB369aNJ598kg8//JAtW7Ywb948JkyYwDXXXAOYp+gHDhzI6NGjmTdvHt9//z19+vQhKSmJHj16WNz7yOR0OlmwYIHXn0okuikT4o9yEQXsdnNf3y5d4K+/oGlT+Pprc7uzmJIliDIhvkKViXK/CG7y5Mk89thj9OvXj127dpGZmcldd93F8OHDPfcZPHgwhw4dol+/fp4LYSxevJjU1FQLey4iIhLBfvkFevSA1avN47vvhmeegaQka/slQgQUwKmpqUyaNMmz7Zk/NpuN7OxssrOzy6xfIiIiUckw4JVXYMAAOHgQqlQxj6++2uqeiXiU+wJYwo/NZqNRo0bYbDaruyJhQpkQf5SLCPT333DnnfDOO+bxhRealzOuVSughysT4itUmSj3u0CUFe0CISIicgy5udCrF/zxB8TFwZNPwkMP+Z3rKxIqUbMLhIQfp9PJ0qVLtYhBPJQJ8Ue5iBAOBwwbBhdcYBa/TZrAV1+Zi99KWfwqE+IrVJlQASxB53K5yM3N1TY24qFMiD/KRQTYuBE6doTRo825v7ffbi56a9PmhJ5OmRBfocqECmAREREpHcMw5/a2agXffAOVK8PcueZit5QUq3snclxaBCciIiKBy8uDe+6BN94wj887D2bPhrp1re2XSCnoDLAEXUxMDFlZWSWu0CfRS5kQf5SLcmj5cvOs7xtvQGwsPPEEfPZZ0IpfZUJ8hSoT2gUiQNoFQkREopbTae7q8Pjj4HZDgwYwZw6cc47VPRPxol0gxDIOh4P58+fjcDis7oqECWVC/FEuyoktW6BzZ8jONovfW26BtWtDUvwqE+IrVJlQASxB53a7WbNmDW632+quSJhQJsQf5aIcePNNaNkSvvwSUlPNub6zZkGI/hKqTIivUGVCi+BERETE24EDcN99MHOmeXzOOZCTAw0bWtsvkSDRGWARERH5n2++gawss/iNiYHHHoMvvlDxKxFFZ4Al6GJjY+nUqROxsbFWd0XChDIh/igXYcblgvHjYfhwc9Fb3brmlIfzziuzLigT4itUmdAuEAHSLhAiIhKxtm0zF7fl5prHN94IL75oXuBCpBzRLhBiGbvdzuzZs7Hb7VZ3RcKEMiH+KBdh4p13zIVuubmQnAyvvWYufrOg+FUmxFeoMqEpEBJ0hmGwceNG9McFKaZMiD/KhcUKC+GBB+Cll8zjNm3MvX2bNLGsS8qE+ApVJnQGWEREJNqsXg1nnWUWvzYbDBlibnVmYfErUpZ0BlhERCRauN0wcSI88gg4HFCrFrz+OlxwgdU9EylTKoAl6OLi4ujWrRtxcYqXmJQJ8Ue5KGN//QW9e8OSJebxNdeYZ4CrVrW2X0dQJsRXqDKhXSACpF0gRESk3FqwAG6/HfbsgYoVYdIk6NvXnP4gEkG0C4RYxm63M3XqVK3iFQ9lQvxRLsrAoUPQvz90724Wv61amfN/77wzLItfZUJ8hSoTKoAl6AzDYPfu3VrFKx7KhPijXITY+vXQti1MnWoeDxoEK1bA6adb269jUCbEV6gyoUk2IiIikcQwYPJkGDwYioogIwNmzYKuXa3umUjYUAEsIiISKXbtgttug48+Mo+vuAJefRWqV7e2XyJhRovgAqRFcIFzu91s2rSJhg0bEhOjWTaiTIh/ykWQLVpk7vKwaxckJsIzz0C/fmE51/dolAnxVdpMBFqvqQAOkApgEREJS0VF5oUsJk0yj5s3hzfeMP8rEmW0C4RYpqioiDFjxlBUVGR1VyRMKBPij3IRBD/+CO3a/a/4vfde+Oabclv8KhPiK1SZUAEsIaEtbMSXMiH+KBcnyDDghRfMyxmvWwfp6eZev5Mnm/v8lmPKhPgKRSa0CE5ERKQ82bMH7rgD3n/fPO7aFWbOhBo1rO2XSDmiM8AiIiLlxSefQIsWZvGbkAATJsDChSp+RUpJi+ACpEVwgXO73ezZs4f09HSt4hVAmRD/lItSsNvhscfgqafM6Q+nnw5z5kBWltU9CyplQnyVNhOB1muaAiFBZ7PZSEtLw1aOtt6R0FImxB/lIkC//go9esCqVebxXXeZZ36TkqztVwgoE+IrVJnQr1cSdHa7nbFjx2ohg3goE+KPcnEchmFexCIryyx+q1SBd981F79FYPELyoSUFKpM6AywiIhIuNm3D+68E95+2zy+4AJ4/XWoVcvafolECJ0BFhERCSeffw4tW5rFb1wcjB0LS5ao+BUJIp0BFhERCQcOBzz+ODz5pDn9oXFjc6Fb27ZW90wk4mgXiABpF4jAGYaB3W4nISFBCxkEUCbEP+XiCJs2mQvdvv7aPL7tNnjuOUhJsbZfZUyZEF+lzYQuhSyWMQyDvLw89LuVFFMmxB/l4h+zZ0OrVmbxm5YGb71lLn6LsuIXlAkpKVSZUAEsQedwOJg2bRoOh8PqrkiYUCbEn6jPRV4e9OwJt9wCBw5Ax47mZY1vvNHqnlkm6jMhJYQqEyqARUREytry5eZZ3zlzIDYWnngCli6FevWs7plIVNAiOBERkbLidMLo0eZiN5cL6tc3i+Bzz7W6ZyJRRQWwhERCQoLVXZAwo0yIP1GVi61boVcvWLbMPO7ZE55/nsK4OFL+WdxTUFBAcnKyhZ20XlRlQgISikxoF4gAaRcIERE5YW+9ZV7COC8PUlNh6lSzGAYKCwtJ+WfBmwpgkZOjXSDEMm63mw0bNuB2u63uioQJZUL8ifRcFBYWUrhjB45bboGbboK8PFxnn83B5cspvOYa8/Z/vrweU4qvSBPpmZDSC1UmVABL0DkcDnJycrSKVzyUCfEn0nPROSWF7TVrEj97Ni7gcaDCN9+QfOaZpKSkeL4yMjI8j8nIyPC67XhfkSbSMyGlF6pMaA6wiIhIMLlc8NRTLAfigd+BXsAX1vZKRI6gAlhERCRY/vgDbr0VPvuMeMB5zTVUfe45Fp5yylEfUlhY6DkLvHPnTs0BFikDKoAl6Gw2G9WqVdNlLMVDmRB/Ii4X774Ld9wB+/ZBcjJMnkxcnz7ElWJ8ycnJUV0AR1wm5KSFKhPaBSJA2gVCRET8KiyEQYNg+nTzuE0bc2/fJk0CfLh2gRAJFu0CIZZxuVysXr0al8tldVckTCgT4k9E5GLNGrPgnT4dbDZ4+GH48suAi18wz/oahoFhGFFf/EZEJiSoQpUJFcASdE6nkwULFuB0Oq3uioQJZUL8Kde5cLvhmWegXTv4+WfIzIQlS2DsWNCFHE5Yuc6EhESoMqE5wCIiIqXx11/Qu7dZ8AJcfTW8/DJUrWppt0QkcDoDLCIiEqgPPoAWLczit2JFeOEFc/Gbil+RckVngCXobDYbjRo10ipe8VAmxJ9ylYtDh2DwYJgyxTxu2RLeeAOaNrW2XxGmXGVCykSoMqFdIAKkXSBERKLU99/DzTeb/wV44AEYMwYSE63tl4iUoF0gxDJOp5OlS5dqEYN4KBPiT9jnwjBg8mRzl4fvv4eMDFi0CCZMUPEbImGfCSlzocqECmAJOpfLRW5urraxEQ9lQvwJ61zs2gVXXgn33w9FRXD55bB+PVxyidU9i2hhnQmxRKgyoQJYRETkSB9/bC50++gj80zvc8+Zi9+qV7e6ZyISJFoEJyIiAuaZ3qFDzSkOAM2amQvdzjzT2n6JSNCpAJagi4mJISsri5gY/YFBTMqE+BNWufjpJ+jRA9auNY/794ennjK3OpMyE1aZkLAQqkxoF4gAaRcIEZEIZBjmZYwfeMDc6iw9HV59Fbp1s7pnInICtAuEWMbhcDB//nwcDofVXZEwoUyIP5bnYu9euPZauPtus/jt0sVc6Kbi1zKWZ0LCTqgyoQJYgs7tdrNmzRrcbrfVXZEwoUyIP5bm4tNPzYVu770H8fHwzDPmFmc1a5Z9X8RDnxXiK1SZ0BxgERGJHnY7DB8O48eb0x9OO81c6JaVZXXPRKQMqQAWEZHo8Ntv5kK3b781j++809zxITnZ2n6JSJnTFAgJutjYWDp16kRsbKzVXZEwoUyIP2WWC8OA114zz/J++y2ccgq88w68+KKK3zCjzwrxFapMaBeIAGkXCBGRcmjfPnOR29y55nHnzvD661C7tqXdEpHQ0C4QYhm73c7s2bOx2+1Wd0XChDIh/oQ8F198Aa1amcVvXByMGQP//a+K3zCmzwrxFapMaA6wBJ1hGGzcuBH9cUGKKRPiT8hy4XTC44/Dk0+C2w2NGsGcOXD22cF9HQk6fVaIr1BlQgWwiIhEjs2boWdP+Oor87h3b5g8GVJTre2XiIQVTYEQEZHIkJMDLVuaxW9amrm92YwZKn5FpASdAZagi4uLo1u3bsTFKV5iUibEn6DlIj8f+veH2bPN4w4dzGK4Xr2T76SUKX1WiK9QZUK7QARIu0CIiIShFSvMvX03b4bYWPMiF0OHmoveRCTqaBcIsYzdbmfq1KlaxSseyoT4c1K5cLlg1Cjo2NEsfuvXh88/NwtgFb/llj4rxFeoMqFPCQk6wzDYvXu3VvGKhzIh/pxwLn7/HXr1Mrc5A7j5Zpg2zZz3K+WaPivEV6gyoTPAIiJSfsyday50++ILSEmBWbPM+b4qfkWkFHQGWEREwl9BAQwYAK++ah63a2cWvo0aWdsvESmXtAguQFoEFzi3282mTZto2LAhMTH6I4MoE+JfwLn49ltzodtvv4HNZi5yGzEC4uPLrrNSJvRZIb5Km4lA6zUVwAFSASwiUsbcbnj6aRg2zLy6W+3a5lZnnTpZ3TMRCVPaBUIsU1RUxJgxYygqKrK6KxImlAnx55i5+PNP6NIFHn7YLH6vvx7WrVPxG+H0WSG+QpUJFcASEtrCRnwpE+KP31y89x60aAGffgpJSfDyy+bitypVyrx/Uvb0WSG+QpEJLYITEZHwcPAgDBoEL75oHrduDXPmwGmnWdsvEYk4EXEG+M8//6RXr15UrVqVpKQkWrVqxapVqzy3G4ZBdnY2mZmZVKxYkc6dO/PDDz9Y2GMREfGydi20afO/4vff/4avvlLxKyIhUe4Xwe3bt4+srCwuuOAC7rnnHqpXr87GjRupX78+jf7ZHmfcuHE8+eSTzJgxg1NPPZVRo0bx+eef88svv5CamhrQ62gRXODcbjd79uwhPT1dq3gFUCbEP7fbzZ5du0ifM4eYRx4Bux1q1jT39r34Yqu7JxbQZ4X4Km0momYXiCFDhvDll1/yRfEVgXwYhkFmZiYDBw7k4YcfBswJ1RkZGYwbN4677roroNdRARw4wzCw2+0kJCRgs9ms7o6EAWVC/DH++gujd29iliwxG7p3h1degfR0azsmltFnhfgqbSYCrdfK/Rzg+fPnc8kll3DDDTeQm5tLrVq16NevH3379gVg8+bN7Nixg65du3oek5iYSKdOnVi+fPlRC+CioiKvFYf5+fkl2mNiYoiPj8fhcOB2uz33jY2NJS4uDrvd7nXpvri4OGJjY0u0x8fHExMTU2KFY3x8PDabrcTk74SEBAzDwOFweLUnJibidru92m02GwkJCbhcLpxOZ4l2p9OJy+XytAdjTEVFRUycOJEhQ4aU+L6W1zEV9z2S3qeyHFNhYSETJ07kgQceoEKFChExpkh8n8p0TEuWYLv9dmJ278aoUAHnU0/hvuMO4hMSsP3zA6/cjSkS36cyHtORPz9iY2MjYkxHtkfK+1SWYzpw4IDn50diYuJxxxTobhHlvgDetGkT06ZNY9CgQQwdOpRvvvmG+++/n8TERG699VZ27NgBQEZGhtfjMjIy2Lp161Gfd8yYMYwcObJE+4QJE6hQoQIAWVlZdO/enYULF7JmzRrPfTp16kTnzp2ZO3cuGzdu9LR369aN1q1b8/LLL7N7925Pe8+ePWncuDETJkzwekPvuece0tLSGDt2rFcfhgwZQl5eHtOmTfO0JSQk8Mgjj7Bp0yZycnI87dWqVaNfv36sW7eOBQsWeNobNWpEr169WLZsGbm5uZ72YI3pyO9XpIwpEt+nsh7TxIkTI25MEHnvUyjHFOtw0GXJEtp98w0AOzIyeOe669izdy+MG1cuxxSJ75OVYyoWSWOKxPeprMb0/PPPA+bPj0DGVHy/4yn3UyASEhJo06YNy5cv97Tdf//9rFy5kq+++orly5fToUMHtm/fTs2aNT336du3L9u2bWPRokV+n9ffGeA6deqwa9cuzyl1/eamM8AaU2Bj0hlgjclut8P33xPXuzcx338PgPO++xhbqRL3Dx5MYmJiuRxTJL5POgOs9ymcxpSfn1+qM8B79+6levXqkT8FombNmpxxxhlebU2bNuWdd94BoEaNGgDs2LHDqwDetWtXibPCR0pMTPR8IB+vPf4ol+NMSEgoVbu/1ztau81m89seExPjtz02NpbY2NgS7XFxccTFlYyBxqQxHa39RMZU3J6YmOjpQ3kfUyS+TyEbk2GQ8NJL8NBDcPgwVK8OM2bguvBCXGPH+v1cDfsxHaO93L5Px2jXmDQmq8fk+zlRmjH5U+6XWHbo0IFffvnFq+3XX3+lXr16ADRo0IAaNWqwpHiRBeaGyrm5ubRv375M+xotEhISGDJkyFHDLNFHmYhiu3ebi9vuvdcsfi+7DNavh8suUy6kBGVCfIUqE+W+AH7ggQdYsWIFo0ePZsOGDcyZM4fp06fTv39/wPxtYODAgYwePZp58+bx/fff06dPH5KSkujRo4fFvY9MhmGQl5fn9acMiW7KRJRassS8otsHH0BCAjz7LHz4Ifzz1zflQnwpE+IrVJko9wVw27ZtmTdvHm+88QbNmzfniSeeYNKkSfTs2dNzn8GDBzNw4ED69etHmzZt+PPPP1m8eHHAewBL6TgcDqZNm1ZibpFEL2UiyhQVmdMdunaFHTvgjDNg5Uq4/344Yhsj5UJ8KRPiK1SZKPdzgAGuvPJKrrzyyqPebrPZyM7OJjs7u+w6JSISjX7+GXr0gOLV3/fcA08/DUlJ1vZLROQI5f4MsIiIhAHDgJdegtatzeK3alV4/32YOlXFr4iEnZMqgBcuXKh5OuKXFjCIL2Uigv39N1x/Pdx5Jxw6ZF7GeP16c/HbcSgX4kuZEF+hyMRJ7QMcExNDZmYmvXr1onfv3jRt2jSYfQsruhSyiIgfS5dCr17w558QHw+jR8OgQRCjPzCKSNkLtF47qU+oZs2asX37dp566imaN2/OOeecwwsvvMD+/ftP5mmlnHO73WzYsMFrA2yJbspEBHI4YOhQuPBCs/g99VRYscJc/BZg8atciC9lQnyFKhMnVQB/9913fPvtt/Tv35+qVavyzTff0L9/f2rWrMlNN93EokWLNEUiCjkcDnJycrSKVzyUiQizYQN06ABjxphzf++4A1avNuf/loJyIb6UCfEVqkyc9N+oWrduzXPPPcf27dt599136d69O263m7lz53LFFVdQu3ZthgwZwk8//RSM/oqIiFUMA2bOhKwsc1uzU06Bt982F78lJ1vdOxGRgAVtklZcXBxXX3018+bNY/v27Tz77LO0atWKv/76yzNFol27dpoiISJSHu3fb25v1qcPFBRAp06wbh1cd53VPRMRKbWQrFKoWrUq9913H9988w1jx44lNjYWwzBYuXIl/fv3JzMzk3/9619s3rw5FC8vFrPZbFSrVg3bERveS3RTJsq5ZcugVSt4802IjYUnn4RPPoE6dU7qaZUL8aVMiK9QZeKkdoE4mh9++IGZM2eSk5PDjh07MAyD9PR0evTowc6dO3n//fc5fPgwSUlJLFy4kPPOOy/YXQg67QIhIlHH6YRRo+CJJ8DthoYNYc4caNfO6p6JiPhVJrtAHGnv3r1MnjyZNm3a0KJFC55++ml27drFpZdeyn/+8x/+/PNPJk2axBtvvMEff/xB//79OXjwIIMHDw5WFyRMuFwuVq9ejcvlsrorEiaUiXJo82ZzmsPIkWbxe+ut5gUuglj8KhfiS5kQX6HKxEkVwE6nk/fee49rrrmGWrVqMXDgQFavXk3jxo158skn+f333/nwww+57rrriI+P9zyuSpUqTJ48mSZNmrBu3bqTHoSEF6fTyYIFC3A6nVZ3RcKEMlHOvPGGOeVh+XKoVMk86ztzpvnvIFIuxJcyIb5ClYm4k3lwZmYme/fuxTAMkpOTufnmm7n99tsDntJQs2ZNNmzYcDJdEBGRYMnPh3vvhddfN4/bt4ecHKhf39JuiYgE20kVwHv27OHcc8/l9ttv5//+7/9ISUkp1eOfffZZ7QghIhIOvv7a3OVh0ybzQhaPPgqPPQZxJ/VjQkQkLJ3UJ9vPP//MqaeeesKPb9my5cm8vIQpm81Go0aNtIpXPJSJMOZywdixMGKE+e969WD2bOjYMeQvrVyIL2VCfIUqEyHZBSISaRcIEYk427ZBr17w+efm8U03wbRpULmypd0SETlRZb4LhEgxp9PJ0qVLtYhBPJSJMPT229CihVn8pqSYi9zmzCnT4le5EF/KhPgKVSZUAEvQuVwucnNztY2NeCgTYaSgAO64A264wby6W9u25vZmt94KZfxnZ+VCfCkT4itUmVABLCISLVatgrPOgldeMYvdoUPhyy+hcWOreyYiUqa0vFdEJNK53fDMMzBsGDgcUKuWudCtc2ereyYiYgmdAZagi4mJISsri5gYxUtMyoSFtm+Hrl1h8GCz+L32Wli//qjFb2FhITabDZvNRmFhYUi7plyIL2VCfIUqE9oFIkDaBUJEyp3334d//Qv27oWkJHj2WfP4GHN9CwsLPXu6FxQUkJycXFa9FRE5adoFQizjcDiYP38+DofD6q5ImFAmytjBg3DPPXD11Wbxm5UFq1ebi9/CaH9V5UJ8KRPiK1SZUAEsQed2u1mzZg1ut9vqrkiYUCbK0Lp10KYNvPCCefzQQ/DVV3Daadb2yw/lQnwpE+IrVJnQIjgRkUhgGPDcc+ZcX7sdatSAWbOgSxereyYiEnZUAIuIlHc7d0KfPrBokXncrZu51Vm1aiXueryFbUfeHsgiOM0RFpHySAWwBF1sbCydOnUiNjbW6q5ImFAmQuijj+C222DXLqhQwdzu7J57jjrXt3iBWyAyMjKOe5+TWUetXIgvZUJ8hSoT2gUiQNoFQkTCyuHD8PDD5rQHgDPPhDfegGbNjvkwW5AXwelHiIiEE+0CIZax2+3Mnj0bu91udVckTCgTQfbDD3D22f8rfu+/H7755rjFL5hbmx3ra+fOnZ777ty587j3PxnKhfhSJsRXqDKhKRASdIZhsHHjRp0ZEg9lIkgMA6ZNgwcfNM8AV6sGM2bA5ZcH/BSlmbObnJwc0jm+yoX4UibEV6gyoQJYRKQ82LPHvIjF/Pnm8aWXwmuvmbs9iIhIqWgKhIhIuPvvf6FFC7P4TUiAiRPhww9V/IqInCCdAZagi4uLo1u3bsTFKV5iUiZOkN0Ow4bB00+bx02bmgvdWra0tl9BolyIL2VCfIUqE9oFIkDaBUJEytQvv8DNN8OaNebx3XebW5wlJYX0ZQsLCz1bpRUUFGifXxEpV7QLhFjGbrczdepUreIVD2WiFAwDXn4ZWrc2i98qVWDePHPxW4iLXzAXvhmGgWEYIS9+lQvxpUyIr1BlQn9jkKAzDIPdu3drFa94KBMB+vtvuPNOeOcd8/jCC83LGdeqZW2/QkS5EF/KhPgKVSZ0BlhEJBwsXWrO7X3nHYiLg/HjYcmSiC1+RUSspDPAIiJWcjggOxvGjDGnPzRpYi50O+ssq3smIhKxtAguQFoEFzi3282mTZto2LAhMTH6I4MoE0e1YQP07GlexQ3MfX4nTYJ/FqFFOuVCfCkT4qu0mQi0XlMBHCAVwCISNIZhzu29914oKIDKlWH6dLjhBqt7JiJSrmkXCLFMUVERY8aMoaioyOquSJhQJo6wfz/06AF9+pjF7/nnw7p1UVn8KhfiS5kQX6HKhApgCQltYSO+lAngyy+hVSt4802IjYVRo+DTT6FuXat7ZhnlQnwpE+IrFJnQIjgRkVBzOs1i94knwO2Ghg0hJwfOOcfqnomIRCUVwCIiobRli7nQbfly8/iWW2DKFNBaAhERy2gRXIC0CC5wbrebPXv2kJ6erlW8AkRxJt54w7yEcX4+pKaaV3Pr2dPqXoWNqM2FHJUyIb5Km4lA6zWdAZags9lspKWlYbPZrO6KhImoy8SBA+YOD7NmmcfnnmtOeWjQwNp+hZmoy4UclzIhvkKVCf16JUFnt9sZO3asFjKIR1Rl4uuvzYVus2ZBTAwMHw6ff67i14+oyoUERJkQX6HKhM4Ai4gEg8sF48aZBa/LZe7skJMDHTta3TMREfGhAlhE5GRt22YubsvNNY9vvBFefNG8wIWIiIQdTYEQEQlQYWEhNpsNm81GYWGh2fj229CypVn8JifDa6+Z+/yq+BURCVvaBSJA2gUicIZhYLfbSUhI0EIGASInE4WFhaSkpABQsHMnyUOHwiuvmDe2aQNz5kCTJhb2sHyJlFxI8CgT4qu0mdClkMUyhmGQl5eHfreSYpGWidZAxY4dzeLXZoMhQ8yrvKn4LZVIy4WcPGVCfIUqEyqAJegcDgfTpk3D4XBY3RUJE+GeicLCwsC+DhzgIeArIOa333BnZnLogw8ofPRRCh2Ooz5O/Av3XEjZUybEV6gyoUVwIhL1iqc1HEtNYCbw1D/H7wJ9t2/n7yuuOO5jdTZLRCS86AywiMhxdAPWA12Ag8CdwHXA31Z2SkRETpjOAEtIJCQkWN0FCTPhnImCggL/Nxw8SMLQocS//DIAjubNaf399/wC7Ny5k+Tk5LLrZIQK51yINZQJ8RWKTGgXiABpFwiRKLN+Pdx8M/z4o3n84IMUDh1KStWqgFk0qwAWEQkv2gVCLON2u9mwYQNut9vqrkiYKFeZMAx49llo29YsfmvUgI8/hqefhsREq3sXUcpVLqRMKBPiK1SZUAEsQedwOMjJydEqXvEoN5nYuROuuAIGDgS7Ha680jwT3LWr1T2LSOUmF1JmlAnxFapMaA6wiAjAwoXQpw/s2mWe6X3mGejXz9zn9x/Jycna0UFEJAKoABaR6Hb4sHkhi2efNY+bN4c33jD/KyIiEUlTICTobDYb1apV02UsxSNsM/HDD9Cu3f+K3/vug2++UfFbRsI2F2IZZUJ8hSoT2gUiQNoFQiSCGAa88AIMGmSeAa5WDV57zZz/KyIi5ZZ2gRDLuFwuVq9ejcvlsrorEibCKhN79sDVV5vzew8fNhe4rV+v4tcCYZULCQvKhPgKVSZUAEvQOZ1OFixYgNPptLorEibCJhP//S+0aAHz50NCAkyYYC5+q1HD2n5FqbDJhYQNZUJ8hSoTWgQnIpHPbodHH4WnnjKPTz/dXOjWqpWl3RIREWuoABaRyPbLL9CjB6xebR7fdZd55jcpydp+iYiIZTQFQoLOZrPRqFEjreIVD0syYRjwyivQurVZ/FapAvPmmYvfVPyGBX1WiC9lQnyFKhPaBSJA2gVCpBzZtw/uvBPefts8vuACeP11qFXL2n6JiEhIaRcIsYzT6WTp0qVaxCAeZZqJ3Fxzodvbb0NcHIwdC0uWqPgNQ/qsEF/KhPgKVSZUAEvQuVwucnNztY2NePjLRGFhITabDZvNRmFh4cm/iMMBw4aZZ3v/+AMaN4bly+HhhyE29uSfX4JOnxXiS5kQX6HKhBbBiUj5t3GjudDtm2/M49tvN6/ulpJibb9ERCQs6QywiJRfhmHO7W3Vyix+09LgrbfMxW8qfkVE5Ch0BliCLiYmhqysLGJi9PtVNPI3ncHhcNC0aVMOHTrkmcd15P0CnQKRnJz8v4O8PLjnHnM/X4DzzoPZs6Fu3RPvvJQpfVaIL2VCfIUqE9oFIkDaBUIkMKHcvsjzcbV8OfTsCVu2mPN7s7PhkUc011dEJMppFwixjMPhYP78+TgcDqu7IpHG6YSRI82zvVu2QIMG8MUX5lXeVPyWO/qsEF/KhPgKVSZUAEvQud1u1qxZg9vttrorYoGCgoISX3v37mXo0KHs3bvX07Zz507PY3bu3On3cUd+Ff74I3TubJ7tdbuhVy9YuxbOPdeqocpJ0meF+FImxFeoMqE5wCISVF7zdP8RFxdHQkICycnJJCYm+n2Mv8d5vPkm3H23Oe83NRWmTTOnQIiIiJwAFcAiEr4OHID77oOZM83jc86BnBxo2NDafomISLmmAliCLjY2lk6dOhGrOZnyD3+ZSE5O5phrcL/5xtzbd+NGiIkxL3Lx2GMQH18GPZayoM8K8aVMiK9QZUK7QARIu0CIlBGXC8aPh+HDzUVvdeqYZ33PO8/qnomISJjTLhBiGbvdzuzZs7Hb7VZ3RcJEwJn44w+4+GIYOtQsfm+8EdatU/EbofRZIb6UCfEVqkyoAJagMwyDjRs3HvvP2xJVAsrEu+9CixawdCkkJ8Orr5qL3045pcz6KWVLnxXiS5kQX6HKRMQVwGPGjMFmszFw4EBPm2EYZGdnk5mZScWKFencuTM//PCDdZ0Ukf8pLIQ774TrroN9+6BNG1izBm67DUJ4UQ0REYleEVUAr1y5kunTp9OiRQuv9vHjxzNhwgSmTJnCypUrqVGjBl26dOHAgQMW9VREAFi9Gs46C156ySx2H34YvvwSmjSxumciIhLBIqYALigooGfPnrz00kuccsSfTA3DYNKkSQwbNoxrr72W5s2bM3PmTA4ePMicOXMs7HHkiouLo1u3bsTFaZMRMZXIhNsNzzxjbmv2yy+QmQn//S+MHQsJCdZ2VsqMPivElzIhvkKViYjZBaJ3795UqVKFiRMn0rlzZ1q1asWkSZPYtGkTjRo1YvXq1WRlZXnuf9VVV1G5cmVmFu8v6qOoqIiioiLPcX5+PnXq1GHXrl2eVYUxMTHEx8fjcDi8rlASGxtLXFwcdrvda85KXFwcsbGxJdrj4+OJiYnxer3idpvNVmLid0JCAoZhlLgsYGJiIm6326vdZrORkJCAy+XC6XSWaHc6nbhcLk+7xqQxhXxMO3aQ0LcvLFkCgKt7d5zTphFTrVr5HVMkvk8ak8akMWlM5XBMe/fupXr16sfdBSIifsV68803Wb16NStXrixx244dOwDIyMjwas/IyGDr1q1Hfc4xY8YwcuTIEu0TJkygQoUKAGRlZdG9e3cWLlzImjVrPPfp1KkTnTt3Zu7cuWzcuNHT3q1bN1q3bs3LL7/M7t27Pe09e/akcePGTJgwwesNveeee0hLS2Ps2LFefRgyZAh5eXlMmzbN05aQkMAjjzzCpk2byMnJ8bRXq1aNfv36sW7dOhYsWOBpb9SoEb169WLZsmXk5uZ62oM1prS0NPr16xdRY4rE96msx3TqL79wzQcfwIEDuBIT+ahLF1ZnZcFLL5XbMUXi+1RWY6pbty7jxo3z+qFX3scUie9TWY8pPj6ehx56iOXLl0fMmCLxfQrXMU2cOJFAlPszwNu2baNNmzYsXryYli1bAnidAV6+fDkdOnRg+/bt1KxZ0/O4vn37sm3bNhYtWuT3eXUG+MTHVFRUxMSJExkyZEiJ72t5HVNx3yPpfSrLMRXu2cPP3bpxdvEvqS1b4nz9dVynnlpuxxSJ71NZj8nhcDB27FgeeOABzyWyy/uYIvF9KssxHfnzIzY2NiLGdGR7pLxPZTmm/Px8Jk6c6Pmc0Bngf6xatYpdu3Zx1llnedpcLheff/45U6ZM4ZdffgHMM8FHFsC7du0qcVb4SImJiZ4P5OO1xx/lylQJR5nLeLR2f693tHabzea3PSYmxm97bGys36uoxMXF+Z1XozFpTEdrL/WYfviBlJtv5uwffzQbHngAxowhLjHR7wdQuRhTJL5PFo0J/H+ulucxReL7pDFpTFaPyfdzojRj8qfcL4K76KKL+O6771i7dq3nq02bNvTs2ZO1a9fSsGFDatSowZJ/5huCualybm4u7du3t7DnIhHOMGDyZDj7bGJ+/JGC5GTsCxbAhAkQ4AeUiIhIKJT7M8Cpqak0b97cqy05OZmqVat62gcOHMjo0aNp0qQJTZo0YfTo0SQlJdGjRw8ruhzx4uPj6dmz51F/+5MosGuXuY/vRx8BYFx+ObuHD6de27YWd0zCiT4rxJcyIb5ClYlyXwAHYvDgwRw6dIh+/fqxb98+2rVrx+LFi0lNTbW6axEpJiaGxo0bW90NscqiRdCnD+zcaZ7pffppbP3700AXtRAf+qwQX8qE+ApVJsr9FAh/li5dyqRJkzzHNpuN7Oxs/vrrLw4fPkxubm6Js8YSPEVFRYwZM6bEZHyJcEVF5vzeyy4zi9/mzWHlSrj3XorsdmVCStBnhfhSJsRXqDIRFWeApez5rs6UCPfjj9CjB6xbZx7fey+MHw8VK3ruokyIP8qF+FImxFcoMhGRZ4BFpIwYBrzwgnk543XrID0dFiwwF78dUfyKiIiEE50BFpETs2cP3HEHvP++edy1K8yYAUdsNygiIhKOyv2FMMpKfn4+aWlpx91YWcDtdrNnzx7S09OJidEfGSLSJ5/ArbfC9u0QHw9jx8LAgXCU91uZEH+UC/GlTIiv0mYi0HpN6ZKgs9lspKWlYdOq/7BUWFiIzWbDZrNRWFhYugfb7fDww9Cli1n8nn46fP01DBp01OIXlAnxT7kQX8qE+ApVJlQAS9DZ7XbGjh2rhQyR5tdfoX17c3GbYcBdd8GqVZCVddyHKhPij3IhvpQJ8RWqTKgAFpFjMwx49VWz0F21CqpUgXffNRe/JSVZ3TsREZFS0yI4kQh1tOkNR7YfbwpEst1unun9z3/MhgsugFmzoHbtoPVTRESkrKkAFolQKSkpx71PRkbGUW87D/i8Th3Ytg3i4mDUKHjoIYiNDWIvRUREyp52gQiQdoEInGEY2O12EhIStJDBQif6vY8DhgNDgViAxo1hzhxo2/aE+6JMiD/KhfhSJsRXaTOhXSDEMoZhkJeXh363slZBQYHfr507d3rus3PnTq/bCr/7jsNt2/IY/xS/t90Ga9acVPELyoT4p1yIL2VCfIUqEyqAJegcDgfTpk3D4XBY3ZWolpycfNQvv/eZN4+k9u2JXbkS0tLgrbfMxW8BTKU4HmVC/FEuxJcyIb5ClQnNARaJdnl50L8/5OSYxx07wuzZUK+etf0SEREJEZ0BFoliMV9/Da1amcVvbCw8/jh89pmKXxERiWg6AywhkZCQYHUX5BhiMBe5VejaFVwuqF/fXOh27rkhe01lQvxRLsSXMiG+QpEJ7QIRIO0CIRFj61bo1QuWLTOPe/aE55835/2KiIiUY9oFQizjdrvZsGEDbrfb6q6Ir7fegpYtzeI3NRVef92c7xvi4leZEH+UC/GlTIivUGVCBbAEncPhICcnR6t4w8mBA+aWZjfdZC56a9cO1q41zwSXAWVC/FEuxJcyIb5ClQkVwCKRbuVKaN0aZswAmw0efRS++AIaNrS6ZyIiIpbQIjiRSOVywVNPwWOPgdMJdeqY0x3OP9/qnomIiFhKBbAEnc1mo1q1arqMpZX++ANuvdXc0gzghhvgxRfhlFMs6Y4yIf4oF+JLmRBfocqEdoEIkHaBkHLj3Xfhjjtg3z5ITobJk6FPH3P6g4iISATTLhBiGZfLxerVq3G5XFZ3JboUFsJdd8F115nF71lnwerV5uI3i4tfZUL8US7ElzIhvkKVCRXAEnROp5MFCxbgdDqt7kr0WLsW2rSB6dPNYnfwYFi+HE491eqeAcqE+KdciC9lQnyFKhMqgEWOobCwEJvNhs1mo7Cw0OrulOR2w4QJcPbZ8PPPkJkJS5bAuHGgqymJiIj4pUVwIuXVX3+Zc3sXLzaPr7oKXn4Z0tMt7ZaIiEi40xlgCTqbzUajRo20ijeUPvgAWrQwi9+KFWHaNJg3L2yLX2VC/FEuxJcyIb5ClQntAhEg7QIRnQoLC0lJSQGgoKCA5ORkazt06JA5v3fKFPO4ZUuYMwfOOMPafomIiIQB7QIhlnE6nSxdulSLGILt++/Nub7Fxe8DD8DXX5eL4leZEH+UC/GlTIivUGVCBbAEncvlIjc3V9vYBIthmEVvmzZmEVy9OixcaC5+S0y0uncBUSbEH+VCfCkT4itUmdAiOJFwtmsX3H47fPiheXzZZfDaa5CRYW2/REREyjGdARYJVx9/bC50+/BD80zvc8+Z/1bxKyIiclJ0BliCLiYmhqysLGJi9PvVCSkqgqFDzSkOYM7xfeMNsxgup5QJ8Ue5EF/KhPgKVSa0C0SAtAtEdCrzXSB++gl69DCv7AbQvz889ZS51ZmIiIgck3aBEMs4HA7mz5+Pw+Gwuivlh2HAiy/CWWeZxW/VqjB/vrn4LQKKX2VC/FEuxJcyIb5ClQkVwBJ0brebNWvW4Ha7re5K+bB3L1x7Ldx9t7nP78UXw/r10K2b1T0LGmVC/FEuxJcyIb5ClQkVwCJW+vRTc27ve+9BfDw8/bS5+C0z0+qeiYiIRCwtgpNyqbCwsMxfJ6ivabcTP2oU8RMnYjMM3E2aUPTaa7hbtTLPAoP1V50TERGJUCqAJehiY2Pp1KkTsbGxIXuN4oVpZSkjSNuPNQbmAG3/OZ4OPPDbbxzs2NHrfpG0PrUsMiHlj3IhvpQJ8RWqTGgXiABpF4jwYrPZrO7CCekDTAZSgL+BvsC7R7mv/tcUEREpnUDrNZ0BlqCz2+3MnTuXG2+8kYSEhJC8RkFBQUie11dhYaHnzO/OnTtPfFrCvn0kDhhA3Ltmues67zwqvPQSs2rXZlawOhvGyiITUv4oF+JLmRBfocqECmAJOsMw2LhxY0jPYFoxPzY5OfnEXveLL6BXL/j9d4iLg8cfJ3bwYJKi6E98ZZEJKX+UC/GlTIivUGVCBbBIqDid8Pjj8OST4HZDo0YwZw6cfbbVPRMREYlqKoBFQmHzZujZE776yjzu3RsmT4bUVGv7JSIiItoHWIIvLi6Obt26ERcXpb9f5eRAy5Zm8VupErzxBsyYEdXFb9RnQvxSLsSXMiG+QpUJ7QIRIO0CEZ0KCws9W64VFBQcew5wfj707w+zZ5vHHTqY/65fP/QdFRERkYDrNZ0BlqCz2+1MnToVu91udVfKzooV0KqVWfDGxEB2NixdquL3H1GZCTku5UJ8KRPiK1SZ0N8YJOgMw2D37t3RsYrX5YIxY8yC1+WCevXMKRAdOljds7ASVZmQgCkX4kuZEF+hyoQKYJET9fvv5vZmX3xhHt90E0ybBpUrW9otEREROTZNgRA5EXPnmgvdvvgCUlJg5kxzizMVvyIiImFPZ4Al6OLj4+nZsyfx8fFWdyX4CgpgwAB49VXz+OyzzcK3USNr+xXmIjoTcsKUC/GlTIivUGVCu0AESLtARKcjd4E4+PnnVPzXv+C338Bmg0ceMef+6oNaREQkLGgXCLFMUVERY8aMoaioyOquBIUN+DdQ4cILzeK3dm349FPzCm8qfgMSaZmQ4FAuxJcyIb5ClQlNgZCQiJQtbGzbt7MYuBjMSxtfdx1Mnw5Vqljcs/InUjIhwaVciC9lQnyFIhM6AyxyNO+9R8VzzuFioBAomjIF/vMfFb8iIiLlnM4Ai/g6eBAefBBeeAEbQFYWyW+8AaedZnXPREREJAi0CC5AWgQXOLfbzZ49e0hPTycmppz9kWHtWujRA376yTx+6CEYNQoSEy3tVnlXrjMhIaNciC9lQnyVNhOB1ms6AyxBZ7PZSEtLw2azWd2VwLnd8OyzMGQI2O1Qs6a5t2+XLlb3LCKUy0xIyCkX4kuZEF+hyoR+vZKgs9vtjB07tvwsZNixAy6/HAYNMovfbt1g3ToVv0FU7jIhZUK5EF/KhPgKVSZUAEt0++gjaNECPv4YKlSA55+H99+HatWs7pmIiIiEiKZASHQ6fBgGD4bJk83jM8+EN96AZs2s7ZeIiIiEnM4AS/T5/nvzEsbFxe/998M336j4FRERiRLaBSJA2gUicIZhYLfbSUhICK+FDIYBU6eaOzscPgzVq8Nrr5nzfyWkwjYTYinlQnwpE+KrtJnQpZDFMoZhkJeXR1j9brV7N3TvDvfeaxa/l14K69er+C0jYZkJsZxyIb6UCfEVqkyoAJagczgcTJs2DYfDYXVXTEuWmAvdPvgAEhJg4kT48EPIyLC6Z1Ej7DIhYUG5EF/KhPgKVSa0CE4iV1ERDBsGzzxjHjdtai50a9nS2n6JiIiIpVQAS2T6+Wfzim5r1pjH99wDTz8NSUnW9ktEREQspykQEhIJCQnWvLBhwEsvQevWZvFbtSq89565+E3Fr6Usy4SENeVCfCkT4isUmdAuEAHSLhDlwN690LcvzJtnHl90EcyaBZmZ1vZLREREyoR2gRDLuN1uNmzYgNvtLrsX/ewzc27vvHkQFwfjx8PixSp+w4QlmZCwp1yIL2VCfIUqEyqAJegcDgc5OTlls4rX4YChQ82zvX/+CU2awIoV8O9/Q4ziHS7KNBNSbigX4kuZEF+hyoQWwUn5tWGDudBt5Urz+F//gkmTICXF0m6JiIhIeNMpMil/DANmzoSsLLP4rVwZ/vMfePllFb8iIiJyXDoDLEFns9moVq1aaC5juX+/uaXZm2+ax+efD7NnQ506wX8tCZqQZkLKLeVCfCkT4itUmdAuEAHSLhBh4MsvoWdP2LoVYmPh8cfh4YfNf4uIiEjU0y4QYhmXy8Xq1atxuVzBeUKnE7KzzbO9W7dCw4ZmMTx0qIrfciLomZCIoFyIL2VCfIUqEyqAJeicTicLFizA6XSe/JNt3gydOsHIkeB2wy23mBe4aNfu5J9bykxQMyERQ7kQX8qE+ApVJlQAS/iaMwdatYLly6FSJcjJMS9soSkoIiIichJUAEv4yc+HW2815/vm5+M65xwa5Odj69mTwsJCq3snIiIi5ZwKYAk6m81Go0aNTmzF5tdfm9ubvf66eSGLESM4vGgRW4LeSylLJ5UJiVjKhfhSJsRXqDKhXSACpF0gQszlgrFjYcQI89/16pnbm3XsSGFhISn/7O9bUFBAcnKyxZ0VERGRcKRdIMQyTqeTpUuXBj5hfds2uPBCePRRs/j9v/+DtWuhY8eQ9lPKTqkzIVFBuRBfyoT4ClUmyn0BPGbMGNq2bUtqairVq1fn6quv5pdffvG6j2EYZGdnk5mZScWKFencuTM//PCDRT2OfC6Xi9zc3MC2LHn7bWjRAj7/HJKTYcYMeOMN8+puEjFKlQmJGsqF+FImxFeoMlHuC+Dc3Fz69+/PihUrWLJkCU6nk65du3otlho/fjwTJkxgypQprFy5kho1atClSxcOHDhgYc+jXGEh3HEH3HCDeXW3tm1h7VoKr7+ewoMHKSws9Pr638MK/X6JiIiIBKrcXwp50aJFXsevvfYa1atXZ9WqVZx//vkYhsGkSZMYNmwY1157LQAzZ84kIyODOXPmcNddd1nR7ei2ahX06AG//go2GwwZYu7zGx9PynEmuWdkZPht11R2ERERCVS5L4B95eXlAVClShUANm/ezI4dO+jatavnPomJiXTq1Inly5cftQAuKiqiqKjIc5yfn1+iPSYmhvj4eBwOB26323Pf2NhY4uLisNvtXoVZXFwcsbGxJdrj4+OJiYnxer3idpvNht1u92pPSEjAMAwcDodXe2JiIm6326vdZrORkJCAy+Xymj9T3O50Or3+rBCMMTkcDlq1alVyTG438ZMnY3v0UWwOB0atWjhefRWjUycS4uIwjni90ioqKgrpmCDy3qeyHJPD4aBFixY4HI6IGVMkvk9WjKk4F5E0pkh8n8pqTMWfFTExMREzpiPbNabSj+nInx+BjMn3e3A0EVUAG4bBoEGD6NixI82bNwdgx44dQMkzhxkZGWzduvWozzVmzBhGjhxZon3ChAlUqFABgKysLLp3787ChQtZs2aN5z6dOnWic+fOzJ07l40bN3rau3XrRuvWrXn55ZfZvXu3p71nz540btyYCRMmeL2h99xzD2lpaYwdO9arD0OGDCEvL49p06Z52hISEnjkkUfYtGkTOTk5nvZq1arRr18/1q1bx4IFCzztjRo1olevXixbtozc3FxPezDHFB8fz5gxY7Db7aTk53PNvHk03LwZgB+bNmVBt24c/uor+Oorz5iGDh3qNaaHHnqI7777jnPPPReAhx56iMzMTO68807Wrl3LRx99BMDYsWPLZEyR+D6V5ZjWr18fcWOCyHufynJMP//8M+vXr4+oMUXi+1TWY4qPj2fp0qURNaZIfJ/KYkyTJ0/Gbrd7PieON6aJEycSiIjaBq1///58+OGHLFu2jNq1awOwfPlyOnTowPbt26lZs6bnvn379mXbtm0lplAU83cGuE6dOuzatcuzrYZ+czv6GeBPPvmEyy+/HLfbTcwHHxB3113Y9u7FSEqCiROx33qrOf0hgDEdOHDA8z3fu3cvKSkp+g27nI3p0KFDLFmyhC5dupCQkBARY4rE96msx+Ryufjggw/o0qUL8fHxETGmSHyfyvoM8JIlS7jyyiux2WwRMaYj2yPlfSrLMRUUFHh+fsTHxx93THv37qV69erH3QYtYs4A33fffcyfP5/PP//cU/wC1KhRAzDPBB9ZAO/ateuo80nBDEFiYmJA7cUf3L4SEhJK1e7v9Y7WbrPZ/LbHxMT4bY+NjSU2NrZEe1xcHHFxJWNwsmNau3Ytl55/PonDhkHxb5hZWdjmzIHTT8ffSI81pmKJiYme1yrrMR3Zh0Dbw/19KhbqMcXHx7N+/Xouv/xyTx/K+5gi8X0q6zE5HA5PLnz7VF7HBJH3PkHZjunITETKmIpF0vtULNRjOvLnx5H9Ks2Y/Cn3u0AYhsG9997Lu+++y6effkqDBg28bm/QoAE1atRgyZIlnja73U5ubi7t27cv6+5Gjeo7dhDfocP/it8HH4SvvoLTT7e2YyIiIhL1yv0Z4P79+zNnzhzef/99UlNTPXN+09LSqFixIjabjYEDBzJ69GiaNGlCkyZNGD16NElJSfTo0cPi3kcgwyB2yhT6vvQSMS4X1KgBs2ZBly5W90xEREQEiIACuHgCd+fOnb3aX3vtNfr06QPA4MGDOXToEP369WPfvn20a9eOxYsXk5qaWsa9jXA7d0KfPsT9M6/afcUVxLz2GlSrZnHHxGqxsbF06tTJ75/PJHopF+JLmRBfocpERC2CC6VAry0dtRYuhD59YNcuqFABnn4a+vXzWuh2ogoLC0lJSQGgoKCA5OTkk35OERERiTyB1mvlfg6wWOzwYRg4EC6/3Cx+mzfHsXw5s9PSsPusLj1RycnJGIaBYRgqfsspu93O7NmzS6zaleimXIgvZUJ8hSoTKoDlxP3wA7RrB88+ax7ffz+sXIn7jDPYuHGjrs4mHoZhKBNSgnIhvpQJ8RWqTKgAltIzDHN3hzZtYP16c47vhx+ahfA/FwkREZHwN2nSJK699loaNGiAzWYrsZ7mWLZs2YLNZvP79eabb3rd94033uD8888nIyODxMREMjMz6datG8uXLy/xvBMnTqRChQolnvPuu+/2ut+nn37K7bffzumnn05ycjK1atXiqquuYtWqVSf0vShLs2bN4qabbuK0004jJiaG+vXrB/zYGTNmHPX7brPZSlwg4kiPPvooNpvNc7GwIw0bNoysrCyqVKlChQoVaNiwIXfeeaffi4Zt2LCBW265hbp161KxYkUaNWrEoEGD2Lt3b8DjsFq5XwQnZWzPHvjXv2D+fPP4kktgxgxztwcRESlXXnjhBZKTk7nwwgu9rvJVGvfdd1+JXZWaNGnidbx37146dOjAgAEDSE9P56+//mLChAmcf/75fPLJJ3Tq1Mnr/ueeey4TJkzwavPdu3/atGns3buXAQMGcMYZZ7B7926eeeYZzjnnHD7++GMuvPDCExpPWXj99dfZsWMHZ599dokLTxzPFVdcwVdffVWiffjw4SxZsoRrrrnG7+PWrl3L008/fdRrIOzfv5+bb76Zpk2bkpqayo8//sioUaOYP38+P/zwA1WrVgVg9+7dnHPOOVSqVIknnniCunXrsmbNGkaMGMFnn33GqlWrvPbvD1uGBCQvL88AjLy8PKu7Yp0lSwyjZk3DAMNISDCMiRMNw+UqcTen02msWrXKcDqdZd9HCUvKhPijXFjPdcRneLNmzYxOnToF/NjNmzcbgPHUU0+d0Gvv37/fiI+PN2655RZPm9PpNGrWrGlcfvnlx338zp07S7QdOHDAyMjIMC666KIT6lO9evWMESNGnNBjS+PI7/sVV1xh1KtX76Ser6CgwEhJSTE6duzo93aHw2G0atXKuP/++41OnToZzZo1C+h5P/roIwMwXnnlFU/bSy+9ZADGf//7X6/7jh492gCM1atXn/hA/Cjt50Sg9Vo5KNHFcnY7DB5s7uX711/QtCl8/bW5+M3Pb3mxsbG0bt1a29iIhzIh/igXR/fzzz9z8803e6YM1K1bl1tvvbXEpW5PlpVn6lJTU6lQoYLXVcRiY2NJSEjAFsAOQtWrVy/RlpKSwhlnnMG2bduC2tdgC/b3/a233qKgoIA77rjD7+1jx47l77//5sknnyzV81b7ZxvTI9+j4iu8paWled23cuXKAFQI8lTIUH1OqACWY/vlFzj3XHjqKfP47rvh22+hVaujPsRutzN16lSt4hUPZUL8US78W7duHW3btmXFihU8/vjjLFy4kDFjxlBUVOT1vapfv36p5o6GytixY0lISCApKYmOHTsyv3iKnB8ulwuHw8GWLVu45557MAyD/v37e2632+0cOHCAzz//nNTUVOLj4znjjDN45plncLlcx+1LXl4eq1evplmzZkEZW3nxyiuvUKlSJW644YYStxVPZZg2bZpnS9FjcTqdHDp0iDVr1jBw4EBOPfVUrr32Ws/tV199NXXr1uXBBx/khx9+oKCggM8//5yxY8fSrVs3mjZtGtSxhepzQnOAxT/DgFdegQED4OBBqFLFPL766gAearB7926t4hUPZUL8US78GzRoEHFxcXzzzTeeM3AAPXv29LrfkWflrJCYmEjfvn3p0qULNWvW5Pfff2fy5MlcddVVvPTSS37PRjZr1oxffvkFgJo1a7Jo0SLOOussz+2GYdCoUSPuuOMOTj/9dPbt28d//vMfHnroIdauXcvrr79+zD7179+fwsJChg0bdtz+G4bht6h2u904nU6vNqu/18fy888/s3z5cu666y6SkpK8bnO73dx+++1ce+21XH755cd9rh07dlCzZk3Pcbt27fjss8+8Cue0tDRWrFjBdddd57WY7oYbbjju+3MiQvU5Eb7vqFjn77/hzjvhnXfM4wsvNC9nXKuWtf0SEYlwBw8eJDc3l3/9619exa8/GzZsCOg5fYu52NjYgKYYHE/NmjWZPn26V9sNN9xAu3btGDJkCH369ClROL7zzjsUFhby+++/88ILL3DZZZcxf/58r90nrrjiCnr37k1iYiIAV111FaeccgpTpkxh0KBBZGVl+e3PY489Rk5ODpMnT/Yqqo9m5syZ3HbbbSXan3jiCZ544gmvtuLiy1/RbHVx/MorrwD4/YVjwoQJ/Pbbb8c8K3+k9PR0Vq5cSVFRET/99BPjx4/nggsuYOnSpZ7CeN++fVx11VUcPHiQnJwc6tSpw/fff88TTzxB9+7d+fDDDy3/ngRCUyDEW24utGxpFr9xcTB+PCxZouJXRKQM7Nu3D5fLRe3atYP2nPHx8V5fM2fODNpz+3ut//u//2Pv3r389ttvJW5v1qwZZ599Ntdffz2LFi2iXr16DBgw4LjP26tXLwBWrFjh9/aRI0cyatQonnzySe69996A+tqtWzdWrlzp9VWzZk369u1bor3YzJkzS3w/reRwOJg1axYtW7akTZs2Xrf9/vvvDB8+nBEjRpCQkMD+/fvZv38/TqcTt9vN/v37OXTokNdj4uLiaNOmDR06dOCOO+7g008/ZdOmTV5bq40bN461a9eyZMkSevTowXnnncc999xDTk4OixcvJicnp0zGfrLCv0SXgJ3UJYMdDsjOhjFjzOkPTZrAnDnmXr+lFB8fT8+ePS3/YJDwoUyIP8pFSVWqVCE2NpY//vgjaM95ZAEH0KBBg6A9tz/FZ0uPt9ArLi6O1q1bM3fuXE/b0TJxrOccOXIk2dnZZGdnM3To0ID7WbVqVc/WXsUSEhLIzMwsUUwWKy6aw8UHH3zArl27eOyxx0rctmnTJg4dOsSAAQP8/pJxyimnMGDAACZNmnTU569duzaZmZn8+uuvnra1a9dSq1Ytr6kSAG3btgXg+++/P8HR+BeqzwkVwAIbN0KPHvDNN+bx7bebF7UIYLK8PzExMTRu3DiIHZTyTpkQf5SLkipWrEinTp34z3/+w5NPPkl6evpJP+fRirlQcDgcvPXWW6Snpx/3vT18+DArVqzwut/RMjFr1iwAzjnnHK/2J554guzsbB599FFGjBgRhBEcm7+i2UqvvPIKFSpUKDE/HKBVq1Z89tlnJdoHDhxIXl4er7322nH/0rBhwwb++OMPunfv7mnLzMzkk08+4c8//6TWEX8dLt6bOJh/vYDQfU6oAI5mhgGzZ0O/flBQAJUrw/Tp4GcVaWkUFRUxYcIEBg0a5JnDJdFNmRB/lAv/JkyYQMeOHT1zaRs3bszOnTuZP38+L774IqmpqQCeoiDQucD+fPvtt2zZsgWA/Px8DMPg7bffBswzevXq1QPMAvT222/n1Vdf5dZbbwXMxXoOh4MOHTpQo0YNtm3bxuTJk1m7di2vvfaa17ZV7du3p3v37jRt2pS0tDS2bNnCtGnT2LhxI/PmzfPcb+bMmTz11FPcf//9NG7cmP379/Of//yHN998kz59+tCyZUvPfZ955hmGDx/OpZdeyhVXXFFieoRvsRxOfvzxR3788UfAXHh28OBBz/f9jDPO4IwzzgAgNzeXiy66iOHDhzN8+HCv59i+fTuLFi3i//7v/zjllFNKvEblypX9XtmvcuXKOJ1Or9vWr1/PAw88wPXXX0/Dhg2JiYnhu+++Y+LEiVStWpWHHnrIc9/+/fuTk5NDly5dGDJkiGcO8KhRo8jIyPBbjJ+MUH1OqACOVnl5cM898MYb5vH558Prr0PdukF5em1rJL6UCfFHuSipZcuWfPPNN4wYMYJHHnmEAwcOUKNGDS688EISEhI89/Nd3HYipkyZUmJOcPFWWq+99hp9+vQBzN0EXC4Xbrfbc7/mzZvz4osvMmfOHPLz80lNTeXss8/m448/pmvXrl7P2b59e9588022bNlCYWEh6enpnHvuuUycOJH27dt77tegQQMOHjzIiBEj2Lt3L/Hx8TRr1oypU6dy1113eT1n8ZXrFi1axKJFi0qMLZx3F5k7dy4jR470aiv+vo8YMYLs7Gzgf4vujvy+F5sxYwYul+uoe/+WRkZGBpmZmTzzzDP89ddfOJ1OateuzZVXXsnQoUOpU6eO575nnXUWK1as4IknnmDYsGHs3r2bWrVq0b17d4YPHx6Uv1r4CsXnhM0I54SEkfz8fNLS0sjLy6NSpUpWd8evgOcAL18OPXvCli0QGwsjR8KQIea/g6CoqIixY8cyZMgQndURQJkQ/5QL8aVMiK/SZiLQek1ngKOJ0wlPPgmPPw5uNzRoYC50C+M/E4mIiIgEm84AB6jcnwHesgV69YIvvzSPb7kFpkyBEIzF7XazZ88e0tPTLb3MpoQPZUL8US7ElzIhvkqbiUDrNaUrGrz5prm375dfQmqqufBt1qyQFL8ANpuNtLS0oGy0LpFBmRB/lAvxpUyIr1BlQgVwGCssLCz1l9djd+zA0asX3Hwz5OfjateOg199ReHVVx/zsSfLbrczduxYLW4RD2VC/FEuxJcyIb5ClQkVwGEsJSWlVF8ZGRmex16ZkcFfNWsSn5ODCxgJJH79NcnNmx/18SIiEl0++OADbr31Vs4880zi4+OPeZbN4XAwcuRI6tevT2JiIqeffjqTJ08O6HU+/fRTbr/9dk4//XSSk5OpVasWV111FatWrTrm4wzD4Pzzz8dms5W4wlthYSE33XQTp512GqmpqSQnJ9OsWTNGjRoV1JM6obJ69WouvvhiUlJSqFy5Mtdeey2bNm0K6LHDhg0jKyuLKlWqUKFCBRo2bMidd97J1q1bS9z30Ucf5corr6RWrVrYbDbPzh7+5OTkkJWVRYUKFUhPT6dHjx5s27atxP3uuOMOmjdvTuXKlalYsSKnnnoq//73v9mzZ0/A47eaCuAIEwMMAb4EGgNbgc5ANuA66qNERCQazZs3jxUrVnDGGWd47bHrT79+/RgzZgz9+/fn448/5pprrmHAgAGMHj36uK8zbdo0tmzZwoABA/joo4949tln2bVrF+eccw6ffvrpUR/3/PPPH3WfY4fDgWEYDBo0iHfeeYf333+f6667jscff5yrrrrquH2y0s8//0znzp2x2+3MnTuXV199lV9//ZXzzjuP3bt3H/fx+/fv5+abb2bmzJksWrSIhx56iA8++IB27dqxd+9er/tOnDiRvXv30r17d69t9HxNnjyZXr160aZNG95//33GjRvH0qVLOe+889i3b5/XfQsLC7nzzjuZM2cOH374IXfccQfTp0+nU6dO5efsvSEBycvLMwAjLy+vzF6zoKCgVF+71683PjUvb2EYYDiuvdYo+OOPgB8fLIcPHzays7ONw4cPB+05pXxTJsQf5cJ6LpfL8+/+/fsbRysLvv/+e8NmsxmjR4/2au/bt69RsWJFY+/evcd8nZ07d5ZoO3DggJGRkWFcdNFFnrYjM7F582YjJSXFePfddw3A6N+/f0BjGjx4sAEYGzduDOj+R+rdu7fRqVOnUj+utG644QYjPT3dq6bYsmWLER8fbwwePPiEnvOjjz4yAOOVV17xaj/yPU5OTjZ69+5d4rGHDx820tLSjG7dunm1L1++3ACMoUOHHvf1p06dagDGJ598ckL9P5rSfk4EWq/pDHAYS05OLtVXxcqVaQIUAEUvvEDc22+TXKtWwI8PloSEBIYMGXLM3zQluigT4k+05yI7Oxubzcb69eu54YYbSEtLo0qVKgwaNAin08kvv/zCpZdeSmpqKvXr12f8+PFB70OgOy289957GIbBbbfd5tV+2223cejQIb8XojhS9erVS7SlpKRwxhlneP2J/chM3HnnnXTp0oVrrrkmoD4Wq1atGgBxceG506vT6eSDDz7guuuu89qloF69elxwwQVeV8YrjaONO5D3+PvvvycvL4/LL7/cq/3cc8+lSpUqvPPOOyf8+icrVJ8TKoAjSZUqXA9kAc5evcCiVbSGYZCXlxfWV+GRsqVMiD/KhenGG2+kZcuWvPPOO/Tt25eJEyfywAMPcPXVV3PFFVcwb948LrzwQh5++GHeffddr8d27ty5THZM+P7776lWrRo1atTwam/RooXn9tLKy8tj9erVNGvWzNNWnImXXnqJb775hilTphz3eQzDwOl0kp+fz6JFi3jmmWe4+eabqRukK5sG28aNGzl06JDne3ekFi1asGHDBg4fPhzQczmdTg4dOsSaNWsYOHAgp556Ktdee22p+1Q8bcHfhSYSExP57bff/PbJ6XRSWFjIl19+yWOPPUbHjh3p0KFDqV//WEL1OaECOMJ8DZz4VeGDw+FwMG3aNBwOh8U9kXChTIg/yoXpzjvv5NFHH+Xiiy9m3LhxtGrViilTpjB69Gjuu+8+Lr74YqZPn061atXIycnxemxsbCyxQbqK57Hs3buXKlWqlGhPTk4mISGhxLzTQPTv35/CwkKGDRvmaXM4HIwbN47Bgwczfvx4MjMzj/s8b731FvHx8aSlpXHZZZdx2WWXMWvWrID64HQ6vb4Mw/AU1L7twVL8vfL3/axSpQqGYZSYc+vPjh07iI+PJykpidatW+N0Ovnss89OaFH7aaedRkxMDF8WXyvgHxs3buSvv/7C7XaX6NOKFSuIj48nJSWFjh070rBhQz766KOg5zFUnxPh+fcBERGRKHHllVd6HTdt2pR169Zx2WWXedri4uJo3LhxiVX+n3zySZn0ETjmmebSnoV+7LHHyMnJYfLkyZx11llet33wwQeceeaZ9O3bN6DnuuSSS1i5ciUHDhzgq6++Yty4cezdu5d58+Yd88//W7ZsoUGDBn5vi4+P9zr+7LPP6Ny5M2AWzUeKjY09obPwJ/v9TE9PZ+XKlRQVFfHTTz8xfvx4LrjgApYuXUrNmjVL1ZcqVarQs2dPZs2aRdu2bbnhhhv4448/uPPOO4mNjcXlcpX4Xp555pmsXLmSgwcPsnbtWsaOHUuXLl349NNPSUpKKtXrW0EFsIiIiIV8zwQmJCSQlJREhQoVSrTn5+eXZdc8qlatytq1a0u0FxYWYrfb/Z7NPJqRI0cyatQonnzyyRJbm7377rts2LCBl19+mby8PK/b7HY7+/fvJzk52atAPeWUU2jTpg0AF1xwAY0aNeKmm27i/fffP+b84czMTFauXFmib9u3b+fFF1/0aj/ttNMA/0XzkcVxIKpWrQrg96z533//jc1mo3Llysd9nri4OM+4O3TowKWXXkqDBg0YO3Yszz77bMD9KTZt2jQMw6Bfv37cfffdxMTEcMstt5CRkcHHH3/s6Xex5ORkz+uff/75tGvXjnPOOYcXX3yRBx54oNSvX9ZUAEtIROuiFjk6ZUL8US7KhzPPPJM333yTHTt2eM0D/u677wBo3rx5QM8zcuRIsrOzyc7OZujQoSVu/+GHH3C73Zx//vklbnvppZd46aWXmDdvHldfffVRX+Pss88G4Ndffz1mXxISEjwFXLGqVaty4MCBEu3F/BXNxcVxoBo1akTFihU937sjfffddzRu3LjELz+BqF27NpmZmccd99EkJyfz+uuv89xzz7Ft2zYyMzNJT0/n9NNPp3379sdd3NamTRtiYmJO+PWPJRSfE5oDLEGXmJjII4884ncyvUQnZUL8US7Kj6uuugqbzcbMmTO92mfMmEHFihW59NJLj/scTzzxBNnZ2Tz66KOMGDHC733uuOMOPvvssxJfAFdffTWfffYZHTt2PObrFN+/cePGgQytVIqL5iO/UlNTS/UccXFxdOvWjXfffZcDBw542n///Xc+++yzE1rEBrBhwwb++OOPkx73KaecQosWLUhPT2f+/Pn88ssvDBgw4LiPy83Nxe12B/37HqrPCZ0BlqBzu91s2rSJhg0bBrzFjkQ2ZUL8US5O3kUXXURubm6JeamB2rp1q+eM5saNGwF4++23Aahfv77nTGizZs3417/+xYgRI4iNjaVt27YsXryY6dOnM2rUKK8pEI8//jiPP/44n3zyCZ06dQLgmWeeYfjw4Vx66aVcccUVrFixwqsf55xzDgB169bF6XT6zUStWrW8phq8+OKLfPHFF3Tt2pU6depQWFjIF198weTJk2nfvn1YXwxj5MiRtG3bliuvvJIhQ4Zw+PBhhg8fTnp6Og8++KDXfePi4ujUqZNnvvf69et54IEHuP766z3fp++++46JEydStWpVHnroIa/H5+bmei6u4XK52Lp1q+c97tSpk2f7snfeeYft27fTtGlTDh8+zNKlS3n22We5++67vb6XH3zwAS+99BLdu3enXr16OBwOvv32WyZNmkTjxo254447gvq9CtnnxIlsShyNrLgQRnmlze3FlzIh/kR7LkaMGGEAxu7du73ae/fubSQnJ5e4f6dOnYxmzZqVaDuZH+WvvfaaAfj98r1ggt1uN0aMGGHUrVvXSEhIME499VTjueeeO+q4PvvssxL9PNpXsaNlAj8Xwvjyyy+NK6+80sjMzDQSEhKMpKQko2XLlsYTTzxhFBYWntD3o6wuhGEYhvHtt98aF110kZGUlGRUqlTJuPrqq40NGzaUuB/g1acdO3YYvXr1Mho1amQkJSUZCQkJRsOGDY27777b+P3330s8/ljf+yPfo3nz5hmtWrUykpOTjYoVKxpt2rQxXnnlFcPtdns9308//WRcf/31Rr169YwKFSoYFSpUME4//XTj3//+93EviHIiQnUhDJthRPkGjAHKz88nLS2NvLw8r42rpaSioiLGjh3LkCFD9KdNAZQJ8U+5EF/KhPgqbSYCrdf0NycRERERiSoqgCXobDYb1apVK5OrE0n5oEyIP8qF+FImxFeoMqEpEAHSFAgRERGR8KYpEGIZl8vF6tWrcblcVndFwoQyIf4oF+JLmRBfocqECmAJOqfTyYIFC054Wx6JPMqE+KNciC9lQnyFKhMqgEVEREQkqqgAFhEREZGoogJYgs5ms9GoUSOt4hUPZUL8US7ElzIhvkKVCe0CESDtAiEiIiIS3rQLhFjG6XSydOlSLWIQD2VC/FEuxJcyIb5ClQkVwBJ0LpeL3NxcbWMjHsqE+KNciC9lQnyFKhMqgEVEREQkqqgAFhEREZGoogJYgi4mJoasrCxiYhQvMSkT4o9yIb6UCfEVqkxoF4gAaRcIERERkfCmXSDEMg6Hg/nz5+NwOKzuioQJZUL8US7ElzIhvkKVCRXAEnRut5s1a9bgdrut7oqECWVC/FEuxJcyIb5ClQkVwCIiIiISVeKs7kB5UTxVOj8/3+KehL+ioiIOHz5Mfn4+iYmJVndHwoAyIf4oF+JLmRBfpc1EcZ12vCVuWgQXoD/++IM6depY3Q0REREROY5t27ZRu3bto96uAjhAbreb7du3k5qais1ms7o7YS0/P586deqwbds27ZghgDIh/ikX4kuZEF+lzYRhGBw4cIDMzMxjbp2mKRABiomJOeZvElJSpUqV9AEmXpQJ8Ue5EF/KhPgqTSbS0tKOex8tghMRERGRqKICWERERESiigpgCbrExERGjBihFbzioUyIP8qF+FImxFeoMqFFcCIiIiISVXQGWERERESiigpgEREREYkqKoBFREREJKqoABYRERGRqKICWIJmzJgxtG3bltTUVKpXr87VV1/NL7/8YnW3JIyMGTMGm83GwIEDre6KWOjPP/+kV69eVK1alaSkJFq1asWqVaus7pZYxOl08uijj9KgQQMqVqxIw4YNefzxx3G73VZ3TcrI559/Trdu3cjMzMRms/Hee+953W4YBtnZ2WRmZlKxYkU6d+7MDz/8cFKvqQJYgiY3N5f+/fuzYsUKlixZgtPppGvXrhQWFlrdNQkDK1euZPr06bRo0cLqroiF9u3bR4cOHYiPj2fhwoX8+OOPPPPMM1SuXNnqrolFxo0bxwsvvMCUKVP46aefGD9+PE899RSTJ0+2umtSRgoLC2nZsiVTpkzxe/v48eOZMGECU6ZMYeXKldSoUYMuXbpw4MCBE35NbYMmIbN7926qV69Obm4u559/vtXdEQsVFBTQunVrpk6dyqhRo2jVqhWTJk2yultigSFDhvDll1/yxRdfWN0VCRNXXnklGRkZvPLKK5626667jqSkJF5//XULeyZWsNlszJs3j6uvvhowz/5mZmYycOBAHn74YQCKiorIyMhg3Lhx3HXXXSf0OjoDLCGTl5cHQJUqVSzuiVitf//+XHHFFVx88cVWd0UsNn/+fNq0acMNN9xA9erVycrK4qWXXrK6W2Khjh078sknn/Drr78CsG7dOpYtW8bll19ucc8kHGzevJkdO3bQtWtXT1tiYiKdOnVi+fLlJ/y8ccHonIgvwzAYNGgQHTt2pHnz5lZ3Ryz05ptvsnr1alauXGl1VyQMbNq0iWnTpjFo0CCGDh3KN998w/33309iYiK33nqr1d0TCzz88MPk5eVx+umnExsbi8vl4sknn+Tmm2+2umsSBnbs2AFARkaGV3tGRgZbt2494edVASwhce+997J+/XqWLVtmdVfEQtu2bWPAgAEsXryYChUqWN0dCQNut5s2bdowevRoALKysvjhhx+YNm2aCuAo9dZbbzF79mzmzJlDs2bNWLt2LQMHDiQzM5PevXtb3T0JEzabzevYMIwSbaWhAliC7r777mP+/Pl8/vnn1K5d2+ruiIVWrVrFrl27OOusszxtLpeLzz//nClTplBUVERsbKyFPZSyVrNmTc444wyvtqZNm/LOO+9Y1COx2r///W+GDBnCTTfdBMCZZ57J1q1bGTNmjApgoUaNGoB5JrhmzZqe9l27dpU4K1wamgMsQWMYBvfeey/vvvsun376KQ0aNLC6S2Kxiy66iO+++461a9d6vtq0aUPPnj1Zu3atit8o1KFDhxLbI/7666/Uq1fPoh6J1Q4ePEhMjHc5Ehsbq23QBIAGDRpQo0YNlixZ4mmz2+3k5ubSvn37E35enQGWoOnfvz9z5szh/fffJzU11TNvJy0tjYoVK1rcO7FCampqiTngycnJVK1aVXPDo9QDDzxA+/btGT16NDfeeCPffPMN06dPZ/r06VZ3TSzSrVs3nnzySerWrUuzZs1Ys2YNEyZM4Pbbb7e6a1JGCgoK2LBhg+d48+bNrF27lipVqlC3bl0GDhzI6NGjadKkCU2aNGH06NEkJSXRo0ePE35NbYMmQXO0uTivvfYaffr0KdvOSNjq3LmztkGLch988AGPPPIIv/32Gw0aNGDQoEH07dvX6m6JRQ4cOMBjjz3GvHnz2LVrF5mZmdx8880MHz6chIQEq7snZWDp0qVccMEFJdp79+7NjBkzMAyDkSNH8uKLL7Jv3z7atWvH888/f1InUlQAi4iIiEhU0RxgEREREYkqKoBFREREJKqoABYRERGRqKICWERERESiigpgEREREYkqKoBFREREJKqoABYRERGRqKICWERERESiigpgEREREYkqKoBFREREJKqoABYRERGRqKICWERERESiigpgEREREYkqKoBFRKLEHXfcgc1mo0uXLhiGUeL24cOHY7PZOPPMMykqKrKghyIiZcNm+PsUFBGRiFPw/+3dr0udbQDG8QtPk4NzdUX/gCMnGAQHhoFitiiiHCwGg+gQQYyLwkw2EWHBaLB52poW4eAWFIZBBH+AOMYQxfGGF1befDzw3p9PvZ9wxS8P98Pz61fq9Xp+/PiRzc3NLC0t/T07OjrK+/fvU6lUcnx8nHq93rmhAG3mDTBAIarVar58+ZJKpZK1tbV8+/YtSfL79+/Mzs7m5eUlnz59Er/A/54ABijI8PBwVldX8/j4mJmZmTw9PeXjx485Pz/PyMhIVlZWOj0RoO1cgQAozPPzc4aGhnJycpLR0dE0m8309PSk1Wqlr6+v0/MA2k4AAxTo+/fvGRwczOPjY5Jkd3c3jUajw6sAXocABijQ09NTBgYGcnZ2ljdv3uTy8jLVarXTswBehTvAAAVaX1/P2dlZurq68vDwkOXl5U5PAng1AhigMF+/fs3nz5/T3d2dZrOZ3t7ebG9v5+DgoNPTAF6FAAYoyM+fP9NoNPLnz59sbGzkw4cP2draSvLvjzJub287vBCg/QQwQEEWFxdzcXGRsbGxLCwsJEmmp6czOTmZm5ubzM/Pd3ghQPv5CA6gEPv7+5mYmMjbt29zenqad+/e/T27v79PrVbL1dVVdnZ2Mjc318GlAO0lgAEKcH19nVqtlru7u+zt7WVqauo/zxweHmZ8fDzVajWtViv9/f2vPxTgFQhgAACK4g4wAABFEcAAABRFAAMAUBQBDABAUQQwAABFEcAAABRFAAMAUBQBDABAUQQwAABFEcAAABRFAAMAUBQBDABAUQQwAABF+QcmMxuMJQ9WVgAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Calculate fitted function values\n",
"yfit = straight_line(final_params, xdata)\n",
"\n",
"fig = plt.figure(figsize = (8, 6))\n",
"plt.title('Fit result')\n",
"plt.xlabel('x', fontsize=16)\n",
"plt.ylabel('y', fontsize=16)\n",
"plt.grid(color = 'grey', linestyle=\"--\")\n",
"\n",
"# Plot data (no line connecting the points!)\n",
"plt.errorbar(xdata, ydata, xerr = xerror, yerr = yerror, fmt='k', linestyle = '', label = \"Data\") \n",
"\n",
"# Plot fit \n",
"plt.plot(xdata, yfit, color = 'r', linestyle = '-', label = \"Fit\")\n",
"\n",
"# Add legend and fit params\n",
"plt.legend(loc = 2, fontsize=16) \n",
"\n",
"text = \"c: {:7.5g} +- {:7.5g}\\n\".format(final_params[0], param_errors[0])\n",
"text += \"m: {:7.5g} +- {:7.5g}\\n\".format(final_params[1], param_errors[1])\n",
"plt.text(0.95, 0.0, text, transform = fig.axes[0].transAxes, ha = \"right\", va = \"bottom\", fontsize=12) \n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The data look to be reasonably well described by a straight line. The value of the statistic $\\chi^2$/NDF and the $\\chi^2$ probability backs this up. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Lecture 6 formative exercise 2\n",
"Fit a straight line to the data in Table 2 (which is the same as Table 1 with two new points added). Do this by first adding the new data points to the arrays, then copying the relevant fitting cells above and pasting them below this cell (use the _Edit_ menu). Unlike above, where we have done the fit in several cells so that we can explain what is happening as we go, bring the code together in a **single** cell. \n",
"\n",
"\n",
"\n",
"|Measurement | $x$ value | Error in $x$ | $y$ value | Error in $y$ |\n",
"|------------|-----------|--------------|-----------|--------------|\n",
"| 1 | 1.50 | 0.21 | 14.3 | 2.1 |\n",
"| 2 | 2.31 | 0.11 | 20.2 | 1.7 |\n",
"| 3 | 2.78 | 0.43 | 30.1 | 3.3 |\n",
"| 4 | 3.58 | 0.13 | 36.5 | 1.1 |\n",
"| 5 | 4.08 | 0.17 | 42.7 | 0.9 |\n",
"| 6 | 4.76 | 0.18 | 47.1 | 1.1 |\n",
"| 7 | 5.62 | 0.15 | 52.9 | 1.5 |\n",
"| 8 | 7.02 | 0.19 | 68.8 | 0.9 |\n",
"| 9 | 8.45 | 0.17 | 85.2 | 1.2 |\n",
"| 10 | 9.65 | 0.11 | 99.4 | 2.9 |\n",
"| 11 | 10.25 | 0.21 | 110.5 | 3.1 |\n",
"| 12 | 11.85 | 0.41 | 122.1 | 3.3 |\n",
"\n",
"**Table 2** *Extended range of data points*\n",
"\n",
"Notes: \n",
"* To extend the data you can either make a new csv file, and read it in like above, or you can use the [np.append](https://numpy.org/doc/stable/reference/generated/numpy.append.html) method. \n",
"* You don't need to copy the definitions of `straight_line`, `straight_line_diff` and `minimise` since they are already defined once you have run the corresponding cells and do not change."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Lecture 6 formative exercise 3\n",
"\n",
"Does the straight line still describe the data? Have the values of the gradient or the intercept changed significantly? Be quantitative in your answer.\n",
"\n",
"Note: two values can be said to be consistent if the difference between them is less than 3 times their combined error (this is called the \"consistency check\" in PHYS106)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Since you will use this code many times, it is useful to package it up into a function (see below) that we can call, which simply takes in the $x$ and $y$ data and their errors, along with the initial parameters. From now on whenever you need to do a fit you can just copy this function and call it with the correct arguments). "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# \n",
"\n",
"def fit(xdata, ydata, xerror, yerror, init_params, \n",
" xlabel = \"x\", ylabel = \"y\", title = \"Fit result\"):\n",
" \"\"\"\n",
" Function to perform least-squares fit for a straight line.\n",
" - xdata array of x data points\n",
" - ydata array of y data points\n",
" - xerr array of x data errors\n",
" - yerr array of y data errors\n",
" - init_params array of function parameters\n",
" - xlabel x-axis label\n",
" - ylabel y-axis label\n",
" - title plot title\n",
" \"\"\" \n",
"\n",
" # Run fit\n",
" result = least_squares(minimise, init_params, args=(xdata, ydata, xerror, yerror))\n",
"\n",
" # Check fit succeeds \n",
" if not result.success or result.status < 1:\n",
" print (\"ERROR: Fit failed with message {}\".format(result.message))\n",
" print (\"Please check the data and inital parameter estimates\")\n",
" return 0,0\n",
" else:\n",
" print (\"Fit succeeded\")\n",
"\n",
" # Get fitted parameters \n",
" final_params = result.x\n",
" c = final_params[0]\n",
" m = final_params[1]\n",
" nparams = len(final_params)\n",
"\n",
" # Calculate chi2\n",
" chi2_array = result.fun ** 2\n",
" chi2 = sum(chi2_array)\n",
" npoints = len(xdata)\n",
" reduced_chi2 = chi2 / (npoints - nparams) \n",
" chi2_prob = stats_chi2.sf(chi2, (npoints - nparams))\n",
"\n",
" # Print chi2\n",
" np.set_printoptions(precision = 3)\n",
" print(\"\\n=== Fit quality ===\")\n",
" print(\"chisq per point = \\n\",chi2_array)\n",
" print(\"chisq = {:7.5g}, ndf = {}, chisq/NDF = {:7.5g}, chisq prob = {:7.5g}\\n\".format(chi2, npoints-nparams, reduced_chi2, chi2_prob))\n",
"\n",
" if reduced_chi2 < 0.25 or reduced_chi2 > 4:\n",
" print(\"WARNING: chi2/ndf suspiciously small or large. Please check the data and initial parameter estimates\")\n",
"\n",
" if chi2_prob < 0.05:\n",
" print(\"WARNING: chi2 probability for given degrees of freedom less than 0.05 . Please check the data and initial parameter estimates\") \n",
" \n",
" # Calculate errors\n",
" jacobian = result.jac\n",
" jacobian2 = np.dot(jacobian.T, jacobian)\n",
" determinant = np.linalg.det(jacobian2)\n",
" \n",
" if determinant < 1E-32:\n",
" print(f\"Matrix singular (determinant = {determinant}, error calculation failed.\")\n",
" param_errors = np.zeros(nparams)\n",
" else:\n",
" covariance = np.linalg.inv(jacobian2)\n",
" param_errors = np.sqrt(covariance.diagonal())\n",
"\n",
" print (\"=== Fitted parameters ===\")\n",
" print (\"c = {:7.5g} +- {:7.5g}\".format(final_params[0], param_errors[0]))\n",
" print (\"m = {:7.5g} +- {:7.5g}\".format(final_params[1], param_errors[1]))\n",
"\n",
" # Calculate fitted function values\n",
" yfit = straight_line(final_params, xdata)\n",
"\n",
" # Visualise result\n",
" fig = plt.figure(figsize = (8, 6))\n",
" plt.title(title)\n",
" plt.xlabel(xlabel, fontsize=16)\n",
" plt.ylabel(ylabel, fontsize=16)\n",
" plt.grid(color = 'grey', linestyle=\"--\")\n",
"\n",
" plt.errorbar(xdata, ydata, xerr = xerror, yerr = yerror, fmt='k', linestyle = '', label = \"Data\") \n",
" plt.plot(xdata, yfit, color = 'r', linestyle = '-', label = \"Fit\")\n",
"\n",
" plt.legend(loc = 2, fontsize=16) \n",
"\n",
" text = \"c: {:7.5g} +- {:7.5g}\\n\".format(final_params[0], param_errors[0])\n",
" text += \"m: {:7.5g} +- {:7.5g}\\n\".format(final_params[1], param_errors[1])\n",
" plt.text(0.95, 0.0, text, transform = fig.axes[0].transAxes, ha = \"right\", va = \"bottom\", fontsize=12) \n",
"\n",
" plt.show()\n",
" \n",
" # Return arrays of parameters and associated errors\n",
" return final_params, param_errors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Lecture 6 summative exercise\n",
"\n",
"So far the data we have been using is imaginary, as we concentrated on the tools to do the fitting. In this exercise you will use them on real data from an experiment on [Boyle's law](https://en.wikipedia.org/wiki/Boyle%27s_law). The data, consisting of a series of measurements of pressure ($P$) [in $Pa = N/m^2$] and volume ($V$) [in $m^3$ ] and their uncertainties ($\\sigma_P$ and $\\sigma_V$) taken at fixed temperature, can be found in the file `boyle.csv` (in the order $V,\\sigma_V, P, \\sigma_P$). Read this data and use the `fit` function to show that it conforms to Boyle's law i.e. \"PV = k\", finding the best-fit proportionality constraint, $k$. Comment on the goodness of fit.\n",
"\n",
"Notes: \n",
"* Think about how to make this into something linear that is easy to fit and don't forget to propagate the errors correctly following what you have leant in PHYS106 and/or the formulae [here](https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulae). All steps should be done in python.\n",
"* You will likely need to plot the data first to find the initial parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Summary \n",
"\n",
"That's the end of Week 3. You should now know to use the `least_squares` function from SciPy, along with NumPy arrays, to fit experimental data, visualising the results using Matplotlib. You will use this technique and the code we have developed in PHYS106 and throughout your labs in later years. We concentrated on simple examples of fitting straight lines, since this is what you will mostly encounter in PHYS106 and many problems can be made linear by a change of variables (as in Ex 4). However, the code we have developed is much more powerful: it can be used to fit any function as long as we can write down the functional form and its derivative. We will come back to fitting later in the course.\n",
"\n",
"Now that our programs are getting more complex, next week we will take a step back and look at good coding practice and debugging errors, before going on to look at more real physics use cases."
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.11.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}