{
"nbformat": 4,
"nbformat_minor": 2,
"metadata": {
"colab": {
"name": "unscented_kalman_filter.ipynb",
"provenance": [],
"collapsed_sections": [],
"toc_visible": true,
"authorship_tag": "ABX9TyNOTGT/taLUTiLelilODZx6",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"source": [
"# Unscented Kalman Filter\n",
"
"
],
"metadata": {
"id": "view-in-github",
"colab_type": "text"
}
},
{
"cell_type": "code",
"execution_count": 1,
"source": [
"import math\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.spatial.transform import Rotation as Rot\n",
"import scipy.linalg"
],
"outputs": [],
"metadata": {
"id": "4RB4SSKW3Q4h"
}
},
{
"cell_type": "markdown",
"source": [
"# Python Robotics\n",
"## States and Control inputs"
],
"metadata": {
"id": "QDt-4Y4q_FhK"
}
},
{
"cell_type": "markdown",
"source": [
"状态量:x坐标、y坐标、偏航角、速度\n",
"> $x$, $y$ are a 2D x-y position, $\\phi$ is orientation, and $v$ is velocity.\n",
"\n",
"$\\textbf{x}_t=[x_t, y_t, \\phi_t, v_t]$\n"
],
"metadata": {
"id": "J0_Gd8Ra4bn9"
}
},
{
"cell_type": "markdown",
"source": [
"控制量:速度、角速度\n",
"> The robot has a speed sensor and a gyro sensor.\n",
"\n",
"$\\textbf{u}_t=[v_t, \\omega_t]$\n"
],
"metadata": {
"id": "Y2qMbqQk4DIn"
}
},
{
"cell_type": "code",
"execution_count": 2,
"source": [
"def calc_input():\n",
" v = 1.0 # [m/s]\n",
" yawRate = 0.1 # [rad/s]\n",
" u = np.array([[v, yawRate]]).T\n",
" return u\n",
"\n",
"u = calc_input()\n",
"u"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1. ],\n",
" [0.1]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 2
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "H8WMFBq63YQu",
"outputId": "be801821-e26b-4f44-c4c8-355a2326515f"
}
},
{
"cell_type": "markdown",
"source": [
"## Motion Model\n",
"The robot model is \n",
"\n",
"$ \\dot{x} = vcos(\\phi)$ \n",
"$ \\dot{y} = vsin((\\phi)$ \n",
"$ \\dot{\\phi} = \\omega$\n",
"\n",
"Given that $\\textbf{x}_t = [x_t, y_t, \\phi_t, v_t]$,$\\textbf{u}_t = [v_t, \\dot{\\phi}]$\n",
"\n",
"the motion model is\n",
"\n",
"$\\textbf{x}_{t+1} = f(x_t, u_t)$ \n",
"$\\qquad = F\\textbf{x}_t+B\\textbf{u}_t$\n",
"\n",
" \n",
"$x_{t+1} = x_t + \\dot{x}\\Delta t$ \n",
"$y_{t+1} = y_t + \\dot{y}\\Delta t$ \n",
"$\\phi_{t+1} = \\phi_t + \\dot{\\phi}\\Delta t$ \n",
"$v_{t+1} = 0 + v_{t+1}$\n",
"\n",
"where\n",
"\n",
"$\\begin{equation*}\n",
"F=\n",
"\\begin{bmatrix}\n",
"1 & 0 & 0 & 0\\\\\n",
"0 & 1 & 0 & 0\\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"0 & 0 & 0 & 0 \\\\\n",
"\\end{bmatrix}\n",
"\\end{equation*}$\n",
"\n",
"$\\begin{equation*}\n",
"B=\n",
"\\begin{bmatrix}\n",
"cos(\\phi)\\Delta t & 0\\\\\n",
"sin(\\phi)\\Delta t & 0\\\\\n",
"0 & \\Delta t\\\\\n",
"1 & 0\\\\\n",
"\\end{bmatrix}\n",
"\\end{equation*}$\n"
],
"metadata": {
"id": "aJfZydyk-4v9"
}
},
{
"cell_type": "code",
"execution_count": 3,
"source": [
"def motion_model(x, u, dt):\n",
" F = np.array([[1.0, 0, 0, 0],\n",
" [0, 1.0, 0, 0],\n",
" [0, 0, 1.0, 0],\n",
" [0, 0, 0, 0]])\n",
"\n",
" B = np.array([[dt * math.cos(x[2]), 0],\n",
" [dt * math.sin(x[2]), 0],\n",
" [0.0, dt],\n",
" [1.0, 0.0]])\n",
" x = F@x + B@u\n",
" return x"
],
"outputs": [],
"metadata": {
"id": "vzHV_JDtAeAH"
}
},
{
"cell_type": "markdown",
"source": [
"## Observation Model\n",
"观测量:x坐标、y坐标\n",
"> The robot has a GNSS sensor.\n",
"\n",
"$\\textbf{z}_t = [x_t, y_t]$\n",
"\n",
"Given that $\\textbf{x}_t = [x_t, y_t, \\phi_t, v_t]$, so the observation model is\n",
"\n",
"$\\textbf{z}_{t} = H\\textbf{x}_t$\n",
"\n",
"where\n",
"\n",
"$\\begin{equation*}\n",
"H=\n",
"\\begin{bmatrix}\n",
"1 & 0 & 0& 0\\\\\n",
"0 & 1 & 0& 0\\\\\n",
"\\end{bmatrix}\n",
"\\end{equation*}$\n"
],
"metadata": {
"id": "v-5e9mAL50Z3"
}
},
{
"cell_type": "code",
"execution_count": 4,
"source": [
"def observation_model(x):\n",
" H = np.array([[1, 0, 0, 0],\n",
" [0, 1, 0, 0]])\n",
" z = H @ x\n",
" return z\n",
"\n",
"x = np.ones((4, 1))\n",
"observation_model(x)"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1.],\n",
" [1.]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 4
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "H0yKhtrz8ksD",
"outputId": "64b15681-eaf3-47e0-eced-14724792bc4a"
}
},
{
"cell_type": "markdown",
"source": [
"给观测量和控制量加上噪音:\n",
"\n",
"$\\textbf{u}_d = \\textbf{u} + \\epsilon$\n",
"\n",
"$\\textbf{z} = H\\textbf{x} + \\delta$"
],
"metadata": {
"id": "gCRc-dLw7wXc"
}
},
{
"cell_type": "code",
"execution_count": 5,
"source": [
"# Simulation parameter\n",
"INPUT_NOISE = np.diag([1.0, np.deg2rad(30.0)]) ** 2\n",
"GPS_NOISE = np.diag([0.5, 0.5]) ** 2\n",
"\n",
"def observation(xTrue, u, dt):\n",
" xTrue = motion_model(xTrue, u, dt)\n",
" # add noise to gps x-y\n",
" z = observation_model(xTrue) + GPS_NOISE @ np.random.randn(2, 1)\n",
" # add noise to input\n",
" ud = u + INPUT_NOISE @ np.random.randn(2, 1)\n",
" return xTrue, z, ud\n",
"\n",
"xTrue = np.zeros((4, 1))\n",
"u = calc_input()\n",
"xTrue, z, ud = observation(xTrue, u, dt=1)\n",
"print('xTrue:\\n', xTrue)\n",
"print('\\nz:\\n', z)\n",
"print('\\nu:\\n', u)\n",
"print('\\nu with noise:\\n', ud)"
],
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"xTrue:\n",
" [[1. ]\n",
" [0. ]\n",
" [0.1]\n",
" [1. ]]\n",
"\n",
"z:\n",
" [[ 1.22539845]\n",
" [-0.11166166]]\n",
"\n",
"u:\n",
" [[1. ]\n",
" [0.1]]\n",
"\n",
"u with noise:\n",
" [[0.3668484 ]\n",
" [0.23745866]]\n"
]
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "sKFRBwjN4CO1",
"outputId": "2b803d58-b51e-4bca-f532-1263cf7833dd"
}
},
{
"cell_type": "markdown",
"source": [
"## Sigma points\n",
"In this case, the functions $f$ and $h$ are not linear. UKF applies UT to performs the approximation with extracted points called ”Sigma points. Sigma points are computed using the matrix $S$ which is defined from the covariance matrix, $\\Sigma_t$. $S_i$ denotes the i-th colum of the matrix $S$\n",
"\n",
"$$S = \\sqrt{\\Sigma_t}$$\n",
"\n",
"We define the sigma points $X_{i,t} \\in X_t$ as follows:\n",
"\n",
"$$\n",
"X_{i, t}=\\left\\{\\begin{array}{ll}\n",
"=\\mu_{t}, & i=0 \\\\\n",
"=\\mu_{t}+\\gamma S_{i}, & i=1, \\ldots, N \\\\\n",
"=\\mu_{t}-\\gamma S_{i-N}, & i=N+1, \\ldots, 2 N\n",
"\\end{array}\\right.\n",
"$$\n",
"\n",
"where $N$ is the dimension of the state, and $\\gamma$ is computed by:\n",
"\n",
"$$\\gamma = \\sqrt{N+\\lambda}$$\n",
"\n",
"$$\\lambda = \\alpha^{2}(N+\\kappa) - N$$\n",
"\n",
"$\\alpha, \\kappa$ are two hyper-parameters (scaling parameters) of UKF.\n"
],
"metadata": {
"id": "kOS42u906buc"
}
},
{
"cell_type": "code",
"execution_count": 6,
"source": [
"def generate_sigma_points(xEst, PEst, gamma):\n",
" sigma = xEst\n",
" Psqrt = scipy.linalg.sqrtm(PEst) # PEst = Psqrt @ Psqrt\n",
" n = len(xEst[:, 0])\n",
" # Positive direction\n",
" for i in range(n):\n",
" sigma = np.hstack((sigma, xEst + gamma*Psqrt[:, i:i+1]))\n",
"\n",
" # Negative direction\n",
" for i in range(n):\n",
" sigma = np.hstack((sigma, xEst - gamma*Psqrt[:, i:i+1]))\n",
" return sigma\n",
"\n",
"N = 4\n",
"ALPHA = 0.001\n",
"KAPPA = 0\n",
"\n",
"lambda_ = ALPHA ** 2 * (N + KAPPA) - N\n",
"gamma = math.sqrt(N + lambda_)\n",
"\n",
"xEst = np.zeros((N, 1))\n",
"xTrue = np.zeros((N, 1))\n",
"PEst = np.eye(N)\n",
"sigma_points = generate_sigma_points(xEst, PEst, gamma)\n",
"sigma_points.shape # (N, 2N+1)"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(4, 9)"
]
},
"metadata": {
"tags": []
},
"execution_count": 6
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "wCr__qnPOS2v",
"outputId": "ff6a1514-cdde-46cd-a616-13c589ae2e3e"
}
},
{
"cell_type": "markdown",
"source": [
"## Predication\n",
"And then, the sigma points are passed through the nonlinear function $f$:\n",
"\n",
"$$\\bar{X} = f(X)$$"
],
"metadata": {
"id": "hgO7MKpjamf-"
}
},
{
"cell_type": "code",
"execution_count": 7,
"source": [
"def predict_sigma_motion(sigma, u, dt):\n",
" \"\"\"\n",
" Sigma Points prediction with motion model\n",
" \"\"\"\n",
" for i in range(sigma.shape[1]):\n",
" sigma[:, i:i+1] = motion_model(sigma[:, i:i+1], u, dt)\n",
" return sigma\n",
"\n",
"u = calc_input()\n",
"perdicted_points = predict_sigma_motion(sigma_points, u, 1)\n",
"perdicted_points.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(4, 9)"
]
},
"metadata": {
"tags": []
},
"execution_count": 7
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "1XZ3hW4KO0gn",
"outputId": "ac92423d-baeb-4f25-a0f0-9b64a9d6fe75"
}
},
{
"cell_type": "markdown",
"source": [
"Here $\\bar{X}$ are transformed sigma points. Finally, the objective parameters $\\mu'$ and $\\Sigma'$ of resultant Gaussian are calculated using the unscented transform on the tranformed sigma points, $\\mu', \\Sigma^\\ =\\ UT(\\bar{X}, \\omega_m, \\omega_c, Q)$:\n",
"\n",
"$$\\mu^{\\prime} = \\sum_{i=0}^{2n} \\omega_i^m \\bar{X}_i$$\n",
"\n",
"$$\\Sigma^{\\prime} = \\sum_{i=0}^{2n} \\omega_i^c \\left(\\bar{X}_i-\\mu^{\\prime}\\right)\\left(\\bar{X}_i-\\mu^{\\prime}\\right)^{T} + Q$$\n",
"\n",
"where $Q$ is Transition model convariance;\n",
"\n",
"and $\\omega^m$ are weights used when the mean is computed:\n",
"\n",
"$$\n",
"w_i^m=\\left\\{\\begin{array}{ll} \n",
"& =\\frac{\\lambda}{N+\\lambda} & i=0 \\\\\n",
"& =\\frac{1}{2(N+\\lambda)}, & i=1, \\ldots, 2 N \\end{array}\\right.\n",
"$$\n",
"\n",
"and $\\omega^c$ are weights used when the covariance of Gaussian is recovered:\n",
"\n",
"$$\n",
"w^c_i=\\left\\{\\begin{aligned}\n",
"&=\\frac{\\lambda}{N+\\lambda}+\\left(1-\\alpha^{2}+\\beta\\right) & i=0 \\\\\n",
"&=\\frac{1}{2(N+\\lambda)}, & i=1, \\ldots, 2 N\n",
"\\end{aligned}\\right.\n",
"$$"
],
"metadata": {
"id": "ipv0ons12ZQI"
}
},
{
"cell_type": "code",
"execution_count": 8,
"source": [
"# UKF Parameter\n",
"ALPHA = 0.001\n",
"BETA = 2\n",
"KAPPA = 0\n",
"\n",
"# calculate weights\n",
"wm = [lambda_ / (lambda_ + N)]\n",
"wc = [(lambda_ / (lambda_ + N)) + (1 - ALPHA**2 + BETA)]\n",
"for i in range(2*N):\n",
" wm.append(1.0 / (2*(N + lambda_)))\n",
" wc.append(1.0 / (2*(N + lambda_)))\n",
"\n",
"wm = np.array([wm])\n",
"wc = np.array([wc])\n",
"wm.shape, wc.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"((1, 9), (1, 9))"
]
},
"metadata": {
"tags": []
},
"execution_count": 8
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "qhSlK7D-19zu",
"outputId": "b3f9da9c-31ff-43a7-f46f-d0bcae01751f"
}
},
{
"cell_type": "markdown",
"source": [
"Recall that:\n",
"\n",
"$\\mu^{\\prime} = \\sum_{i=0}^{2n} \\omega_{m}^{[i]} \\bar{X}^{[i]}$\n",
"\n",
"$\\Sigma^{\\prime} = \\sum_{i=0}^{2n} \\omega_{c}^{[i]}\\left(\\bar{X}^{[i]}-\\mu^{\\prime}\\right)\\left(\\bar{X}^{[i]}-\\mu^{\\prime}\\right)^{T}$\n"
],
"metadata": {
"id": "yqOVCk-k7j8U"
}
},
{
"cell_type": "code",
"execution_count": 9,
"source": [
"# Covariance for UKF simulation\n",
"Q = np.diag([\n",
" 0.1, # variance of location on x-axis\n",
" 0.1, # variance of location on y-axis\n",
" np.deg2rad(1.0), # variance of yaw angle\n",
" 1.0 # variance of velocity\n",
"]) ** 2 # predict state covariance\n",
"\n",
"def calc_sigma_covariance(x, sigma, wc, Pi):\n",
" nSigma = sigma.shape[1]\n",
" d = sigma - x[0:sigma.shape[0]]\n",
" P = Pi\n",
" for i in range(nSigma):\n",
" P = P + wc[0, i] * d[:, i:i+1] @ d[:, i:i+1].T\n",
" return P\n",
"\n",
"xPred = (wm @ perdicted_points.T).T\n",
"PPred = calc_sigma_covariance(xPred, perdicted_points, wc, Q)\n",
"xPred.shape, PPred.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"((4, 1), (4, 4))"
]
},
"metadata": {
"tags": []
},
"execution_count": 9
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "CszGeICbPftV",
"outputId": "8e763d77-e75c-4c1a-8568-5c3d76b50d91"
}
},
{
"cell_type": "markdown",
"source": [
"## Update\n",
"Kalman filters perform the update in measurement space. Thus we must convert the sigma points of the prior into measurements using the measurement function $h$:\n",
"\n",
"$$\n",
"Z = h(\\bar{X})\n",
"$$\n",
"\n",
"Then we compute the mean and covariance of these points using the unscented transform, $\\mu_z, P_z = UT(Z, w_m, w_c, R)$.\n",
"> The $z$ subscript denotes that these are the mean and covariance of the measurement sigma points.)\n",
"\n",
"$$\\begin{aligned}\n",
"\\mu_z &= \\sum_{i=0}^{2n} w^m_i Z_i \\\\\n",
"P_z &= \\sum_{i=0}^{2n} w^c_i{(Z_i-\\mu_z)(Z_i-\\mu_z)^\\mathsf T} + R\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"where $R$ is Measurement model covariance."
],
"metadata": {
"id": "fVcxOMCk8qoN"
}
},
{
"cell_type": "code",
"execution_count": 10,
"source": [
"# Update\n",
"R = np.diag([1.0, 1.0]) ** 2 # Observation x,y position covariance\n",
"\n",
"zSigmaPoints = observation_model(perdicted_points)\n",
"mu_z = (wm @ zSigmaPoints.T).T\n",
"P_z = calc_sigma_covariance(mu_z, zSigmaPoints, wc, R)\n",
"\n",
"mu_z.shape, P_z.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"((2, 1), (2, 2))"
]
},
"metadata": {
"tags": []
},
"execution_count": 10
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "USc2o11lxOY-",
"outputId": "098d0936-ac4c-4fd5-e942-19db0ffa5445"
}
},
{
"cell_type": "markdown",
"source": [
"Next we compute the residual and Kalman gain. The residual of the measurement $z$ is trivial to compute:\n",
"\n",
"$$y = \\mathbf z_m - \\mu_z$$\n",
"\n",
"To compute the Kalman gain we first compute the [cross covariance](https://en.wikipedia.org/wiki/Cross-covariance) of the state and the measurements, which is defined as: \n",
"\n",
"$$P_{xz} = \\sum_{i=0}^{2n} w^c_i(\\bar{X}_i - \\mu')(Z_i - \\mu_z)^\\mathsf T$$\n",
"\n",
"And then the Kalman gain is defined as\n",
"\n",
"$$K = P_{xz}P^{-1}_z$$\n",
"\n",
"If you think of the inverse as a *kind of* matrix reciprocal, you can see that the Kalman gain is a simple ratio which computes:\n",
"\n",
"$$K \\approx \\frac{P_{xz}}{P_z} \n",
"\\approx \\frac{\\text{belief in state}}{\\text{belief in measurement}}$$\n"
],
"metadata": {
"id": "gvLXSjvE1wn1"
}
},
{
"cell_type": "code",
"execution_count": 11,
"source": [
"def calc_pxz(sigma, x, z_sigma, zb, wc):\n",
" nSigma = sigma.shape[1]\n",
" dx = sigma - x\n",
" dz = z_sigma - zb[0:2]\n",
" P = np.zeros((dx.shape[0], dz.shape[0]))\n",
"\n",
" for i in range(nSigma):\n",
" P = P + wc[0, i] * dx[:, i:i+1] @ dz[:, i:i+1].T\n",
" return P\n",
"\n",
"y = z - mu_z\n",
"Pxz = calc_pxz(perdicted_points, xPred, zSigmaPoints, mu_z, wc)\n",
"K = Pxz @ np.linalg.inv(P_z)\n",
"y.shape, Pxz.shape, K.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"((2, 1), (4, 2), (4, 2))"
]
},
"metadata": {
"tags": []
},
"execution_count": 11
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "lEzxEdnPzBkc",
"outputId": "d48f1760-5141-41f5-d93f-d266931ee639"
}
},
{
"cell_type": "markdown",
"source": [
"Finally, we compute the new state estimate using the residual and Kalman gain:\n",
"\n",
"$$\\mathbf x = \\mu' + Ky$$\n",
"\n",
"and the new covariance is computed as:\n",
"\n",
"$$\n",
"\\mathbf P = \\Sigma' - KP_zK^{\\mathsf T}\n",
"$$"
],
"metadata": {
"id": "PTLt_vJa4HF9"
}
},
{
"cell_type": "code",
"execution_count": 12,
"source": [
"xEst = xPred + K @ y\n",
"PEst = PPred - K @ P_z @ K.T\n",
"\n",
"xEst.shape, PEst.shape"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"((4, 1), (4, 4))"
]
},
"metadata": {
"tags": []
},
"execution_count": 12
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uOkhb0nBDCY3",
"outputId": "20f82800-8936-4b21-c259-84732f5eebe2"
}
},
{
"cell_type": "code",
"execution_count": 13,
"source": [
"xTrue, xEst"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(array([[0.],\n",
" [0.],\n",
" [0.],\n",
" [0.]]), array([[ 0.93523919],\n",
" [-0.07444109],\n",
" [ 0.06277945],\n",
" [ 1. ]]))"
]
},
"metadata": {
"tags": []
},
"execution_count": 13
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "11FVdDRPl-Rx",
"outputId": "577dad73-dd26-4304-dcd4-3fc4f42682e7"
}
},
{
"cell_type": "markdown",
"source": [
"## Sum-up\n",
"This table compares the equations of the linear KF and UKF equations.\n",
"\n",
"$$\\begin{array}{l|l}\n",
"\\textrm{Kalman Filter} & \\textrm{Unscented Kalman Filter} \\\\\n",
"\\hline \n",
"& \\boldsymbol{\\mathcal Y} = f(\\boldsymbol\\chi) \\\\\n",
"\\mathbf{\\bar x} = \\mathbf{Fx} & \n",
"\\mathbf{\\bar x} = \\sum w^m\\boldsymbol{\\mathcal Y} \\\\\n",
"\\mathbf{\\bar P} = \\mathbf{FPF}^\\mathsf T+\\mathbf Q & \n",
"\\mathbf{\\bar P} = \\sum w^c({\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})(\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})^\\mathsf T}+\\mathbf Q \\\\\n",
"\\hline \n",
"& \\boldsymbol{\\mathcal Z} = h(\\boldsymbol{\\mathcal{Y}}) \\\\\n",
"& \\boldsymbol\\mu_z = \\sum w^m\\boldsymbol{\\mathcal{Z}} \\\\\n",
"\\mathbf y = \\mathbf z - \\mathbf{Hx} &\n",
"\\mathbf y = \\mathbf z - \\boldsymbol\\mu_z \\\\\n",
"\\mathbf S = \\mathbf{H\\bar PH}^\\mathsf{T} + \\mathbf R & \n",
"\\mathbf P_z = \\sum w^c{(\\boldsymbol{\\mathcal Z}-\\boldsymbol\\mu_z)(\\boldsymbol{\\mathcal{Z}}-\\boldsymbol\\mu_z)^\\mathsf{T}} + \\mathbf R \\\\ \n",
"\\mathbf K = \\mathbf{\\bar PH}^\\mathsf T \\mathbf S^{-1} &\n",
"\\mathbf K = \\left[\\sum w^c(\\boldsymbol{\\mathcal Y}-\\bar{\\mathbf x})(\\boldsymbol{\\mathcal{Z}}-\\boldsymbol\\mu_z)^\\mathsf{T}\\right] \\mathbf P_z^{-1} \\\\\n",
"\\mathbf x = \\mathbf{\\bar x} + \\mathbf{Ky} & \\mathbf x = \\mathbf{\\bar x} + \\mathbf{Ky}\\\\\n",
"\\mathbf P = (\\mathbf{I}-\\mathbf{KH})\\mathbf{\\bar P} & \\mathbf P = \\bar{\\mathbf P} - \\mathbf{KP_z}\\mathbf{K}^\\mathsf{T}\n",
"\\end{array}$$"
],
"metadata": {
"id": "vREtoA624SXC"
}
},
{
"cell_type": "code",
"execution_count": 14,
"source": [
"def ukf_estimation(xEst, PEst, z, u, wm, wc, gamma, dt):\n",
" # Predict\n",
" sigma = generate_sigma_points(xEst, PEst, gamma)\n",
" perdicted_points = predict_sigma_motion(sigma, u, dt)\n",
" xPred = (wm @ sigma.T).T\n",
" PPred = calc_sigma_covariance(xPred, sigma, wc, Q)\n",
"\n",
" # Update\n",
" zSigmaPoints = observation_model(sigma)\n",
" mu_z = (wm @ zSigmaPoints.T).T\n",
" P_z = calc_sigma_covariance(mu_z, zSigmaPoints, wc, R)\n",
" Pxz = calc_pxz(perdicted_points, xPred, zSigmaPoints, mu_z, wc)\n",
" K = Pxz @ np.linalg.inv(P_z)\n",
" y = z - mu_z\n",
" xEst = xPred + K @ y\n",
" PEst = PPred - K @ P_z @ K.T\n",
"\n",
" return xEst, PEst\n",
"\n",
"def setup_ukf(nx, ALPHA, BETA, KAPPA):\n",
" lamb = ALPHA ** 2 * (nx + KAPPA) - nx\n",
" # calculate weights\n",
" wm = [lamb / (lamb + nx)]\n",
" wc = [(lamb / (lamb + nx)) + (1 - ALPHA ** 2 + BETA)]\n",
" for i in range(2 * nx):\n",
" wm.append(1.0 / (2 * (nx + lamb)))\n",
" wc.append(1.0 / (2 * (nx + lamb)))\n",
" gamma = math.sqrt(nx + lamb)\n",
"\n",
" wm = np.array([wm])\n",
" wc = np.array([wc])\n",
"\n",
" return wm, wc, gamma"
],
"outputs": [],
"metadata": {
"id": "jVuggEWkPuFt"
}
},
{
"cell_type": "code",
"execution_count": 15,
"source": [
"def plot_covariance_ellipse(xEst, PEst): # pragma: no cover\n",
" Pxy = PEst[0:2, 0:2]\n",
" eigval, eigvec = np.linalg.eig(Pxy)\n",
"\n",
" if eigval[0] >= eigval[1]:\n",
" bigind = 0\n",
" smallind = 1\n",
" else:\n",
" bigind = 1\n",
" smallind = 0\n",
"\n",
" t = np.arange(0, 2 * math.pi + 0.1, 0.1)\n",
" a = math.sqrt(eigval[bigind])\n",
" b = math.sqrt(eigval[smallind])\n",
" x = [a * math.cos(it) for it in t]\n",
" y = [b * math.sin(it) for it in t]\n",
" angle = math.atan2(eigvec[1, bigind], eigvec[0, bigind])\n",
" rot = Rot.from_euler('z', angle).as_matrix()[0:2, 0:2]\n",
" fx = rot @ np.array([x, y])\n",
" px = np.array(fx[0, :] + xEst[0, 0]).flatten()\n",
" py = np.array(fx[1, :] + xEst[1, 0]).flatten()\n",
" plt.plot(px, py, \"--r\")"
],
"outputs": [],
"metadata": {
"id": "5l11IMebtopi"
}
},
{
"cell_type": "code",
"execution_count": 16,
"source": [
"DT = 0.1 # time tick [s]\n",
"SIM_TIME = 50.0 # simulation time [s]\n",
"\n",
"# Covariance for UKF simulation\n",
"Q = np.diag([\n",
" 0.1, # variance of location on x-axis\n",
" 0.1, # variance of location on y-axis\n",
" np.deg2rad(1.0), # variance of yaw angle\n",
" 1.0 # variance of velocity\n",
"]) ** 2 # predict state covariance\n",
"R = np.diag([1.0, 1.0]) ** 2 # Observation x,y position covariance\n",
"\n",
"# UKF Parameter\n",
"ALPHA = 0.001\n",
"BETA = 2\n",
"KAPPA = 0\n",
"\n",
"def main():\n",
" nx = 4 # State Vector [x y yaw v]'\n",
" xEst = np.zeros((nx, 1))\n",
" xTrue = np.zeros((nx, 1))\n",
" PEst = np.eye(nx)\n",
" xDR = np.zeros((nx, 1)) # Dead reckoning\n",
"\n",
" wm, wc, gamma = setup_ukf(nx, ALPHA, BETA, KAPPA)\n",
"\n",
" # history\n",
" hxEst = xEst\n",
" hxTrue = xTrue\n",
" hxDR = xTrue\n",
" hz = np.zeros((2, 1))\n",
"\n",
" plt.figure(figsize=(12,9))\n",
" time = 0.0\n",
" while SIM_TIME >= time:\n",
" time += DT\n",
" u = calc_input()\n",
"\n",
" xTrue, z, ud = observation(xTrue, u, DT)\n",
" xDR = motion_model(xDR, ud, DT)\n",
"\n",
" xEst, PEst = ukf_estimation(xEst, PEst, z, ud, wm, wc, gamma, dt=DT)\n",
"\n",
" # store data history\n",
" hxEst = np.hstack((hxEst, xEst))\n",
" hxDR = np.hstack((hxDR, xDR))\n",
" hxTrue = np.hstack((hxTrue, xTrue))\n",
" hz = np.hstack((hz, z))\n",
" if (int(time*10) % 10) == 0:\n",
" plot_covariance_ellipse(xEst, PEst)\n",
"\n",
" plt.plot(hxEst[0, :], hxEst[1, :], color='k', lw=2)\n",
" plt.plot(hxTrue[0, :], hxTrue[1, :], color='b', lw=2)\n",
" plt.plot(hxDR[0, :], hxDR[1, :], color='g', lw=2)\n",
"\n",
" plt.axis('equal')\n",
" plt.title(\"UKF Robot localization\")\n",
" plt.show()\n",
"\n",
"main()"
],
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
""
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAIYCAYAAABkAIS4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1yV5fvA8c/NFEFURBG34MKtOHKbM3NmOcptWq5Ks2/zV5ptS7Ms09xluUfuvffGLW5NVJChIojAuX9/3IiomDjgAOd6v17nJeeZ13MyvM79XM91K601QgghhBBC2AI7awcghBBCCCFEWpHkVwghhBBC2AxJfoUQQgghhM2Q5FcIIYQQQtgMSX6FEEIIIYTNkORXCCGEEELYDEl+hRDiKSmlziqlGqXSsbsrpTanxrGTnGO9UqpXws+dlFIrU+EcHyulJjzr4wohxOOS5FcIkW4ppbRSqth9y4YqpaYl/FxfKfVvknVOSql5SqktSin3hG1jlVKRSV7vP+RcZ5VS0QnbXFZKTVFKuaXuFSZ/jdaktf5La93kaY5x/3+XhON+rbXu9XTRCSHE05PkVwiRKSilnIF5QA6gidb6esKqmVprtySv4f9xmJZaazegIlAJ+Ch1oxZCCJHWJPkVQmR4SqmswCLAAWiutb75NMfTWl8GVmCS4DvnaKWUOqyUikgoE/C7b7eqSqkjSqlwpdRkpVSWJPv2VkqdVEqFKaUWKqXyJSzfmLBJQMKIc4cUXGtNpdQupdS1hD9rJlnnkXDuoIQ4FiQsz6mUWqyUCklYvlgpVeAhx08ss1BKvX/fqHmsUmpKwroeSqmjSqkbSqnTSqk3E5a7AsuAfEn2y5d0xP5Rn2fCKPx7SqkDCdc5M+nnKYQQT0OSXyFERueMSbZuAa211tFPe8CExLAZcDLhfQlgOjAQyA0sBRYppZyS7NYJaAr4AiWA/0vYtwHwDdAe8AbOATMAtNZ1E/atkDAqPfMRcXkAS4CfgVzASGCJUipXwiZ/AlmBMkAe4MeE5XbAZKAwUAiIBn551OegtR5+Z8Qc8ANCgDsxBgMtAHegB/CjUqpywhePZkBQktH2oPuuIyWfZ3vgBaAoUB7o/qh4hRAiJST5FUJkdNmAGsBUrXVMMuvbJ4wu3nnl+49jLVBK3QAuYJK7IQnLOwBLtNartNaxwA+AC1Azyb6/aK0vaK3DgK+AVxOWdwImaa33JsT3EVBDKVXkCa61OXBCa/2n1jpOaz0dOAa0VEp5Y5LOPlrrcK11rNZ6A4DWOlRrPVdrHaW1vpEQX72UnlQp5QIsAH7SWi9LOOYSrfUpbWwAVgJ1UnjIlHyeP2utgxI+z0UkGYUXQoinIcmvECI9iwcc71vmCMQmeX8V6AhMVUo1TeYYs7TWOZK8gpLZ5o42WutsQH2gFOCZsDwfZsQWAK21BZMg50+y74UkP59L2Ce5fSOB0Pv2Tal7jpXkXPmBgkCY1jr8/p2UUlmVUuOUUueUUteBjUAOpZR9Cs87ETiutf4uyTGbKaW2J5RyRAAvcvfzeqzreMjneTnJz1FAqj98KISwDZL8CiHSs/NAkfuWFeW+BFBrPQ/oDcxRSj3/tCdNGMmcghmRBAjClAwAoJRSmGTzYpLdCib5uVDCPsnt64opWUi6b0rdc6wk57qISR49lFI5ktlvMFASqK61dgfulFuoR51QKfUhpozj9STLnIG5mM/HS2udA1O6cOd4+nGu4yGfpxBCpApJfoUQ6dlM4P+UUgWUUnbK9NJtCcy5f8OEEoABwD9KqVrP4NyjgMZKqQrALKC5UqqhUsoRk0zGAFuTbN8/IU4P4BPu1sZOB3oopSomJI1fAzu01mcT1l8BfFIY01KghFLqNaWUQ8IDcqWBxVrrS5ja5zEJD7g5KqXuJLnZMHW+EQnxDUn26PdRSjUD3gZeuq+W2glTax0CxCVsl7Q92hUgl1Iq+0MOnZLPUwghUoUkv0KI9GwYJiHaDIQDw4FOWutDyW2stZ6KSaSWKKWqPc2JtdYhwB/AZ1rr40BnYDSmzKIlpi3a7SS7/I2pez0NnAK+TDjOauBTzEjpJcwDcR2T7DcUU7IRoZRq/4iYQjEPmQ3GlE68D7TQWl9N2KQLpiTkGKZmeWDC8lGYmtqrwHZgeQo/hg6YB9KOJuncMDahbvhtTBIbDrwGLEwS5zFM0n86uTrrFH6eQgiRKpTWj7o7JYQQQgghROYgI79CCCGEEMJmSPIrhBBCCCFshiS/QgghhBDCZkjyK4QQQgghbMYjk1+lVEGl1Dpl5qw/rJR6J2H5UKXURaXU/oTXi6kfrhBCCCGEEE/ukd0eEqbM9NZa71VKZQP2AG0w865Haq1/+M8DJOHp6amLFCnyFOEKIYQQQgjxaHv27Lmqtc59/3KHR+2Y0Dj9UsLPN5RSR3myaTkpUqQIu3fvfpJdhRBCCCGESDGl1P3TwQOPWfOrlCoCVAJ2JCwaoJQ6oJSapJTK+ZB93lBK7VZK7Q4JCXmc0wkhhBBCCPFMpTj5VUq5YWYoGqi1vg78hpmpqCJmZHhEcvtprX/XWlfRWlfJnfuBkWchhBBCCCHSTIqS34S51+cCf2mt5wFora9oreO11hZgPPBUU4kKIYQQQgiR2lLS7UEBE4GjWuuRSZZ7J9nsJeDQsw9PCCGEEEKIZ+eRD7wBtYAuwEGl1P6EZR8DryqlKgIaOAu8mSoRCiGEEEII8YykpNvDZkAls2rpsw9HCCGEEEKI1CMzvAkhhBBCCJshya8QQgghhLAZkvwKIYQQQgibIcmvEEIIIYSwGZL8CiGEEEIImyHJrxBCCCGEsBmS/AohhBBCCJshya8QQgghhLAZkvwKIYQQQgibIcmvEEIIIYSwGZL8CiGEEEIImyHJrxBCCCGEsBkO1g5ACCGE7YiKjSIwNJBzEedo4tuELA5ZuBZzjejYaJRSONg54Oroiouji7VDFUJkUpL8CiGEeGLxlniOhx7ndPhprkReIfhmMFdumj+DbwYTfiucGzE3iLwdyaXIS4917LqF6+Lh4oFHFg/zp4sHOV1yJv7s4eJBziw58XLzIqtj1lS6QiFEZiPJrxBCiBSJiYvhUPAh9l3ex95Le9l3eR8BlwOIjotOlfNtPLfxifarX6Q+/ar0o0LeCvjm9MXezv4ZRyaEyMgk+RVCCJGsyNuRbD6/mWkHpvHXwb8eul3h7IXxy+2Hl6sXXq5e5HHNQx7XPHi5eeHh4kE2p2y4Obnh7uxOvI7H2d4Zezt77JU9dsoOpRRaa2ItsUTFRrHr4i48s3oScSuCsOgwwqLDCL8VnuzPoVGhXLl5hdvxtxPjWX92PevPrgcgq2NWyuYpSwWvCuaVtwLl8pTD3dkdpVRqf4RCiHRIaa3T7GRVqlTRu3fvTrPzCSGESDmLtrD30l4WBy5m9enV7Li4gzhL3D3beLt5U79IfSp7V6ZS3kpU8q6Eh4uHlSI2tNaERocydf9UVp9ZjZ+nH4GhgQRcCeDf6//+575jXhxDQ5+GFPcoLsmwEJmMUmqP1rrKA8sl+RVCCNsVHRvNmjNrWHR8EYsCF91Tl2un7KiSrwoNijTA18OXJr5NKJS9kBWjfXyhUaEcDD5IwOUAAq6Y195Lex/YroB7ARr5NKJ58eY09W1KNudsVohWCPEsSfIrhBACgMuRl1kSuISFgQtZdWrVPTW7BdwL0LJES14o9gL1Ctcje5bsVow0dcRZ4ph1eBaXIy+z4+IO1p5Zy9Woq4nrneydaFi0IR3LdqR1ydaZ8jMQwhZI8iuEEDZKa82h4EMsPL6QRYGL2HFxxz3r/b39aVWyFS1LtKRi3oo2d/vfoi0cuHKAFSdXsChwEVsvbEVj/m10tHMk1hILmBKJvlX7WjNUIcRjkORXCCFsSFRsFJvObWLpiaUsDFzI2Yizieuc7Z1p5NOIliVa0qJEC/K757deoOlQ8M1g5h2dx4xDM9hwbsMD60+9fQqfnD5WiEwI8Tgk+RVCiEzuTPgZ5h+bz7KTy9h0bhMx8TGJ6/K45qFF8Ra0LNmSxj6NcXVytWKkGUfQjSDK/VaOsOiwB9Zd/d9VcmXNZYWohBAp8bDkV1qdCSFEBhYYGsicI3OYe3TuAw9y+Xv708S3Ca1KtqJa/mrYKZnR/nHly5aP0PdDibfEk/v73ITfCk9c5/m9J61LtmZBxwVWjFAI8bgk+RVCiAzmcuRlZhyawbQD09hzaU/icjcnN5oXb07rkq1p5NOI3K65rRhl5mJvZ0/YB6avsOf3nonL/zn+D+evnc9wXTCEsGWS/AohRAYQeTuS+UfnM+3gNFafXo1FWwBwd3bnpVIv8bLfyzT2bUwWhyxWjjRzy5U1F/GfxfPFhi8YumEoAINXDmZK6ylSSiJEBiHJrxBCpFOx8bGsOr2KaQem8c/xf4iKjQJMB4KWJVrSuXxnWpRoIQlvGrNTdgypP4RK3pXoPK8zc47MITA0kAUdFlA0Z1FrhyeEeAR54E0IIdIRrTW7gnYx7cA0ZhyaQUhUSOK62oVq07lcZ9qVaWf1WdWEcezqMVrPaE1gaCAeLh7MfGUmjXwaWTssIQTS7UEIIdI1rTVrzqxh2IZhbDq/KXF5Kc9SdCnfhdfKvUaRHEWsF6B4qGu3rtFpXieWnFiCnbLj+8bfM+i5QTbXL1mI9EaSXyGESIe01qw4tYJhG4ax7d9tAOTMkpPuFbvTuXxnKuWtJElUBmDRFj5b9xlfbfoKgF6VejGm+Rgc7R2tHJkQtktanQkhRDqitWblqZV8uu5TdgXtAiCXSy7eq/ke/av2J5tzNitHKB6HnbLjywZfUsGrAl0XdGXCvgmcjjjNnHZzyOmS09rhCSGSkORXCCHS2KZzm/hk7SeJ5Q15XPPwv5r/o0+VPrg5uVk5OvE02pVpR+EchWk1vRVrz6yl3pR6rOm6RtrOCZGOSPIrhBBp5FDwIcr9Vi7xvYeLBx/U+oAB1QaQ1TGrFSMTz1K1/NXY2XsnTac15WDwQRr80YA1XdeQxzWPtUMTQgAy3Y8QQqSy6zHXeW/le/ckvkPqDeH026d5v9b7kvhmQoWyF2Jdt3X4efpxKPgQDaY2IPhmsLXDEkIgya8QQqSq+UfnU/KXkozYNgKFokWJFpx46wRD6w8le5bs1g5PpKK8bnlZ120dpXOX5nDIYZ6f+jxXIq9YOywhbJ4kv0IIkQquRl3l1bmv0nZWWy5HXqZ6/urs7L2TRa8uophHMWuHJ9KIl5sXa7uupXTu0hwJOcLzU5/ncuRla4clhE2T5FcIIZ6x+UfnU2ZMGWYcmkFWx6z8/MLPbH19K1XyPdBxR9gALzcv1nVbR5ncZTh69agkwEJYmSS/QgjxjATdCKL97Pa0ndWW4JvB1CtcjwN9DvBW9bewU/Lr1pblcc3Dum7rKJunLMeuHsN7hDerTq2ydlhC2CT5bSyEEE8p3hLPLzt/we9XP2YfmU1Wx6yMbjaatd3W4uvha+3wRDqR2zU3a7uuTXzfZFoTzl87b8WIhLBNkvwKIcRT2HdpHzUm1uCtZW9xPeY6LUu05Ei/IwyoNkBGe8UD7k+A602px5nwM1aMSAjbI7+ZhRDiCUTFRjF4xWCqjK/CrqBd5M+Wn3nt5/FPx38onKOwtcMT6djzRZ8n/INwquevztmIs9SbUo+TYSetHZYQNkOSXyGEeExbzm+h4tiKjNw+EoCB1QdytP9RXvJ7CaWUlaMTGUGOLDlY2WUlNQvW5ML1C1QeV5l9l/ZZOywhbIIkv0IIkULRsdEMXjGYOpPrcCLsBGVyl2FHrx38+MKPZHPOZu3wRAbj7uzOis4rqFOoDjdu36Dy75UJDA20dlhCZHqS/AohRAoE3wym7pS6jNw+Ejtlx8e1P2bPG3ukfZl4Km5ObizrtCzxfdNpTbl4/aIVIxIi85PkVwghHqHr/K54/eDF7qDdFM1RlO29tvNVw69wdnC2dmgiE3B1cuX6h9eplr8aZyPO0mRaE0KjQq0dlhCZliS/QgjxH+pMrsOfB/5MfL/t9W0y2iueuWzO2VjWaRllcpfhSMgRmkxrQsStCGuHJUSmJMmvEEI8RL0p9dh8fnPi+4gPIvBy87JiRCIz83DxYGWXlfjm9GXvpb28MO0Frsdct3ZYQmQ6kvwKIUQyvt70NRvPbUx8f+uTW2TPkt2KEQlbkC9bPtZ2W0uRHEXYcXEHL/71IpG3I60dlhCZiiS/QghxH79f/fhk7ScAFHAvQPxn8VLfK9JMoeyFWNt1LQXcC7DlwhZaTm9JVGyUtcMSItOQ5FcIIZI4fvU4x64eS3x/fuB5malNpLmiOYuyrts6vN28WX92Pa1ntObm7ZvWDkuITEF+owshRILrMddpM7MNAPbKnvjP4mXSCmE1xTyKsbbbWrxcvVh9ejX1ptQj6EaQtcMSIsOT5FcIIQCLttB5XmeOXT1GmdxlCP8gXEZ8hdWV8izF+u7r8cnpw55Le6g2vhp7L+21dlhCZGjym10IIYDP13/OosBF5MiSg386/iMztol0o5RnKXb02kGdQnW4eOMidSbXYf7R+dYOS4gMS5JfIYTNm3d0HsM2DsNO2THzlZn4evhaOyQh7uGZ1ZNVXVbRrUI3omKjaDurLd9t/g6ttbVDEyLDkeRXCGHTDgUfouv8rgB81+g7mvg2sXJEQiTP2cGZya0n823DbwH4cM2H9FzYk9vxt60cmRAZiyS/QgibFR4dTpsZbbgZe5PXyr3G4BqDrR2SEP9JKcUHtT9gbvu5uDi4MGX/FBr/2ZirUVetHZoQGYYkv0IIm6S1pufCnpwKP0WlvJUY33K8dHYQGUZbv7Zs6rGJfNnysfHcRp6b8Nw9LfqEEA8nya8Qwib9vONnFhxbQHbn7MxpP4esjlmtHZIQj8U/nz87e+2ksndlToWf4rkJz7H69GprhyVEuifJrxDC5uy8uJP/rfofAJNaT8Inp4+VIxLiyeR3z8/G7htp69eWazHXeGHaC4zdPdbaYQmRrknyK4SwKeHR4XSY04FYSyxvV3ubtn5trR2SEE/F1cmV2e1m81Htj4jX8fRd0peBywcSb4m3dmhCpEuS/AohbMadOt+zEWepmq8qwxsPt3ZIQjwTdsqOrxt+zZTWU3C0c+SnHT/RakYrrsdct3ZoQqQ7kvwKIWzGTzt+SqzznfnKTJwdnK0dkhDPVLeK3VjTdQ25XHKx9MRSak2qxdmIs9YOS4h0RZJfIYRN2HlxJ++veh+Aya0nUzRnUStHJETqqFO4Djt67aCUZykOBR+i+oTqbLuwzdphCZFuSPIrhMj0wqPDaT+7PbGWWN6p/g4v+b1k7ZCESFW+Hr5se30bjX0aE3wzmOenPs9fB/6ydlhCpAuS/AohMjWtNZ3mdeLctXNS5ytsSo4sOVjy2hL6VulLTHwMned35sPVH8qDcMLmSfIrhMjU7IbZsezkMgBmvjITJ3snK0ckRNpxtHdkTPMx/Prir9gre77b8h1tZraRB+GETZPkVwiRaR0NOZr48yulX5E6X2Gz+lXtx4rOK8iZJSeLAxdTY2INToadtHZYQliFJL9CiEwpzhJH1wVdASjvVZ7Z7WZbOSIhrKuhT0N29d5F6dylORJyhKrjq7Li5AprhyVEmpPkVwiRKX27+Vt2B+2mUPZCbOqxydrhCJEu+Hr4sv317bQu2ZqIWxG8+PeLfL/le7TW1g5NiDQjya8QItPZd2kfn2/4HIBJrSbh7uxu5YiESD+yOWdjXod5DKk3BIu28P7q9+k8vzPRsdHWDk2INCHJrxAiU4mJi6Hrgq7EWeIYUHUADX0aWjskIdIdO2XH0PpDmdd+Hq6Orvx98G9qT67N+WvnrR2aEKlOkl8hRKYydP1QDgUfophHMb5t9K21wxEiXXvJ7yW299qOT04f9l7aS5Xfq7DpnJQJicxNkl8hRKax7cI2hm8djp2yY2qbqbg6uVo7JCHSvbJ5yrKr9y4a+zQmJCqEBn804Nedv1o7LCFSjSS/QohMISo2im4LumHRFv5X83/ULFjT2iEJkWF4uHiwtNNS3n3uXVMytGwA6nPFxesXrR2aEM+cJL9CiEzhw9UfciLsBGVyl+Hz+p9bOxwhMhwHOwdGNB3B+JbjE5cV+LEAlyMvWzEqIZ49SX6FEBne2jNrGb1zNA52Dvzx0h84OzinbQAWS9qeT4hU1KtyL5r6Nk187z3C24rRCPHsSfIrhMjQIm9H0vOfngB8WvdTKntXTtsAPv0U7O3B0xNeeQW2bEnb8wuRCpZ3Xk7Pij0T3y8JXGLFaIR4tiT5FUJkaJ+t+4xz185RMW9FPqr9UeqfMCYGuneHgADzvkkTGDIEWreGjRuhdm147z0ZDRYZ3sTWExNLiF6d+yqHgw9bOSIhng1JfoUQGdaui7v4acdP2Ck7JrScgKO9Y8p31hrCwuDSJYiPT/k+nTrB1Kmwd69ZVqcODB0KEyfC2bPQrx/Exj5eHEKkU5/W/ZQOZTpw4/YNWk5vScjNEGuHJMRTk+RXCJEhxcbH0mtRLyzawrvPvYt/Pv/HO4DFAvnzQ758kD07dOgAhw799z5Tp8LcufDDD9Cjx4Prs2aFX36BUaPALgW/Xrt2NdvlzQudO8O+fY93DUKkMqUUk1tPpmq+qpyJOMPLs17mdvxta4clxFOR5FcIkSH9sPUHDlw5QNEcRfn8+RR2d9i/3ySZFoup0x07FsaMgS5dYMUKqFwZZs5Mfl+t4auvoGpVGDTo4edQyrzmzoVf7+uVGhUFvXqZ0WYwCfdnn5nSiSVLoEoV+FYm5hDpi4ujCws6LiB/tvxsOr+Jvov7ouWOhcjAHKwdgBBCPK4ToSf4fINJeMe1GEdWx6yP3unsWWjUyIzOhoSAlxd063Z3/ZdfwscfQ82H9Ac+dAhOnoQJE4izWNi8cSNaaypUqICHhwcWi4Xg4GBu3rxJ4cKFcZg+3ZRG9O9v9rdYoH17WLoUmjWDl1+G5s3NC+DaNVMy4fyMO1VcvgxXr0KRIuDm9myPLWxGvmz5+KfjP9SZXIdJ+ydRJV8V+lbta+2whHgikvwKITKcQSsGERMfQ9cKXWns2/i/N960CX7+GbZuhchImDPHJL73y5ULxo1Da82G9esJDw/HxcUl8ZXVYsFlwgQC7OzolC0bt27dSty1SJEinD17NvG9r68vH5YsyZkzZ3AcMoS333kHj1mzzOjuL7+YxPd+2bPDtGlm1PhpnTsHBQuakoopU+Cjj8zPDRqYkeY6dVJ2HK1N0l+mTMrKOESm5p/PnwmtJtBpXifeWf4OFfJWkMlkRIakHnXrQilVEPgD8AI08LvW+iellAcwEygCnAXaa63D/+tYVapU0bt3734GYQshbNXyk8tp9lczsjll48RbJ/BySyaRPXYMHB3B1xf27IF27eDMGbPO3h4GDiT+iy84deECERERia/jx4+zdOFCtqfw91T16tU5cOAA0dHR/7ldvnz52KI1RYoVgw0bHp3gjh9vtunVK0Vx3GPBAvNQ3q+/mq4Up0/Drl1w8KBJhC9ehK+/Ngnx/c6ehf/7P6hY0XSsCAszXwry5zej4n36SBIsGLR8EKN2jMLbzZs9b+zBO5v0ARbpk1Jqj9a6yv3LUzLyGwcM1lrvVUplA/YopVYB3YE1WutvlVIfAh8CHzzLoIUQIqnY+FgGrTD1tp/W/TT5xDcgAOrXh+eeg2XLwN/fJJGffAK7d8O4cTBrFi8FBLBo9epkz+MC1CtYkPhSpYiOjjavmzeJDgsjODKSm1FR9OvXj19//ZX4+HgCAgKYPn06Fy5coFy5cqxfv558R49SNCaGJUWKsHv3bsYAw0eMSNnI7rRpEB39+Mnvpk0m0ff3NyUeAD4+5tWhg0lge/c2Ncda3xvLsmVmm/h4k/wCuLiYh/wmTTLlGxs3wl9/mS8QwmYNbzycfZf3seHcBtrNbsfabmtxsneydlhCpNgjR34f2EGpf4BfEl71tdaXlFLewHqtdcn/2ldGfoUQT+On7T8xcMVAinkU41DfQw/O5BYTA+XKmQfLNm82da4AR47AgQPQsaN5HxnJa2+8wfTp0wGoWLEixYoVI3/+/NSqVYsGX35JrgIFTJnCHWFhkCcP+u13iPpiBKGhEBpqFl+7ZnLVW7cS/oyMI3r4aG4VL8fJIqWYNWsGTg52ZM+Rg9uxMYAFe3s77O0VWbKAl5crBQrkoGnTWhQq5IHb7yNxW7+YHHvXkicPuLqmIGeOi4PSpU1t8Z49powiORbL3Yfy7ti+HerVM+UN8+dD4cL37qM1fPedGS0eOtT0NRY27UrkFfx/9+fijYu8Ve0tfm72s7VDEuIBDxv5fazkVylVBNgIlAXOa61zJCxXQPid9/ft8wbwBkChQoX8z5079yTxCyFs3NWoqxQfXZyIWxEs7LiQliVbPrjRxIlmtHTZMnjhBa5evcqePXuIjo7Gzc2NS5cucenSJU6cOMGkSZOwJExEsWzZMl544QWuX4cLF+BC8z6cz1GeCy37mfcXIDgYwk6FERqdlRiypOm1OznF4elpIV8+R7y8FF5epmy5cGEoWtTk+IUPLMK5XSuYNw9eeum/D2ixmO4Wbm5QvTqULWuS5927wcPj4fv162cS5DsP8QmbtuPfHdSdUpfb8bf5o80fdKnQxdohCXGPp05+lVJuwAbgK631PKVURNJkVykVrrXO+V/HkJFfIcST6rekH7/t/o3GPo1Z0XkFKrmh0GbNTI3rsWOEhYeTK1euZI7kBBQHSlKIUhQp+wrRLhUJDFRcu5ayWLKoW+TK40AuLwc8PCBHDlMh4OIUR5bAg7hUKUMWdydcXEzzhhthoVzefxAn7/w45/QAFLGxFmJj47h27TaHD58jIOAk4HbfKyfmcQuXR8aklCanQwg5C8bhmTuSAgWi8PGJpnjx2zz3XG7KlvW7+5lpDblzmyT599+JmzePvpMn88/OnQwaNIiBAwfi4vLocwoxfs943lj8BlkcsrD99e1UyFvB2iEJkeipkl+llCOwGFihtR6ZsOw4UvYghHhSN26YB7CWLjUPo2XJAg0bwogRD2x68MpBKo6riEIR0CeAMnnKPIIDwd0AACAASURBVLBNaGgowb17c9XDg7P16nH48GG++24yUDnx5ZalOjdj8qF18g9tubhAQe9YCobup2CNghSsmpeCBU3jhLx5zbNfuY5vJWu75qas4PRpE39cnKk1njzZdFqYNOneSTACA6FkSfjxRxg4MNlzX7hwgX379nH08GGOfvstR11difP2plKlyoSF3WLHjjMEBcUBeTAJsTfmeeMiQFGgIA9/jOMWDg4nyZ//GiVLxuLrG0nZBZ/gVcCZ423aMGTIEOLi4hK3LliwIF999RWdOnXC7v4H3EJCTHlJgQIPOZewNb0W9mLivomUzFWSPW/swdXJ1dohCQE8RfKbUNIwFQjTWg9Msvx7IDTJA28eWuv3/+tYkvwKIQAz8lirFmzbZupUS5c2xbJlypjaUoDbt8HJCa01Df9oyLqz6xhQdQCjXxyN1prjx4+zcuVKVq5cya5duwgOjgFqAtW5m/Dmf+DUdnYaHx9FSZ/blFr5MyXbVaDUW40pUQLy5Ekohb3ze/FhhbYXLpga4saNTd/g+HizbZ06ppVYw4YP7lOjBgQFweHDj+63Gxtr6pbvq9uNiIjg4sWLHD58mHXr1rFp0yZOnz6Ni4sL+fMXxjNXBbK6luHGjdyEhHgSFubNtWsFuHUrz0NOdAbYmfDaBewFbiaurVSpEj/88AMNGjS4u0vx4maij7///u9rEDYjKjaKauOrcTjkMD0q9mBS60nWDkkI4OHJL1rr/3wBtTEtzg4A+xNeLwK5gDXACWA1Jvn9z2P5+/trIYTQWmu9dat5JWfmTK0rVND6+nU978g8zVC0x3ce+vDpw3r9+vW6adOmGrw1tNPws4Z9GuK1yVrvvrJkidGlSwfrBg0C9JRBe/RByuhbTVtpHR6u9Zo1ZqNFi8w5T5zQumdPsy6l4uK0PnhQ6z17tA4L++9tN2/W2s5O6+bNtY6KenC9xaL1hAlaR0Sk/PxJjRqltaur1tevP7AqIsKi//77tH7ttTW6SpVN2jv3Ee3AjQc+L4jTHh7ndZMmh3WuXD015NCAbt68uT58+LDW8fFaZ8um9YABTxajyLQOXTmks3yZRTMU/deBv6wdjhBaa62B3Tq53Da5han1kuRXCBsXHa312LEmafwvK1ZobWeno/v21j4/+WiGol8b1Uvb23fQMF7DyQcSN0dHi65RKlS/Zz9CT//qlA4MNLnaPcaN09reXutWrbQeNkzrLFm0/vVXrbt319rZWescObTesiXVLl+PHWuCffFF8/7iRa1Pn9Z69myta9c267744smOvW2b2X/EiEdv+913Os4xiz64LkRPzDlYv1lsta5USWsHh3s/U6Us2s5un4YRWqnWulb5uroH6GFt2+rx48fr0qVL6549e+rTp08/WcwiUxm3e5xmKDrb19n0ydCT1g5HiIcmv4/d6uxpSNmDEDZu6lQz8cKGDVC37n9uGtunHx0PLWZe4wuoEF/0b4fBcre1maurhVq17KhTx1QbVKsGLrevkVi/sGmTeRLtfnv3mp5kxYrdnenN3d20QfvsMzOhQ2patQquXzezvPn5mQk5APLlg88/h9dff7JZ3rSGF180vXg3bIAqD97pu2fbo0dNuckrr5jZ786fJzrWgR07YN06WL/edEC7fTvpjnHAFswjIEuAowAopWjRogUDBgygUaNGD9YJC5ugtabDnA7MPjKbKvmqsKXnFun/K6zqmbQ6e1qS/Aph49q3N3W+58/fk+BZLBauXr3K6dOXWLQolnnzNMfPeKLfqQJZw2DaMjjZkPz5z9G/vy+NGysqVgSH5J7vWrUKmjc3s7v99pvpX3vnXOHh8MMPpstBxYpw8qSZCa5IEetM3DBzpqntLVHCtBxL9oIew6VLULOmeSht8mRo08ZcH5hr/f57066sQpIn8hctglatzLTL97Uwi4oyCfD69bBi4Q32BGQhHsckW5wGluDgsIK4uNVADMWLF6d///50796d7A/rNSwyrYhbEVQcW5Fz184xuMZgfmjyg7VDEjbsiWt+n+VLyh6EsHH+/lo3a5b49tSpU3rgwE+1u/sbGmZpiLx7273e55qhaNc3/PSQoXv1lSvRKT/P2rVaFyhgDlShgtYNG2pdrZrWTk5m2UcfpcLFpRNBQVo3aaL1smXmMyhaVOu8ee/Uhpja4KQsFrN9587JH89iMS+tdfjizXrmtNu6a1etPT3vLZFwcrqls2adp6GlBiedNWtWPXDgQH358uXkjxsfr/Wff2pdt67Wu3ebZatWaf3881rPmZN4TpHxbD2/Vdt/bq8Zit54dqO1wxE2DCl7EEJYXfXqaHd35vUZzLBhBzhwoBTQBJJMGuHufpzS/vsIeL430ZZINlT7jbrN+jz+uaKizAxtCxfCqVOmK0OlStC5870jn5lVdDR8/bVpI+fsDOXLmxKH5Mo6IiJMBwoHh7ulH56eppXbzJmmPVvnzvfsEh8Pu3aZj3jxYti//+46e/ubxMfPAWbh4rKFt956k/fee4/cuXObDUJDTSzr10OpUqYcplo1MznJW2+Z/16tWpmOEq7SNisjGrJuCMM2DqOUZyn2v7n/wdkYhUgDUvYghHg60dEmwwkJMe2u/PxSvGtUVBQbN26mWbNvgK7AK0C2hLUWKlSIpFu3bLz8sqJQIfh4zcd8s/kbGvs0ZmWXlalwMSJZN26YOuB//727rFw5+Pjju1NDP8Tp0zBrlsmVkybCEAbMxcnpb7p0Kcagfn0p8+abcPAgjBljasCT1gjHx8PPP8N770GLFrBgwZPVQAuruhV3iwpjKxAYGsjQekMZUl+mxBZpT5JfIcSTuX4dhg6F8eMhMtIs+/BD+OYbc8c7IgJy3ju5Y0REBKNHj2br1q3s3XuT4OAmQBegcOI2hQtfpF+/7HTp4oa39919g28GU/SnokTFRrG91hSqN+qW6pcoktDaTNRx/Tp4e5uZ4B5TYODdRPjQoaRrjgGTqMsffPh/vWn6+ecPfzhu1CgYNMiM3LdMZiprke5tOLuB+lPr42TvRECfAEp5lrJ2SMLGSPIrhHh88fHmdvT+/dCp093b5l5eZoav+fPNQ1LLlkGFCsTGxjJu3DiGDPmasLB6QD+gTpIDnsXZeTY7N/aifIFo0+HgPu+ueJcft/9Iy0DFQvUaTJuWVlcrUsGRI/DHHzBpUhwhIXce6IsDllCw4Cq++qo2r73WDvv7HziMjzdlEc8/f+/IsMhQ7sz+VrdwXdZ1W4edkv+WIu1I8iuEeDIbN4LFAvXrP7ju4EF48UV0bCxLvv+ed4ZO4vTpRkAvzBS8piVZu3bQrZsddesm5DHNm5uEeuxYc2s74bb2v2FnKTi6KAD7Fheg4vL9Zk5hkeHFxZnvSGPH3mb5cnssljvJbhCenrP48MNcODiEUbNmTapWrWrVWMWzExYdht+vfgTfDGZCywm8Xvl1a4ckbIgkv0KIxxMXl6LWW9um/U3fLlMJoC/QEjBJTblymv79FZ06JTOb74ED0KGD6XFbsKBpO9aqFepi78RNdK9/U7/nrkhbO3bA5MlcGfAFU5bkZNSoSC5fvtOLOQaYDvzM7t3j8ff3N33W5s41D+Fly3bvsbQ2X5q0Bg8P84BjyZLQtq15aE6+NKUb0w9O57V5rwFw5p0zFMlRxLoBCZvxsORX7j8IIZLXvDn07ZvsKovFwoIFSylZ8mtqdilPACuANtjbQ4cO8WzaBAEBijffTCbxBdN5ICDA3A+vVg3OnePIhb2Jqz+u9ZEkvpnRpUswbhxet87xwQcOBAXlYPnyOCpWPA84At2BvTRvno1Roy6y6fdJXP3hhwe/hM2bZ+4YWCwmAR4yBN55x9QnDxtmOkisX5/mlyeS17FsR7zdTGF/0Z+KWjkaIWTkVwjxMNmzm/ZWv/5KREQEO3bsQCnF4cPn+e67cK5ceRUoAIBb1jDevv0zb89ohNfLtZ/odB3mdGDW4VkA6CFp93tJpKHTp83kIyNHmofZ7lvVs+deNmzwBe5MjnEa+I5ixbZSp05V2rRpQ4v4eOxeftlMCrJ06QMPW3LwIPTqZab969jxv2e6E2lm07lN1J1iZnUc33I8vSr3snJEwhZI2YMQIuW0NsW5n37KqW7dKFasGJAHeAfzEJu5Ve3ldZWhQ13p2d0ZJyee+MGkPUF7qDK+Ck72ThwfcFxui2Zmzz0HV66YNhDJ9PBds2YHgwbt5cjhJsRbfBOWBgE/AL9T2v4WHxYqRMd9+3B82Axyt2+bdnw5csC+ffLAXDrh/7s/ey/tpVyecuzvs18efhOpTsoehBAppxR4enL8yBFatHgdGAmcBT4GclCixBXmz48lKMiTPn1ccMpi98QJhtaawSsHA/BO9Xck8c3svvvOtFLr1cvUld+nYcPqHJjbiBj36sws+B7ly8UD+YCR2KlzHIn/iK5nwilRsSJjxowhOjr6wXM4OcGXX5ra8o0bU/2SRMps7rGZAu4FOBh8kOkHp1s7HGHDJPkVQtxDa82SJUuoEZedUnOrcOzYEmAQ4EKbNpqtW+H4cS/atHG8m+8uXw4vvQQXLz72+RYFLmLDuQ14uHjwcZ2Pn+WliPSoXj0YPRrs7c0EGlOmmGQ4LMw0CAYoXBj7Vs1pv6E/+wPsWbwYatQAi84FfIGd3QXOnu1C//4f4Ovry7hx47h8+fK953npJfOFTGp/0w0XRxe+eP4LAD5Z+wm34m5ZOSJhqyT5FSIz27ULPvkE2reHGTP+c9P4+HhmzJhBuXJ1aNFiD9sj9gIfAq40bXqLfftg/nxFjRr37ag1fPuteTI/T57HCs+iLXy4+kMAhtQbQo4sOR6xh8gU+veHnj3NHYYePaBIEdOdoWJF07XByclMeVy0KEqZZy+3bIF1i2/SqM4tLBY3YBgODue4dKkdffq8jbe3N0oplFLMmzfPPGmZKxfcSYq1Ni8wB6tf38wuN3OmKZMQaaJL+S6Uy1OOc9fO8cvOX6wdjrBRkvwKkRmFhpqRr2rVYPhw01khPt6sO3MGunY1I21ATEwM48ePp0SJ8rz66n4OH14EfAa406hRHDt3wvLZsVSsmMx5tIYvvoANG8wscI6OjxXmouOLOHr1KIWyF6JPlT5Pc8Uio2nQwNT97toF48bBjz+a7h8PoRTUb+7Kqo1ZWL/ejATHxXkAP+HsfA7TKcK02Xv55ZeZPGkSuLubBzf//ReaNoU//zQHc3c3nSKWLjUPxZUrB3v2pPIFCwB7O3uGNx4OwFebviIsOszKEQmbpLVOs5e/v78WQqQyi0XrunW1dnLS+ptvtI6IuHf97NlaOzpqXaqUXjtnjvbx8dXQQcMZfWd4rF69eL15c8L2p09rnTu31gMHan38uDl+fLzWoaFat21rdujSxSx/TDUn1tQMRY/aNurpr1tkfmvWaP3WW1pbLNpi0XrhQq3LltX67rDuYQ0tNaABvXrVKq2PHtW6QAGt3dy0/uOPe48XH6/1P/9oXbCg1p6eWt+4YZ3rsjEWi0U3nNpQMxQ9eMVga4cjMjFgt04mH5VuD0JkRocPmyfqGzRIdvW1pUv5X6tWjI+vAvwImFqGsmU1I0cqGjdOuvE1eO89mDTJjJa5uJhR5E2bzK3rTp3g/fdNDedj2HJ+C7Un1yZnlpycH3QeN6fkGgILkcSff5q7FnPmwMsvA+av4vTpZh6MM2fubLgCGEj27JfYWrAgpS9cMHcnKlRI/rhBQWbkd9cuc6AUTO4ins7eS3vx/90fgBNvnaCYRzErRyQyI+n2IIQtKVMm2cQ3IiICpRQ5mvdlfPwfwHagBl5emvHjYf/++xJfMLeNx4+Hs2fNdMR9+8LAgaamMiAAPvrosRNfgOFbza3P/lX7S+IrUqZjR5PAvvGG+buH+avXubOZLHDUKMiRQwNNURzk2rWh1Dt0gRcKF+avQ4ceftx8+cwXuy++gAUL0uZabFxl78rULFgTgOKji1s5GmFrJPkVIjOZM8fUNoaGPrAqJCSEnDnzYB5iOwa8hhO3+LjEHE6cUPTq9YgctmBBePNNGDHCtKsqXfqJkl6AoyFHWXh8Ic72zrxV/a0nOoawQY6OZrpjFxdT9Pvll4mrnCLDeKfAXAILN+FNxmL+eRvIVU6w4kB1OnfuxrZt2wDYtWsX7777Ln/++ScX73QoadHCTJO8YkXaX5eN+qzuZ4k/LwlcYsVIhK2R5FeIzOT4cVi58oHJA5YvX07Rot2B/cA3gAsdOmgCA2L46tjLZMuWtmH+sPUHAHpU7EEe18frECFsnK+vKU9o3hyOHLm7vFAheOUVcocdZ+z0HOzdZ0ft2hrIDfwO7KRTp5Hcvn2btm3b8uOPP9K1a1eKFi3KhAkTzBe54sVN2zWRJpoWa5r4c4vpLawYibA1UtgkREajtWnVtGyZmcL1hRfMsuDguxNNWCyJm//883TeeScGMCMrefNe588/3WnUSHF3Gtm0E3QjiD8P/IlCMbjm4DQ/v8gEvL1h9ux7J8m4czeiTh1wcKAisHGjYvZsGDzYwr//VubMmRmUK7eMf/81HQZKlSrFsWPH6N27NwEBAYyMjMQxb17rXJONqpS3Evsu7wMg3hKPvd2T3U0S4nHIyK8QGcmZM2aSgDp1TAuzO3WMhw+bsoTt2837gAC0hokT4d13m2HaQMUAQ9i9O5ZGjRKON38+vP46xMam2SX8tP0nYi2xvFz6ZXnIRTydpA+m9e8Pzz9/zzLVpTPt/27DsWN2tG//LwCBgS2AQ0AThg0bxsSJE3F0dOSXX36h6fHjXC0u9adpadvr2xJ/3ntprxUjEbZEkl8hMopTp0yd48GDMGYMhIebLgxgmvl36QILF4JSnP9pPo0aaXr1gvj4HMBqWrf+lJs3PyB//lx3jzllCqxZ89j9eZ/UtVvXGLtnLADv13w/Tc4pbJi3NyxdimvYBWbOLMDMmedwczsJFAVWMHfuS7Rs2ZP169fj5eXFurg4qs6axbFjx6wduc1wdnCmb5W+APx98G8rRyNshSS/QmQUixebvk5bt5qOC25JOiR4e8PEidz+fTyjdQ9Kz/6MtWsVcBU7u6789tspFiwYTtasWe/uc/iwOWaHDml2Cb/v+Z3rMdepX6Q+VfNXTbPzChvVv7+p5e3VC2Jjad/eh9BQX9599zIuLpqZMx3w84OgdbnYvWkT/v7+nD1/Hj8/P3x8fPjrr7+sfQU2oWelngD8fehvYuPT7i6UsF2S/AqRUbzzDpw+DX5+ya6eOnU1zm/k420mctPiBswnR/ZarF7dgz593rx344sXTZ9UD4+7o8epLCYuhlE7RgEy6ivSSJEiMHq0eQi0cWO4dAknJ8WIEXk5eFDRqH4soaHQ7v9K8nGLUBYu3ECjhJqgM2fO0LlzZwYPHkxc0tribdtg507z89Wr5q5LkyYwbdrdWRTFY/H39sfP04/gm8GsPLXS2uEIGyDJrxAZwZ1/VJNpy3Dy5Cn8/UfRvbs/8CIQBnSiZ8spHL99geenTYNVq0zZxN695lhvvmkS4HnzIHfuNLmE6YemE3QjiLJ5yvJCsRfS5JxC0KuXmTb5zBm4c+dj+HB8u9Zi5dZsjKEvLvYx/Bn4HM8958rgwUsYMWIEbdu2xcHBgZEjR9KsWTNCg4Kgd2+oWRN+MN1KcHaGdu1MD+wuXaB2bbh0yWqXmlEppehaoSsAfxx4+BTXQjwrMsObEBnBiy+aTg6LFycu0lrTr99HjB1bCbhTurAYeIMdOxZQrXhxM1vVxIkQHX33WNu3Q86cpiNEqVJpEr7WGv/f/dl3eR+TW0+me8XuaXJeIRLdvg1OTubngQNh/36oXBk6dybQrTJdutwd0B00CL7+Gnbu3Mgrr7xCSEgInk5OvH/7Nv0GDcL1iy/ubSeoNfz9t/lSWaSIOVDSEiPxSBeuXaDwqMI42TtxafAlcrrktHZIIhOQGd6EyMi0hsuX71k0fPh6xo59A+iAg0M0w4b9y7p12YiP/5dqYWFm+uHRo00LtFWrYPJkMwmGry+UKJFmiS/A9n+3s+/yPnK55KJj2Y5pdl4hEt1JfMFMBbd+PYwcCZUrU6KE6R74+eemRPjHH8HfHzw86rJnzx5q+/py9fZt3geKTpvG2D//vPfYSplpvhcsgB497u0/LFKkYPaCNCjagJj4GGYfmW3tcEQmJ8mvEBlBiRJw9CjcusXt23F07LibDz+sDfhQsGAwR4648OmnBahfvx52dnbwv//B77+bfd3coFEj6N7d1Pl6eqZ5+L/s+gWA3pV7k8UhS5qfX4hHcXAwN0q2bYOSJU3+Wq0arFlTkI0+PiwrW5Zq1aoREhJC3759qVOnDu3bt+f69et3D9KokSmvaNTIjDSLx5JY+hAgpQ8idUnyK0RG0KoVOiqKaQOG4eGxjZkzqwCOeHlN59ixXNzTmvTUKdP/t2ZNa0V7j8uRl5l9eDZ2yo4+VfpYOxwh/lPVqrBnD3TrZqqFevSAnvlXUnfdDrZv384333wDwObNm5k9ezZ+fn5YkkwqQ/36cO0aBARY5wIysLZ+bcnqmJUtF7ZwKuyUtcMRmZgkv0JkALF16lDHrRldJr7FzZt1sLML4803FxMY2JysWZPMiGSxwODBkCWLGelNB8bvGU+sJZZWJVtROEdha4cjxCO5upoW2JMmmf+VpkyB6s9n5fhxxfvvv0/jxo0Ttw0KCuLVV18lMjLSLChY0PwZHHxnA/jyS4iJMe+nTIEWLUwt/p1lAgA3JzdalmgJwPKTy60cjcjMJPkVIp2Li4un2YtL2BL5D+CNr++/BB7PytixLXB3d79347ffhn/+MU/r5MtnlXiTio2PTZzUYkDVAVaORojH06MH7PxtDyVdznPoEFSpAjNm2LF06VImTJhA9uxmevBZs2ZRrVo1MzlGSIjZOUcOGDHC1Nh/9pkZTgbzIFxgoOlCUa4cHDhgpatLn54v8jwAG89vtHIkIjOT5FeIdEprzU8/jcPR8Q/WrHkJcKRjx385dqwAvmummqlcP//cPLTTpw+0bWtKHT77zDzNng4sPbGUoBtBlPIsRYOiDawdjhCPrVx5xa7oMrxW4zQ3b5rn2j74wIFu3V4nIiKCo0eP4ufnx9GjRylfvjzjR4wgztnZdH947z144QU4efJuGVL79nD8OCxfDlFRpkzixAmrXmN6UrdwXQA2ndtEWnajErZFkl8h0qFTp07RvHkfBg6sDPQAoujadTnTpxfAwQHT9iw4GIYONWUOf/9tnjhv2dIkxEpZ9wISzDw8E4DuFbqj0klMQjyWihXJ5pOHabfaMeYXCw4O5vtmixYQEQGlSpVi586dtGjRgtjYWN5Yv55q7u4cGjPGfAmdNw98fO49plLQtKnpyOLqavpvCwBKeZbCM6snlyIvcSpc6n5F6pDkVwhrOH8eBgyA/PnNv6BgGvH37cuC337Dz68Py5Z9AVTF0/MGy5dfZ+rUJBND9O5tpieOjobwcPOAzdy5yU6CYS1RsVEsPL4QgHZl2lk5GiGekJ0dfP01at9e+h59mzWrLHh6wooVUL26GcR1c3Nj+vDhlEmYfXFfSAiV7Ox4/do1fh0zhqCgoOSPXbSoGRWuX9/MFidQSiWO/m48J6UPInVI8itEWlu0CMqUgfHjoVYtiE2Yy/7yZWZMnEjbfjuJjV0K5KFOnVscO5aNpk3zJn+sLFlMbWE6HFVdcXIFN2NvUiVfFXxy+jx6ByHSqw4dzB2WX3+lLhvZtQvKl9cEBkL1qvEs6zYDt5o12V28OEeOHKFPnz5YgEmTJzNgwABKly7N5MmTuXXr1oPHvn4d8uaFv/5K88tKr+oWkuRXpC5JfoVISzt2mF67pUqZIaNZsxKnF15cugyvxn6GZjLgyKDaO1m7Ngu5clk35Ce1MNCM+rYt1dbKkQjxDHz/PaxeDfXrU6QIbPHtRlvmcu2GPS3+aMfIvMNx/mEEfg4O/NalCwH79jFgwACqV6/OtWvX6NmzJy4uLvzxx309bD09zZfXOw/KCeoUrgNI8itSjyS/QqSlcePMKM/KlWYa1AQxMfD22zmB/0MpC2P8JzJyaw0cThy1WqhPI94Sz5LAJQC0KtnKytEI8QwoBQ0bJr51a9WA2Z8EMKTFHizYM/hYb94dUwzLrDlQqxZlixdn9OjRbNu2jT+TzAjXrVs3tm7deve4ISFmBkcPj7S8mnStglcF3J3dORNxhvPXzls7HJEJSfIrRFqaNMk0v895d9760FCoXz+WM2dqApEMHLiWvmvbmQdl0nAK4mdpx8UdhESF4JPTh9K5S1s7HCGeve7dsftyGEMX+fP33+DoaGZNfnXRq9zC2dzZwdSwdu7cmfPnz1O0aFEA6tWrx4oVK8xxlpgvidSqZY2rSJfs7eypX6Q+AKtOrbJuMCJTkuRXiLSWJPE9fRpq1LCwfbsjcJGSJXvzxRc1wN0dWrdOl7W8KXHnQbdWJVpJlweR6b36qulc5u4Os3YUoSkrCf/l3hreggULcvz4cV5++WXi4uKYPHmyGfEdOdL0+61WzUrRp09NfZsCsOLUCitHIjIjSX6FSCt//QV168Lt24CZgbh2bThxwg4IwMurNWvW/ICrq6vZ/tw5GDQIjma80ofE5FdKHoSNaNDAdC7Llw82Upc6E7txYcqae7ZxdHTkiy++AGDhwoVcCQ6G6dPh998z7Bfd1HIn+V19ejXxlngrRyMyG0l+hUgrUVHmX8cLF9i5E+rVg0uXANbh6NiQJUvGkT9//rvbnzlj7qOajTKME6EnOHr1KDmy5KB2odrWDkeINFO+PGzbBn4l4zlMWWq8U9V8d00yWYOfjw+tq1QhOjqab7/9Fv6fvfuOjqrq+jj+vTOpQAJJ6L33KlV6EZEmCgKCBXtBxfrYO6JiAUQRUUEsrw0UsYCIqCDSUamCFOk9nfRk7vvHSSEgUlLuTPL7rMWayb1ni+oA8QAAIABJREFU5m6exXPdOfecvZs2hQ4dnAvaS9UJr0PtsNpEJ0ez5sAap8ORIkbJr0hhyezw9POkdfTqBVFR4O8/H+jHSy89TuvWrXOPz+r6VK1a4caZR9/8/Q0A/er1w9/t73A0IoWrenX4bbmbLh3T2R8XSrdusO6xz01N34YNoUwZnl5jkrmpU6eevgawcHHtiwFY9M+iM4wUOTdKfkUKS5MmfNPgAfq+0Y/jx6F06W9JS7uUSy7pzpgxY3KPtW348EOoUQPq1nUm3vP05V9fAma9r0hxFBYG3y/0o08fU8yhx+TLWNPwajM1fMcdtJw3jyGDB5OSksITTzzhdLheK6vZxdI9Sx2ORIoaJb8iheSTT2Dw9vGkEESl4JnExl5K48b1+fjjj3G5Tvq/4ltvmSUS//ufT60FPBh/kGV7lxHoDqRfvX5OhyPimBIlYO5cs281OiGQXsvG8tvdn8Mrr0Dfvox97jkCAgKYMWMGc+bMcTpcr5S1bGrZ3mV4bI/D0UhRouRXpBB8/DFcfTWkZ7jo2+w7DiZdT6VKFZk/axZhpUvnDIyMhPXrTS3gyy+H225zLujzMGfLHGxsLq5zMSGB3tNqWcQJgYEwaxYMG2YauV18Mfz0kznXqFEjxo8fD8BNN91EZGSkg5F6p6qhValYqiKxKbHsjN7pdDhShCj5FSlgn30G11wDHg8MGbKeH7dcDsDzzz9P9ZdeggoVzHrgli1N0jt0qJku+uILcLsdjv7cfPHXFwAMaTTE4UhEvIO/v/nld9Qos+e1f3/T4wZgzJgxdO/enaioKN5++21nA/VClmXRupLZC7H2wFqHo5GiRMmvSAH64gu46iqT+PbsuZQvvmhBWloa99xzD9deey0MGAADB5pnpFWrwoMPmg+5XD613AEgMS2RJbuX4LJcKnEmcgK32/S3ufVWSE6Gyy6DxYvB5XLx8MMPAzDljTdIe/VVU+UF4LffTP20CRMgNtbB6J2VnfweVPIr+cfP6QBEiqqvvoIrr4SMDAgLm8pPP43GsixeffVV7r33XjPoiivMnyJgzYE1pHvSaVGhBWHBYWf+gEgx4nLB1KnmfvDuu+b33oULoXfv3jSoUoWt+/cz6IEH+K50aaybbjIZc3Q03H8/vPSSqRN+Qnvl4qJ1ZSW/kv808ytSAL75BoYNs0lPBxhPdPRoqlevzuzZs3MS3yJmxb4VAFxY9UKHIxHxTpZl9rKOHAnHj0PfvrDuue94dv9+AOYDn2Y1uenQAf74A1atgrJloV8/swm2mMma+f394O/YJ9RLFskLJb8i+eznn+GKK2zS0izgVdzux5k0aRI7duxg8ODBTodXYLKS3w5VVbBf5HTcbnj/fRg8GGJi4OKnOtC0ww28PWUKAP/73/9ISEjI+UDbtibprVULfvnFmaAdVDmkMhVKViAmOYZ/Yv5xOhwpIpT8iuSjP/6AQYNsUlMt4A0qVHiFX375mbvvvhs/v6K9ykjJr8jZ8fMzpQ/7dkvkGOXotfNtul10G61bt2b//v288MILuT8QFgZr18KFF8LKlc4E7RDLsnKWPmjTm+QTJb8i+WTHDvMYMz7eAj6jbNmxrFmzms6di36L35jkGA4eP0iwXzD1I+o7HY6I1wsIgC/ml6BHDzh0xE3fvi6eeWYqAK+88go7d55U2isoyEwXf/ihA9E6S5veJL8p+RU5V7YNs2dD164wfToAh/84QJ+Guzl8GFyun4FrmTp1ClWrVnU21kLyT7R5HFkrrBaWj1WpEHFEbCzBG1bx9ScJtG4NO3fCM8+0ZfjwG0hJSeGuu+7KvcbV7TYVYQ4fdi5mhyj5lfym5FfkXKSmwogRphbvoUMQGkpcHPS9JoId6TVowBo8nkuBVPr37+90tIUmay1e7bDaDkci4iN27ID27Sm1fCHffWeW9K5eDceOTaF06QjmzZvHfffdlzPetk2v5DJlnIvZIS0rtgRg05FNDkciRYWSX5FzMWaM6Vrx/PPw11+kXTaUwYPhj02B1K0L3Ya9DxynORD0ySdOR1todsXsAqBm6ZqOxiHiM2rWNK8bNlChAsyfDxERsGhREE2bLgFg0qRJzJo1y4xbswaOHTPrfouZaqWrEegO5ODxgySkJpz5AyJnoORX5Gz99RdMm2bqbj7yCLbLzR13wKJFpknb88+v4Z1ZU3C73cxo3x5r2TIzW1MMnLjsQUTOQng4dOoE770HKSk0aABff22W9v72W2PgUQCGDRtGbGwszJ0LoaEwpPh1T3RZruynStujtjscjRQFSn5FzlajRqacwxNPAPDaa/DOO+Y/Vm+/fYgxYwZi2zYPPvggrZcuNZXsi8n6112xuwCoVUbJr8hZe+op09Ft9GjweOjY0bRCNreNccBIAG699VZ49llT6qx0aQcDdk69iHoAbIva5nAkUhQo+RU5Fy1bQunSzJtnJoABpkxJ4Ikn+nDo0CF69OjB008/bWoZ2TZs2mSKeRZx2cseytR0NA4Rn9K7t/lleuVKSEkB4PJBHl57zZwO8PuAIFc7PvvsM2Z98QW0auVgsM6qG1YX0Myv5A8lvyJna+ZMmDqVjRtN22KPB+688ygPP1yL9evXU79+fWbPnk1AQIAZv3o1NG0KS5Y4GnZhUPIrcp6eecbcI4KDISoKSpfmzpeqc4vrXVLT3QTZXwIVueOOOzh69Kj5zN69MGdOzncsWWJ+0S7Csmd+IzXzK3mn5FfkbC1YwJGx0xg40CY+HoYPhwMHbufo0aOUL1+eb7/9lvDw8Jzx5cub18hIZ+ItJIlpiRxPPU6QXxBlgorfTnSRPLEss/4XIDkZbrwRq1dPXh+zjS6NjhFjVyE09EeOHo3lqqFDSb/6arNZbsQIU30G4IUXzC/affvC7t2O/VUKUt3wzJnfaM38St4V7ZZTIvkovftFDPv0VnZh0bYtjBy5kEGDviAwMJA//viDypUr5/7AwYPmNSKi8IMtRMcSjwEQERyhGr8ieVG5MkyaBEAAMPsItGkDe/c2IShgBgsXX40/sOG662j68MPg728+N2OGaX7x3HPQrp2ZCW7QwLG/RkGoF66ZX8k/mvkVOUuPbbmGxXSnov8xXhm/hWuvHQrAI488cmriC+axpMsFHYp2u9/IRDOzHVGiaCf5IoWtfHlT5CE42CY59SpgDADNZs40yW3WL5uVKsGDD8KqVWavwRVXQFqaY3EXhKqhVQlwB6jcmeQLJb8iZ2HOHHhpUgBul4f30y7n5sEXEhsby+WXX84TmdUfcomNNWXRhg3LWf5QREUmZSa/wUp+RfJbq1bw3nsmybWsCYBpl34w68nSiRo2NFVmevXK3kBXVLhdbuqE1QG06U3yTsmvyBls2wbXXWfej3/JxbwO6fwdE0OTJk34cPp0XCc+6s/IMMsdSpeGBQtgwgRHYi5MmvkVKVjDh5uJXdt2ExQ0ByjLZZddRnx8/KmDL70UBgyAefMKPc6ClrXpTcmv5JWSX5H/kJhoniDGxcHgwdCu3VImr1yJ2+3mgw8+oOTYsaYv6YAB0K8fVK0KPXtCerpZ7lCpktN/hQJ34ppfEclnR47AU0/x3A076dQJkpPLEhw8i1WrVhMaGsrEiRNP/cyECTB+fOHHWsCyyp2p1q/klZJfkSzbt5uZW4CffsK+515GD9jD+vVQp04GNWs+Q9++l2DbNo888ggXXHCB2Y3Spg0cOACHD0OPHqZ0UTHa+JW17KFsibIORyJSBB0+DM8+i//6tXzyidk/m5TUHXgQgPvuu4/Dhw/n/kxqas5muCKkRpkaAOyP2+9wJOLrlPyKJCbCbbdB/fqweLE5tmcPH02N5/2fqxNsJWEl9GXChKdJSEhgyJAhPP7442bcyJEwezb8/jusXWvaMw0bBm63c3+fQpY1CxMSEOJwJCJFUNaegb17qVbNFHUAcLtfIGv975tvvpkzPiXF3IuaNi3cOAtB+ZLmf4sjiUccjkR8nZJfKd5SU6F/f3j7bbjnnuz/YPzT7TruCHwHgAtLPsL2QwupGx7O8uXLmT17NoGBgU5G7VU+Wv8RADuidzgciUgRVKGC2cj2+edg2/TtCw8/DBkZFmXL/giUZerUqSQlJZnx06ebrpLDhzsadkEoV6IcAEcSlPxK3ij5leJt/Hj45Rd4/32zTq58edLT4eqrIT7eomXLHfx0/DX8XS4+u/FGOrRv73TEXis8OPzMg0Tk3N17r2mB/NZbAIwdC507w7FjgZQuPYujR4/yWlZP5KVLzb6Diy5yMOCCkT3zq+RX8kjJrxRfKSkwcSJcdhlcc0324eefh2XLICIimT//bAfAq5MmccFLLxWrtbxn67Euj1EmqAx3tbvL6VBEiqYbbzQbaj0eAPz8zAqr0qUhNrY7cB3PPvMMO3fuhA8+MLPERfBelZX8Hk046nAk4uss27YL7WJt2rSx16xZU2jXE/lP6elmlqR8eWjcGIDly6FLF7Pv7dZbZzFt2jBq1apl/qNi26Z2b0QEDB3qcPDexWN7cFn6XVqkwKSl5Wxie/llWLqUj+IHcc3PN+BHHOk055I+DZk3f36R7bSY7knHf6w/FhZpT6ThdhWfvRVyfizLWmvbdpuTj+u/VlJ8+flB9+7ZiW98vFnukJEBdet+xbRpwwC4//77zXjLMksj/u//HArYeynxFSlgJ1ZvSE6G7du5atOjDAlZQDqhuF0f8f2ChXz//ffOxVjA/Fx+RARHYGNnV5kROR/6L5YUX7Gx5tnhP/8Apoj8zp1Qs2Y027fnbBa59tprcz5Tu7YpayYi4pQnnoBNm7AOH+KtnX2oUAEyPJ2Be3nooYfIyCrZWARp3a/kByW/UnwlJ8NVV8H77/PLL2Yvib8/XHXVD0AqALfeeishISeU8IqNhZIlHQlXRORkZcuaAg/GODZssPnyyy+dDKlAZSW/h44fcjgS8WVKfqX4qlAB+vQh8c2Z3HSDmSl57DGwrI0AlCxZkpdeeiln/LFjpn5mu3ZORCsi8q/694dbbgEIBGZw0023snv37pwBGRnZT7h8Xe2w2gBsi1SXNzl/Sn6leBs3jiePjWHHP26aNvFw330p2bMmb775JqGhoTlj166FgAAYNcqhYEVE/t0rr0CVKh6gLXFx11CzZk0iV62C66+H8HCzZCsx0QyeORNef91sovMxDSIaALA1cqvDkYgv83M6ABEnrcpozURa4SKDGTWf49VXbTZv3ky9evUYmlXRISkJliyBPn1g3z4oU8bZoEVEThISAlOnurj0UoBxwBwGdujAshIlTNfJTp3AlTnftXy5aezz0UfwzTc5XeR8QMOyDQHYcmyLw5GIL9PMrxRbKSlwww3gsV3c128rQYOr8fzzzwPwbocOBE+ebJ4l1qxpamzu2KHEV0S81sCBcMUVNlAKeJPlts22H36AGTNMreCgIDNw2jSYNQs2bDBrJlJTnQz7nDQoa2Z+lfxKXij5lWLrlVdg0yaoWxee+qwhN02bRlpaGrf36kXXjz4yPUQ//9zMmPz8M9Sp43TIIiL/afJki9KlAQYAVzDy7ruJi4s7deAVV5iZ3zVrYMqUQo7y/NUJq4Ofy489sXtITEt0OhzxUUp+pVjaswfGjTPvp02DTZtWs2rVKipWrMiLX35ppoUTEiA6Gr78Erp2dTZgEZGzUKlUPOPHm/du95usWbOdAQMGkJj4L4ni4MFmuvjfkmMv5e/2p254XWxsNh/d7HQ44qOU/EqxdN99ZinvsGHQsyf8+uuvAAwYMMBscvP3hxIlimSLUBEpwvr04eZVN9OpE2RklKNkyUn8+uuvNGnShE8++eTU8e+/D507mxuij2hbuS0AK/etdDgS8VVnTH4ty5phWdYRK6v+kzn2tGVZ+y3L+jPzT7+CDVMk/yxcCF98YXLbV1+FjIwMZs6cCUCXLl2cDU5EJC/27sVlZ/DWW+B2Q1LStdSqdSm7du1i5MiRrF69Ovf4WbPgoosg0nc6pnWo2gGAFftXOByJ+KqzmfmdCVzyL8cn2rbdMvPPvPwNS6RgpKbCXSOOAvBElRlU/eB5ul14IZs2bQKgW7duToYnIpI3JUtCTAxNm8Ltt4PHY1GjxlfcdtvtALzwwgu5x+/caVq9lyvnQLDn58KqFwKwfO9yhyMRX3XG5Ne27SVAVCHEIlIwdu82ZcoOH+a112BrZDnql9rPvYFvkvzYY/yWORNy1513UqNGDYeDFRHJg7ZtYfFiSE7mmWdMid9ffrFo2/Y5AgMDmTNnDps3Z66VTUszm3q7doXAQGfjPgfNKjQDYEf0DvbH7Xc4GvFFeVnze6dlWeszl0WEnW6QZVm3WJa1xrKsNUePHs3D5UTOw8GD0KULrFzJwRW7efZZc/i2Z2MJ3LCGb15/HYAWERFMnjzZwUBFRPLBjTdCVBS88Qbh4WTf8557LpxRo24FTpj9HTfOdH677z6Hgj0/fq6cFgXN32ruYCTiq843+Z0K1AFaAgeBV0830Lbtt23bbmPbdptyPvRYRYqIO+4wbYl//pln5rfj+HGAudx3XxMOHTrEhz/8AEDv664zm9s8HiejFRHJm+7dTcZ71VUA3HorNGlictzSpZ/Gz8+Pjz/+mC2bNsH8+XDttabWr4+KStKDaTl355X82rZ92LbtDNu2PcA7QLv8DUskH+zcCXPmwEMPsbVEK959F1wuD/AQACEhIXzzzTcAhIWHw+rV0LixmS0WEfFVTzwBlSpBejp+wwbzWudZALz5eghXNu6Fx+PhqWefhblzTZL80kumrvmRIw4HfvYe7Pig0yGIDzuv5NeyrEon/Hg5sPF0Y0Uc8/PP5nXECB55BDIyoEKF7wDTE37+/PkABAUFMWrUKChdGrZuhc8+cyhgEZF8tHcvbN9Or2nDuJS5JCT74V4/gEA/Pz7//HNWtWhBap8+8NBDMHGi2fgGpiPcO++AbTsb/3+46YKbAKgSUsXhSMQXnU2ps0+A5UADy7L2WZZ1I/CSZVkbLMtaD/QA7i3gOEXOnsdj2nYGB8OkSSxbX4o5cyAoyMPBg7dkDxs6dCgAd911F1WqVIH69aF6ddPxSETE19WqBevXw4EDPP9hdSzL5mP/OxjZyySO7Y8coYTLRXiZMnw9a5bZHQdm4uCWW+Dmm712KVjtsNoEuAPYH7+f46nHnQ5HfMzZVHsYYdt2Jdu2/W3brmrb9nTbtq+xbbuZbdvNbdu+1LZtPScW7xEXBxdcAFddhX3PPTw4dCcADat8Bhyia9euBJ6ws/mhhx7K+Wzp0hAfX8gBi4gUoEqVaHJ1K66+2iItzeKI37jsUxkeD9ExMUycODFn/Pvvw2OPwfTpZkbYC7ldbuqG1wVgW+Q2h6MRX6MOb1I0JCXBa6+Zx3RlysC338Ly5Xx96bv8RmdKuyLZvPtOAN555x327t3LkCFD+Pjjj4mIiDDfkZZmdoVUq+bgX0REpGA8c8t+/P08zJsXBjQGYPTo0QCsWrUKO2uZg8sFY8ea1sfPPOO17Y/rR9QH4O/Ivx2ORHyNkl8pGm6/He69F37/3fzcpw+etu15cqHp2BbreZrU9Cj69etH/YoVKVeuHLNnz2bEiBE53xETA5dcAv3UsFBEip5aC97iloy3sG2Lli2/48knn8x+8pWYmMi0adNyBluW2QQXHw8//eRQxP+tfriSXzk/Sn7F961aZR7TPfootG6dffjrbyzWJ9UH9gFvA3Bv27ZmHdzJNad37IC//zatPpX8ikhRtHIljzWbS3Aw/PlnTfr3f4ZSpUpln7799tu5/vrrc2aA27WD/fvhssscCvi/Zc/8Rin5lXOj5Fd838yZpqXnCWt3bTunuDu8SA1S+axCBS6qUcMUgP/4Y3MqKcmsa2vfHoYOJbMQsIhI0RMfT6XyHu6+2/z41FMQHh7O8OHDs4fMnDmTDz/80Pxg2+b+mNkF09to2YOcLyW/4vvWrIELL4SQkOxD330Hf/wBbvdh4F3uvuMOhoWFQWiomfmdOxeaN4ewMLjpJqhXD375BU6YBRERKVKqVoVt23jgfpuSJeH778198tNPPyUuLo6ZM2cC8OCDDxIfH2/urU8+aZ6MeaETk1/bi8uyifdR8iu+r2JFaNMm+0fbhiefTAcgI+MFSpcOYtB998HGjebxXcWKZgNH7dqmA9yiRbBsmSl1JiJSVF1yCezeTcTO1dxqOh2T1ek4JCSEa6+9lvbt23P48GFefvllePVV81Stb1/nYv4P5UuWJzQwlJjkGI4lHnM6HPEhVmH+ttSmTRt7jWqoSgH7/vuse/VhGjbsx++/LyU4ONictG2oUQO6dIH/+z8nwxQRKVxJSWbzWr9+7D9gUbu2KXLz11/QoIEZ8ttvv9G5c2dCAgM5lJJCibFj4fHHnY37P7R9py1rDqxh6fVL6VS9k9PhiJexLGutbdttTj6umV8pGjIyIDoaMBV6jJdp3LhmTuIL5i7fvj1ccUWhhygi4qjgYOjfHyyLKn9+x3U9dmPbMH58zpBOnTrRtl494lNS+Lp1a3jkEefiPQta9yvnQ8mv+D6PBzp3hoEDWb44lWXLIDg4CXiLWrVq5YxLTTXre3/6CXr3dixcERFHeTzw9NM8uKAnLjL4cGY6ewaMzm5occ3ttwMwYu1afv36a7PmNynJyYhPS+XO5Hwo+RXf53LB3XfDb78xcdhyAAIC3gMSaNeunRlz6BBceiksXw5Tp2pjm4gUXy4X/PYbdf5vLMOrLyfd9uOV3zqYJj/AlTVqZA/tOngw1K1rNgSDqZYTE+NA0P9O5c7kfCj5laLhyivZ9fIsvjjSGTdpxMaOo3HjxgwZMsTU/61Tx/Srf/ddGDbM6WhFRJwVEAAjR/LwN50BmJF2LbFjJ8OaNZQbMiTX0KOTJ0PbtuaHCROgYUNzP/UCWcmvWhzLuVDyK0XG6wevwIObGqHzgQPcdNNNuN1uCAyEESNMtYcbb3Q6TBERr9G8OfToAQkJ8N57mMo5c+awdNGi7DG/VqkCZcuaH4YMgfBws6t4xQpngj5B9dLVAdgfv9/hSMSXKPmVIiEuDt55x7xPKmnedOqUufP3qafMjG+9eg5FJyLivcaMMa9vvJyIxwNcdhmtOnTIPr948eKcwa1awZIlULkyjBoF6emFG+xJIkpE4O/yJyopiuT0ZEdjEd+h5FeKhBkzTAv6li3jOHjwWypVqkTrE1odi4jIvxvYci812MWOAyWYP98c27dvX/b5yZMn5/5A2bLw8sumJfyCBYUY6alclotKIZUAOBh/0NFYxHco+RWf5/HA66+b99WqzQbgqquuMkseRETkP7nnf8sdTAFy7qX1TnpSdkpPgAEDzOuSJQUd3hlVDqkMwIH4Aw5HIr5Cya/4vJ9+gp07oVo1iIkxPel79erlcFQiIj5ixw5uDPqY4GCbBQtg61awLIs///wze0hGRkbuzwQGQmwsvPhiIQd7KiW/cq6U/IrPe/tt83rjjTYbN64DoFmzZg5GJCLiQ4KDCU87zNUjPABMMZPAtGjRgrCwMABi/q28WWgoWFZhRXlalUuZ5Feb3uRsKfkVn3b0KHz1lSlb2bv3XqKjo4mIiKBy5cpOhyYi4hvatoXq1Rk9+BAAH30EyZl7x8pmVnk4fPhw7s+89hrcfDOkpBRmpP9KM79yrpT8ik97/33Tm75fP9i+/RcAOnbsiOUFsxEiIj6hf3/Yvp2W/avQqpXpFP/11+ZU3bp1Afj77xOaSKxda9oeHz1qlj84TMmvnCslv+KbPB7sBT/wzvNHAOi86zFGjRoFQJcuXZyMTETEt7jd5vFZcjLXl5kDmAo6AA0aNABg69atYNvw+efQsyeULw9vveVUxLko+ZVzpeRXfFNyMktGvsXf0eWp7HeY6VsmZJ/q2bOng4GJiPiotWsZuWIMAaTwwwIPeyfPoX7mqR07dsDtt8Pw4aZm+pIlULGio+FmqVCqAgCHEw6fYaSIoeRXfMuOHaa2WYkSvNdtJgAdLo1kW2Zx8y9mzDD1fbdsgfXrHQxURMTHdOpExF9LGVR7AzYuPrh7DWUmTQIgPj4eBg403YRWrIDq1R0ONkdoYCgAx1OPOxyJ+Aolv+I7oqOhY0d44AGSkuDLH80NLzR0LgDjxo1j8PXXm0dzo0bB5ZdDYqKTEYuI+JYaNbhhShsA3qv2JMHPjgUgMTHRrA2+6Sbw83MywlOUCigFKPmVs6fkV3zHlClw5Ahccw3z5pmObq1bQ2zsauCEouyWBePHm+K/06Y5GLCIiO/p3RuqVIEdewPZFdIHgLi4OIejOr2QgBBAya+cPSW/4ju++go6d4ZWrfj4Y3No5EjYtWsXADVr1swZ2707tGsHn31W2FGKiPg0t9vcWwHWrWsIwIYNG07t8uYlAtwB+Ln8SM1IJTUj1elwxAco+RXfsX07tGxJTAx8952Z4B02zOaff/4BTkp+wSyR2LCh8OMUEfFxw4aZ1++/L0W5chWJiorKvtd6G8uyspc+JKQmOByN+AIlv+I7/PwgLY05c0xd9U6d0qhWzUVMTAwVK1bMLsaeLSgITm7JKSIiZ9S6NdSqBYcOWdSqdQ0A69atcziq09O6XzkXSn7Fd0yfDo88kr3koVmznFndbt26ndrYYvBg2LevEAMUESkazJM18/748X4AHDjgvXV0s5LendE7HY5EfIGSX/EdgwZxrGQNfvrJxs8Pmjffnn1q9OjRucdu3Ajt28ObbxZykCIiRUNW8rt7dxvA5dXJb0xyDABxKd67MU+8h5Jf8SnzpvyDx2PRo91xUlIOAtCnTx+6du2aM+joUbjySggPhzvucChSERHf1qoV1KkDCQmlgK4cPHjQ6ZBOq3mF5gCUK1nO4UjEFyj5FZ/y9XKzrvfSNU8R9c03ALRv396cTE+H//s/uOAC0wzj888hIsKpUEVEfNqJSx9gKPv373cynP+UtdHaL9KSAAAgAElEQVQtLCjM4UjEFyj5FZ+RnAzfLzX1HC/teIzIRYsAiJg7N2fQ/febGd9ffzX950VE5LxddlnWu/78+ec6ry13Fp0cDUBYsJJfOTMlv+Izfv4ZEhKgZUuo/vP7HL7oIgDC62d2n/fzM0nvH39AmzYORioiUjS0aQPly9tADY4ciWDPnj1Oh3QKj+3JXvOrmV85G0p+xWdkTfAOGmR2HX+7dCkAFzz9dM6gevXApX/WIiL5weWCvn2zKun0Z8WKFY7G82/iU+Lx2B5K+pfE3+3vdDjiA5QliE/weCBziS+XXgrTp08nOTmZwYMH07hxY2eDExEpwvr3z37HypUrnQzlX2nJg5wrJb/iE9avhwMHTL/5Vq1g/vz5AIwaNcrhyEREirbevcHl8gCdWLp0k9PhnCIqKQqAiGBtcJazo+RXfELm3jYuvhiio6NYuXIl/v7+9OjRw9nARESKuDKRO7iw4nbAj99XlybJzw+uv97psLJlJb/hweEORyK+Qsmv+IQffzSvvXrBQw89hMfjoWvXroSEhDgbmIhIUZaWBr16MejQDAAy6MeyoUNzykBERZkqOykpjoWo5FfOlZJf8W6ffUbqgMEsWZAIQMzdzXj33XcJCgri5Zdfdjg4EZEizt8fZs+m9/f3Zx7ozk+1a5udx2DK8EyYAEOHms0ZDlDyK+dKya94t7lzWbHWn0S7BI3DD/BSkunbPmnSJFq1amVmJUREJH999x08+iikpkKbNjTvVY5SpVKBmsybtzln3JAhMHmy2ZHsUDv5Y4nHACW/cvaU/Ir32b0bsjoJvfUWi27+FIByzY6wJzGRpk2bcvPNN5u1EG3bQkyMg8GKiBRBY8fC7NmmzRum5FnnzqbBxZYtFXKPvfNO6N4dxo93ZPZ3X9w+AKqGVi30a4tvUvIr3sW2zUaKvn3N+9BQflxkbr5//jkBgLFjx+JyucDthk2b4O67nYxYRKRoiYyElSvhuuvMsodMvXr5AZCS0iF3pzfLglGjYN8++OuvQg4W9sbtBaBaaLVCv7b4JiW/4l1WrzZryG66CSyLxERYtQosy0Ns7Fzatm3LoKy1Zj16wH33wYcfws6dzsYtIlJUZHVxa9Qo1+FevdwA2HY3jh49mvszPXvCa6+Z9vKFbE+sibd66eqFfm3xTUp+xbv88IOZRbjmGgDWroX0dKhY8SgQR8+ePbEsK2f8bbeZGeJ585yJV0SkqClZ0rzGxuY63Lw5+PkdB2oxadJXuT9TsSKMGQOVKhVOjCfYG5s581taM79ydpT8infZvx8iIiDMdOpZvtwc9nh+A6BUqVK5x9esCSVKaOZXRCS/1Klj7q3p6bkOu91wwQXHAZg+/W88J67vHT0aOnQwkxGFKC4ljtiUWIL9gtXkQs6akl/xLqGhEBeXXcUhq4384cNzARg5cmTu8ZZldhhffnlhRikiUnS53bB9u1l+dpJBg8oDcORILRZldR9aswbefx9at87eIFdYTpz1tQr52uK7lPyKdxk1ytxI/fyw7ZyZX1hOixYtqF27du7xa9aYZPmCCwo7UhGRostt1vcybRosXZp9uH37rLShLVOnToVly6B/f6hcGZ59ttDDXHNgDZCTBIucDSW/4l0aN4ZmzcDjYc/2VA4dgpIlk4FtNGvWLPfYjAx44AF4+mknIhURKdri4+GVV6BrVxg+HL78kjYhWzNPtmDOnG/5uVMnKFUKFi40S9YK2XVzrwMgKT2p0K8tvkvJr3gfjweGD2f5iMkAlCu3A4CmTZvmjElJgVtugcWLzc05a4OGiIjkj5AQ+P13ePBBsxl5yBBKt29Ig4qxQCDQnJ5Axpo1UL++IyEOrD8QgO41uztyffFNSn7F+7hc0L07q9aax26uuIUANK1Xz5xfuhRatIAZM+DJJ00tShERyX8hIfDii3D4sFni8PHHtG2bddK8+SurKZET4QWGAHBdi+sci0F8j5Jf8U533sm61jcAEBX1IwBN95kuPqSnQ0CAab/5zDOFvsFCRKTYCQiACy+EESNoe1HpzIMm+W3cuLFjYanMmZwPP6cDEPk3tg3rdpkbbAzrKRUURPWB5vEW3brBunVKekVEHHDizG/t2rVNx02HqLubnA/N/IpXOnTIdNiEGGAvN9xyC1atWuakZSnxFRFxSIsWYFk20JDYWOc2mnlsD/vjzJKLqqFVHYtDfI+SX/FK69dnv6Nx48ZMnDjRyXBERCRTiRJQs6YN+BMXVw67kBtbZDl8/DBpnjTKlihLsH+wIzGIb1LyK14pJ/ndwPbt20lOTnYyHBEROUGTJiZ9SEurx44dOxyJQUse5Hwp+RXv4vFAQgIb1puZBLd7I2lpaSQlqYajiIi3aNIk+x1r1qxxJAZtdpPzpeRXnJecDK+/bvrCBwVBqVKs/z8z9ZuR8Sf9+/cnwoHi6SIi8u9yCjw0di751cyvnCclv+K8jAx4/nlTwuyee7BfHM82d8PMk1u4+/rrYds2R0MUEZEcOcmvczO/e2L3AEp+5dyp1Jk4JzISwsNNd7Z166B8eQAOHYTEhwEiKV8+gF5z5sAdd8CqVVBNNzkREac1apT1rj6LF/8Gt90GZcuaWsC9e5u6wAUse+ZXyx7kHGnmV5xh2zB0KNx8s/k5M/EFyNk7sYOUlBSshx+G48dNO2MREXFcyZJQJSQaCABqMOfTT00nuEGDICqqUGLIXvOrmV85R0p+xRl//gk//wzNmp1yavt28+rnt5vY2Fh2BAXBE0/A99/Dhg2FHKiIiGSzbYiNBaBOs1KZB2swJiQET1wcLF8OFSuacT/+WKChaOZXzpeSX3HGokXmdcSIXIfT0tJYsMBkv+npWwD4448/4NprzYAffii0EEVE5CTvvGO6XBw8SM26/gCUKtWMffv2cTA6Oqf923ffmeUPn35aIGGkZaRxMP4gFhZVQqoUyDWk6FLyK844eNBUSj9huQPAs88+y6efrgLAz28P9957LwMGDDAzCWXKmM+JiEjhS001T+Fq1oSKFalZ0xwODW0OwD///JMz9pJLTCL8v/9BWlq+h3Ig/gA2NhVLVcTf7Z/v3y9Fm5JfcUZYGCQmQnx8rsMLFiwA6gLQsWN5JkyYQFBQkKkEceQIvPKKA8GKiAi//mruw/ffD5aVnfz6+9cDTkp+/fzg0Udh3z5YsiTfQ9GSB8kLJb/ijN69zeYIV+5/goGBgYC5mf3vf8NyTnz3HVSoYKpCiIhI4fvrL/Parh1AdvKbkVEVgF27duUe362bef3zz3wPRZvdJC9U6kyc0b69+QNmY4RlARAXlwCUx7Js+vRpYc4nJcFjj0Fo6InFJUVExEFZyW9CQjngpJlfgOBg85qSku/XVoMLyQvN/Iqz5s2DXr2y1/LGxgYCbsLCMvD3x9QCvuwy2LwZpk7FHBQRkULXpQtMmWLqnAFVq5p5i5iYkoD71OQ3Lg62boXRo/M9FLU2lrzQzK84KzkZVqyABg3gxhuJjzSF0StWtE2tyPr1zbrg6dOhb1+HgxURKcZatDB/Mvn7Q0QEHDtmARH88ssvREdHExYWZgY8+ih8/DEcOpTvoRxKMN9ZOaRyvn+3FH2a+RVnDR5s1vH26wdTppBw3NSNrFTJMt3fbrgBVq+G6693OFAREcG2zezvCy8AOQV7QkPNRuUxY8aYA198YSYtbrvNLFnLZ0cSjgBQoWSFfP9uKfqU/Irz6tUztSAjI8lwm0dYFSpkmHMvv5xrpkFERBy2YoWZ1R0xgnKlkgAYMeJuAKIOHYJHHoFhw0yr43HjCiSEw8cPA1C+ZPkzjBQ5lZJf8R4hIVSs1AoAj+eAw8GIiMgpLAtmzoRnn4U5cyi/6hsASsWYJWvuI0dMJZ9rrjFNibI2veWz7JnfUpr5lXOn5Fe8StmypppDdPQWhyMREZF/5XabZhc7d1K+o1nuEH3cJLnu6tVNG/qZM6FUqf/4kvMXnxJPdHI0Ae4AwoPDC+QaUrQp+RWvEppZtmbfvo0ORyIiIv+pcmXKXXwBADGB5t7tDgyEpk0L9LIbjmwAoHG5xrgspTFy7vSvRrxGVFQUGzaY8jWpqUcdjkZERM6knCnxS2ysWfbg51fwRaTWH14PQPMKzQv8WlI0KfkVr/HDDz8QHW02uo0c2d/haERE5ExKlzavycmmBntaWlqBX3PdIdPps0UFbYaW86M6v+I1mjdvDiQDMGBAV2eDERGRM8rsdwGY9b0HMxsWFaT1RzTzK3mjmV/xDhs30mjaNNyYaYRjA7vCc885HJSIiPyXrD1tHk8JAPbt21eg1/PYnuxlD5r5lfOl5Fecl5ICF12ENW0abqsMAOvKBUDFiuZ8air89ZeDAYqIyL/JmvlNj04B4OD+/XgmTIA1a0xDjHy2K2YXx1OPU7FURcqVLJfv3y/Fg5JfcV5gIMyZAwcOkG6ZTkAbW9SCm24y5198ETp0MJ3gRETEa2Qlv0l7oogA0j0ejtx/P7RtWyAt6bXeV/KDkl9xzrRppkuQxwMXXsjvu3bj8ZhNE1dfPSxn3HXXmWdrV14JhbCZQkREzsK771IyORKA4+HVCch8Wpe+dq25v48ale+XVKUHyQ9KfsUZHo9pe7l0qekYBPz002IALCuDPn0uzhlbvTq8+SZs2QKzZzsRrYiInGjBArj5Zkr93zQAElL8iIqOBiC8YUO45RYYMcKMfe89WL06Xy677rBmfiXvlPyKMzZuhL17zdKGzOR31ao/AHC7PaSkpOQeP3CgWQP87beFHamIiJzItuGhh6B+fQIfuQ+AlBSblJQUAgICCD6xpXFSEjz2GDzwQL5cWjO/kh/OmPxaljXDsqwjlmVtPOFYuGVZCy3L2pb5GlawYUqRs3u3eW3UKPtQ5crVAUhPT+a2227LPd7lMl2DCngnsYiInMHff5s9GGPGYAUHAbn3tmVkZOT8EBwM994LS5bAP//k6bLxKfHsiN6Bv8ufhmUb5um7pHg7m5nfmcAlJx17GFhk23Y9YFHmzyJnr4Qpi0N8fPahfv0GZb5LY+bMmad+Zv58WLy4wEMTEZH/sHmzee3QIevBHbZtUb16dVJTU1m4cGHu8b16mdfff8/TZTceMXNwjcs1xt/tn6fvkuLtjMmvbdtLgKiTDg8C3s98/z5wWT7HJUVdy5ZQu7aZ0c3UvPkFme/SAXIvfYiONq2EXnqpEIMUEZFTuFxQty6UKXNC8gtVq1YFYOPGjbnHh5oqPiQk5Omy2fV9K2q9r+TN+a75rWDbdlYbl0NAhdMNtCzrFsuy1liWtebo0aPneTkpciIiYPt26N49+1BGRtY/R5P8rj5xg8S4cZCYCJec/BBCREQK1aBBsG0b1KmTK/kNCjJLIBqdsJwNgKgos/yhWrU8XTZrs1vz8lrvK3mT5w1vtm3bwGkrWdu2/bZt221s225TrpwKUssJLAsyMuCJJ2DfvuxJ4BIlTMugr7/+2txR33gDXn0Vbr0VmuumJyLiFY4exbI9ANi2zfr1mZvRTr5Pf/klTJkCXbrk6XKa+ZX8cr7J72HLsioBZL4eyb+QpFjZvBkmTYKmTfGfMB4AP7eZPZg7dy6MHg133QUDBsBrrzkZqYiIZNm82cz8jn8RMPMUx44dIzQ0lGonzvB++SWMHw87doCf33lf7sS2xqr0IHl1vsnv10BW9epRwNz8CUeKnWbNzCaIrl0JeGUcAHZ8MgB///03iRdcYJLer74yneBERMR5jRrBoEHY48x9m8wZ4Li4OD755BNT4mz8eBg+HNq3h8cfz9PldsfsJj41ngolK1C+ZPm8Ri/F3Bl/DbMs6xOgO1DWsqx9wFPAi8DnlmXdCOwGhp3+G0TOoF49+PprArbvg3qQ7A6GzEo56cOH52yWEBER72BZMHMmyeH1YDL4p8SSnHnqqquuYmSlSnDwIAweDDNmQOZ64PO14cgGQLO+kj/OmPzatj3iNKd65XMsUsz51zI7hdM9Of8scxVLFxER7+F2k/jAkyb5DbLJyn4rVapkavu2bw9du+bLpTYcNslvs/LN8uX7pHg7/wU4IvnM7TZ/MjIswI1lefB4PE6HJSIip5GUZF7tYLKT3ylTpsDll+frdbJmfptVUPIreaf2xuI9bJsAf5PsNq1ZH9u2Ty2WLiIiXiMr+Q0IyOnq1qdPn3y/Tnbyq5lfyQdKfsV5GRkwdSrUq0dQcgwAobsOAxAYGAixsaAZYBERr5OV/EZEmK6dZYODKfHMM/D22zlt7PMoJT2Frce24rJcNC7XOF++U4o3Jb/ivG+/NSXNKlSgdIRZibMjvDYAdevWhVGjzCO05OT/+hYRESlkSQlmYsLlMvfnNsnJpnzlrbdCrVpw9dV5nrzYcmwLGXYGdcPrEuyvfSCSd0p+xXmDBsHKlbB0KWWqmcoOh6My8Pf3p1rVqtCjB3zzDdx4o8OBiohItshIEm+7FwA/P9OZs0TWRMXWrfDwwyYBduUt1dCSB8lv2vAmzvnuOyhVCrp1g3btAChTJutkGWrWrImfvz/cfbdZ+vDUU3DDDdBLhUZERBx3441E/RMGQIkSKQD4+fmZMmj168Pzz+eM/ftvqFzZ3PPPkSo9SH7TzK84Z9w4073NzumOfWLyW7du3ZyxDz0E5cubdWQiIuKslSth7lwi+14NQMmSZtmD2+0+deyRI5DVsOg8rD9iOrup0oPkFyW/4ozUVHPzHDjQzBJkOm3yGxgIF10E27YVapgiIvIvZs+GwEAiG3UGwN8/FoDQf2tKVL48dOxoml2cB838Sn5T8ivOiIszmyAqVMh1OCgoa1NbGQYNGpT7M5Mnw+LFhROfiIic3tat0LAhx+JN23nbPgbAtGnTaNWqFaNGjSImJiZn/MUXw86dEBV1TpeJTopmf/x+gv2CqR1WO9/Cl+JNya84o0wZCAiAPXtyHQ4ISAAgJKQGvU5e23vNNXDhhYUVoYiInM5dd8HEiURGmh+DgxOyT/3555988MEHbN68OWd8mFkbTHz8OV0ma7Nbk/JNcLv+ZUmFyHlQ8ivO8PMzbS/Xrs11uGTJ4wD4+9fIPf7AAVi0CHr3LqwIRUTkdHr3hh49OGYmfLn00k40b94815C2bdvm/LBjh2nhWb78OV1m3aF1ADQv3/wMI0XOnqo9iHPeeguqVs11qFQps27MtivlHExPNzUjwcw2iIiI8xYuJHJHGyCM5s0rM3/+fKpUqQLANddcg7+/vxmXlgaXXAJNm0LwudXp/fPQnwC0qtQqPyOXYk4zv+KcOnXMRrboaFMSJzGRwEDzDC09PXMtcFISXHaZaYQxYQLU1povERGv8NprHNhj6vtWqEB2O/qBAwfywQcf5Ix77DGz5rdNm3O+xB+H/gCgZcWWeY9XJJOSX3He55+bm2ONGvjPNqVwUpMyyz4EBZm1wW+8AXfc4WCQIiJyoqQpMzhGOfxJpeKs11nx228AdO5sKkBw5IhpTvTyy3D99ab27zlISU9h45GNWFi0qNAiv8OXYkzLHsR5t94KjRvD5MnY384D0khJL0PKsXgCy4bAF1/kKocmIiLO25ti1u9WDY7Edc8Ysqqwx8TEmPv6e+9BRgY8+iiMHXvO37/hyAbSPGk0KtuIkMCQfIxcijvN/Ip36NIFZs2izcIFwEEAfllrXpX4ioh4nz1bEgGo3q4inm7d8GQe/2v2bChZ0nTn3LzZNDQ6jxbHK/etBKBN5XNfLiHyX5T8ilfp1Lkz5cubNWR33z2etLQ0hyMSEZFTTJvGnisfBKB6dQvXnXdmn2q3fbtZqla69DkvdTjRon8WAdC1Rte8xSpyEiW/4nU6d64GwNatqaxcudLhaEREJJfp0+G229hbpQMA1asDV1xB+/btAejw6adwxRVmQ/N5PrnL8GTw866fAehVq9cZRoucG635Fa/TqFFmeRzqUa1aNUdjERGRExw+bJYzXHQRe6qOhO0m+c3IyGD9+vUAtLjoIhg6NOcz6emmtvs5+OPQH8Qkx1A7rDa1wmrl599ARDO/4n3q1ct+R/A51oQUEZEC9MEHkJAAb7zBzl0mhahRA/7++2+SkpKoUaMG4eHhZsbXsmD+fGjYEI4fP6fL/LjzR0CzvlIwlPyK1zkx+Z08ebKToYiIyImioqBXL2jQgC1bzKGGDU1LY4CWLU+qx1uqlOnu9u2353SZrPW+Sn6lICj5Fa9Tt27Wu3q8+eZU4s+xF7yIiBSQF16AH38kJgYOHTIN26pVg9WrVwPQqtVJndg6djRrf3///awvkZyezNI9SwHoWatnvoUukkXJr3idcuUgNNQGShMd7eKtOnVy2huDueOKiIgzbJutW2wAGjQA8DB79mwAevY8KVl1u83sb0LCWX/9sr3LSE5PpkWFFpQrWS6fghbJoeRXvI7112bqpW3O/Kker8TFkda8ufnx2DHT4viee0y/eBERKTwrVkCtWmz5bgdgljwsXryYvXv3UrNmTTp16pR7/NGjEBmZWRLi7PT6wCx1KBNUJt/CFjmRkl/xPm++SQOPWUxWtmx3jqSksOHCC825gAC44QZ47TW49lqwbQcDFREpZpo0gYQEtsxYBpjk97777gPgqquuwnVyM4vjx+Gyy6Bfv7P6+nRPes6lyjXJn5hFTqLkV7zPxIm0eKA3ABERPQBYvny5ORcaaoqnv/ACfPopfPihU1GKiBQ/ISEwZQp/HQgFoGGtlOzNbh06dMg9dt06qFoV5syBZs3O6uu/2fpN9vvX+72ePzGLnETJr3iP5cth61bw96dFF3NjTU83v/mvWLEi99gHH4TWreHFFws7ShGR4m3YMDaX6wZAo8cGZx9u0yazDfGePfDYY9C2LTz66Dl99eRVpsLP5Esm47KUokjBUJML8R4vvAB//QXbtpFVLefQoQrAvyS/Lhc8/LDZQZyWBv7+iIhIwYuNhW1Hwwjw91C/RyVc/+fC4/EQMXo0LFli1vgCjBgBjzxy1t+7/vB6ftn1CyEBIYxqOaqAohdR8ive5O+/IbNMToUK5s/hw24CAxuyffsWjh07RtmyZXPGd+xoiqgfPw5hYQ4FLSJSvPzxh3lt3sJFyuuv4vlwOqVKlcK/USNz427QAAYMOLFu5Vl5faVZ5nBdy+sIDQzN77BFsumZgngPlwsyMrJ/bNHCvNaqdTlAdimdbD/9ZPrHHzxYWBGKiBR7WSV7W7eG6OhoAMLCwmDcOJg61VTjOcfENzIxko82fATAne3uzNd4RU6m5Fe8R5MmpoxOZgKctfQhNrYmALfffjv79+/PGT9/PpQpA/XrF3KgIiLF19q15vWU5DcP3v39XZLTk+lbty/1I3RPl4Kl5Fe8x5VXwoEDsMyU0MlKfsuWvSR7SP/+/cnIyDBTD599BtddB35avSMiUliykt8LLoCEzOYVJUuWPO/vS/ekM2X1FADGtB+T5/hEzkTJr3iPyy+HhQuhSxcAskr77t9ThU/7D8ANrFu3jr9//BEGDjSzvo8/7ly8IiLFTHy82Z7h7w9Nm+bM+MbExJz3d87dMpe9cXupH1Gfi+tcnF+hipyWkl/xHi4XXHSReT93LjWev5VKJWKIinXT8oddlM+s6BBStapJfCMjzWzx0aMOBi0iUnz8/rvpLdSsGQQGQkhICADHjh077+/MKm92V7u7VN5MCoX+lYl32rABa/q7dEz8EYDf0toQmdnOOKJ2bbPdeMoUWLoUevY00xEiIlKglv5qumpmPZl79oYbAKh/nnsv1h1ax5LdS0x5sxYqbyaFQ8mveKfHH4elS+nYOhWAZa3HEFyiBAA7d+40bY5Hj4avv4bNm01BdRERKTg//cSS8abbZteusHLpUqYvWoQLmJaSAr/9ds5f+dyvzwFwQ6sbCAkMyc9oRU5Lya94rwsvpOMbIwFYntSKeo0bA/DKK6/kjOndG264Ad5+21ReFxGR/PfOO6T36sOyhOaA2ZqxftMmAMoEB9Pk4EHo3BmeesqsizgLaw+sZfbm2QT5BfG/jv8rsNBFTqbkV7zPsWMmod2wgVatzLqyzZshMtIDwMyZM0lKSsoZf++9ZvzOnQ4FLCJShK1ZA7fdxp8dR3PcLkXdulCpEjTMnJAoW62aaU1/3XXw2muwY8cZvzLDk0Gbd0w75Dvb3kmV0CoF+TcQyUXJr3iff/6B996DXbsIDISsdvE1alyVPeTll1/OGR8ebgqrn9wCWURE8u6ZZ6BsWX7t/yKQXZCH0FDThc2yLChZEmbMgE2bzqrBRdC4oOz3D3d+OP9jFvkPSn7F+2RWdSAlBTD72QDCwoZmD/n2229zxmctd8hDnUkRETmNp5+GTz9lyepgwKz3Bfj1118BaJXZlh7LgipVID0dFi8+7dfFJseS7knP/jmiRESBhC1yOkp+xfvUr282tGVunujd2xzesqUae/fu5YorrmDixIk54380FSFo3bqQAxURKQZat8bTrQeZuW72zO/3338PwCWXXJJ7/EsvmVmL09T+LTO+TPb7qAej8j1ckTNR8ivep0QJGDTIPEI7eJAOHaBUKfjrL7CsqsyaNYtOnTqZsfHx8PLLZm1E5vozERHJJ+np8NVX/P7xFiIjoXp1qF0bli9fzoIFC7Asi4svPqkxRcuW4PGYzRon2R+3P9fPYcF5a4sscj6U/Ip3GjsWKlaEY8fw94fu3c3hrEnebLt3mwR4wgTzyE1ERPKPywU33MD8V0xlh759ITk5iWHDhpGens69995LpUqVcn8mKHM9b2rqKV9XdWLV7Pcpj6cUWNgi/0XJr3inBg3MrEGzZpCaykUllwGw8PsM2LPHbIhLTzf9Nbdvz3kOJyIi+cflgiuu4Pv1lQGT/K5fv559+/ZRu3ZtXnzxxVM/s2GDea1VK9fhg/EHs9+XLVGWAHdAgYUt8l+U/Ir3crvN65w59P7sJgB+/PQodo0aprTZokXmfJgem4mIFJSo2x5lhUmmHwgAACAASURBVN0OfyuNnl3SSM2c0a1UqRL+WRuUs9g2zJwJLVpAjRq5TlWeUDn7/a67dxVw1CKn5+d0ACJnNHw4jRo3oXLXBA7EVGTj45/RbEh9c3MVEZEC9eP2mniAbvbi/2/vzuNsLP8/jr+uM5sZy9hlyb4UZR1LyE72LIOkRWRJKOWLfigJlaUihKRCRYTsy7ShslP2fcnYl8HsM+fcvz/uSSaGLDNnzHk/H4/zOOdc93Xf87nvx3Hm45rr/lxkfPlTjjZpAoCfn9/1nY2BOXPgzJlEzU6X8+rrfJnykd5X1XnEfZT8yn3BPPoIfQbaMx2yPdcW8tx6HxERuXvLltnPDVulJ6pZMwYnLCffskwZOHECIiLs6jzLl8PMmXad33/V+p29c/bV1wd733oRDJHkpORX7ht9+7o7AhERz+Jy2TktQMO3KjNqwQqOHDnCo3ny0P3DD+HaspPFitk3IRcpkugYsc5Y3vzpTQCmNpuqub7idkp+RURE5IbWrYNTp+wSZ7lyneHdd98F4OOJE/F2OCA0FPz9oXRpeyqa4/pbiaZsnsLBiwd5KPtDPF/2+ZQ+BZHrKPkVERGRG5o7134ODoZ1634nOjoagKB69f7TqppXYq4w9JehALxb9128HUo7xP1U7UFERESuY1nw3Xf26+Bg2LBhw9VtEydO/E/HGP3baM5GnqXqg1V5ssSTyRGmyG1T8isiIiLX2bTJLqueJw9UrgwzZsy4um3p0qW33P9U+CnG/D4GgPfrvY/RQkSSSij5FRERkev8PeWhdWt7Km+lSpWubnO5XLfcf+gvQ4mIi6B5ieZUz189ucIUuW1KfkVERCSRf095iI+PZ3vCym2NGzfm22+/ven++87vY8rmKTiMg3frvpvc4YrcFs08FxERkUS2boWDByFnTqhWDWbMmMm+ffsoUqQICxYsuH5lt38Z+ONAnJaTF8u9SMkcJVMoapH/RsmviIiI2OLiYM8epk9/FIC2hTcR/+5y3p40CYAhQ4bcMvFdf3w9c3fNxd/bnyG1hiR3xCK3TcmviIiIwNq10LEjcdFOvo49BBie50s+GDyeI0DJTJloX7fuTQ9hWRb9QvoB8GqVV8mbKW+yhy1yuzTnV0RExNP9+CPUrQvGsOz5WZw9ayhZEnLO7scwf38AxkZF4VW5Mhw4kORhlu5fyuqjq8nqn5X+1fqnVPQit0XJr4iIiCe7cgU6dICiRWH9er7cUxmA55+Hfv3+R2RUFMHBwdRbtw6yZ4ewsBsexulyMuCHAQAMenwQgekCU+wURG6Hpj2IiIh4spkz7TWM583jvJWVRYvs0ma1ah2nf//Z+Pv7M2bMGHuN482bIYl6veM3jGfHmR0UzFyQHhV7pPBJiPx3GvkVERHxZC+8YE97eOwxZs2y73mrXx927/4BgAYNGpA/f367rzF2ojxqFDidVw9x/PJxBv00CIBxDcfh5+2X4qch8l8p+RUREfFk6dJB7dpYFkybZjd17Ag//fQTADVr1kzcf9Uq6NcPdu262tR7WW/CY8Np+VBLmpVolkKBi9wZJb8iIiKebMUKmDyZ9ethyxbIlg0ee+w0c+bMAeyR30QKF7afT54EYNHeRczfM58MvhkY12hcSkYuckeU/IqIiHiykBDo3ZsJH8YC0LkzjBjxJpGRkTRv3pxSpUol7n/xov2cIQPhseH0XNYTgGG1h5EvU76UjFzkjij5FRER8WQdOnAmNpBvv3NgDJQu/RtTpkzBx8eHESNGXN9/+XLw9YXSpRny8xCOXTpG+dzl6VmpZ8rHLnIHlPyKiIh4srJl+azMx8Q6vWkYdIrBg58BYPDgwdeP+sbF2dMk2rdnW/gBPlr3EQ7jYErTKXg5vNwQvMjtU6kzERERD+Z0wqTzbQDw3dSJw9ZhypQpw4ABAxJ3XLcOSpWCzZtxRkXSbdGTOC0nvSv1pkKeCm6IXOTOKPkVERHxYIsXw7HjDvLlieD7Eyvw8vLi888/x2fMGLtDeDisWQOrV0O3bjBpEpP3zmRD6AbyZszLO3Xece8JiNwmJb8iIiIe7MMP7efL4aMBFwMGDKRcmTLQtCmcOGGveFGqFIwYAT17ciTsCANC7FHhcY3Gkckvk/uCF7kDSn5FREQ81Pr18MsvAGFcvjyGqlWrMnjwYDvhDQ2FqCjw8QFvO11wupw892UTrsReodXDrWj5UEu3xi9yJ3TDm4iIiId67bVTCa8+oUePZ/npp5/w87tmdTZ//6uJL8Do30az5tgaHsjwAJObTsYksdSxSGqm5FdERMQDTZv2K7/9lhOIpnPnCMaPH4+vr2+S/bee3Mrgnwbb+zafRvaA7CkUqci9peRXRETEw8TGxtKz51HAQcmSm/j003duOoobFRfFM/OfIc4VR4+gHjQq1ijlghW5xzTnV0REJC1zuew5vAAHDkBMDFN/3k9UVDDgYt68qrecvvDGD2+w6+wuSmQrwagGo5I/ZpFkpJFfERGRtMjphHHjoGBBOH/ebps6FR55hOGvHAZ8KVfuECVK3DwVWLp/KWPXj8Xb4c3MVjMJ8AlI9tBFkpOSXxERkbQmNhZat4ZXXoFixexavQDdu7PytZGccHYDYCzvw6VLSR7mwIUDdJjXAYChtYYSlCco2UMXSW5KfkVERNKaQYPg++9h7FgICYECBez2ggV5acFDQABFcqzl8e1fwJtv3vAQV2Ku0Gp2K8Kiw3iyxJP0r94/xcIXSU6a8ysiIpKWhIbaK1d07gy9eyfa9M03azl0qAHgYsqsEpDpdyhT5rpDxLviaTe3HdvPbKdEthJMbzkdh9F4maQNSn5FRETSEm9ve+T32WcTNYeGhtK9+19AdR59dBd16pQEctgbY2PBywu8vLAsi97LerPswDKy+Wdj8dOLtYqbpCn6b5yIiEhakisXvPUWFC58tSk0NJRixRpy+XIw4GTGjCL/9N+8GXLnhg0bAHhn9Tt8sukT/Lz8+P6p7ymatWgKn4BI8tLIr4iISFpy4QIcPgwVKlxtypAhA1FRfQEfypX7gzLXTnUIDLT32bePjx2beOvnt3AYBzNbzaRa/mopH79IMtPIr4iISFryxRcQFATHj19tWrToKPAMEMe0aYUS97csAL6K+J3ey+05wlOaTiG4ZHDKxCuSwpT8ioiIpCVNmtjPn30GgMtl8eqrLsCLcuXWU7bsv+bvbtzID4Wg47mpAIyqP4rO5TunYMAiKeuupj0YY44AVwAnEG9ZlgoAioiIuFOJEtCyJbz/PgQH88qkUM6fb4AxF5k+vUjivk4n2z8dTqv2hnjLyeuPvU7fqn3dE7dICrkXc35rW5Z17h4cR0RERO6F8eOhQgW+qVab8ZfWAtChw0EeeeSaMar4eI5fCaVx3VNcdlq0KdmGkfVHuilgkZSjG95ERETSmjx5mNC5Mz2HXwGKkznzGaa+5g/790NMDKxfz+WJH9Kk6RmOOy5QPX911fIVj3G3ya8FrDTGWMBky7Km/LuDMaYr0BUgf/78d/njRERE5FYmT55Mz+ETgQMAfPllDvyerwfbtwMQ54DgLgH86YikRLYSfP/U96TzTufGiEVSzt0mv9Utywo1xuQEVhlj9liWtfraDgkJ8RSAoKAg6y5/noiIiNzE/Pnz6dGjB/ARkJW6daFZMwNZJ8Lhw1je3nSNnsWqYwvJmT4nyzosI6t/VneHLZJi7urvG5ZlhSY8nwHmA5XuRVAiIiJy+/r06UOrVq1wucphTE+8vOCDD8AYoHp1ePZZ3s69ly+OLSTAJ4AlTy+hUJZCtzyuSFpyx8mvMSa9MSbj36+BBsCOexWYiIiI/HcbNmzgo48+wsfHn9y5F2JZhj59oHTpf/p8vvVz3v7lbRzGwezg2QTlUZEm8Tx3M/KbC1hrjPkD2AAssSxr+b0JS0RERG7HqFGjAHj88VmcPJmH/PlhyJB/tq88uJKui7sCMKHxBJoWb+qGKEXc747n/FqWdQgoc8uOIiIikqwOHDjAd999h7d3Qdavt5Pa8eMhfXp7+7ZT22j9bWviXfEMqDaA7kHd3RitiHup1JmIiMj95PJl8PEBf384fBiWLGHMr79iWRZ5887l6FEHLVpAs2Z29yNhR2j8VWPCY8Np/0h7htcd7t74RdxMBf1ERETuBy4XjBwJefPC3Ll2265dnOnViy9mzQKacfRoBTJkgHHj7M3nIs/xxMwnOBl+kloFa/H5k5+rlq94PP0LEBERSe1cLujYEfr3h7p1oVw5u71xY8b37k00WfFjMgDvdA/lwQchPDacJl83Yd/5fZTJVYYF7Rbg5+3nvnMQSSWU/IqIiKR2kybBjBnw9tswfz488ggA4RERjJ8xA/iYGHJTzXcjvXb3IM4ZR/C3wWwI3UDBzAVZ1mEZgekC3XsOIqmE5vyKiIikZrGxMHQo1K4NgwcnFO21TZs2jYsX6wBPExBg8cWiXJjHvqHTwk6sOLiCHAE5WPnMSnJnzO2++EVSGSW/IiIiqZm3N8ybB5kyJUp8AUaN+hJYnvDaULROfvqu7MvMP2eS3ic9SzsspVi2Ym4IWiT1UvIrIiKSmjkcULXqdc2rVoVw/PggIAcVKlyke/csjPltDGN+H4O3E+bVHKtFLERuQHN+RUREUrNjx2DqVLh48WrT2rVrad78W6AlPj5RzJmTiXl75tJ3VV8AvlwADdKVclPAIqmbkl8REZHU7MgR6NIFfvkFgGXLlvH44y8QHT0GgIkT/Tjts5Fn5z8LwHuZWvP0diBbNjcFLJK6KfkVERFJzapUgRw57IoPlkW/foOBWUBGGjcOp06rozw560mi46PpUr4L/X6MhgcegCJF3B25SKqk5FdERCQ1iImB6Gj7dXy8XdsXwNcXBgwgfMUK/hwyhN27XwAqkDPnFSZOi6fpN004E3GG+oXrM8E/GLN4Cbz0kj1XWESuo38ZIiIi7nTqFHTrZk9TCAmx20JC4MEHYdgwzv/1F0+tX09GoMzQbTidL2NMHPMX+tFpRWt2n9tNqRylmON4Cp8LYdCvH/zvf249JZHUTNUeRERE3GXbNnjiCQgLg2efheLF7fbcuaFMGZYPHkynt9/mZHw8Xl6FsJyf4wIGDrrCZxva8eOFH8nlCmDJl3EEbuhsr/62atV1JdFE5B8a+RUREXGH8+ehcWNIlw62brUrOvyd/JYpw3edO9MIOBkfz2O+GSjzyJ+4yErjxhBQfwrTLoTgHweLpsdTwDcnfPYZLF+uxFfkFjTyKyIi4g4//ghXrsCaNVCyZKJNMTEx9O1rly1747nnCJtZhU/+yEDevNBi0Gy6rnwDg+Grp2ZT8Z1gJbwit0EjvyIiIu7Qpg2cPg1lyyZqdjqddO/enSNHjlCyZEkKVpvGJ66X8PWF1yct5eUQu6TZqPqjaPlIGyW+IrdJI78iIiLuEhCQ6K3T6aRTp05Mnz6dgIAAevf+ml69vADo/dEPvLGtFXGuOF6r8hqvPfaaOyIWue9p5FdERCSlHTlij/guXZqoeezYsUyfPp306dMzY8Yqhg4tQ2wstOy9moknGhDjjOGloJcY3WA0RiO+IndEya+IiEhKy5IF/vgDtmy52rRnzx7efPNNAKZPn8UHH1TlxAko23gjIbmaEuntoqNXBcY3Hq/EV+QuKPkVERFJaYGBUKECzJ0LlsX+/fupU6cOERERBAe3YcmSpvz6K+Qs/QdHHn+CK3FXaLcDppYehMPoV7fI3dC/IBEREXd4+WX44w8Ojx5NnTp1OHnyJLVq1aJUqZlMmwZ++XYT91R9wmIu0uJCLmb8GIhXvQbujlrkvqcb3kRERNzhuedwfv45hfv1u9r0zDPLePFFX8h6gIDudbkYe5aG6R5l1oTt+Iwcc90NciJy+5T8ioiIJLdjx2D8eHvZ4qpV7ddeXpQ8dOhql76vL+Sl7r6Q/gxZe9flQvxJakfmYl6Hb/GLmwuvvurGExBJOzTtQUREJDnNmAEPPQQffgjZskHp0gCc2rqVfaGhAAT4l2fq5EbEuZzkbVOVCxyj8nFY+GUs/gePwqBB4NCvbJF7QSO/IiIiyeX77+G556B2bfj8cyhQ4OqmMF/fhFe5yBYzj79c3hTp2JWDBQ+S22RifrNPyDCuJfj7uyd2kTRKya+IiEhyKVsW+vSBESMgXbpEm374+WcgCxkCfuWvyAIULv8eBwt+io/Dh7kdl5H7wapuCVkkrVPyKyIiklwKFIAPPriuOTw8nPHjvwSWER5ZhIJlfudkozcAGNtwLFWV+IokG00gEhERSQ4rVsCCBdc1h4WFUa9eU/bseQ+oTL5i53A914EoH3ih7At0D+qe8rGKeBAlvyIiIsnhyy+hb99ETS6Xi4cfLs369a8DdUgXcI5CfZ/h2JXDBOUJYmKTiVq9TSSZKfkVERFJDlmywPnzYFlXmzZv3sapU+8DzYDzlH6lI2tOriB7QHa+KzOCdM91gn373BayiCdQ8isiIpIcKlSAsDDYuhWAY8dCqVHjINAeuEy9Ya+zwW8JDuNgdvBs8h86D998AzExbg1bJK1T8isiIpIcWrSwy5SNHk14eDTlyu0kOroNDkcU7SYPIyT+SwB6V+pNnUJ1YO1auyJEsWJuDlwkbVPyKyIikhyyZoXXXyd20x88UupPLlxogDFXGD/7D2afHHW12wdPfADHj9tzhFu1uq4kmojcW0p+RUREkklYz9coFf81R49VAi4xcso6hh1rfXX7X33+wjid9iixZcHQoe4LVsRDKPkVERG5xw4dOsTHH08hd8HNHDhcBrhIn4ovMG5/E05cOUGNDKWILPkV+TbuBW9v6NwZ5s6FIkXcHbpImqdFLkRERP6L2Fi7bm9ICBgDkyfb7U4neHld7TZ79mzat++GZc0F6uHldZFhE1Yz/uzPhDrjqHoMFn21E/+YDlClCtStCy+95J5zEvFASn5FRERuZf166NABDh60S5jVqmW3WxbUrg01asBbb7FgyRKefroPlvUjUJ5MmSL5YvFJem7swQnnRarnr87SF6eTsfs5e8S3eHF3npWIR1LyKyIicjMbN9oJ7gMPwOLF0KgROBJmDcbGQsGCMHw4c5csof2OKFyu1UBRiha1+PCbfbz4c0NOR5ymRoEaLHl6CRl8M0DeQu48IxGPpjm/IiIiN1OoEDzzDPz+OzRp8k/iC+DnR9xnn/G/Bg1os80QH/8LUJQKFSze//Znnl5Vg9MRp6lbqC5Ln15qJ74i4lYa+RUREbmZ7NlhypQbbjp58iRt27Zl7Vo/4BcgI/WqXOa58Stpv6QDsc5Y2pZqy/QW0/Hz9kvRsEXkxjTyKyIikpRvv4WFC2+4ac2aNZQvX561a0sAy4CMPJ1/Lc1eGcbzi9sS64ylZ8WefNP6GyW+IqmIRn5FRESSMnUqXLoEzZtfbYqLi2PMmDEMHPgWLtd7QB8AXuvr5IMMj/P1XrvfsNrD+L/H/w9jjBsCF5GkKPkVERFJiq8vxMRcfRsSEkKvXr3Ys+cEsABohI+PxfsTTrPYrwMctvu9VfMtBtYY6JaQReTmNO1BREQkKQ8/DLt347x0if79+1O/fn327InFx2cz0Ihs2eDdb0N471IZfjz8IwDD1vgypNYQt4YtIknTyK+IiEhSmjfn/OjRPFWlCiF79uBw1MXPbyFRUQGUKgU72xj6/mF3rVmgJt98cYXcUfrVKpKaaeRXREQkCVvTpyfIz4+QPXtJn/4dYBVRUQE0ahrLzjaJ5/KGFHuH3L9sgdat3ROsiPwn+u+piIh4hrNnYc0aOHnSns5Qp85Nu8+cOZMuXboQHeNPYOAqLl16HIBe/3eCzYXbwvF/+rqeO4SpUwfy5dNSxSKpnJJfERFJ2y5dgn79YNo0iI+327p3t5Nfy4L9+xMtM3zmzBly5cqV8K4cGTKs4NKlHGTJ7KJ/oWZ86P0zp49HkicgF1+UfpP61Z+3V34LC4NVqyBjxpQ/RxH5zzTtQURE0rbgYPjsMzvhXbcOTp2CESPsbYsWQalS8MUXAEyfPp2iRYsm7NgFb+8NhIfnoHwFi5cnD2dg82WcdkRS+zBsffs09Ru+DPPmQbt2dhIdFOSWUxSR/04jvyIikrZNmgTHjkHt2tdvq1EDatcm7oUX6Dt/PuMWLgQyA18D7YmPh2d7nORslRcZtnspGOhfvjfDytXGu/o5CAyEatXsY2XPnoInJSJ3SsmviIikTfHx4O0NRYrYjxtwZcrEmtdfZ8DataxbuBAvr7pkyrSAixczkD49dBr1HV9f6cb5Q+fJki4LU5tPpdXDrVL4RETkXtK0BxERSZveegsqVwan84abf/vtNypXrkythg1ZF+UkPe/jcq7k4sUMVKh2mQYTn+fjM8GcjzpPgyIN2P7SdiW+ImmARn5FRCRt2rEDIiPByytRc2hoKP369ePrr78GIEuWqvj6zuH06Tw4jIuOg3/j56zPsPnwYfy9/RlVfxQ9KvbQMsUiaYSSXxERSZtiYiAgIFFTeHg41apV4+jRo/j6+lOjxnesXduQixcNBQvHUeutd/jyyHBcYS7K5y7PV62+4qHsD7npBEQkOSj5FRGRtClPHtiyxS5nZgxRUVH06tWLo0ePUrx4K/z9vyYkxA+Alp0P8FfFZ/ji8HoMhv7V+jO09lB8vXzdfBIicq8p+RURkbSpRg34/HNYv55N3t60bduWw4dDMeZtDh0aRHy8gzx5LdoM+4KpJ3oRcSqCfNG+zOi2gloFa7k7ehFJJrrhTURE7g+WBfPn28sHFy9uV3AIDrbbb6RVK2jZknmbNlGjRg0OH86Bn99OLOtN4uMddHzpAkHvt2Hs0U5ExEXQbpeDPyOfV+IrksZp5FdERO4Pzz8PM2ZA3rxQtapdxixfPvj7RrSjR6FAAQBOnDjBjBkzuPzww4zoNRgYCfQgJsZBsWLQfeQPfHDoeUIPhJLRNyMTTpTlmblrMdtfddvpiUjKUPIrIiL3h7597aS3S5frKjiwYgU8+SR/ffYZZXv35sKFC4ABOgH7gBx4GSev9rxIfMPhvL5xDABV81Zh5paCFPpkFgwcCCVLpvBJiUhKU/IrIiKp24ULkDUrlC5tP26kYkUoXpwVXbtyITISqAiMByoBUPPBQ7wS04i34/fzx0YLL8vwZtHO/F+78XhPeRz69YOhQ1PqjETEjZT8iohI6lajhj3iO2VK0n2yZoW5c6lQqjYwFngRgFy54njvg2h25ZhMu3WHiHNZFAn3ZebK9FTpWQl8/GD1akiXLkVORUTcT8mviIikXmfOwM6d0LHjTbsdO3aebq+cZUX8biATEEv3lyIIemEp/7f6f5zcfxKDoUv5LoxpMIaMozL+s7MSXxGPouRXRERSr7/+sp+LFUvUHB4ezsiRIzl69ASbN5dh585goBoAGVnCB2P/4osMM5m09FcAKuWtxPhG46mYt2JKRi8iqZCSXxERSb18fOznuDgADh48yMSJExk/fgKxsU2B4UAJAAID9/JSn72cK7uIrts+w7pokTN9Tt6t+y4dy3bEYVTdU0SU/IqISGpWuDA4HGxZsYL358xhzpy5WFZ94BegMgC5c19hwKAwYsosYsTa4YRtC8Pb4U3vSr15s+abBKYLdOspiEjqouRXRERSJcuy+GHdOt4vUICQqVOBZhizHggCIFu2ON5+x+BbcQ5D17zF8ZDjANQPz8nYp2fycLn67gteRFIt/Q1IRERSnZCQECpWrEj9+g0IOVweh+MPYCGWFUTOnPD++xZjVy1igqs0XZd05vjl45TJVYZlJ+uyYuIVHi5e1d2nICKplEZ+RUQk+UVH2yXF9u8HX18oWxYqVABH4jGYqKgo3njjDcaOnQi0wctrBk7nw7hckDvdBfp3uUDpXmep83VVWGjvUyhzIYbVGcZTZ3Li6NEA+vSB9OlT/hxF5L6g5FdERG7u8mWIiIAcOewlhW/X4sXwwgtw7lzi9tWr4fHHr77dv38/LVp0ZteuasBhIC9OJzz4IPTvcZmiS4IYeeUIr35tXd2nxUMtmN16Fr6ffQGvdYGHH4YhQ+7kLEXEQ2jag4iIJO3FFyEwEPLkgUyZoHVr2LLl9o4RFAR168LSpXDyJBw5AvPmQfXq9vbt2xk3dgWlSq1m167lwLtAXkqWhClTLD5evpSvcjSkYb3D/FjQIlM0/N+vXpz5oQLzg+fi6+0HEybYI8k//AAZM94kGBHxdBr5FRGRxGJj7RFehwNatLArLmTNai828c03sGOH/fi7DFlSdu+26/M+8ADMmpV4W4ECOJ2wePwBXn31OEdodHVTvXrxvN7XQfiD8xmxZjhb52wFIJt/Nl6t8io9vR4jc8wS2LcPvLzsnVauhJw5r5tGISLyb8ayrFv3ukeCgoKsTZs2pdjPExGRO9CtG8THw9SpYEzibZcu2dMXihS5+TFcLihaFOrUsY9zjaNHYdo0+3H8+N+tkVQN2s1HnxVnU8xMxm0Yx55zewB4IMMD9H2sL92CupHBN8O9OUcRSfOMMZstywr6d7tGfkVE5B+bNsGUKfC//12f+II9BSIwEKKioF8/6NIFSpe+vt/GjXD4MNSuDdiDyYsWwaef2oO0/4y77AemUjLrZCo1KkH9pXu5FHMJgPyB+elfrT+dynUinbeWIBaRe0PJr4iI/GPSJHvO7ODBN+8XFQVffQWhofb83X/bsgULWJe+Ht/0htmz4cwZe5Ovr4scOdYQeuJNKLwaKsPuYrDLbIAYqPZgNXpV6kWrh1vh43WLqRUiIrdJya+IiCfo3BkiI+2bzDp0gMyZb9wvJASeeOLWN41lzQpdu8Lo0XD+PGTLBtgjutu3w9ezyzCDQ5xomevqLj4+e/H1nU5E3ERC84RhWhksezd8LS+e3uqk18frKfdgpXtxxiIiN6Q7A0RE0qJz52Do0H/mF1gW/P479Oxp38D27xvQAOLi7Am5LQDwIgAADVFJREFUpUpdt2nHjh1EREQkbmzeHJxOrNVr2LEDhg2DRx6xKFMG3v+lKicoBBwHxgBBxGV5iIiaIzCvX4ImYGWzyJcpHyPqjOCvlw8wbdp5yuWreK+vhIhIIhr5FRFJ7eLiYOFCu2pCtWp2IjtlCrRqZdfe/bcrV+wbzfbtg3btoEQJ++4ysMuU9eoF7dvbx2nfPvG+770HNWsmalqyZAlNmzalePHiLFu2DB8fH8LCIvh1XXp+4CNWP1eVM+F/9zbAOWAuWTMs4sm22QmsmZWfrzjZllDm18KiRoEa9KrUixYPtcDboV9FIpJyVO1BRCQ1++UX6NQJDh2ypxlMngwHDtglxAIC4O234fXXE9+c1r27fWfZ8uVQv/71x4yNhY8+gh49IEPi6glxcXGMHDkSh8NB8+bNOX36NM2aNSMyMhLIAdQDmgGNgGunTpwBFgNzyZ59G0M+Gcgu/13M3D6TyzGXAQj0C+TZ0s/SpUIXSuf6101yfftC5crQps3dXjERESDpag9KfkVEUqvvvrNHbgsXhjFjoHHjf+ra7twJAwfC99/b83k//dROgI8etfv37Aljx976Z/z9O8AYIiIiCG7ViuUrVyZszAjUBOqSKVNLLl8ukGjXdOkOki3bb+TLt4VGjbJRonxR5u+Zz4EMB9hy5p+FMKr4FaVbw4G0LdWWAJ+A62P46Sd7pPrDD+HVV2/7MomI3IhKnYmI3E927rRvTKtcGZYts1dXu1apUjB/vl2VYfhwePRReOUVexEKl8seDb6FyK1bCXj2WfjsMy4WL069eu3YssUbf4aRPnMrzoUV4+9fE5cvg7+/xUMPnSM4OD3t2gVQJG9eLlhNmL87ltk7Z/PO1iE4LSdEJozyPtqBLhPXUzpkOxQOgLI3SHz374dnnrHrBnfteg8unIjIzSn5FRFJjfLlgwED7KkJ/058/2YMvPMOhIVBoUIAxFepgtfo0Zj8+ZM8tGVZvPRSbyZP/o2c1CRTo7OERuclKsoe8Y0CosLAy8uibNloGjXyo149Q5UqBj+/HIRFh7Fg1yx6jenDqpzhxBsXAN4ObxoXbUy7Uu0ILhlsj/JWvQhNmtgj2F99BS+8AE2b2ivIDR5sT79Il86eohFwg+RYROQe07QHEZH70Lhx41i2bBmNGjUid+7c5MiRg++++45JkyZRqFAhWrduTc2aNfn99985fdrCsh4hMrIIly4V4Pffwzl3Lhfg/6+jRlKxooOaOQ9Ra8n/qPFMATJ+Po54B+w5t4f1x9ezYO8CVhxYQZwrDgAvHNQtUo+2JdvS8uGWZPXPen2wsbEwahR88AFER9s35Dkcdkk1f397ekaBAtfvJyJyF5Jlzq8xpiEwFvACplqW9d7N+iv5FRFP9ttvv7Fjxw66dOmCudHqaX9buBDOnyf+2Wfx8vK6ru/KlSt54okn/rWTP1A44VEk4VECKA3k4kby5LlCmbi1+MauIaKik1dfrUvDRnXZe34vmyYMZPP6BWwql5Ot6a8QFR91dT+HZah12KJd5mq0HD6PHBly/rcLEBsLu3ZB2bL2+7g48NEiFiKSPO558muM8QL2AfWxCzluBNpblrUrqX2U/IqIp5ozZw5t27YF4PDhwxQsWDDpzq1asXjDBjrFxuLnl47g4M7Urt2BDBmKMmfOGiZNWgA8AOTmgQcqExaWjejoG4y4JvDziyFLlhOcOrUC2EGNGlkY9GZzSlbMQ2jnNuyPOcmmHi3YfHIzW05uISIu4rpjFPbPQ4Wtp6l10EnrAz7kemWgPW3BoXLxIpI6JUfy+xgwxLKsJxLevwFgWda7Se2j5FdEPFGZHr35c/sOcDnA8iJH9txkDsxO7gcexGF8cMY7iIn2IjrSi+goL84fOMX5WD9wZQQrI//cnmGBsW747OXlJKvPQSj9A8X3PkBg49L4ZbmMX2AYLp9LhMWEcfbKWS5GX+TI5SM3jbdAYAEq5KlAUO4ggvIEUT53ebIdPQNTp9pzi4OD7ZrDIiKpWHJUe8gL/HXN++NA5bs4nohImrN3717+zD4J6sVdbTub8Nh/bcfAa16XvP2f40w4JsDZ/EDcHLv07pmb71fugXIUzFyQ8rnLE5QniAq5K5Aj/Q0Wzng4m11uTUTkPpfs1R6MMV2BrgD5b3L3sYhIWlSiRAkKnuiA04ojT95cREaFExF1mdOnQ4mIvAzGRZasmXggT1a8fZx4+7rYtm09liOWho3q4pvOgcMLXC4Xp0+dxmDImzcvxhgM5p/ndetYZx3nSCYnzY/5UzI6E5m7vUJgQBYyp8tMoF8gmdNlJltANrL5ZyNzusx4bd1mz7/9u3awiIgH0LQHERE3iI+P55NPPmHQoEFcvnz5uu1VjeHXI0fgvwwauFx2abTHH4fZs+Hbb+Hpp+G336BSpRvvs2mTvVRynz72ksYiImlMcsz59ca+4a0uEIp9w9vTlmXtTGofJb8iIomdOnWK0aNHc/z4caKjo4mJicGKi6Nfz57UadHivx/o6FGIiICSCXMmdu60F8KwLHuFisCEeRWWBYsWQceOdv3gjRshxw2mOYiI3OeSq9RZY+Aj7FJn0yzLGn6z/kp+RURug8sFhw5B0aJJ9wkLs5PYpKoufP89tG8P1avbNXbPnIG9e+0V4RYssJdCFhFJg5JKfu+qRo1lWUstyypuWVaRWyW+IiJym/r0gfLlYcWKG28/eBCqVoWnnkr6GMWLQ6dOcP487NtnV2uYMsWe9qDEV0Q8kFZ4ExFJrUJDoWFD2LEDWrWylwguWBAuXIClS+3SY35+9ghuzZrujlZEJFVJjlJnIiKSnPLmhfXr7RvSJkyAefMgfXp7bq+Pjz3iO3w4PPiguyMVEblvaORXROR+EBcHf/xhz/H19YVy5SBjRndHJSKSamnkV0TkfubjA0HXfYeLiMht0qLsIiIiIuIxlPyKiIiIiMdQ8isiIiIiHkPJr4iIiIh4DCW/IiIiIuIxlPyKiIiIiMdQ8isiIiIiHkPJr4iIiIh4DCW/IiIiIuIxlPyKiIiIiMdQ8isiIiIiHkPJr4iIiIh4DCW/IiIiIuIxlPyKiIiIiMdQ8isiIiIiHkPJr4iIiIh4DCW/IiIiIuIxlPyKiIiIiMdQ8isiIiIiHkPJr4iIiIh4DGNZVsr9MGPOAkdT7AfeG9mBc+4OwoPoeqcsXe+Uo2udsnS9U46udcrS9f7vCliWlePfjSma/N6PjDGbLMsKcnccnkLXO2XpeqccXeuUpeudcnStU5au993TtAcRERER8RhKfkVERETEYyj5vbUp7g7Aw+h6pyxd75Sja52ydL1Tjq51ytL1vkua8ysiIiIiHkMjvyIiIiLiMZT8JsEY08YYs9MY4zLGBF3TXtAYE2WM2ZbwmOTOONOKpK53wrY3jDEHjDF7jTFPuCvGtMgYM8QYE3rN57mxu2NKi4wxDRM+vweMMQPcHU9aZow5YozZnvB53uTueNIaY8w0Y8wZY8yOa9qyGmNWGWP2JzxncWeMaUUS11rf2feAkt+k7QBaAatvsO2gZVllEx7dUziutOqG19sYUxJ4CigFNAQmGmO8Uj68NO3Daz7PS90dTFqT8HmdADQCSgLtEz7XknxqJ3yeVQ7q3vsC+7v4WgOAHyzLKgb8kPBe7t4XXH+tQd/Zd03JbxIsy9ptWdZed8fhKW5yvZ8EZlmWFWNZ1mHgAFApZaMTuSuVgAOWZR2yLCsWmIX9uRa571iWtRq48K/mJ4EvE15/CbRI0aDSqCSutdwDSn7vTCFjzFZjzC/GmMfdHUwalxf465r3xxPa5N7paYz5M+FPbPpz5b2nz3DKsoCVxpjNxpiu7g7GQ+SyLOtkwutTQC53BuMB9J19lzw6+TXGhBhjdtzgcbNRmZNAfsuyygGvAV8bYzKlTMT3tzu83nKXbnHdPwGKAGWxP9tj3BqsyN2rbllWeexpJi8bY2q4OyBPYtklpFRGKvnoO/se8HZ3AO5kWVa9O9gnBohJeL3ZGHMQKA7oxopbuJPrDYQCD17zPl9Cm/xH//W6G2M+BRYnczieSJ/hFGRZVmjC8xljzHzsaSc3undD7p3TxpjclmWdNMbkBs64O6C0yrKs03+/1nf2nfPokd87YYzJ8fcNV8aYwkAx4JB7o0rTFgJPGWP8jDGFsK/3BjfHlGYk/KL6W0vsGw/l3toIFDPGFDLG+GLfwLnQzTGlScaY9MaYjH+/Bhqgz3RKWAg8n/D6eeB7N8aSpuk7+97w6JHfmzHGtAQ+BnIAS4wx2yzLegKoAQw1xsQBLqC7ZVmakH6XkrrelmXtNMZ8C+wC4oGXLctyujPWNGakMaYs9p8pjwDd3BtO2mNZVrwxpiewAvACplmWtdPNYaVVuYD5xhiwf799bVnWcveGlLYYY74BagHZjTHHgbeA94BvjTGdgaNAW/dFmHYkca1r6Tv77mmFNxERERHxGJr2ICIiIiIeQ8mviIiIiHgMJb8iIiIi4jGU/IqIiIiIx1DyKyIiIiIeQ8mviIiIiHgMJb8iIiIi4jGU/IqIiIiIx/h/M1GC8wQUAqAAAAAASUVORK5CYII="
},
"metadata": {
"tags": [],
"needs_background": "light"
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 553
},
"id": "d7NL-b2K-Pgu",
"outputId": "b0f133fa-aee8-4981-d488-74c64073a6a7"
}
},
{
"cell_type": "markdown",
"source": [
"## cholesky VS sqrtm\n"
],
"metadata": {
"id": "8t2dD3WzBRKr"
}
},
{
"cell_type": "markdown",
"source": [
"### 正定矩阵 & 半正定矩阵\n",
"对于任意长度的非零向量$X$,若:\n",
"- $X^TAX > 0$,则矩阵$A$为正定矩阵\n",
"- $X^TAX \\ge 0$,则矩阵$A$为半正定矩阵\n",
"\n",
"记转换后的向量$M = AX$,则$X^TAX = X^TM = \\cos(\\theta)\\cdot ||M||\\cdot||X||$,那么意味着:\n",
"- 对于正定矩阵,$\\cos(\\theta) > 0$,即$\\theta < 90^{\\circ}$\n",
"- 对于半正定举证,$\\cos(\\theta) \\ge 0$,即$\\theta \\le 90^{\\circ}$\n",
"\n",
"矩阵变换$AX$表示向量$X$会沿着该矩阵特征向量的方向进行变换(缩放),缩放比例由特征值$\\lambda$决定。从特征值角度来看:\n",
"- 正定矩阵所有特征值都大于0\n",
"- 半正定矩阵所有特征值大于等于0\n",
"\n",
"那么正定矩阵小于90度的含义是变换后的向量$M$是沿着原向量$X$的正方向进行缩放的(即$M$投影回原向量时方向不变)。\n",
"\n",
"参考资料:\n",
"- [矩阵的特征值是虚数,它的几何解释是什么?](https://www.zhihu.com/question/319276709/answer/1195620832)\n",
"- [如何理解正定矩阵和半正定矩阵](https://www.cnblogs.com/marsggbo/p/11461155.html)\n",
"- 为什么协方差矩阵要是半正定的?[浅谈「正定矩阵」和「半正定矩阵」](https://zhuanlan.zhihu.com/p/44860862)"
],
"metadata": {
"id": "WtEgINSRHxEn"
}
},
{
"cell_type": "markdown",
"source": [
"### cholesky分解\n",
"Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解$A = L L^T$。它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的。"
],
"metadata": {
"id": "jMaQ4-LQHrj9"
}
},
{
"cell_type": "code",
"execution_count": 17,
"source": [
"A = np.array([[1, 0, 0],\n",
" [1, 1, 1],\n",
" [2, 1, 1]])\n",
"\n",
"np.linalg.eig(A)[0]"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([2., 0., 1.])"
]
},
"metadata": {
"tags": []
},
"execution_count": 17
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "9wMJ4CH5BQoM",
"outputId": "1ae2d99d-c9fd-4931-90a5-4d190b666707"
}
},
{
"cell_type": "markdown",
"source": [
"由于$A$的特征值有0(半正定),所以不能用Cholesky分解,这里用`sqrtm`:"
],
"metadata": {
"id": "Ro-QTUDyIR9I"
}
},
{
"cell_type": "code",
"execution_count": 18,
"source": [
"scipy.linalg.sqrtm(A) @ scipy.linalg.sqrtm(A)"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[1., 0., 0.],\n",
" [1., 1., 1.],\n",
" [2., 1., 1.]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 18
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "LUdtZi2pIRN3",
"outputId": "060e32fb-bf29-4215-db87-a98494d632ed"
}
},
{
"cell_type": "markdown",
"source": [
"再看下对称正定矩阵:\n"
],
"metadata": {
"id": "DvtqHfFAIr5r"
}
},
{
"cell_type": "code",
"execution_count": 19,
"source": [
"B = np.array([[6, -3, 1],\n",
" [-3, 2, 0],\n",
" [1, 0, 4]])\n",
"\n",
"np.linalg.eig(B)[0]"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([7.81113862, 0.33192769, 3.85693369])"
]
},
"metadata": {
"tags": []
},
"execution_count": 19
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "JjWdFxbeI-fK",
"outputId": "ac06b60e-cb26-43f4-e692-aa1cdc9aef19"
}
},
{
"cell_type": "markdown",
"source": [
"其特征值都大于0,那么可以用`cholesky`分解:"
],
"metadata": {
"id": "hwiwdfpBJCCR"
}
},
{
"cell_type": "code",
"execution_count": 20,
"source": [
"scipy.linalg.cholesky(B).T @ scipy.linalg.cholesky(B)"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[ 6.00000000e+00, -3.00000000e+00, 1.00000000e+00],\n",
" [-3.00000000e+00, 2.00000000e+00, -4.26642159e-17],\n",
" [ 1.00000000e+00, -4.26642159e-17, 4.00000000e+00]])"
]
},
"metadata": {
"tags": []
},
"execution_count": 20
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "J_DEMf1BJBUS",
"outputId": "0083aa5a-79a0-4bd2-d2de-adde80dc7aba"
}
},
{
"cell_type": "markdown",
"source": [
"# filterpy\n",
"[Kalman-and-Bayesian-Filters-in-Python](https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python)中实现的[filterpy](https://filterpy.readthedocs.io/en/latest/kalman/UnscentedKalmanFilter.html)。\n",
"\n",
"```shell\n",
"!pip install filterpy\n",
"```"
],
"metadata": {
"id": "lYuJgE4NQnmb"
}
},
{
"cell_type": "code",
"execution_count": 21,
"source": [
"def motion_model(x, dt, u):\n",
" F = np.array([[1.0, 0, 0, 0],\n",
" [0, 1.0, 0, 0],\n",
" [0, 0, 1.0, 0],\n",
" [0, 0, 0, 0]])\n",
"\n",
" B = np.array([[dt * math.cos(x[2]), 0],\n",
" [dt * math.sin(x[2]), 0],\n",
" [0.0, dt],\n",
" [1.0, 0.0]])\n",
" x = F@x + B@u\n",
" return x\n",
"\n",
"def observation_model(x):\n",
" H = np.array([[1, 0, 0, 0],\n",
" [0, 1, 0, 0]])\n",
" z = H @ x\n",
" return z\n",
"\n",
"def observation(xTrue, dt, u):\n",
" # Simulation parameter\n",
" INPUT_NOISE = np.diag([1.0, np.deg2rad(30.0)]) ** 2\n",
" GPS_NOISE = np.diag([0.5, 0.5]) ** 2\n",
"\n",
" xTrue = motion_model(xTrue, dt, u)\n",
" # add noise to gps x-y\n",
" z = observation_model(xTrue) + GPS_NOISE @ np.random.randn(2)\n",
" # add noise to input\n",
" ud = u + INPUT_NOISE @ np.random.randn(2, 1)\n",
" return xTrue, z, ud"
],
"outputs": [],
"metadata": {
"id": "6uar_19gRFtT"
}
},
{
"cell_type": "code",
"execution_count": 22,
"source": [
"from filterpy.kalman import MerweScaledSigmaPoints\n",
"from filterpy.kalman import UnscentedKalmanFilter as UKF\n",
"\n",
"import numpy as np\n",
"\n",
"dt = 0.1\n",
"sigmas = MerweScaledSigmaPoints(4, alpha=.001, beta=2., kappa=0)\n",
"ukf = UKF(dim_x=4, dim_z=2, fx=motion_model,\n",
" hx=observation_model, dt=dt, points=sigmas)\n",
"ukf.x = np.zeros((4))\n",
"\n",
"# Covariance for UKF simulation\n",
"ukf.Q = np.diag([\n",
" 0.1, # variance of location on x-axis\n",
" 0.1, # variance of location on y-axis\n",
" np.deg2rad(1.0), # variance of yaw angle\n",
" 1.0 # variance of velocity\n",
"]) ** 2 # predict state covariance\n",
"ukf.R = np.diag([1.0, 1.0]) ** 2 # Observation x,y position covariance\n",
"\n",
"\n",
"xTrue = np.zeros((4))\n",
"u = np.array([1, 0.1])\n",
"uxs = []\n",
"xs = []\n",
"for i in range(500):\n",
" ukf.predict(u=u)\n",
" xTrue, z, _ = observation(xTrue, dt, u)\n",
" ukf.update(z)\n",
" uxs.append(ukf.x.copy())\n",
" xs.append(xTrue)\n",
"uxs = np.array(uxs).T\n",
"xs = np.array(xs).T\n",
"\n",
"plt.figure(figsize=(12, 9))\n",
"plt.plot(uxs[0], uxs[1], color='k', lw=2, label='Estimated')\n",
"plt.plot(xs[0], xs[1], color='g', lw=2, label='True')\n",
"plt.axis('equal')\n",
"plt.title(\"UKF Robot localization\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"print(f'UKF standard deviation {np.std(uxs - xs):.3f} meters')"
],
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
""
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAIYCAYAAAB0YMpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3wU1cLG8d/ZdEgINXQIXRCEUAQpNpQmCkpHURSkWVFRsIGKFQRFQVEEC12lCXhBQbp0gpSAdAi9phBS97x/bOSNGCTSJuX53s9+sjs7M/vs3nvDk9kzZ4y1FhERERER+TuX0wFERERERDIjFWURERERkXSoKIuIiIiIpENFWUREREQkHSrKIiIiIiLpUFEWEREREUmHirKIyHVijNlrjLnrGu27qzFm2bXYd5rXWGSM6Z56/0FjzPxr8BovG2PGXO39iohcDhVlEcnyjDHWGFP+gmWDjDHjU+/fboyJTPOcrzFmmjFmuTEmT+q6ScaY2DS3Fy/yWnuNMedS1zlijPnaGBN4bd9h+u/RSdbaCdbaJleyjwv/e0nd7zvW2u5Xlk5E5OpQURaRHMUY4wdMA/ICTay10alPTbHWBqa5ffAvu7nXWhsI1ADCgAHXNrWIiDhBRVlEcgxjTC7gJ8AbuMdae/ZK9metPQLMw1OY/3qN+4wxW4wxZ1KHKlS+YLM6xpitxpjTxphxxhj/NNs+bozZaYw5ZYyZZYwplrp8SeoqG1OPZHfIwHutb4xZY4yJSv1ZP81z+VNf+1Bqjhmpy/MZY2YbY46nLp9tjClxkf2fH+phjHnxgqPxScaYr1Ofe9QYE2GMiTHG7DbG9Exdnhv4GSiWZrtiab8JuNTnmXp0/wVjzB+p73NK2s9TRORKqSiLSE7hh6eYxQOtrLXnrnSHqSWyObAz9XFFYBLwLFAImAv8ZIzxTbPZg0BToBxQEXg1dds7gXeB9kBRYB8wGcBae2vqttVTj3ZPuUSu/MAcYARQABgGzDHGFEhd5TsgF3AjEAIMT13uAsYBpYFSwDng00t9DtbaD/46Eg9UBo4Df2U8BrQE8gCPAsONMTVT/0hpDhxKcxT/0AXvIyOfZ3ugGVAGuAnoeqm8IiIZpaIsIjlFEHAL8I21NiGd59unHrX861bsX/Y1wxgTAxzAUwQHpi7vAMyx1v5irU0ChgIBQP00235qrT1grT0FvA10Sl3+IDDWWrs+Nd8A4BZjTOhlvNd7gB3W2u+stcnW2knANuBeY0xRPAW1l7X2tLU2yVq7GMBae9Ja+6O1Ns5aG5Oa77aMvqgxJgCYAXxsrf05dZ9zrLW7rMdiYD7QKIO7zMjnOcJaeyj18/yJNEf3RUSulIqyiGQHKYDPBct8gKQ0j08AHYFvjDFN09nHVGtt3jS3Q+ms85fW1tog4HbgBqBg6vJieI4EA2CtdeMp08XTbHsgzf19qdukt20scPKCbTPqb/tK81rFgZLAKWvt6Qs3MsbkMsaMNsbsM8ZEA0uAvMYYrwy+7lfAdmvt+2n22dwYszJ1OMkZoAX//3n9p/dxkc/zSJr7ccA1P7FSRHIOFWURyQ72A6EXLCvDBWXRWjsNeBz4wRhzx5W+aOoR0q/xHOkEOIRn2AIAxhiDp5geTLNZyTT3S6Vuk962ufEMm0i7bUb9bV9pXusgnqKZ3xiTN53tngcqAXWttXmAv4Z8mEu9oDGmP56hJN3SLPMDfsTz+RS21ubFM3zir/3Z//I+LvJ5iohcMyrKIpIdTAFeNcaUMMa4jGeu4nuBHy5cMXUYwpPATGNMg6vw2h8BdxtjqgNTgXuMMY2NMT54imcCsCLN+k+k5swPvML/j+WdBDxqjKmRWjDfAVZZa/emPn8UKJvBTHOBisaYzsYY79ST/6oAs621h/GM1R6VevKejzHmr0IchGdc8pnUfAPT3fsFjDHNgaeB+y8Y++2LZ2z4cSA5db20U8odBQoYY4IvsuuMfJ4iIteMirKIZAdv4ilPy4DTwAfAg9bazemtbK39Bk/pmmOMuflKXthaexz4FnjdWrsdeAj4BM9Qj3vxTCWXmGaTiXjG6e4GdgGDU/fzK/AaniOwh/Gc7NcxzXaD8AwbOWOMaX+JTCfxnED3PJ7hGy8CLa21J1JX6YJnWMo2PGOsn01d/hGeMcAngJXA/zL4MXTAc7JdRJoZLD5PHef8NJ7CexroDMxKk3Mbnj8Qdqc3LjyDn6eIyDVjrL3UN18iIiIiIjmPjiiLiIiIiKRDRVlEREREJB0qyiIiIiIi6VBRFhERERFJh4qyiIiIiEg6vJ0OkJ6CBQva0NBQp2OIiIiISDa2bt26E9baQhd7PlMW5dDQUNauXet0DBERERHJxowx+/7teQ29EBERERFJh4qyiIiIiEg6VJRFRERERNKRKccoi4iIiORUSUlJREZGEh8f73SUbMPf358SJUrg4+Pzn7ZTURYRERHJRCIjIwkKCiI0NBRjjNNxsjxrLSdPniQyMpIyZcr8p2019EJEREQkE4mPj6dAgQIqyVeJMYYCBQpc1hF6FWURERGRTEYl+eq63M9TRVlERERE/sbLy4saNWqcv7333nsXXXfGjBls3br1/OPXX3+dX3/99YoznDlzhlGjRv3n7QYNGsTQoUOv+PVBY5RFRERE5AIBAQGEh4dnaN0ZM2bQsmVLqlSpAsCbb755VTL8VZT79OlzVfZ3OXREWUREREQypH///lSpUoWbbrqJF154gRUrVjBr1iz69etHjRo12LVrF127duWHH34APFdbHjBgADVq1KB27dqsX7+epk2bUq5cOT7//HMAYmNjady4MTVr1qRatWrMnDnz/Gvt2rWLGjVq0K9fPwCGDBlCnTp1uOmmmxg4cOD5XG+//TYVK1akYcOGbN++/aq9Xx1RFhEREcmkrtVYZWvtvz5/7tw5atSocf7xgAEDuOuuu5g+fTrbtm3DGMOZM2fImzcv9913Hy1btqRt27bp7qtUqVKEh4fTt29funbtyvLly4mPj6dq1ar06tULf39/pk+fTp48eThx4gT16tXjvvvu47333mPz5s3nj2zPnz+fHTt2sHr1aqy13HfffSxZsoTcuXMzefJkwsPDSU5OpmbNmtSqVeuqfE4qyiIiIiLyN+kNvUhOTsbf359u3brRsmVLWrZsmaF93XfffQBUq1aN2NhYgoKCCAoKws/PjzNnzpA7d25efvlllixZgsvl4uDBgxw9evQf+5k/fz7z588nLCwM8ByJ3rFjBzExMdx///3kypXrb693NWjohYiIiEgmZa29JrfL4e3tzerVq2nbti2zZ8+mWbNmGdrOz88PAJfLdf7+X4+Tk5OZMGECx48fZ926dYSHh1O4cOF0p3Kz1jJgwADCw8MJDw9n586ddOvW7bLeS0apKIuIiIjIJcXGxhIVFUWLFi0YPnw4GzduBCAoKIiYmJjL3m9UVBQhISH4+Pjw22+/sW/fvnT327RpU8aOHUtsbCwABw8e5NixY9x6663MmDGDc+fOERMTw08//XQF7/LvNPRCRERERP7mwjHKzZo145lnnqFVq1bEx8djrWXYsGEAdOzYkccff5wRI0acP4nvv3jwwQe59957qVatGrVr1+aGG24AoECBAjRo0ICqVavSvHlzhgwZQkREBLfccgsAgYGBjB8/npo1a9KhQweqV69OSEgIderUuQqfgIe53MPv11Lt2rXt2rVrnY4hIiIict1FRERQuXJlp2NkO+l9rsaYddba2hfbRkMvRERERETSoaIsIiIiIpKOSxZlY0xJY8xvxpitxpgtxphnUpfnN8b8YozZkfoz30W2fyR1nR3GmEeu9hsQEcnOjhw5wubNm4mMjCQlJcXpOCIiOUpGTuZLBp631q43xgQB64wxvwBdgQXW2veMMf2B/sBLaTc0xuQHBgK1AZu67Sxr7emr+SZERLKSuLg4IiIi2L59O7t37yY4OJjSpUtTqFAhgoODiYqKYunSpcycOZMVK1ac365MmTKMHDmSpk2bsnTpUnbv3k358uWpX78+Xl5eDr4jEZHs6ZJF2Vp7GDicej/GGBMBFAdaAbenrvYNsIgLijLQFPjFWnsKILVgNwMmXYXsIiJZSkxMDMOGDWPo0KHnpze6FD9/P0LLhXLs9DH2HN9Di04tKFikICfOnAAvwAXlK5Snb9++hNUI+9tVvHy9fPH39sff258A74Dz9/28/XAZjbwTEbmU/zQ9nDEmFAgDVgGFU0s0wBGgcDqbFAcOpHkcmbpMRCRbO3z4MF9//TUREREkJiaSmJjIsuXLOB57HAKhdO3SFC5bmMAigZxOPM3xuOOctWdJIAG3rxvfIF+snyXexrPdvf1v+z7Bib893slOnlj/BKzPWDaDIcgviLz+ec/fgv2CyeuflwIBBSgSWIQigUUoHFj4/P2CuQri7dKMoiKSs2T4t54xJhD4EXjWWhud9qiFtdYaY65onjljTA+gB3iuCS4iklWcPHmSo0ePUqJECeavmM/IySNZsmkJ7iA3BOO5FQIeB3w82+xL/Q8AvkDg3/cZTzykDkn29fIlwDsAP28/fIwPLreL4NzBniPDuDh8+DAHDx78x9W2vP29CcwbiPE2uL3cnEs6R6I7EettiU6IJjohmv1R+zP0Hl3GRUG/goTmC6V8wfKEBodSOm9pSgeXpmy+spTJV0ZFWiSbOHnyJI0bNwY850l4eXlRqFAhAFavXo2vr6+T8a6rDP1WM8b44CnJE6y101IXHzXGFLXWHjbGFAWOpbPpQf5/eAZACTxDNP7BWvsF8AV45lHOUHoRkevg4MGDrFmzhsjISA4dOsSxk8c45j7G/vj97IraRaxvLBTAc/MDyqTe0pHXP6/naG3uwoTkDqFw7sIUDvTcLxBQwHN01z+YYL/g8z/9vP3S39kFGdeuXUuePHk4evQo7733Hhs3buQMZ/65sknN6Q+1G9Wmx1M9KFSyEGfiz3Ai7gSb9mzi5yU/c/zccUyQwZXHRYpfCsfij3Hs8DFWH179z32mQGBiICVzleSuGncRVjKMSgUrUTWkKnn88mT4sxYR5xUoUIDw8HAABg0aRGBgIC+88ML555OTk/H2zhl/GF/yXRrPoeOvgAhr7bA0T80CHgHeS/05M53N5wHvpJkRowkw4IoSi4hcQ0eOHGHSpElMmDCBw4cP4wpyEZkc6RlcVgTPz8J4xgenw8QbivsVp3b52lQuWplSwaUoHVyaUsGlKBVciiC/oGuSu3jx4hQv/v8j2zp06MBvv/3Gli1b2L9/P/v27aNQoUL06NGDfPnyMXHiRD744APW/rSWdbPX8eCDD+Lv78+qVavYvHkz1lp8fX1JTEwkhRTwgtCbQjkUd4jEgETPUfK8qbf8QDDEBsQSYSOI2BABG/4/W5m8ZbixwI2U9C1JSd+SuA+5ObLtCGdjz9KiRQuaN29O7ty5r8nnIiJXR9euXfH392fDhg00aNCAPHny/K1AV61aldmzZxMaGsr48eMZMWIEiYmJ1K1bl1GjRmXZE44z8udAA6ALsMkYE5667GU8BXmqMaYbsA9oD2CMqQ30stZ2t9aeMsa8BaxJ3e7Nv07sExFx2sGDB1m6dCnLli1j1apVHDpxiMPmMLaohXLArXgK4QUMhhDvEEJzhRJWMow65epQuWBlyucrT8HcBf92Qp1TjDHceeed3Hnnnek+379/f3r06MHgwYP59NNPGT9+/PnnfH196dmzJ2+++SbGGHbu3EmRIkUoXrw4SUlJbN68mfj4eAICAggICMDf359Tsaf43+r/8e2cb9l2YhsUAP+S/iTnS2bPmT3sObPn7wFyA2dg3Gfj8H3DlzrF6lCjQg3KlSuHl5cXLpcLLy8vqlWrxs0335yjvuoVScu8cW1+n9iB//3L+8jISFasWIGXlxeDBg1Kd52IiAimTJnC8uXL8fHxoU+fPkyYMIGHH374ChM7IyOzXizD80Vdehqns/5aoHuax2OBsZcbUETkajt06BAdO3ZkafhSKI3nVgfPOOILJoPI5Z2LmkVrUqNIDW4qfBM3Fb6JG0NuJNA38J87zmLy58/PsGHD6N27N2PHjqVQoULUq1ePsLAwAgICzq9Xq1at8/d9fHwICwv7x75KU5qwG8N48eEXGTVqFB988AGR0yI9n2cBoAjkuyEfyQWTSS6QzDn/c54/RspBIoksZznLTy+HFXgOvewDjgMWcuXKxbPPPsugQYPw8fG5th+KiFxUu3btLnlkeMGCBaxbt446deoAcO7cOUJCQq5HvGsiZwwwEREBVmxdwSc/fcLMjTM5F3buH3/qextvqoZUpV6JetQtUZebi99MpQKV8HJlza8MM6pChQq8++67V2VfXl5ePPXUU/Ts2ZMpU6YwY8YMatasSZcuXf52ovbB6IOsOriKVZGrWLpnKRuObiA+XzzkA6p61vFL8cP3sC8xm2J4Z+w7/G/e/+jZoyeJiYkUL16cpk2bkitXrquSWySzupwjv9dK2iFS3t7euN3u84/j4+MBsNbyyCOPXLXfKU5TURaRbCs2MZbf9vzG18u+Zu72ucTn9vwip5LnR6BPII1KN+LW0rfSqFQjahatSYBPwMV3KBnm6+tLly5d6NKlS7rPF89TnAfyPMADlR8AIMWdwtbjW/k98ncW71vM4r2LORhzkIQSCZ7TwIH1cevpOb8n7AR2QW6bm1tvvZXg4GDuuOMOHn30UR1xFrlOQkNDmT17NgDr169nzx7P8KrGjRvTqlUr+vbtS0hICKdOnSImJobSpUs7GfeyqSiLSJa3YMECfvjhB/wD/EnKm8T+gP3sce1he9x2ktxJnpVyA/FQNKEod1W4iyfueYJaJWppSrNMwsvlRbXC1ahWuBo9avXAWsvu07tZtHcRi/ct5rc9vxFJJFTDcwPOHjnLzzt/hpUweepkhgwZwqBBg+jQoUOOOSNfxClt2rTh22+/5cYbb6Ru3bpUrFgRgCpVqjB48GCaNGmC2+3Gx8eHkSNHZtmibC6cdzMzqF27tl27dq3TMUQkk4uKiuK5F55j7K9jPUeJK+GZgeEvbjyTVO6CpuWbMmXYFIKD0jk7TzI9ay1/nvyTebvmMW/XPH7b8xvnks+df96V4MK9zQ3boFRSKR576DFuu+026tat+7fx1iJZQUREBJUrV3Y6RraT3udqjFlnra19sW1UlEUkS3G73cQnx9P/y/58ufxL4kvGQ5phqrlsLkqeK4n/AX9Orz9N5I5IWrduzeTJk/W1fDYSnxzPsv3LmLdzHrN3zPbMtPGXZGA3sA389/rzdPenGTBgAHnz5nUqrsh/oqJ8bagoi0i2ZK1l7ry59P2kLzt8d0BlPBfMSBUaGEq7m9rRqlIr6pWo97eT71JSUrLs/J2ScX+e/JOZ22Yyfdt0VkauxJL6b5sbz3jmvbl5pc0r9O3dF39/f0ezilyKivK1oaIsItmCtZZjx46xNWIrP236iclbJnM472FI8w26zwkf7q90P290eIMbCt3gXFjJdI7GHuWnP39iWsQ05u+aT4pNvRZ4Mvgd8OOBCg/QsVZHqlWqRqlSpfSHlGQ6KsrXhoqyiGRZsbGxTJo0ifnz57Ng/QJOlz4N1fFc+S1VYVOYx295nFblWxFWKkwFRy7pZNxJfoz4kZGLR/JH9B//f1WABGALeG3yoqxPWSqUr0CFChWoWrUqrVu3pmDBgk7GlhwuIiKCG264IVNcvCi7sNaybds2FWURyVqOHDnCxx9/zGdjPyOqeBTUwHMBkFS5k3JzS+5beKPdG9QvX9+xnJL1HYo+xMCpA5m5eybH/Y7//xMngXBgIxDtmR+2UaNGlCtXjuLFi+NyucidOzePPfYY+fLlcyi95CR79uwhKCiIAgUKqCxfBdZaTp48SUxMDGXKlPnbcyrKIpIpRUdHM2TIEIaOH0p8tXjPlF+pVykO8Aqg3Y3t6FqjK7eF3obLuP51XyL/1bYT2/g6/Gu+Cf+GI2ePAJ5Lkxc4XYBT807h3u6GC/55rFy5MnPnziU0NPT6B5YcJSkpicjIyPMX8ZAr5+/vT4kSJf5xUreKsohkKomJiYz4bARv/PAGsTfEnr+YBECjUo14tMajtK3SliC/IOdCSo6R7E5m/q75jAsfx6zts0hMSQQgxDeEW3xvoeyZsgSaQKZNm8aWLVsoXLgwkydP5vbbb3c2uIhcFSrKIpJpjJ02lucmP0dU2ajzJ+YFegfSrVY3etbqSeVCOnlFnHMi7gRfh3/N52s/Z9fpXQD4uHxoW6Utj1R5hCFPD2HBrwtwuVz069eP119/XZfQFsniVJRFxBGJiYlMnjyZdevWERETwXr/9ZwsdBJSR1GUz1We/o3706laJ3L5qGxI5uG2bn7Z9Quj1o5i9p+zcVs3AHWK1aH4/uLMeGcGuKFMmTK8/PLLdOrUidy5czucWkQuh4qyiFxX1lomT55Mv5f6cTDoINwClEx9MgVq+9fmk4c+oV6pek7GFMmQ/VH7Gb12NKPXjebkuZMAFPEvglltODznMCRAcHAwn332GZ06dXI4rYj8VyrKInJdWGtZtmwZb7z9BgtOLYD6QOoEAblduWlZpCX97+xPjXI1HM0pcjnikuL4Jvwbhq8czo5TOwAIMAHk3Z2Xwz8ehlgYMmQIzz//vGYpEMlCVJRF5Jr7/fffeaLvE2zw2uA5gpx6Hl6F/BXoW68vD1d/mNy++mpasj63dTP7z9l8+PuHLNm3BABvvElelQwr4OlHnmbYsGGa41ski1BRFpFrIjExkenTp/Pd998x5/gcqAekDjWuWrAqA+8YyP033P+3y0mLZCdrDq7h3WXvMn3bdM+CFCAcWhdszY9jfsTl0rSGIpmdirKIXDUpKSls2bKF6dOnM+qrURwre8xzBNnf83y9YvV47fbXaF6+ub5+lhxj87HNvLP0HaZsnoIbN7ihYkJFZvSdQeXCmslFJDNTURaRK2at5bXXXuOjjz7ibOJZqAM05PwR5AZFGzC4yWBuK32bCrLkWH+e/JOnpjzF/KPzPbO7uOGB0AcY0WYExfMUdzqeiKRDRVlELkt0dDR9+vQhPDycggULsnjZYqgFXrd7kZIrBYCGJRvyduO3ubX0rQ6nFck8psyfwuPfPU5M2Rhwga/x5al6T9G/YX8K5irodDwRSUNFWUT+s0OHDtGiRQs2btzoWVAVaMz5WSxqFa3F4DsH07RcUx1BFklHXFwcj774KFOPTYUbPcuCfIN4/pbneaH+Czq5VSSTuFRR1pkGInKe2+3myy+/pFq1amzcuJGS9Utyw4c3QFsgH9xY6EamtZ/GmsfX0Kx8M5VkkYvIlSsXkz+ZzBMhT8BocO12EZMYw6DFg6j4aUW+Dv/6/IVMRCTz0hFlEQHgwIEDPPzwwyxatAjyQ0jnEI4VPAZAkcAivHXHWzxa41HNYiHyH7jdbvr06cPo0aOhNAS0DuBcvnMAhBUJ48MmH3JHmTscTimSc2nohYhc0sqVK2nevDlnzp0hoHkAiTUSSSGFAO8A+tXvR78G/Qj0DXQ6pkiWNX/+fJ588kl27NxB4C2B+Lbw5VTyKQDuq3QfQ+8eSoUCFRxOKZLzqCiLyL9yu93Uql2L8ORw/O71I8E3AYOha42uvHXHWzpbX+QqiYqKolOnTvz888/gA96NvKEhJLuS8fXy5cX6LzKg0QBy+eRyOqpIjqExyiJyUQkJCTz11lOEVw+HNpDgm0C9EvVY22MtY1uNVUkWuYqCg4P56aefGD58OHVq1CF5YTLJw5JxbXKRmJLI4KWDqTKyCjO3zSQzHsQSyYlUlEVyqEW/L6LQw4UY5R4FoRDkCmJcq3Esf2w5NYvWdDqeSLbk5eXFs88+y+rVq9m+fTvtmrfD/aMbvgKv417si9pH6ymtaTmpJbtO7XI6rkiOp6IskgNN3TiVu6bfRUyVGADuKXgP+/vtp2uNrriMfi2IXA8VK1Zk6tSprFq1itvL3U7KZykwF4iHuTvmUmVkFd5Z+g5JKUlORxXJsfQvokgOcjjmMHeNvosOMzqQkjuFgJMBrOm+htlPzCavf16n44nkSDfffDMLFy7k5zk/0zioMa5RLgiHRHciryx8hTpf1mHtIZ23I+IEFWWRbO7s2bP0H9Cfmo/XpOT7JVlwZAEkgs8vPvzW5Tdql7zoOQwicp0YY2jWrBm//vorR3cd5eE8D8M3wGnYeHQjN395M+3HtOd07Gmno4rkKCrKItmY2+2mTbc2vH/4fTaU2ECKTwreu73pkdSDnRN3UrdOXacjisgFChYsyNdff02/tv1gFLACrNvy/cHvKfhaQdr2a0tiYqLTMUVyBE0PJ5JNJSUl0fzV5izwWgB+npP1eof25pVWr5AnTx6n44lIBmzdupWVK1fy8x8/M9vMJj5vPAANvBsw/8X5mkpO5AppHmWRHGjllpU0G9WMqJAoABrlb8S0btMomKugw8lE5HIlpSTRbkQ7Zp6eCV5QzK8Y0x6aRt0S+mZI5HJpHmWRHMRay6eLPqXBhAZEhUThSnAxoOIAFj+5WCVZJIvz8fJhUq9JVFtVDY7BoYRD1BtTj6enP01iioZiiFwLKsoi2UR0QjRtJ7blqcVP4fZzE3w8mM19NvNOp3cwxjgdT0SugoCAANbPWc/QikPxWesDwCd/fEL598qz7cQ2h9OJZD8qyiLZwKrIVVQaXolpO6dBIpTYUILdb+6mconKTkcTkavM29ub5595nr2j93LHvjvgNBxIPsBNn97E2A1jdVU/katIRVkkC3NbN09MeoJbxtzCkYQjcBhqrKrBui/WkT9/fqfjicg1VKxYMRaOW8hn1T/DbDIkmSS6zepGh6kdiE6IdjqeSLagoiySRc3/fT4Fni3AqD9HYY3FtcrFW6XfYu28tYSEhDgdT0Suk16P9mLaQ9Pwnu0NifD9tu8p9kYx3v/ufVJSUpyOJ5KlqSiLZDFJSUl0eb0LTac15Uz+M5g4Q9uEtuwdvZdXB7yKl5eX0xFF5Dpr3bo1C4ctJGx1GByGs35n6b+jPyH3hfDRRx8RExPjdESRLEnTw4lkIX/++SeNX29MZKVIcEGJpBL80vsXbih+g9PRRCST2KVmmtYAACAASURBVHNgD13Hd2VJ4hLPgs2QZ1Eeej3Wi+eee47ChQs7G1AkE9H0cCLZQHJyMkM/HUqVN6sQWdlTkh8s9SB739yrkiwif1OmZBkWD1jM1DZTCXAFQFWI7hDNB2M/IDQ0lGeffZbVq1drWIZIBqgoi2RyS5Ys4Ybbb6Dfjn6kVEjBJ8WHCfdOYPyj4/FyaZiFiKSvXdV2bOi9gSqFqkAh8OrtRXzZeD7++GPq1q1LoUKFaNeuHUuXLnU6qkimpaEXIpnY1q1bqd65Osktk8EHSvuXZsHjCyiXv5zT0UQki4hNjOXxnx5n8ubJAFQ/W53oadHs2bUHAJfLxQcffMBzzz2nOdclx9HQC5EsKu5cHI0HNyb5fk9JfrDqg0Q8F6GSLCL/SaBvIBMfmMjHzT7G2+XNxtwbqfxmZTZs3cALL7yA2+3mhRdeoGPHjsTGxjodVyRTUVEWyYS27txK8b7FOVLpCLjh7Vvf5rsHviPAJ8DpaCKSBRljeLru0/zS5RfyB+Rn7o65dF7YmZ4DejJt2jSCgoKYOnUqNWvWZMqUKbjdbqcji2QKKsoimcyCdQuoPqI6Z4qewSQaRjUcxct3vKyvREXkit0eejtrHl9DlUJViDgRQd0xdclbIy+rV6+mcuXK7Nixg44dO3LvvffqCn8iZKAoG2PGGmOOGWM2p1k2xRgTnnrba4wJv8i2e40xm1LX06BjkUuYsnwKTb5vQnKBZALiAljZbSW97+7tdCwRyUbK5ivL791+554K93Dq3CmajG/C4tjFhIeHM3r0aPLkycPcuXNZs2aN01FFHJeRI8pfA83SLrDWdrDW1rDW1gB+BKb9y/Z3pK570YHSIgLjfh9Hp/91wh3gJt/pfOx4aQc3l73Z6Vgikg3l8cvDzI4z6Ve/H8nuZHrN6cXrS16n++PdefzxxwEYMmSIwylFnHfJomytXQKcSu854/kuuD0w6SrnEslRPlj8AY/NewzrbQmJDGHPW3sonr+407FEJBvzcnnxwd0f8NV9X+FlvHh/+ft0md6F3k/1JiAggB9++EFTx0mOd6VjlBsBR621Oy7yvAXmG2PWGWN6XOFriWQ7buvm6TlP89Kil8BA/vD8bHxrI8FBwU5HE5Ec4rGwx5jdebZndoxNE+m+qDtPv/g0AH379tWJfZKjXWlR7sS/H01uaK2tCTQHnjDG3HqxFY0xPYwxa40xa48fP36FsUQyv/jkeNp/355P1n4CKZBvUT7WDV9HkSJFnI4mIjlMs/LNWNJ1CUUCi7Bo7yJmFZxF4YqFWbduHZMnT3Y6nohjLrsoG2O8gQeAKRdbx1p7MPXnMWA6cNEBl9baL6y1ta21tQsVKnS5sUSyhJiEGO6ZeA8/RvwI8ZBrWi6WjVxGaGio09FEJIcKKxrGym4rqVywMhEnI4h/MB4Kwvz5852OJuKYKzmifBewzVobmd6Txpjcxpigv+4DTYDN6a0rkpOcjDtJ428bs3DPQogF1zcupn04jSpVqjgdTURyuNJ5S7P8seU0LNWQKBsFj8KqA6ucjiXimIxMDzcJ+B2oZIyJNMZ0S32qIxcMuzDGFDPGzE19WBhYZozZCKwG5lhr/3f1ootkPQejD3LruFtZc2gNnAa+ghGvjKBp06ZORxMRASBfQD7mPTSPu0rfBblhW91t3PnoncTFxTkdTeS6M5lxQvHatWvbtWs17bJkLztO7uCOcXdw8OxBOAZ8B892e5bhw4c7HU1E5B8SUxJp+GFD1pxbA0nQ0dWRSW9qkivJXowx6/5tCmNdmU/kOth6fCt1PqvjKcmRUGh2IX4Y94NKsohkWr5evvz+wu/cV+w+8IHJTGbcmnFOxxK5rlSURa6xzcc2U/+L+kSlRMFu6JzcmW0bttGmTRuno4mI/Csvlxczus+gyM4i4AXd53Zn4qaJTscSuW5UlEWuoT+O/kHDLxsSlRwFO6F/6f5MGDeB/PnzOx1NRCRDjDF0KdoFfgM3brpM78L4P8Y7HUvkulBRFrlGwo+Ec/vXt3tK8g7ok78P77zxjtOxRET+s06dOuG30s9Tlq2bh6c/zHcbv3M6lsg1p6Iscg1sOLyBO7+5k9Pxp+FPaHqmKZ8M/wTPVd9FRLKWsLAwli9fTui+UFgIFssjMx7h6/CvHU4mcm2pKItcZZuPbebu7+72lORt4D3Nm0+Gf4LLpf+7iUjWVatWLdatW0fHYh1hgacsPzrjUQ3DkGxN/3KLXEU7Tu7g7u/u5uS5kwQfDYbvoe9TfalQoYLT0URErlj+/PmZNGkSE3pNgAWAga4zujI9YrrT0USuCRVlkatk35l9NP62MUdij+B9wJuoL6MoGlKUV1991eloIiJXVefOnelZpScshhSbQocfOjBv5zynY4lcdSrKIlfB4ZjDNP62MQeiD2AOGJK/S+buO+5m2bJl5MmTx+l4IiJX3XvvvUep3aVgJSS5k7h/yv0s2bfE6VgiV5WKssgVOhl3kru+u4tdp3fhd8oPO97S67FezJs3j7JlyzodT0TkmsibNy8TJ0yE/4HXRi/OJZ+j5cSWrDm4xuloIleNirLIFYhLiuPeSfey9fhWAmICSBiTQOWylRk2bJhmuBCRbK9BgwbUqVOHlBkpNMrbiJjEGJpPaM72E9udjiZyVagoi1ymZHcyHX/oyO+Rv0MUnPviHEWDizJlyhQCAgKcjicicl20bt0aLBz57Aj1C9bn5LmTNJvQjMMxh52OJnLFVJRFLoO1lt6ze/PTnz9BHLgmuni227Ns27aNatWqOR1PROS66d69O5UqVWLH9h2seG4FhZMLs/fMXlpMbEF0QrTT8USuiIqyyGUYtGgQYzaMgSRgEkwaMYnhw4frxD0RyXFCQkIIDw9nwIABeKV4cXTYUbyjvAk/Ek6bqW1ITEl0OqLIZVNRFvmPvlj3BW8ueRPcwPcwasAo2rdv73QsERHH+Pv7884777BmzRrCKoWR/HUyxMKvu3/lsZmP4bZupyOKXBYVZZH/4Nfdv9J7dm/Pg9nw5kNv0rt3b2dDiYhkEmFhYaxatYpnH3kWJgCJMGHTBAb+NtDpaCKXRUVZJIMijkfwwOQHcOOGpfB0w6d1MRERkQv4+Pjw4Ycf0rJWS5gKWBi8dDAT/pjgdDSR/0xFWSQDTsSd4J6J9xCTFANboW3+tgwfPlxTwImIpMPlcjF69GgCDgbAz55lj816jBUHVjgbTOQ/UlEWuYSE5ATun3w/e87sgUNQdmNZxn41FpdL//cREbmYYsWK8dxzz8FqyLMtD4kpiTT/pjm7T+12OppIhulfepF/Ya2l5+yeLDuwDKIh18xcTJsyjaCgIKejiYhkei+++CKlS5cmemo07ILolGjqflxX08ZJlqGiLPIvPln9Cd9s/AYSwUw2fD/2e6pXr+50LBGRLCFPnjyEh4cza8Ys+pfrDyfghOsEt310m2bCkCxBRVnkIhbvXUzf//X1PJgJo14bRYsWLZwNJSKSxeTNm5d7772Xdwe+y5s3vAnnIDw+nLYft3U6msglqSiLpGPH0R3c8809nhkulsE7D75Dr169nI4lIpKlvfbEa3QN6goWpp+ZTr/R/ZyOJPKvVJRFLrBz706qDq7KWc7CLnipzkv079/f6VgiItnCuFfG0cS3CRgYumcoE37WtHGSeakoi1zg7o/uJrFgIj5nfZjRZQbvvfOepoETEbmKfh7wM+WSy0EAdP+lOzHxMU5HEkmXirJIGr0+78XefHshCWZ1nkWru1s5HUlEJNtxGRfLX1iOV5QX8cHxtPy8JdZap2OJ/IOKskiq7+Z9x+jI0QA8mPdBmtVo5nAiEZHsq3BwYV4p+wokwpKoJfSfrCFukvmoKIsAsQmx9JjfA3ygUnwlvn3uW6cjiYhke6/2fJWww2EAfLD5A8bNHudwIpG/U1GWHC8xMZGwV8KIzxOP12kvlg1YpqvuiYhcBz4+Pqz5cg0V4iqALzyz7BnikuKcjiVyntqA5GjWWm5+/GZ2Bu2EJPj01k8pmKeg07FERHIMLy8vfn7yZzgBMQEx9JqpqTgl81BRlhxt7IyxbCyxEYBXar5Crwf0C1pE5HorV7IcdffXhST4bst3TNo0yelIIoCKsuRQ1lo+/+Jz+vzaB3zhJq+beOv+t5yOJSKSY418fSSuXzy1pPuM7uyP2u9wIhEVZcmhBg4cSO+JvUkMScQnzodZPWZprmQREQfVqlWL99u9D9sgzh1Hq3GtcFu307Ekh1NRlhxn4sSJvDXuLbgNDIa5PeZSOqS007FERHK8559/nkfyPwKxEB4VTu9vejsdSXI4FWXJccaOHwsPAC54/pbnuavcXU5HEhERwBjD2E/Gclv0bQB8sesLmjzUhH379jmcTHIqFWXJcTYW2Qj5oXxQeQbfOdjpOCIikobL5eLXkb9Sx1UHvOGXPL9wQ9UbmDJlitPRJAdSUZYcZfKayZwofQKSYXDYYPy8/ZyOJCIiF/D29mbhSwsJzRMKhSH+5ng6duzI0KFDnY4mOYyKsuQYn4/9nM4TOwPgvdybBhUaOJxIREQuJtA3kPFtxmMwuG51QVHo168fL7/8MjExMU7HkxxCRVlyBGstfWf3xea1BMcFs/7T9ZQoUcLpWCIi8i8alGrAUzc/hRs3pZ4sBV7w7rvvEhoayvz5852OJzmAirLkCDM2zCC+ejykwKJnF1GtSjWnI4mISAa83fhtQvOGsj9pP4+NeYz69etz6tQp2rVrx/bt252OJ9mcirJke3FJcTwx/wkASu4tSY2iNRxOJCIiGRXoG8iYe8cAMD5yPF9M/4I2bdoQHR1Nnz59HE4n2Z2KsmR7by1+i8MJh+EoPFLuEafjiIjIf9S4bGO6h3UnMSWR7j91Z+SokQAsX76cpKQkh9NJdqaiLNnalmNbGLpiKFjgJ+jxWA+nI4mIyGUY2mQoxYKKsTJyJTP2z6BChQokJCSwadMmp6NJNqaiLNmWtZbec3qTbJNhLXS+tTMlS5Z0OpaIiFyGYP9gPmr6EQD9F/TnpltuAmDevHlOxpJsTkVZsq1vNn7D0v1LIRb8l/szbNgwpyOJiMgVaFulLU3KNeFM/BnO1DkDwJAhQzh16pTDySS7UlGWbOlk3ElemP+C58E8qFu9LoULF3Y2lIiIXBFjDCNbjMTPy48FJxdQ84GanD59msGDdZVVuTYuWZSNMWONMceMMZvTLBtkjDlojAlPvbW4yLbNjDHbjTE7jTH9r2ZwkX/T/9f+nDx3klIppWAT1K1b1+lIIiJyFZTPX54BDQcAcKbBGfCCESNGsGHDBoeTSXaUkSPKXwPN0lk+3FpbI/U298InjTFewEigOVAF6GSMqXIlYUUyYt2hdXy14St8XD4kz0oGoEWLdP+WExGRLOilhi9RPn95dsfspuHzDUlJSeHRRx/VDBhy1V2yKFtrlwCXM/jnZmCntXa3tTYRmAy0uoz9iGSYtZZn/vcMFkv7Uu05tPEQxYoVo2HDhk5HExGRq8Tf258RzUYA8EfePyhVuRQbN27k/fffdziZZDdXMkb5SWPMH6lDM/Kl83xx4ECax5Gpy0SumalbprL8wHLy+eZj4aCFALRv3x4vLy+Hk4mIyNXUvEJzmpdvTnRiNDc945kB48033yQqKsrhZJKdXG5R/gwoB9QADgMfXmkQY0wPY8xaY8za48ePX+nuJAeKS4qj3y/9AIieEc3hvYdp0KABr776qsPJRETkWviwyYd4GS/mHp2Lq6gLay1ut9vpWJKNXFZRttYetdamWGvdwJd4hllc6CCQdtLaEqnLLrbPL6y1ta21tQsVKnQ5sSSHG7piKAeiD+Bz0oeUNSk89thjLFy4kAIFCjgdTUREroHKhSrTp04f3NaN+243dza+k3z50vuSW+TyXFZRNsYUTfPwfmBzOqutASoYY8oYY3yBjsCsy3k9kUs5FHOI95d7xqYlzUqidq3afP755/j6+jqcTERErqWBtw0ktys3lIXY4rFOx5FsJiPTw00CfgcqGWMijTHdgA+MMZuMMX8AdwB9U9ctZoyZC2CtTQaeBOYBEcBUa+2Wa/Q+JId7Y9EbxCXFQQS4Drj47rvv8PHxcTqWiIhcYwVyFaBLyS4ArAlew7ETxxxOJNmJsdY6neEfateubdeuXet0DMkitp/Yzo2jbsRicX/iprhfcSIjI52OJSIi10l0bDSF3ihEYmAiFSIqsOW7LTpYIhlijFlnra19sed1ZT7J8l5Z+AopNoU7894JJ6Bs2bJORxIRkesoT2AePmzpmVdgR/EdzP3lH5d3ELksKsqSpa0+uJofI37E39ufXGtyAXDfffc5nEpERK63Prf2oUBSAcgDYzePdTqOZBMqypJlWWt56deXAHiy1pMsmL4AgDZt2jgZS0REHOAyLh7I8wAAc6LmcPj0YYcTSXagoixZ1oI9C1i0dxH5/PNR42wNzp49S61atShTpozT0URExAFvdHkD38O+pPim0ODFBuzfv9/pSJLFqShLlmStZdCiQQD0q9+PaROnAdCuXTsHU4mIiJOKFi3KuM7jANgTsoeKNSry2muvERuraePk8qgoS5a0YM8Clh9YTv6A/FQ7V41p06bh7+9Pp06dnI4mIiIO6nxrZ24rfhv4QkJYAoMHD6ZSpUps3pzeJR9E/p2KsmQ51lreWPwGAM/VfY5+z3guW/3qq69SqlQpJ6OJiEgm8F6z9wDIdUcuwuqHcejQIdq3b09cXJzDySSrUVGWLGfhnoUs27+M/AH5KbS3ENu2baNs2bL069fP6WgiIpIJ1CtRjyblmhCXEkfTgU2pXLkyERERPPPMM05HkyxGRVmylLRHk/vW7cuQwUMAeO2113S5ahEROe/1W18HYNT6UYz+ZjR+fn6MGTOG0aNHO5xMshIVZclSlu5fytL9S8nnn498O/Kxc+dOypcvz0MPPeR0NBERyUQalGpA4zKNiU6IZsHZBXzxxRcAPPXUU+zZs8fhdJJVqChLlvLB8g8AaFOiDS8965lDeeDAgXh7ezsZS0REMqHXb/McVf5k9Se06diGDh06kJSUxLhx4xxOJlmFirJkGZuPbWbOjjn4e/nz40s/cvbsWR5++GE6d+7sdDQREcmEGpVqRL0S9Th17hTjwsfRq1cvAD7++GMOHDjgcDrJClSUJcsYssIzHrlCbAVOHzxN06ZNGTNmDC6X/mcsIiL/ZIyhX33Pid4f/v4hDRo1oHXr1kRHR58vzSL/Rg1DsoQDUQeYuGkiLuNi/xTPlZbeffddfHx8HE4mIiKZWatKrSifvzx7z+xlWsQ0Pv/8c3x9fZk7dy5nzpxxOp5kcirKkiV8tPIjkt3J1PStSdTeKKpXr05YWJjTsUREJJPzcnnx/C3PA55vJkNCQqhVqxYAq1atcjKaZAEqypLpRSdE88V6z9nKEV9FAPDKK684GUlERLKQR6o/QqFchVh3eB2L9i6ifv36AKxYscLhZJLZqShLpvdN+DfEJsZSMLYgZ3ee5d5776Vt27ZOxxIRkSwiwCeAJ+o8AcCI1SPOF+UffviB+Ph4J6NJJqeiLJma27r5dM2nAJyYc4LAwEBGjhyJMcbhZCIikpX0rN0TH5cPs7bPonK9ypQrV46tW7fqqq7yr1SUJVP7Zdcv/HnyT4gC713ejB8/npIlSzodS0REspgigUVod2M73NbNt1u/5fvvvwdgzJgxxMXFOZxOMisVZcnUhi4d6rmzBr6f8j2tWrVyNpCIiGRZT9Z5EoAv139J5WqVqV27NvHx8fz6668OJ5PMSkVZMq1dp3bx675fIRmaFGpC69atnY4kIiJZWL0S9ahZtCYnz51kyuYp5893eeutt3C73Q6nk8xIRVkyrY+WfQQG2AxD3xjqdBwREcnijDHnjyp/svoTnnzySYoWLcratWuZOHGiw+kkM1JRlkwpMSWRCVsmAFAhqgLVqlVzOJGIiGQHHat2JH9AftYdXsf26O288847APTv35+jR486nE4yGxVlyZR+2v4TpxNPw1EI9Ql1Oo6IiGQTAT4BPFTtIQC+Wv8VDz/8MPXq1ePgwYM0a9ZMV+uTv1FRlkxp1MpRnjvroVHDRs6GERGRbKVbzW4ATNg0gYSUBGbOnEnFihUJDw+nS5cuDqeTzERFWTKdPaf2sHD/QkiGmwNupn///k5HEhGRbOSmwjdRp1gdohKimBYxjZCQEH744QcAlixZ4nA6yUxUlCXTeXj4w2DAb48f0ydMx8fHx+lIIiKSzXQL8xxV/mrDVwDExMQAUL58eccySeajoiyZyr79+1gWuwyAt9u8TbFixRxOJCIi2VHHqh0J8A7gt72/sevUrvNXfNWVXyUtFWXJVEb/PBryQkBCAH1b93U6joiIZFPB/sG0u7EdAN9u/JZKlSoBsHnzZs6dO+dkNMlEVJQl00hMTGTc2nEA3BJ0Cy6j/3mKiMi10+Umz4l7kzZPIl++fNSsWZOEhASNU5bz1EQk03ji6Sc4UuAIAK+1es3hNCIikt3dEXoHhXMXZsepHaw9tJamTZsCMH/+fIeTSWahoiyZwrFjxxizeAwEQIWgCtxe5XanI4mISDbn5fKiY9WOAEzcNJEmTZoAMG/ePCdjSSaioiyZwuzZsyH14nuP133c2TAiIpJjdK7WGYDJWyZTt15dcufOzZYtWzhy5IjDySQzUFEWx6WkpDByzEioBAZDp2qdnI4kIiI5RJ1idSifvzxHYo+w4tAKbrnlFgCWLVvmcDLJDFSUxXEjRoxgfex68Ib6xetTIk8JpyOJiEgOYYyhc1XPUeVJmyfRqJHnarBLly51MpZkEirK4qj169d7rrxX2fO4fbX2zgYSEZEcp/2Nnn97Zm6fSf2G9QEVZfFQURbHREdH06FDBxJtIt6VvQG4/4b7HU4lIiI5TZVCVaiQvwIn4k6QXCwZHx8fNm7cSHR0tNPRxGEqyuKYV155hZ07d1K6cWmSTTI3F7+ZksElnY4lIiI5jDGGByo/AMDcPXOpVasWbreb33//3eFk4jQVZXHEpk2bGDVqFC6Xixvb3ghAm8ptHE4lIiI51V9Fefq26dx2+20AfPvtt05GkkxARVkc8dZbb+F2u+nZpyfLjnvOLP7rl5SIiMj1VrtYbUrkKUFkdCQN2jbA19eXSZMmcfToUaejiYNUlMURcXFxAARWDSQ6IZpqIdUon7+8w6lERCSnchnX+fNkVpxeQbVq1bDWsnfvXmeDiaNUlMURd999NwCzImYB0LJiSyfjiIiI0KpSKwDm7JhD2bJlAdi6dauTkcRhKsriiPLlPUePd7l2AdCiQgsn44iIiNCwVENy++Rm07FNlLzRc3L5li1bHE4lTlJRlusuPj6e3r17Q35IDk4mn38+6pWo53QsERHJ4fy8/bir7F0ArDyxEoAyZco4GUkcpqIs193ChQs5cOAAIQ1DAGhavineLm+HU4mIiEDz8s0B2JLoOZIcFhbmZBxxmIqyXHdz5swBIDAsEIAW5TXsQkREMofmFTxFOapAFHhBlSpVHE4kTrpkUTbGjDXGHDPGbE6zbIgxZpsx5g9jzHRj/q+9+w6zqjzUNn6/U5ih16F3QQRURFEDIoIiomCviTWaT83RJH5JjOn1nC8mJyYaNeZoomiMLTkSiSKxN6yAiAgKSu8zgNSBKfv9/tgjQR0Upuw15f5d11yz1tpr7/0Ma4b9zDvvXiu02cN9l4QQ3g4hzA4hzKjJ4Kq/XnjhBciFFdkrCATG9xufdCRJkgDo2bon+7fZH/Kg1eBWtGlTacVRI7E3I8qTgE82mSeBA2OMBwMLgO99xv3HxBgPiTEOq1pENSRlZWUsWLAAekBJqoShXYZS0Lwg6ViSJO0ytGV6ukWLIS0STqKkfW5RjjG+AGz4xLYnYoxlFauvAt1rIZsaoFWrVlFSUkLzA5sDcGzvYxNOJEnSx3XZ0QWAnV13JpxESauJOcqXAo/v4bYIPBFCmBlCuLwGnkv13MqVK9MLvdOfju1jUZYk1S05q3IgBRvyN7Bl55ak4yhB1SrKIYQfAGXAX/ewy8gY46HAicBVIYRRn/FYl4cQZoQQZhQWFlYnluqwW2+9FfJge9vt5GTlMLLnyKQjSZL0MV3bd4VVEEPkucXPJR1HCapyUQ4hXAJMBM6PMcbK9okxrqz4vA6YDByxp8eLMd4eYxwWYxxWUOCc1Ybo1Vdf5a9//Ss5++UQQ+SIbkfQMq9l0rEkSfqYK664gtYbWgNww8M3JJxGSapSUQ4hjAe+A5wSY9y+h32ahxBafrQMjAPmVravGoff/e53AAw9I/0mCecnS5LqombNmvGdc74DwIvLX2Tp0qUJJ1JS9ub0cPcDrwADQggrQgiXAbcALYEnK0799seKfbuGEKZW3LUT8FII4S3gdeCxGOO0WvkqVOelUimeeeYZAHZ2Sr854pjexyQZSZKkPbrmjGsIqUCqU4rv/fSzTu6lhuxzL4cWY/xiJZv/vId9VwEnVSwvAoZUK50ajBdeeIGioiJ67teTeRvnkRWyOLLbkUnHkiSpUs1ymzG43WDmfjiXl5e/nHQcJcQr8ykj7r33XgCO+eIxlKXKOKjjQc5PliTVaUf3ORqAlWElqVQq4TRKgkVZta68vJx//OMfAHQY2gGAET1GJBlJkqTPNXq/0QCUdS7j9ddfTzaMEmFRVq17++23Wb9+PT179uSDnR8AFmVJUt03vPvw9EJ3eHjyw8mGUSIsyqp1bdq0AWDturVMXzYdsChLkuq+Hq17UJBXAPkw6dFJlJaWJh1JGWZRVq3r3bs3J510Ejub7mT9jvV0bN6RPm36JB1LkqTPdUzf9BmaCvMKmTx5csJplGkWZWXE1772NeiSXh7WdRghhGQDSZK0F47sXnGGpq5w8803JxtGGWdRVkaMGjVqV1Ee2nlosmEkSdpLH71mZXXL4qWXu7XvwwAAIABJREFUXmL+/PkJJ1ImWZSVESGEXUX50C6HJhtGkqS9NLRLuiiHzgECvPXWWwknUiZZlJURJSUlFmVJUr3Trmk7erXuRXlWOXSAtWvXJh1JGWRRVka8Ou9VaA5ZO7Po1bpX0nEkSdprH40q0xny8/OTDaOMsigrI55+52kA2pW08418kqR6Zdd7a7pAXl5esmGUURZlZcRrS14DoHfT3skGkSRpH+0qyp1h0aJFyYZRRlmUlRFvr3kbgJH7j0w4iSRJ+2Zwx8HphQ7w5ptvJhtGGWVRVq3bsmULG7M3AnDK8FMSTiNJ0r7p1boXedl50ApmvD0j6TjKIIuyal1+03woSC+33Nky2TCSJO2j7KxsDuhwAABrytZQWFiYcCJlikVZtW7l1pWQC2yBGS/6m7gkqf4ZWDAwvVDguZQbE4uyat1Tbz2VXiiEPn36JBtGkqQqGNjh30V59uzZyYZRxliUVesmPToJgP5t+3PCCSckG0aSpCrYVZS96EijYlFWrVtYtBCA4w87PuEkkiRVTXvapxfawfnnn59sGGWMRVm1bn1qPQDDBw5POIkkSVUz+9n0dIus9lkcPOTghNMoUyzKqn1t05+6Ne2WbA5JkqrokEGHwDZIZaVYvnF50nGUIRZl1aryVDmpVikAWpZ7ajhJUv00ZswY8nfkA3D3I3cnnEaZYlFWrVr24TJidoQt0KtLr6TjSJJUJSEEBnUdBMC9j92bcBplikVZteqpGelTw+UX51NQUJBwGkmSqm70kNEALCxc6KWsGwmLsmrV9LnTAejUtFPCSSRJqp4BHQekF9rCSy+9lGwYZYRFWbVqc9wMQPPy5gknkSSperq1rHhTekvIzs5ONowywqKsWpXbPheADcs2JJxEkqTq6ZDXIb3QAo444ohkwygjLMqqVSV5JQCsfX8tpaWlCaeRJKnqlsxdAkB222wOO+ywZMMoIyzKqlXvLHsHgI75HcnJyUk4jSRJ1bANSEF5XjmlKQd/GgOLsmrVkvVLAPi/l/1fQgjJhpEkqRp27tgJW9PLa7auSTaMMsKirFoTY6Q0P/0b9/iR4xNOI0lS9WzYsGFXUV69ZXWyYZQRFmXVmm2l2yAbKIV1K9clHUeSpCpLpVLcfvvtu4ry2m1rkw2kjLAoq9as3Vzxn8h2yMvLSzaMJEnV8PDDDzN//nyaZTUDYGPxxoQTKRMsyqo1D/3zIQDyYh4jR45MOI0kSVX3j3/8A4ChA4YCsKHY0542BhZl1YoYI5MenARA7469ycryW02SVH/17dsXgLg9AhblxsL2oloxa9YsFqxYAMDA3gMTTiNJUvWMGDECgFWLVgEW5cbCoqxa0alTJ3Japs+bXLalLOE0kiRVT4zpkeTNazYDsGGHRbkxsCirVnTv3p0xJ44B4O033k44jSRJVRdj5Nvf/jYAE4+bCDii3FhYlFVrDjrsIABWLl7Jzp07E04jSVLVlJeXM2/ePAAu/uLFAGwr2ZZkJGWIRVm1p+KK1WXFZbz55pvJZpEkqYpycnJo1ix9WrgmoQkA20u3JxlJGWJRVq3ZVlrx23YJLF26NNkwkiRVQ+vWrQEo3Z6+4qxFuXGwKKvW7CrKpbB8+fJkw0iSVA1du3YF4MPCDwGLcmNhUVat2TV/qxTWrFmTbBhJkqph8ODBAEy6YxJgUW4sLMqqNes2rksvlMLYsWOTDSNJUjX87Gc/o3nz5kz53ymARbmxsCir1nyw+AMAThh7AuPHj084jSRJVde7d2++/vWvQ8WlAYrLipMNpIywKKvWFO9M/ycy7vhxCSeRJKn6tm7dCvHf6x9dhEQN114V5RDCnSGEdSGEubttaxdCeDKEsLDic9s93Pfiin0WhhAurqngqttijOwsSZ87uUO7DgmnkSSp+j56v00gAFAey5OMowzY2xHlScAn/3b+XeDpGGN/4OmK9Y8JIbQDfgIcCRwB/GRPhVoNywMPPMCOkh0AdO/aPeE0kiRVX9u26QqTVVGfUjGVZBxlwF4V5RjjC8Anr9V4KnB3xfLdwGmV3PUE4MkY44YY40bgST5duNXAbNq0iW984xtU/MJNs6bNkg0kSVINKCgoAHYbUU45otzQVWeOcqcY4+qK5TVAp0r26QbsfgLdFRXb1IC99957FBYWktMkfWm+nKychBNJklR9WVnp2vRRUXZEueGrkTfzxfRs9mrNaA8hXB5CmBFCmFFYWFgTsZSQQw45hDZt2lBWVpZ0FEmSakxJSQngHOXGpDpFeW0IoQtAxed1leyzEuix23r3im2fEmO8PcY4LMY47KM/bah+atKkCSeffDJU/KLtn6YkSQ3BnDlz0gsVUwuzgicPa+iqc4SnAB+dxeJi4JFK9vkXMC6E0LbiTXzjKrapgevTp8+uolyaKk02jCRJ1VRUVMSTTz5JdnY2MaT/iJ6blZtwKtW2vT093P3AK8CAEMKKEMJlwPXA8SGEhcDYinVCCMNCCH8CiDFuAH4BvFHx8fOKbWrgOnXqtKsol6WcgiFJqt+eeOIJysrKGD1m9K4BoNxsi3JDt1fvsooxfnEPNx1Xyb4zgK/stn4ncGeV0qne6tixo0VZktRgtGnTBoDtO9KXrs4O2U69aAQ8wqoVLVu2hIqpyaXlTr2QJNVvX/jCF8jNzeWV114BHE1uLCzKqhV5eXnOUZYkNRjt2rXjBz/4wa7m5PzkxsGirFoxc+ZMqOjHxaXFyYaRJKkGXHfddeS3yAe8RkBjYVFWrXj00UchfbpJtpZsTTaMJEk1ID8/n8OPOhyA7LLshNMoEyzKqhVr1qyxKEuSGpye/XoCkJNyRLkxsCirVgwdOtSiLElqcDbt2ARAs+xmCSdRJliUVSuaNGmyqyhvKdmSbBhJkmrIpuJ0UW6e2zzhJMoEi7Jq3KJFi7j77rvJKkt/ezmiLElqKFYVrQKgXfN2CSdRJliUVeNSqfR54VrltwJg887NScaRJKlGxBhZWbQSgC7tuyScRplgUVaN69u3L82aNePDVR8CsL54fcKJJEmqvrKyMnawA4CClgUJp1EmWJRV47Kyshg3bhykr/JJ0faiZANJklQDcnNzaVaQfhNfXllewmmUCRZl1Yrf//73/y7K2yzKkqT6b86cOWyveHHr3r57wmmUCRZl1YomTZr8uygXW5QlSfXfrbfeChVnhevRrkeyYZQRFmXVimXLlsFOIJV+M19JeUnSkSRJqpb33nsPKs4KV9DMOcqNgUVZtaKgoIBAcJ6yJKnB2LRp064R5YLmFuXGwKKsWtG7d28uvPBC2JZeX7t1bbKBJEmqCRVFuUOzDsnmUEZYlFVrfv7znxO2BgBmL56dcBpJkqrnxFNPhKaQFbNo37R90nGUARZl1ZpevXrRq20vAP753D8TTiNJUvUcd9pxAMQPI5s3ezGtxsCirFo1dL+hAMxdNjfhJJIkVVP6grPETZF333032SzKCIuyalWb7DYAlDcrTziJJEnVs3zz8vTCZigv93WtMbAoq1blbM8BoLRpacJJJEmqnhWbV6QXNkEIIdkwygiLsmpV05KmAGzL2ZZwEkmSqmdR0SIAwtbAQQcdlHAaZYJFWbVqQKcBAHwYPqQ85Z+pJEn115xlcwDYv+P+tGjRIuE0ygSLsmrVfj32gy2QCilWblmZdBxJkqps2dZlAAzuMjjhJMoUi7Jq1fr162FjevmDDR8kG0aSpCoqLS9lffl6iFBa6PtuGguLsmrVzJkzYUN6+f0N7ycbRpKkKlry4RJSpGATbFq/Kek4yhCLsmrVmDFjdhXlN5e+mWwYSZKqaNdgzwYYNGhQsmGUMRZl1aqJEydySK9DAHjlvVcSTiNJUtXsXpSPPfbYZMMoYyzKqnWH9z0cgJU7fTOfJKl+mvbGtPTCBhg9enSiWZQ5FmXVuoO7HAzAetZTWu4bICRJ9cvtt9/O1DemAnDayNMoKChIOJEyxaKsWnfMiGNgY/oUcQvWL0g6jiRJe23atGl89atfhY7p9d9993fJBlJGWZRV6wYPHkzuxlwAnp//fMJpJEnaO+vXr+e8884j1TQFLaBFkxb0atMr6VjKIIuyal1WVhY9m/YE4Om3n044jSRJe2f69Ols2rSJA0YdAMCBHQ8khJBwKmWSRVkZ0a9VPwAWfrgw4SSSJO2defPmAdD5oM4AHNTxoCTjKAEWZWVE/9b9AVhRuiLhJJIk7Z01a9YAsL3ldiA9oqzGxaKsjBjcaTCUw8awkS07tyQdR5Kkz9WtWzcAFm1fBFiUGyOLsjKie+fusBYI8OYar9AnSar75s+fD9mwsclGAA7tcmjCiZRpFmVlRIcOHWBVennGqhnJhpEkaS8sXboUCqA8lNO/XX/a5LdJOpIyzKKsjFi+fDmsTi/PXD0z2TCSJH2OGCOLFy+Grun1YV2HJRtIibAoKyNee+21XSPKM1dZlCVJddvLL7/M4sWLabpfU8Ci3FhZlJURCxcuhHWQE3J4b/17bN65OelIkiTt0WOPPQZAiwEtAItyY2VRVkZs374dymG/5vsB8PrK1xNOJEnSng0aNAhyoCiriEBgaOehSUdSAizKyoiPrmTUPy99PuXpy6YnGUeSpM90+OGHQ1eIWZFBHQbRMq9l0pGUAIuyMuKoo44CoOSDEgCmL7coS5LqptLSUi677DLomV4f1XtUsoGUGIuyMmL8+PEAvDklfQ7lV1a8QlmqLMlIkiRV6k9/+hPTp08nb/88AI7ueXTCiZQUi7Iy4ogjjuDggw+mcFEhHbI7sLVkK2+vfTvpWJIkfUr79u0hQFmX9IDOyJ4jE06kpFS5KIcQBoQQZu/2sTmEcM0n9hkdQti02z4/rn5k1UchBL73ve8BUPpBKeD0C0lS3TRhwgSa9mxKeW45XZt3pUfrHklHUkKqXJRjjO/FGA+JMR4CHAZsByZXsuuLH+0XY/x5VZ9P9d+ZZ55J27Zt2fT2JgBeWPpCwokkSfq05s2b02tULwB6fjRRWY1STU29OA74IMa4tIYeTw1Qbm4up556KixOrz+75FlSMZVsKEmSPqGwsJAFJQsAGD9wfMJplKSaKsrnAffv4bbhIYS3QgiPhxAG19DzqZ764he/CBuATVC0vYg5a+ckHUmSpI+59Q+3kuqVHsg5f/j5CadRkqpdlEMITYBTgL9VcvMsoFeMcQhwM/CPz3icy0MIM0IIMwoLC6sbS3XU8ccfz1VXXQWL0ut/feWvyQaSJGk35eXl/PEff4Rm0DmvM/u13S/pSEpQTYwonwjMijGu/eQNMcbNMcatFctTgdwQQofKHiTGeHuMcViMcVhBQUENxFJdFELgpptuYli79KVA/+eJ/yHGmHAqSZLSnnnmGdY2T1eaCYMm7LpglhqnmijKX2QP0y5CCJ1DxXdYCOGIiudbXwPPqXosOzubB3/1IABb2m3h9ZlezlqSlLzi4mK+//3vQ9/0+vF9j082kBJXraIcQmgOHA88vNu2K0MIV1asngXMDSG8BfweOC86fCigb0Ff2pa2hSZw15N3JR1HkiT+4z/+gxmzZ0Dv9PpxfY9LNI+Sl1OdO8cYtwHtP7Htj7st3wLcUp3nUMO1f87+vMZrzN05N+kokqRGrqioiHvuuYfsftmUZ5cztPNQOjSrdLaoGhGvzKfEDModBMDCuDDhJJKkxm7q1KmkUim6jekGwEn9T0o4keoCi7ISMyB/AJTAdraztWRr0nEkSY3YP//5TwC2ddsGwMT9JyYZR3WERVmJaZHfAv4AFxReQIsmLZKOI0lqxFasWAEFsD61noJmBRzR7YikI6kOqNYcZak68vPzySvOIyv4+5okKVnt2rXb1Yom7D/B1yYBFmUl6LLLLuOyyy5LOoYkSbRv3x46p5dP3v/kZMOozvDXJUmS1OjtzN4J3SGHHM+frF0sypIkqVF75JFHmPzuZMiCIzsdScu8lklHUh1hUZYkSY3WLbfcwmmnnUZp/1IALj3y0oQTqS5xjrIkSWq0Jk+eDM0g9A1kZ2Vz2gGnJR1JdYgjypIkqdFq3bo1HAAxRI7rcxztmrZLOpLqEEeUJUlSo9W2bVuoOJX/OYPPSTaM6hxHlCVJUqMUY+SVt16BPpCN0y70aRZlSZLUKM2aNYv5WfMhC47r67QLfZpFWZIkNUr33XcfHJxevviQi5MNozrJoixJkhqlh59/GHpA0+ymTrtQpSzKkiSp0bnnnntY0moJAGcPPptmuc2SDaQ6yaIsSZIalZdeeonLvnLZrmkXlxxySaJ5VHdZlCVJUqOxevVqTj/9dMq6lkFb6NGqB8f0PibpWKqjLMqSJKnRuP766ykqKqLziZ0BuODgC8gK1iFVzu8MSZLUKNx111384Q9/gHzY2HUjAJcOvTThVKrLLMqSJKnBe/7557n00kspKyvjhGtPYGdqJ2P7jqVfu35JR1MdZlGWJEkN3k9/+lMArvvudazsshKAyw+9PMFEqg8sypIkqUFbt24dzz33HM2aNWPMhWOYu24uHZt35NQDTk06muo4i7IkSWrQPvzwQwC6du3K/QvuB+DSQy6lSXaTJGOpHrAoS5KkBm3p0qUA5LfL58F3HgTgK4d+JclIqicsypIkqcFKpVJce+21ALQ7vh07ynZwUv+T2K/dfgknU31gUZYkSQ3WtGnTeOutt+jeszsL2y4E4Jojr0k4leoLi7IkSWqwbr31VgBGXjGS1VtXM6hgEGP7jk04leoLi7IkSWqQFi1axOOPP05eXh4L2i4A4BtHfoMQQsLJVF9YlCVJUoP0hz/8gRgjx150LLPWzaJd03ZccPAFScdSPWJRliRJDc68efO4+eabAdh52E4A/s+h/4dmuc2SjKV6xqIsSZIanKuuuoqSkhLOuvIsnl3zLE2ym/D1I7+edCzVMxZlSZLU4Lz55psAZI/KJhK5eMjFdG3ZNeFUqm8sypIkqUEpKysjxgit4OH3HyYrZPGdo76TdCzVQxZlSZLUoDzzzDNs3ryZNie1oTRVylmDzqJfu35Jx1I9ZFGWJEkNylNPPQVNYduAbQB896jvJpxI9ZVFWZIkNShvvvkmDIfSUMr4fuMZ2mVo0pFUT1mUJUlSgxFjZOb8mXBkev3Ho36cbCDVaxZlSZLUYKxYsYKNAzdCHpzY70SG9xiedCTVYxZlSZLUYDzz2jNwRHr5Z6N/lmwY1XsWZUmS1GDc/ObN0AT2K9uPw7sdnnQc1XMWZUmS1CCs2ryKWVmzAPje8O8lnEYNgUVZkiQ1CF/7368RcyL5i/L58vgvJx1HDUBO0gEkSZKqa37hfCYvnQwRzmx7JllZjgWq+vwukiRJ9d5Vk68ihgiz4Iozr0g6jhoIi7IkSarX/vzkn3l29bNQAsdlH8dRRx2VdCQ1ENUuyiGEJSGEt0MIs0MIMyq5PYQQfh9CeD+EMCeEcGh1n1OSJAng9ddf54r/TY8g9y/qz2MPPua0C9WYmpqjPCbGWLSH204E+ld8HAncxq7r5UiSJFXNW2+9xej/GE35yeXklebxym9eIS8vL+lYakAy8SvXqcA9Me1VoE0IoUsGnleSJDVg3/j2Nyg+uhiAG06+gfYt2yecSA1NTRTlCDwRQpgZQri8ktu7Act3W19RsU2SJKlK/v73v/N86fPQBg7scCBXHn5l0pHUANXE1IuRMcaVIYSOwJMhhHdjjC/s64NUlOzLAXr27FkDsSRJUkO0YsUKLvvWZXBhev3WibeSnZWdbCg1SNUeUY4xrqz4vA6YzK4rrO+yEuix23r3im2ffJzbY4zDYozDCgoKqhtLkiQ1UDfeeCObj9wMuXDe4PMY1WtU0pHUQFWrKIcQmocQWn60DIwD5n5itynARRVnv/gCsCnGuLo6zytJkhqnGCP3vHgPDIb8rHx+ffyvk46kBqy6Uy86AZNDCB891n0xxmkhhCsBYox/BKYCJwHvA9sBrykpSZKq5KVXXqJwWCEAPxj1A3q07vE595CqrlpFOca4CBhSyfY/7rYcgauq8zySJEkA3/z7N6EjtC5vzbeP+nbScdTAeUZuSZJUL7w07yVmNE9f2+x3x/2O/Jz8hBOpobMoS5KkOi+VSnHWpLMgB3p92IsvH+NMTtU+i7IkSarzvvL7r7C2+VrYDvddfF/ScdRIWJQlSVKdtn77ev5S+BcALutxGSMOGZFwIjUWFmVJklSnXfTARZQ1KYPF8JsLfpN0HDUiNXFlPkmSpBqXSqW49s/XMnXVVCiBE8tOpE2bNknHUiNiUZYkSXVOcXExo8aPYsYRM6AFHLjmQB7+08NJx1IjY1GWJEl1zk033cSMjumSvH+T/Zn5x5k0yW2SdCw1Ms5RliRJdcpbb73Fz/72MzgwfZnqx6983JKsRFiUJUlSnbFkyRKOP/14dhy3A4Abxt9A37Z9E06lxsqiLEmS6oTt27dz0oSTKBxRCM1hTK8xXDnsyqRjqRGzKEuSpDrhO9/5DvNbzof+0Da/LX854y9kBauKkuN3nyRJStx9993HrX+/FY5Pr9956p10a9Ut2VBq9CzKkiQpUfPmzePSKy6FM4EcuPzQyzntgNOSjiVZlCVJUrImTZrEzlE7oSMMaD+A357w26QjSYBFWZIkJWzK+1PgCMgJOdx/5v00b9I86UgSYFGWJEkJmrt6Lu8d8B4A/zXqvxjaZWjCiaR/syhLkqREbC/dzmn3nQZ50GpZK6495tqkI0kfY1GWJEmJuHrq1Xyw9QNYD6dln0YIIelI0sdYlCVJUsbd9eZd3DX7LkJ5gIdg9PDRSUeSPsWiLEmSMuqNlW9w+ZTLAYj/jHTN7srEiRMTTiV9mkVZkiRlzKJ1ixh922jKKIMZcM7+5zBnzhwKCgqSjiZ9ikVZkiRlxHMvPsdBvziI7bnbyVqRxaRzJ/Hggw/Svn37pKNJlcpJOoAkSWr4rrnmGm764CYYBjnbcnjyiicZPWx00rGkz+SIsiRJqlVr1qzhpukVJTnm8PTlT1uSVS9YlCVJUq26/qHr4cT08l1n3MWofqOSDSTtJYuyJEmqFfPmzeO4Lx7HTatvgmwYnTuaCw6+IOlY0l6zKEuSpBo1f/58JkyYwOAjB/NM52cgHwYxiH9d+6+ko0n7xDfzSZKkGlNYWMi4ceNYsXYF4cuB2CYytONQpn9lOk1ymyQdT9onjihLkqQaUVRUxMSJE1mxcgXtvtKO2DXSu01vpl00jaa5TZOOJ+0zR5QlSVK1LV++nLFjx7JgwQJantOSDZ020DqvNVO/NJWOzTsmHU+qEkeUJUlStaRSKS688EIWLFhA53M6s2XQFppkN2HyuZMZWDAw6XhSlVmUJUlStdxyyy08//zztBzdkjWD1pAVsrj/zPsZ02dM0tGkarEoS5KkKlu4cCHf/e53YRBsHb0VgNsm3MYZA89IOJlUfRZlSZJUJVu3buWss86iuHMxWWdnEYn855j/5PLDLk86mlQjLMqSJGmfpVIpLrjgAuZsnEP4UiAVUnz9iK/z/aO/n3Q0qcZYlCVJ0j77wQ9+wCMzHoELIOZGLhpyEb8b/ztCCElHk2qMRVmSJO2Tv/zlL1x/9/VwAZAH5x14HneecidZwVqhhsXvaEmStNdmz57Npd+7FC4EmsLpB5zOPafdQ3ZWdtLRpBpnUZYkSXvtt3/5LWVfKoNmMHH/iTxw1gPkZucmHUuqFV6ZT5Ik7ZWn5j7FfU3ug3wY2moofzv7bzTJbpJ0LKnWOKIsSZI+122P3Ma4v46jPL+cJsua8NiFj5Gfk590LKlWWZQlSdJnemzOY1z12lXE/EjBhgLm/XgeXTp0STqWVOuceiFJkipVWlrK7x/5Pde9dR0xL9KxsCPLfreMvNy8pKNJGeGIsiRJ+pgdO3Zw44030vGojnx79rcpzyknd34ur3z7FUuyGhVHlCVJ0i7vvvsu48aNY3nb5XAqkA1D41Ae+q+H6Nu7b9LxpIyq8ohyCKFHCOHZEMK8EMI7IYRvVLLP6BDCphDC7IqPH1cvriRJqi3r169n4sSJLO+2HM4AsuG6o65j5k9m0m+/fknHkzKuOiPKZcC3YoyzQggtgZkhhCdjjPM+sd+LMcaJ1XgeSZJUy5555hm++a1v8kHvD+Do9LYbxt3AN4d/M9lgUoKqXJRjjKuB1RXLW0II84FuwCeLsiRJqqNijFxwwQXc98B9MBE4FHJCDneeeicXDrkw6XhSomrkzXwhhN7AUOC1Sm4eHkJ4K4TweAhhcE08nyRJqhkvvPAC9/3vfWRdlAWHQtOcpjzyxUcsyRI18Ga+EEIL4H+Ba2KMmz9x8yygV4xxawjhJOAfQP89PM7lwOUAPXv2rG4sSZK0B2VlZbz88sssWLCA39/ze7gMUh1TdGzekSnnTeHI7kcmHVGqE0KMsep3DiEXeBT4V4zxt3ux/xJgWIyx6LP2GzZsWJwxY0aVc0mSpMpt2bKFCRMm8OKLL0JX4EtACxjQdgDTLppG7za9E04oZU4IYWaMcdiebq/OWS8C8Gdg/p5Kcgihc8V+hBCOqHi+9VV9TkmSVHUbN27kpJNO4sUXX6T1F1qT/ZVsaAGHtz+cVy9/1ZIsfUJ1pl4cBVwIvB1CmF2x7ftAT4AY4x+Bs4CvhhDKgGLgvFidIWxJklQlL774IhdccAHLli+j1cmt2HTYJgAuG3oZt024jdzs3IQTSnVPdc568RIQPmefW4BbqvockiSpesrKyvjFL37Bf/7nf5LKSdH2irZs7LyRQOCXx/2S7xz1HSr++CvpE7wynyRJDdTixYs5//zzeeWVV6AddLi6A0VZRbTOa839Z97Pif1PTDqiVKdZlCVJaoAmT57MJZdcwubNm2l/RHtKTimhqKyIAzocwCPnPcL+7fdPOqJU51mUJUlqYNatW8f5559P8Y5iBl4xkPe6vkeqLMXJ+5/MvWfcS6u8VklHlOqFGrngiCRJqjtuvPFU+mqXAAARAElEQVRGikMxBd8oYH6X+aRiih+N+hH/OO8flmRpHziiLElSAxFjZOrUqdz08E1wBRS2LqR90/bce8a9jO83Pul4Ur1jUZYkqQHYtm0bE0+eyHPFz8E5QDYM7z6cB896kB6teyQdT6qXLMqSJDUA1/zgGp7r/BwMSK9/bdjXuGH8DZ4fWaoGi7IkSfXUjh07uP3227nj2TuY228uDICWuS2ZdPokzhh4RtLxpHrPoixJUj304YcfcvLpJ/NS3kswPL2tX04/nrnqGadaSDXEoixJUj2zatUqjjn7GN4/+H3oDNlk88OjfsiPjv0R2VnZSceTGgyLsiRJ9cjK1Ss55OpDKDq2CHKgd8vePHTuQxze7fCko0kNjudRliSpHti6dSvX//l6+v2yH0VD0iX5wkEX8vbVb1uSpVriiLIkSXXY1q1b+a9f/hc3vHwDpSNLoT3kbMvh3nPu5dxh5yYdT2rQLMqSJNVBO3bs4J577uFHN/+IdV9YB6PT249qehQPXPUA3Tt0TzSf1BhYlCVJqiNeffVVHnzwQYqKinjimSdYN3AdnA5kQ0FeAXefeTcn9j8x6ZhSo2FRliSpDpg/fz5jxoxhx44dsB9wFtAOAoErDruCX479JW3y2yQdU2pULMqSJCWkuLiYW2+9laVLl/LUU0+xI3sH3b7WjZXtVwJwUMeD+J+J/8PwHsMTTio1ThZlSZISsHnzZk455RSef/759DmoDodwemBl3kqa5jTlJ8f8hG8O/6aXoJYSZFGWJCnD1q9fz4knnsgbb7xB+8Pbk3NyDmtTa4lETux3IrecdAt92/ZNOqbU6FmUJUnKoNWrVzNu3DjmrphLsy83Y32v9ZCC/drux43jb2RC/wmEEJKOKQkvOCJJUkaUlZXx29/+lsNGHMbcTnMJVwe299pO89zm/PK4X/LOf7zDxP0nWpKlOsQRZUmSMuAHP/4Bv37m13AO0AwikfMPOp9fjf0V3Vp1SzqepEpYlCVJqkU7S3Zy7aRruXnbzVBxCuSjex7Nfx//3xzZ/chkw0n6TBZlSZJqwfbt27ni11fwYOGDlHYshXbQrrwdk86f5BQLqZ6wKEuSVINijEx9dyrn33E+m1pvgo6QvT2bs9qfxaRrJpHfJD/piJL2kkVZkqQa8sLSF/jh0z/kxeUvQmsIxYGvDPwKvz3vt7TIa5F0PEn7yKIsSVI1xBh5dsmz/L8X/x9PL346vbEYWr7dkud+/RyHDj402YCSqsyiLElSFaRiikfefYTrp1/P6ytfB6BJqgklz5fQal4rpj8znQMHH5hwSknVYVGWJGkflJSXcN/b9/Gr6b/i3aJ3AWif355ea3sx67ZZ5MU8pvxrCgceaEmW6juLsiRJe2FD8QbumHkHt7xxCys2rwCgc35neq3sxey7ZjNr6yzy8/N5ZMojHHPMMQmnlVQTLMqSJO3BkiVLmL5gOk9vfZoH5j9AcVkxAB1SHQgvB9Y8s4Y1qTUAnHDCCfzkJz9h+PDhSUaWVIMsypIkfUJ5qpwbH7uR6x6+jvLe5bu2tylqw4fTPqTogyKI0LlzZy644AKuuOIK+vXrl2BiSbXBoixJUoWVm1dyx8w7uPXlWykqK4LeEEoDYU4g9UqKD4s+JC8vj/MuOo8vf/nLjBw5kuzs7KRjS6olFmVJUoNUUlJCTk4OWVlZld6+Y8cOXn/9dea8PYcnFz/JWzlvsSx/GTHE9A4bYPCOwTz7m2dpndea+fPns2jRIkaOHElBQUEGvxJJSbEoS5IahI0bN/L444/z1FNPMWPGDObNm0eLFi0YNWoURx11FI8//jjLli1j7NixLFy4kOnvTad0YCkcDLSpeJBy4F3oWdiTa8+6lisuv4Lc3FwAhgwZwpAhQ5L68iQlIMQYk87wKcOGDYszZsxIOoYkqQ4rKyvjnXfe4emnn2bKlCm89NJLlJf/ez5xCIFPvcY1Aw4EhgDd/r25XWjHMc2PYUSzERw+6HBGjRpFCCETX4akBIUQZsYYh+3pdkeUJUn1yoYNG/jyl7/Mk08+SXFx8a7tOTk5HHvssUyYMIHhw4czZMgQCgsLmfbsNB6a/RAr267k/fA+5TFdplvktuCcwedw4ZALGdVrFFmh8ikakhovi7Ikqd6IMfLVr36VKVOmANC3b19GjBjBhAkTGD9+PG3apOdQbN65mcnvTeZv8/7GtBXT2Nl2JwDZZHNS/5O48OALOWXAKTTLbZbY1yKp7rMoS5LqtLVr1/Lggw8yZ84cZs+ezcyZM2nevDlvvPEGAwcO3LXfum3ruHv23Tz87sP86/1/sbM8XY4DgZE9R3L2oLM5Z/A5dG7ROakvRVI9Y1GWJNVZM2bM4OSTT2bNmjW7tuXm5nLbbbdxwAEHMGftHB5d8Cj/XPBPXlvxGpH0nORA4OieR3P2oLM5c9CZdG3ZNakvQVI9ZlGWJNVJkydP5vzzz6e4uJgRI0Zw3nnn0WdgHza22cgra1/hhzf9kGWblu3av0l2E8b0HsPJ+5/MGQPPoEvLLgmml9QQWJQlSXVKjJHf/OY3XHfddcSsyPgrxzPkjCHcu/ReZrw8g1RM7dq3U/NOTOg/gYn7T+T4/Y6nRZMWCSaX1NBYlCVJdcasd2bxrRu/xXMfPAcXQW6fXKYxjWkvTwMgNyuXkT1HMrbPWE7odwLDug7zbBWSao1FWZKUiFRM8cGGD3hj1Rs8/PrDPLXgKTblb4LupD+AUko5uNPBjO0zlrF9x3J0r6MdNZaUMRZlSVKtS8UUizcuZubqmbyx8g0enfUoC7YuIJX772kUNAPKoWNJRyYePJFThpzCiB4jKGju5aIlJcOiLEmqMTFGVm9dzdx1c5m7bi7vrHuHuYXpz9tKt31851xgM7AKmm5sytlHnM1PL/8pfbr3SSK6JH1KtYpyCGE8cBOQDfwpxnj9J27PA+4BDgPWA+fGGJdU5zklSdVXWlrKE088wapVq7jkkkvIzc3d6/umYoq1W9fy/ob3//2x8d/Lm3durvR+rbNbs+39bZQtLaP55ub84qu/4OCRB7NlyxbGjRtHs2Ze/ENS3VLlohxCyAZuBY4HVgBvhBCmxBjn7bbbZcDGGGO/EMJ5wK+Ac6sTWJJUPXfffTfXXnsthYWFANx1111ce+21vLfwPeYvnc/qTatp2r4p488cT0mTElZsXsGKLSvSnzevYOXmlZSmSvf4+G3z23JgxwN3fXSkI3f9+i6m/n0qABMmTOD2v95O166e21hS3RZijFW7YwjDgZ/GGE+oWP8eQIzxl7vt86+KfV4JIeQAa4CC+DlPOmzYsDhjxowq5ZKkxijGyO3P3k52m2zKUmWUlpdSmir92OdtpdvYsG0Dd99/N6mcFM3aNKM0q5TSrNL0/OCme/987Zu2p0/rPsQNkVlPzSIWRVqVteLc48/lorMuokWLFhQXF3PXXXcxadIkSktLadWqFTfeeCOXXHIJIYRa+7eQpL0VQpgZYxy2p9urM/WiG7B8t/UVwJF72ifGWBZC2AS0B4qq8bySpE+YNGkSV7585a6zRXymwelP29n+8e0Rmmc1p01uG3Zu2EnRsiLYBp2admLE4BEM7DaQQd0HsXDmQqY9PI0ZM2ZQXl4OQJ8+fVi8eDF3zLuDO26642MPG0LgnHPO4b//+7/p2bNnDXy1kpQZdebNfCGEy4HLAf8jlaR9dPXVV8MwoBBGHDmCvr37UtCugNysXMpKypg3dx6b129m8buLWb1kNRd/6WIuOvciWjZpScu8lnRo1oG2+W3Jzsre9ZhTp07l6quvZvHixUz++2QmM/ljz5mdnc2wYcP41re+xbnnnsusWbO46667eOmll4gxEkJg2LBhXHvttQwYMCDD/yKSVH1OvZCkBqBDhw6sX7/+Y9sOP/xwBgwYwN///nd27Nixa3tOTg4LFy6kd+/en/u427dv5+GHH2bevHksXLiQJUuWcNBBB3Hqqady7LHH0rJly5r+UiQpYz5v6kV1inIOsAA4DlgJvAF8Kcb4zm77XAUcFGO8suLNfGfEGM/5vMe2KEvS3osx8vrrrzN69GhatGjBqFGjeOKJJ9i6deuufcaPH88JJ5xAu3btOOiggxg6dGiCiSWpbqi1OcoVc46vBv5F+vRwd8YY3wkh/ByYEWOcAvwZ+EsI4X1gA3BeVZ9PklS5EAJHHnkkxcXFu7YVFxfz6KOP8v7773PGGWc49UGSqqDKI8q1yRFlSZIk1bbPG1HOymQYSZIkqb6wKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVCDHGpDN8SgihEFi626YOQFFCcVR1Hrf6yeNWP3nc6iePW/3kcaufKjtuvWKMBXu6Q50syp8UQpgRYxyWdA7tG49b/eRxq588bvWTx61+8rjVT1U5bk69kCRJkiphUZYkSZIqUV+K8u1JB1CVeNzqJ49b/eRxq588bvWTx61+2ufjVi/mKEuSJEmZVl9GlCVJkqSMqtNFOYRwdgjhnRBCKoQwbLftvUMIxSGE2RUff0wypz5uT8et4rbvhRDeDyG8F0I4IamM+mwhhJ+GEFbu9jN2UtKZtGchhPEVP1PvhxC+m3Qe7Z0QwpIQwtsVP2Mzks6jyoUQ7gwhrAshzN1tW7sQwpMhhIUVn9smmVGftofjts+vbXW6KANzgTOAFyq57YMY4yEVH1dmOJc+W6XHLYQwCDgPGAyMB/4QQsjOfDztpd/t9jM2NekwqlzFz9CtwInAIOCLFT9rqh/GVPyMeaqxumsS6des3X0XeDrG2B94umJddcskPn3cYB9f2+p0UY4xzo8xvpd0Du2bzzhupwIPxBh3xhgXA+8DR2Q2ndTgHAG8H2NcFGMsAR4g/bMmqQbEGF8ANnxi86nA3RXLdwOnZTSUPtcejts+q9NF+XP0CSG8GUJ4PoRwdNJhtFe6Act3W19RsU1109UhhDkVf77yz4p1lz9X9VcEngghzAwhXJ50GO2TTjHG1RXLa4BOSYbRPtmn17bEi3II4akQwtxKPj5rRGQ10DPGOBT4JnBfCKFVZhILqnzcVId8zjG8DdgPOIT0z9sNiYaVGqaRMcZDSU+buSqEMCrpQNp3MX36ME8hVj/s82tbTm0n+jwxxrFVuM9OYGfF8swQwgfA/oBvhsiQqhw3YCXQY7f17hXblIC9PYYhhDuAR2s5jqrOn6t6Ksa4suLzuhDCZNLTaCp7T47qnrUhhC4xxtUhhC7AuqQD6fPFGNd+tLy3r22JjyhXRQih4KM3gYUQ+gL9gUXJptJemAKcF0LICyH0IX3cXk84kypR8R//R04n/QZN1U1vAP1DCH1CCE1Iv2F2SsKZ9DlCCM1DCC0/WgbG4c9ZfTIFuLhi+WLgkQSzaC9V5bUt8RHlzxJCOB24GSgAHgshzI4xngCMAn4eQigFUsCVMcZqT9hWzdjTcYsxvhNCeAiYB5QBV8UYy5PMqj36dQjhENJ/TlwCXJFsHO1JjLEshHA18C8gG7gzxvhOwrH0+ToBk0MIkH4tvi/GOC3ZSKpMCOF+YDTQIYSwAvgJcD3wUAjhMmApcE5yCVWZPRy30fv62uaV+SRJkqRK1MupF5IkSVJtsyhLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlbAoS5IkSZWwKEuSJEmVsChLkiRJlfj/FlISk1eNyQYAAAAASUVORK5CYII="
},
"metadata": {
"tags": [],
"needs_background": "light"
}
},
{
"output_type": "stream",
"name": "stdout",
"text": [
"UKF standard deviation 0.050 meters\n"
]
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 571
},
"id": "VaHsojrmSGZX",
"outputId": "adf62250-c843-4649-d2e4-5766905c022e"
}
}
]
}