{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "[Table of Contents](table_of_contents.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Topic 22. Invariant Subspaces and Differential Equations\n", "Author: Bryan Redd, bryan.d.redd@gmail.com\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Invariant subspaces provide one way to both analytically and intuitively analyze the behavior of linear operators. Intuitively, when an element of an invariant subspace is operated on by its linear operator, the result of the operation will also be in the same subspace. That is why these spaces are \"invariant.\" Invariant subspaces allow us to decompose a linear operator into simpler components. This can be leveraged for efficiency reasons, when computing all of the information associated with the operator might be unnecessary or impractical.\n", "\n", "When working with matrices, the subspace spanned by the eigenvectors form the invariant subspaces of the matrix. Summing the projections onto each invariant subspace and scaling by the appropriate eigenvalue results in the same operation as left multiplying by the matrix itself.\n", "\n", "Differential and difference equations involving matrices can be analyzed using invariant subspaces." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Explanation of the theory\n", "\n", "### Definition of Invariant Subspace\n", "\n", "Let $A$ be a linear operator. If $S \\subset R(A)$ is such that $x \\in S$ means that $Ax \\in S$, then $S$ is said to be an __invariant subspace__ for $A$.\n", "\n", "If we know that an element of a given invariant subspace goes into the linear operator, the output will also be an element of that same subspace.\n", "\n", "### Analytical Example\n", "\n", "If $S$ is an invariant subspace of $A$, then $S$ is an invariant subspace of $e^{At} = I + At + \\frac{1}{2!} A^2 t^2 + \\frac{1}{3!} A^3 t^3 + \\ldots $.\n", "\n", "We can prove this if we let $x \\in S$ where $S$ is an invariant subspace of $A$. Then by definition of an invariant subspace, we know that $Ax \\in S$. Because $S$ is a vector space, we can multiply by a scalar $t$ and know that $Atx \\in S$. Operate over $A$ and multiply by the scalar $t/2$ to show that $\\frac{1}{2} A^2 t^2 x \\in S$. Repeat this process indefinitely to show that every term of $e^{At}x = Ix + Atx + \\frac{1}{2!} A^2 t^2 x + \\frac{1}{3!} A^3 t^3 x + \\ldots \\in S$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple Numerical Examples\n", "\n", "One way of visualizing invariant subspaces for matrices is using a phase diagram. Consider the differential equation $\\dot{x}(t) = A x(t)$. It can be shown that the solution to this equation is $e^{At}x_0$ where $x(0)=x_0$ is the initial condition and $e^{At}$ is the matrix exponential as defined above. In the Python example below, we plot the solutions to a differential equation with a two-by-two matrix for various intial conditions in red. The initial condition is shown in blue. Note their general tendency to approach the eigenvectors of the matrix $A$ for large values of t.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2020-12-07T11:17:45.154626\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.1, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEICAYAAAB8uBDgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAACurUlEQVR4nOxdZ3gUZRc9k00vQIDQe+8t1NB7ExAEBBQRRUSETykCggqCSrEAig1RqaJiQRRpofcmIErvvZMQ0nf3fj9O3sxmszMpJBAk53nmSXZnduo7973l3Hs1EUE2spGNbNwr3B70CWQjG9n4byBbmGQjG9nIEGQLk2xkIxsZgmxhko1sZCNDkC1MspGNbGQIsoVJNrKRjQxBtjD5j0HTtA2apg140OeREjRNW6FpWr8HfR7ZyDhkC5OHEJqmndE0LVrTtLuapl3VNO1bTdP8H/R5KWiaJpqmRSac301N09Zqmvak4zYi0l5E5j2oc8xGxiNbmDy86CQi/gBqAagD4I0HfD7OqJ5wfuUBzAUwS9O08Zl9UE3T3DP7GNlwjWxh8pBDRC4CWAGgisPXxTVN26ppWoSmaas1TcurVmiatkTTtCuapoVrmrZJ07TKDus6aJp2KOF3FzVNG+mw7jFN0/Zrmhamado2TdOqpfL8bojIAgAvAXhd07Q8CftLNMc0TSutadq6BC3mhqZpizRNy+Vw7Fqapu1LOK8lmqb9oGnaOwnrmmmadkHTtNGapl0B8K2maYGapv2hadp1TdNuJ/xfxGF/GzRNeyfhOu5qmva7pml5Eo57R9O03ZqmlUjLc8hGtjB56KFpWlEAHQDsc/i6D4D+APIB8AQw0mHdCgBlE9b9BWCRw7qvAbwoIgGgcFqXcIxaAL4B8CKAPAC+BLBM0zSvNJzqbwDcAdR1dRkAJgMoBKAigKIAJiQc2xPAr6B2kxvAYgBdnX5fIGFdcQADwXH9bcLnYgCiAcxy+k0vAH0BFAZQGsD2hN/kBnAYQKZrUf85iEj28pAtAM4AuAsgDMBZAJ8B8ElYtwHAGw7bDgaw0mA/uQAIgJwJn8+BAiOH03afA5jk9N1RAE0N9isAyrj4/gqApxzOc4DB7x8HsC/h/yYALgLQHNZvAfBOwv/NAMQB8Da5XzUA3Hb4vAHAOIfPHwJY4fC5E4D9D/o5P2xLtmby8OJxEcklIsVFZLCIRDusu+LwfxQAfwDQNM2iadoUTdNOapp2BxRKAKDMoCdALeespmkbNU1rkPB9cQAjEkycME3TwkDtoVBqT1bTNA8AQQBuuViXT9O07xNMqzsAFjqcUyEAFyXhLU/AeaddXBeRGIf9+Wqa9qWmaWcT9rcJQC5N0ywOv7nq8H+0i89ZxqH9sCBbmDxa6AOgC4BWAHICKJHwvQYAIrJbRLqAJtBSAD8mrD8P4N0E4aUWXxFZnIZjdwFgBbDLxbrJoDZTTURyAHhanROAywAKa5qmOWxf1On3zqnvI0DHb72E/TVxvM5sZA6yhcmjhQAAsQBuAvAF8J5aoWmap6ZpT2mallNE4gHcAWBLWP0VgEGaptXTCD9N0zpqmhaQ0gE1TcutadpTAD4FMFVEbhqc110AYZqmFQbwmsO67QnnMUTTNHdN07rAtd/FeX/RCfvLjWz/x31BtjB5tDAf9LFcBHAIwA6n9X0BnEkwDQaBGgJEZA+AF0An5m0AJwA8m8KxDmiadjdh2wEAhonIWwbbvg2GuMMBLAfwi1ohInEAugF4HvQRPQ3gD1AoGmEGAB8ANxKucWUK55qNDICW1BTNRjayPjRN2wngCxH59kGfSzZ0ZGsm2cjy0DStqaZpBRLMnH4AqiFb28hyyBBhomnaN5qmXdM07R+H73JrmrZG07TjCX8DDX7bTtO0o5qmndA0bUxGnE82/nMoD+AAaAaNANBdRC4/2FPKhjMyxMzRNK0J6ECbLyJVEr6bBuCWiExJEBKBIjLa6XcWAMcAtAZwAcBuAL1F5NA9n1Q2spGN+4oM0UxEZBOS8we6AFCJXPNAIpIz6gI4ISKnEhxt3yf8LhvZyMZDhsxMisqvVFERuaxpWj4X2xRGUgLSBQD1XO1M07SBIFUafn5+wRUqVLi3szt7FrhxA6hYEThxAnB3B7y8gPBwwNMTsFoBGyOjxzw8EBkfj4qaBm9PT66/excoXx5wcwOOHOF3FSoAFgsgApw6BYSFAcWKAUFBSY9ts3H9nTtA/vxAkSLJz88R168DFy7w/6JFgbx5zbd3hgiv9coVIC4O8PEBChYEAl1anhmDyEjey/BwHIkqikj4ww93UcH7LODnx8XXl4uWTf/IMhAB4uOx9+DBGyISlPIPkvw2wyjeJQD84/A5zGn9bRe/6QFgjsPnvgA+SelYwcHBcs+4dk0kZ06Rjh1FfvpJBBCZMEEkIECkbl0RTRMpUEDE3V0uWiyS18dHqgESrWkiXbqIlCwpUqiQyOXLIqGhIu7uIq1aicTFcf+xsSKdOnG/s2YlP358vMjLL3N9ly4iERHm53vqlEjz5ty+bVuR8+fTfs1xcSLffCNStiz3U7GiyIIFPJdMROiSWxJS4aaEPv2tSPv2Innz8viAiKenSO3aIi+9xHM7eFDEas3U83mkERsrcuKEyNq1Il9/LfLWWyLPPCPStKlI8eIiFosISYB7JK0yIK0/MNxRcmFyFEDBhP8LAjjq4jcNAKxy+Pw6gNdTOlaGCBMRkcmTeQvWrxdp146C5N13+d1jj/Gvj4+Ip6f8UaiQAJAhmsbvR48W8fUVadiQD+jbb/l9//4idrv+4Lp04fczZ7o+h48/FnFzE6lRI2UBYbNRMPn6UhB++61+rLTAahVZvFikShWeW6lSIp98InLrVtr3lR7Y7SKnT4v8+KPIa69RSAYE6ALGz0+kcWOR4cNF5s4V2bBB5MyZTBd6/wlER4scPSqyapXI7NkiY8eKPPUUx2nhwpwk1X0G+LlIEZFGjbjduHEis2dnOWHyPoAxCf+PATDNxW/cAZwCUBLMbj0AoHJKx8owYRIVxRtcr57IsWMiXl4ivXuL1K8vkiePSPnyIoGBiTd9WKVKAkB+9fQU8fAQmTSJ6wYN4v7eeoufJ03SjxEbK/L44/x++nTX5/Hnn3yZChUS2bMn5fM+cYIvG0DN6uLF9F2/zSaydCk1MUDE21ukTx/OWjZb+vaZXthsIocPi8yfLzJ0KJ+Bt3fSgW+xiJQoIdKsmcizz1KTfBSEjd0ucveuyNmzIn/9RU34hx9EPv2Uk9qTT3IMFyiQ9H453rOmTUX69eMY/eYbkXXrRE6e5Ph0gfQIk4yK5iwGszfzgglT46HndhQDs1F7iMgtTdMKJZg2HRJ+2wFkLFoAfCMi76Z0vNq1a8uePXvu+bwBAF9/DQwYAPz0E7BvH/Duu8CCBUD//kDbtsCffwLFiwPnzyPWZkPDwoVx6uJFHPDxQdG8eYHHHwc++QT48kvghReAfv34+7lz+T8AxMcDvXoBv/wCfPQRMGxY8vP45x/gsceAa9eAhQuBbt3Mz9tuBz7+GHj9dfprxo8Hhgzh/+nBX3/xXixaRF9HqVK8B88+m7JPJ7MQH0/f1pkzrpdLl/jKKFgs9CmVKMGlWDH6hQICzBcvr/vnt4mPB27dAm7e5OL4v6vP6rtYA8Kvhwevs0QJjlPnv4UL0x+YRmiatldEaqfpNxkhTO43MlSYWK1A9er8u2MHUKUKH0DLlsCUKUCPHsCSJUCOHEB8PE74+KDm3buoEReH9RYL3Fu0oBN27Vou9esDHTsC69cDv/3G/wEOot69gZ9/Bj74ABgxIvm5XL1K4bRjB/D228Abb3DfZjh2DHjlFWDlSjqEp08H2rdP//2IjqbQ+/prXoObG4Xq888DnTqlX1hlBmJjgfPnUy9sjODubixo/Py4D6s15cVmM18fGUmnu9l55MmjL7lzm3/OkwfIl49CNIORLUzSi2XLgC5dgM8/B/z9gb59gdmzgalT9cEYGcmX3d0dC6tWRd99+/CWpuFtEWDMGODXXzmL7NnDh968OXDoELBuHQUMQIHy1FMUTtOmAa+9lvxcoqOBgQOpnbRty78pRW9EqEENGwYcPw506EANqHz5e7svJ08C335LLeviRUal+valYKlU6d72fT8QHw9ERKR9uXtX/z8ykgLV3d14sVhSt97X11w4+PtnmchWeoRJhvlM7ueSYT4TBbudDqj8+UXCw3X7c9ky2p39+vFvhQp0lgLSr3ZtcQNkva8vnVhz5tApWr067durV0VKlxbJnVvk0CH9WPHxtHEBOntdOVDtdpEvvmCko0gRkW3bUncdsbEiH3wgkiMHo0vDhoncvn3v98dqFVm+XKRbN+4XoD9pxAja3iqClY3/DPAgHbD3c8lwYSIisnUrb8fEiSI7dvD/MWPokFWOWU2jY9bPTyJ8faVcgQJSCJDrOXLw+/nzKWy6d6dAOHGCAqpoUZELF/RjxcfTcw7whTSKyOzdyxC0u7vIRx+lPnJz9arICy/wfPPmFfnyy4wLt169yghUmzYUdgCFV48edIZevZoxx8nGA0W2MLlXdO0q4u/PF+KZZ/iybNvGSEvLloz8FC/O2+blJX9VqyaeFos8Bojdy4sRiClTdKEkQu97QADDsI6hV5tNZMgQbvvss8aRiNu39WhQt24iYWGpv56//tKjPtWri6xcmb5QshEiIkR+/VXk+edFChaUxFBjvXq8/r/+ytjjZSPzYbOJnDyZLUzuGYcPU7MYMYLhVj8/vsjTp/NWjR3LvzVrJsbrZ7ZvLwBkRgLZR4YMEenbl///8gv3u3YtBVOjRgxHK9jtDG8q4lp0tOvzsttpvlgsNJ327Uv9NdntDCMWK8bjBAdTAGR06Ndmoyb19tt6qBkQKVRIQtt/ICFlrkro12fuf8g5G65htYocPy7y228i770n8vTTHNc+PvLASWv3c8k0YSJCQeDjQ+1EEdhWrhSpVo3mypNP0uwoVEgkIEDsFot0qltXPAH5K2dObj93Ll8oPz+RAwe43x9+oADq3Dm5j+Hjj/m7pk3pszHC5s08rpeXyGefpW3Wj40V+eorEtQAakrffZd5bNMrV0iq695dQty2CyASgi3U0po1I1ntxx9JXsvWXjIHNhvN661bRRYtEnnnHZrXNWok5/AULUpm9fDhInPmZAuTDMHRo9RORo6kplCiBG/+5s28XUOHigQF0RkLiPj7y/XSpaVwnjxSFpA7QUFkqK5bxxe/WDG+WCJkrwJ8oM4z9MKFFFK1apn7Ha5eFWndmvtp3Vrk3Lm0XV98PI9VsSL3UbYsadUG5KWMQOgqq4TUiJTQ11aKDB4sUqeO7m8B6Ndp107kzTfp9L58OdPO5T8Fq5VEtk2b6K+bOJEmZ8uWImXKJL3HjkKjXTtq319/Tf+giwksW5hkFJ56igLh6lW+eAD/PvccX/ipU/ldo0aJD2lD167ipmnyDMCoTrlyZGb6+NCXokwYpe28+GLyGXn5cm5frhwZnUaw2ch+9PWl8/Obb9I+u9tsIj//TNUWoNCbNSupGZaZiIkR2b2bGtZzz4lUrZoYKROAUawWLbhu4kS+LJs2UXg+Krk7cXHMyVq3jlre+PH0rzVrxklORdYclwIFON6efJLs2M8+I8P60CFGGVOJ9AiTbJ6JKxw5Qh7Fa68BkycDwcHMAN68GahWjYuXF7BtG5ArF3kJYWGY8NRTeHvRIsz38EBfu53clT59gO7dSVhbtIg8gtdfJyFuxAjg/feTcgu2biUT1tcX+P13oFYt4/M8dYos1U2byC2ZPZuEu7RABFixgszfbduAAgXI5O3fHyhZMj13L/2IjCQLefduMnJPndLJZ45wd0/KdHVcihcn69XfP2XC3/1EXFzaWa/XrpHprKBpQKFCSa/V8W+xYoC3d4acbjZpLSPx1FNksJ4+zQHeti3ZpT4+wKBBpLKPHs2XfetWIHduWL290aJAAfz111/YFxSEstevU1jExwNjxwITJwJvvskX+H//A2bNAiZMIBXeEQcPUqDcuEEB9Pjjxudpt3M/Y8ZQwH38MfD002knP4kAGzeSqLdqFT+3bEmCWteuGTZI04WYGODcOWOm62WDomt+filT6Z0XN7d7Y7nGxQG3bycXDHfvGl+fl5drIluBAkkFR9Gi942BnC1MMhKHDwOVKwOjRlGLaN2aQuXYMb5kt25xBn/zTQqaVasAiwXnO3VC9dBQlLx7F9tKloTXuXNAaCiZpPPnAz/8APTsSSHw/PNkl374ITB8eNLjX7lCIbJrF4//2mvmAuL4cWoTW7cCnTszV6hAgfRd+7lzPK9vv+XLGhhI4fr880CNGunbZ2bCUdicO0ctMjXs1ogIICrq3o7t4ZGU5erhwfuVWlp87txZsqZLNgM2o9GrFyMy168zm1eFhzdu5P9vvcUoT8GC9DnkySMCyNLXXhMAMiwhNCr58zNDs2FDetF37eL+rVaSvQASy5wRFaWzZZ97LmUnqdUq8uGHPEauXIwS3Qs71WYTWbOG98HLi+dRqxb9NRnBrM0KsFrJ3Tl/nn6Ff/+lE/7kSTo3L16k7+zmTToqIyP5HP7jIW5kO2AzGP/+y3DuuHH83Ls3HaQXL4r07Mn/ly3jNqpuSYECIoGBMiQ4WADIH15edJTWr88wXYkS3EZFYWJjWUZA0xhSdobNppc2aNaMgzolHDnCQk2ASKVKIqtX3/u9uHmTwql6de7X25skuq+/zo6+/JcQHS2yZ0+2MMkUPP4482siIzlbeXiQqn7mDF+o3r0ZLtY0Mmg1TcTTU6ItFqnu4yN5NU0uFCnCWz14sMg//zACU6WKzmaNimKYV9MYmXGFhQsZ6itbljNnSrDbSUgqXZrH7tyZJKV7hd1OLe3ll8kIVlGE2rVJwNu9+z8/a/8nEBVF8uPChdS2u3Th2EqIqGULk8yAMmm++IKfhw4lE/XkSV1jWLWK5ky1agxpJhSpOQKIr8UizQCxqhl97lyaDu7uFCDKDImKYr6LShp0hS1byHHJlYuCIjWIiSHF39+fgnDUKHNiXFpgt4vs308yVIMGehWvAgVolv38s8idOxlzrGykD1FRTGtYsEDk9dc5qZQpkzQM7+5O3tETT3BM//BDtjDJFNjtpKBXqMAZ9+JF+g+ef55x+8KFuX7xYknkjwBkmlos8o2mCQCZBIhUrkxtZu9evczjc8/pHJHoaBKKAJbcc4XTp3k8gCzS1PpELl1iSUmAPpyvv854vsa1a+SDPPkkuTaAiIeHhNZ6TUKKnJXQd3dwm2xkPCIjOa7mz2eCaqdO1EodyzS6u9Ps7d6dnJUffqCm7MIXly1MMguKuLZ8OT8PHcoHc+oUacoAX85WrfgSdetG7cXPT+yentLH318sgGz29uasXby4yI0bumbzzjv6saKjWXTZURtyRnQ0CzAr4pxjRnJK2LWL/htVRmDu3MwpIRAXR9Lea69JiM9fOp0e4PV3707y37p1Gacp/Zdht9MRvGsX0xDef5+mZseOnLgchYaHByeuHj1oev74I/1/aWA5ZwuTzEJsLDWQVq34+cIF+i8GDOBDrl+f6/ft44Ps3p3mSIkSIoCEA1La31+KaprcLF2av23ThtR2lRS4cKF+vJgYDhKADEYjLFrEaFNQEE2n1MJu5wCrVo3HKFFC5PPPedxMQGioSEjdeAn9cD9fgp49WVrBkblZvjyTzWbOpGl5/vyjw3QVodZ76RKz1BcvZrHzF1+kplqhQmICXpIlVy46xHv2ZILlkiWMSGXA5JDlhAnY1nG/w3IHwKtO2zQD2z6qbd5Kab/3XZiI6JXsVeLekCHUTk6f1v0qkycz8qPCxiqaYrHIbkA83NykKyD2pk25bswYCqpmzSiENmzQjxcTo7fKmDHD+LwOHeIxNI0DKi3OT7td5PffWTJAhbE/+ihNtOt7wvXrIitWkC7fubNexsBRLS9V6uGm1UdFUTAeOEAtbMkSapyTJnEyatWKjk9XeTR58tCkfeIJJuB9/DF9ZQcOpK0URTqQHmFy30hrCa1ALwKoJyJnHb5vBmCkiDyW2n3dF9KaM27dIgPxySeBb75hU6zSpVk0evZs1kfdvJns1caNSUQqUwZYvZrs0fh4fAhgZFQUPgUwuGVL1oz98UegVSugYUPSxjdvBqpW5THj4kjD/+UXkuPefts1uSkykqzchQtJrvv227TR6kV4Lu++C2zYwDKRw4cDgwcDOXNmwM1LAy5eZHFtV4WknZmu7u6kkCsKfVrZrl5e6We6RkSkTIu/dYtlOI2QP79xIejixZkS8ICQpUlrANoA2Ori+2YA/kjLvh6IZiLC0K6nJ/0d6rO7O8PEBw/SQz58uF7ucdw4+lAqVxYBxGaxSPu8ecVL0+SAjw+T7Hx9Rf7+m/soVIiLY5JffDydvQD9JEazsd1O4puPD485b176Uvu3bNF9NgEBDIPv2JE1ygQ49oT58ktGJ3r3FgkJoXZWrBhbk7hKgMvMxWIRyZePEZGGDall9e/PzPP33uO5/vQT+zMdOEAz2ah2TRYBsrhm8g2Av0RkltP3zQD8DLYGvQRqKf+a7euBaCYAcOAA6eSffMK2EufPU/t4/nngs8+A555jLs2xY1y/YQMrzI8ZA9SrB+zejWt2O6r7+yNXVBT2lC8Pv9u3qcXs3s1ZuXFj0uC3bNELSYtwH9OmUTOaP984R+PECban2LqV2tKXX7IVaFqxdy+vc8kSUs4rV+Z19u2b9vak9xsirFyfUuHo2FidDp/a4tBqG39/nRafI0eWo8OnCSLU/PbuZUH0PXugrV2bNTUTsMHWDbD/sPO6HAD8E/7vAOC4wT4GAtgDYE+xYsUyWA6nATVrklKu8Nxz1AZu3KAd7+1Np+qxY/SD9O/PmTMwkHU7AgIk1NNTNE2T5wGShTw86JC1WukP8PKiH8PZdzFtGmfCtm3N/RpWK30f3t487qJF6dcswsM5s6rqacrBvGLFw+GzyEZS2O1ME/j5Z5LV2rQhKdMxElSrVtZzwCYeBOgCYHUqtz0DIK/ZNg/MzBHRq6Lt38/PBw/y87vv8vOoUXSG7t9PNVfTRL7/ng+pWTNu6+MjYwsWFACyGNDNmNGjuY9ff6XJ1KFDcs/8119zXf36KVPrjxwhmQwgO1cVaUovDh4UefXVxBwkKVqUA3LHjmzWa1aE3U6TaulSFp5q356RP0cHd40adAR/8QXZywkRvawsTL4H0N9gXQHo2ct1we5/mtn+HqgwuXGDfpNXXtG/a92akYjYWBaNzpmTL29YGG3pkBA9ytO8uYjFInGANChSRAIsFjnp46Mn9C1ezH1++SU/9+uXXKv45ReeQ7ly1IDMYLUyHOvlxRnos8/uvY1mTAxDy23b6kzKoCCe648/ZnqkIRsucP06o4GffUb+SfPmHHuOfp1q1agpf/qpyM6dpn6bLClMAPgCuAkgp8N3gwAMSvh/CIB/wT7DOwCEpLTPBypMRKjm582rk4BWrOCtXLCAn8eP17WXOXP4/9y5DAGWKMGHHBgop93cJJefn9R1d5dY5bzz8dH7DU+cyN+++mpygbJ5M88hMJCOvZRw6JCuGVWtyjBlRuDGDZpRvXvrfZnd3TmYP/iA2lFWcN7+V3D1Kp/3rFkMADRtmlTbAJj7Vb8+TfCPPyZ3JTIyTYfJksIkM5YHLkyWL+et+/lnfrbb6cmvWZP/37rFB9qtGzWDWrWYs6N+98QTokhHPxUoIADkNYDaSbFi3PbyZe7r1Ve57ZtvJj+Pkyd5XA8P4wRBR9jtjCokkOmkWzeyeDMK8fEUcmPGMJExYXCHFnxaQoKOSejQpZwRM4kc95+B3U6TdN06kU8+YRSvSRNOHo5CI2dOmrEDBtBHtmoVOS0ZILyzhcn9Qnw8zZpOnfTvZs/m7VTEszfflESSmypGPX48+/F4eJDharGIaJq8WKmSAJCVAOnPvr40jWJiODAGDODvp05Nfi63b+sFpkePTp3vIiqKFH5fX5o/Y8eyB05G48wZkc8+k5DAQ5KETu/hQTLWoEH0Af39972bXg8T7Hb6u/76iybrRx/RbO7ShT6MXLmSC42QEIbpp09nSYkLFzJV48sWJvcTr77KF1FFVaKiOHN07crPN29SO3niCX7u3p3U94MHOTiaNKG5kz+/RAFSpXhxyefhIZe9vOjjAPQkQKuVZgRAe9cZcXF8MZWjNbUM1gsXSGEHKBw//jhTCkqHhoqENLBL6KIr1IxGjyarNUcO/YXx9WWe0YABFHQLF5LzcuHCw+fctduZ0LhrFxmvKo/msceosQUEJBUWAMdG5cqcZAYPJut5zRomlj4AMzE9wiS7bGN6sX490KIFG5arGq2jRrFO7IULZDe++SbwzjvA33+TbVmpEvDii+RsvPwy8OqrwIwZQP78+NdmQ527d9HQbseqYsXg1q0beSUzZgCvvMI6st27s8n6vHnAM88kPR8R8kKGDQMqVAB++gmoWDF117J9O899yxYgXz4Wun7pJbJEMxN2O8tN7t6dyG/AiRNsEO8IT08yXV0VjnbFbPXxyTjeR3x86opAO3+Oi0u6nxw5kp+/4/+5cz94rsqlSxwDW7dC+/jj7Bqw9w3x8UBQENCtG+n1AOvGVqoEfPABX8hbtzhYOnQAvv+e9PSvviLl/umnSVKrXZuUe7sdX9WogYF79mCKpmF09+48xrJlwPLlQLt2rHX62GMUZN99RwKbM9auZUX8yEjS/Pv0Sf01bdpE4bdmDV/UV14Bhg7lQL+fiI42LyB95Yr5793ckgsYX18Kr9TS5a1WktrMCkF7eLiu6ZonD9MZHIVGrlwZcmsyDHY7cOhQovDAli28twDg4wMtOjpbmNxX9OnDYtGXL5MVCQD163MAHjzImWbUKBaMPnWK2kmZMhQuI0dy2wEDmJ8TGAg5cwZP1quHX3bvxha7HfU/+ABYsIAV8rdvp6CKjATat+cAWLiQuTvOuHQJ6NWLeT6DBlFbSkt1+V27mKezbBmZni+/TMGSHiZtZiA6muzj8HAyWe/cSZntGhWVlN2aGqarhweFg1FBaH//B69NpBZ37rB9yLZtHDvbtrHwNkAtumFDoFEj/q1ZE5qnZ9ZkwGb0kiV8JiJsrwkw9KbwxRf8ThWNPneOjtYRI/hZZRPv3MkUc4uF2b6ASMmScjtHDilepIiU8PGR2+7uzBItUIAp+6qwUEQEQ4JubgzLukJ8PAl0AKNJJ0+m/foOHGCESdN4np06kQCVGfVPspFxiIig03/6dDaUq1Ahab2TSpXozJ07V+TECZc+GWQ7YO8zbt0ip+L11/XvwsJIY3/pJf27Xr3obLxzh0tQEHkYN2+SSNa8OSMyvr4i3t6yvX59sVgs0sPXV+xFitB77+1NB6UKq969S96Im5vOb3GF335jdCBnTjo10+PMO36cTtOEcpRSoAA/p6YWbTYyF5GR7CU8cyYjhaochRIchQsz8XDiRFITVJJqCsgWJg8CzZvTQ++Ip57iC6wiIzt38laruiQzZkhiGFn1H/7iCwqTBH7G5G7dBIDMtlhYIEeVhXzmGV0g3L3L42sas4SNcPo0Q4sAC2Snl1YfF0fNpHNnaiqASOPGLEF561b69pmN1CMqSmT7do6ZZ5/lWHGs5VqgACNGEyawTs2lS+k+VLYweRBQvYMd+8isWSNJSG0ifJlLlmSYNyqKD755c5ojlSuzCNB77/F3VaqIzctLWmma+Li7yz8AZxZlDk2apO83MpJhVk1z3XtHwZFWnycP84XuBZcusVB12bKSSNdu0oRcmH//zWa93gusVmqDv/3GMfH006TCKwGu0hc6dCCf6bffGELOQGQLkwcBJTgcyybGxVEzefZZ/buffkoqYKZP5+eNG/V9vPMOSUtBQSL+/nLZy0vyublJ5Rw5JAogw9FVmcfISA4sgIPP7EU+dEjPAO7R494LPNvtnC3HjtV76gBk2Q4ZwlSDLF6744FBCY2lSzkpPfUUn7+3t34fATKi27Vjftcvv9APl8nCOluYPAjcvs3bqLKGFXr3plBQafpWqz4oRKid5M9PrUKE7Ec/P758bm6JeTQrAQEgg1QJgxMnuM7TM2mZx7g4kT59eC7Dh5sTveLjWWLS05P7/OqrjCsncO4cTbZOnRLrloZ6dZCQHAcltMcXrIh+6tSjpblERlKI//orx0mfPq6FRtGiHB/Dh5MZvH37/Su2feMGJ7zBg0UqVMgWJg8M5cvTj+AI5ePYulX/btw4CgpVTf6jj7jN5s0UEp6eLEeg8nEaNBBxd5fXEgTKT97erHNy+TJzcgIDRQ4f1vdvs1EbUNnGKVHUDx5kciHAwb1xY4bcjkRERYn8+aeEFDhJOr22TX9x8uThi/PGG1TTz59/+JiuCnfu8F7+/jv9GSNHkvFcp07yJDyA+Vft2zPC9803LOFwvyv0h4fzfIcNS6pR+vmJtGuXzYB9YHjmGRLPLl/WeQdhYSS1jRwJTJ7M706cAMqWBd57D3j9dXIfSpUCqldn4/Phw4GZM4EdO4AnniDR6vJlxNlsaBQTg+NeXtgfFYXigwezkXm9eoCfHzko+fPzGCLApEnA+PFsYL54MfdjBBE2U3/tNTJ3e/Yk87Z48Qy7PWvXAm+9BUx8Mx4tg/7WGa+7dwP//ksCGUCmqyMz1HkpUICEtPuB2FhjluvVqyR4qTq1t24l/a2XV3KGa4kS5BhVrJj5zGJXiIjgONmwAVi3jvffZuO5NmwING9ORnedOoCHR7pqwGYLk4zArFlkip47x6LTCi1bcuD984/+XZMmZHAePUrBM2UKBcvff5MUVro00KwZSy9260aW6w8/4JS7O2oCqJwjBzbeugWPBQuA8uW5baVKHCCOg/TTT3lOtWqRfFaokPk1REVRiEybRgHz2msUbpnN3IyKAvbt4/U7M12vXUu6rYcHzyc1xaI1zZzd6siAvXMnOSU+MtL4nH19XVPi1d98+e6f0DPChQs6u3XrVpYctdtJxqtXTxceDRq4JDRmC5MHhZ07yWZ1zNMBmFczbBhfDDXTf/sta8Vu2cIZwbnq/eTJwNixpLZPmQJs3MiHv3kzvo+PR28AY7288K6bGzWY8+eBLl04MP74I2lt2N9/J0M2Vy7+X7Nmytdy7hwwejTp/zlysJbtq69Sy7rfiIpKWqX+7Fmd9Wq03L1LYegMxXh1Zr5aLBQ+zuxWs8++vlmL+WqzccJSgmPLFj5HgJpr/focaw0bAiEhqap6n6Wr02fkkuV8Jlev0t78+OOk3//1F79X1dNEyE7082N2rIKqen/lCp11hQrRX3L8OEO5nTqR9JY7tzyvaaIBEurjw/aPt2/rrUZ7907ud9i/n449P7/U9ydW5969O0POvr60rTM4/JgpsNl4jyMiGEWKj//vOXvv3mWtk0mTWO3OMfu6UCFG6WbOZJGtdJZ2QLYD9gHBZmONDlXDVSEujhGNYcOSft+nD6MoKoJy9CgfxVtv8fNXX/Hzr7/q3JIRI0QAuevrKxXc3KQAIFctFpKUbDZyPgDXVdkuXaIzUNO4XVperkOHGI62WCjwBg1i/ZFs3B/YbCzNuXgxn22dOnorD00jcW3QILKgMzBKli1MHiSKFye5yBkNG3JxhIr0bNmif9epEwVMVBRnk7JlGWGJimLX+nLl2P3N21sOAOIFSHt3d7EBFDh2OwvsuApTi1Dj6dmT67t0SUqySw1OnhQZOFDvPFe3Lkly2X2CMw52O+/zDz+wKX2LFnoDeICh5KZNyelZvjzjWccXLjC69OST2cLkgaJhQ/I/nDFsGLUTx+S427c5uzhqMmvX8nGoPJt58/h56VKG8ACyHf38RAoXlk8TwsUf5sjBGWr5cs5iTz0lLk0uEQ7W6dN57JIl9VqzacG1awxpJzQWE19fhqE3bfrvmROZCbudaQ5LlrDMZatWeg1dgEK7Th1qHXPm0FzN6ATL6GjmfY0YkaTMprCUaNYTJmDrioNgH+FkJwhAA/AxgBMA/gZQK6V9Zklh0rMntQlnfP89b/NffyX9vkULJmUp2Gx8wRWJLT6eGknNmlzXti1nqQTKvT1vXunq6SkegOwOCuK648c54Lp04TGN6sJu20YCnacnG5anRwjY7eRHvPCCXjmsaFEO/j/+yJSKbQ8tIiNF9u5ln+TRo/ksVbsQgMK9Vi1qfl9+yW1VsfKMhN1OXtL06eT4qGbonp4iLVuyL9OBAyJ2e5YWJoZ9cMDGWysShEp9ADtT2meWFCbDhnGWdsbJk7zNs2cn/V7R6U+c0L9T/pHTp/lZOVZ/+42+C4uFAy4kRMTfX24CUtRikdKAhAcEUFu4c4eZxa1bkyD344+uz/f6dQ5qgIWlr15N/7Xfvct09scfp+akVPKOHdl64ezZ9O/7YUJkJLW9efMoNB57jDlXjlm8Hh7Ms3n+eQryXbsyN93g9m0yW194gWQ5dR7lyokMHUrB76LM58MqTL4E0Nvh81EABc32mSWFyRtv8OV1htVKITB2bNLvT5zg7f/kE/27M2c48CZM4Of4eEZsgoN1n4ibG/MzPDxEKlaUzYC4AdLH21vsmsaas3Y7B0ijRpz1jKI4NhsT8xStfsmSe78PMTHMIfrf//giARKKFhLiuVtCG0+gA3jt2oe3t05MDDXANWtofowaRaFRsmRyoVGlCjXWt9/mvT10KPNrwcTG0uScMIGTjkoOzJGD9YG/+EKfrEyQVYXJaQB/AdgLYKCL9X8AaOTweS2A2i62yxrtQY0wZgxfSlcoUoR+BUfY7cwc7ts36fetWjFJToV4VTOuDRvocMudmy0dVZGlEiVkoo+PAJBvixSRJA7YsDA6Sj08GBkywj//iNSuzd8++SS1loxAglodUuIi6fRee/SXDWAawlNP0XT77juaX5cuPVhafWQkTYGVK/niqeboDRow7OooMJSJULUq79vEidQCDh++fwWk4uNZ4mLyZGqjynTRND7TceOYrpHG80mPMMl00pqmaYVE5JKmafkArAEwVEQ2OaxfDmCyiGxJ+LwWwCgR2Wu0zyxHWgPIGP3sM9fMyfr1SYxasybp948/Dhw5wkVh0SLWh1WktuhoktoaNQKWLmUJxuHDWRd2+HAgIgK2a9fQyssLuyIj8VeVKij/778kqXXsSJJX27ZsSv3jj0DXrq7P32ol+3XCBNZ/nTqVaQIZwORMpNNPBFrWuKkXj969m8ulS0l/4EhHV8WW/f2N2a5eXikzXWNikjJcjajy0dFJz8XdncWsHc/H8W+RItzmfsFuZ0nQdeu4bNpEBi8AVKlC8mKLFmRaBwamff82G7BzJ7SGDbM2aQ3ABAAjnb77b5g5r75KVdIVunVjYp4zXNVCCQtLHul54w3ONMeOUc0uXpwOu3Xr+PtGjeQCIHksFqmuaRJdvjzP5dAh/j48nB3e3N05c5rhwAFuC3BmcyxJmVm4e5c1UJYvZyuP114j8apOHbYD8fJKqg3c6+Luzv1WrEhTsHNnts0cOZIz/MKFnM3Pn3/wzdmV0/TTT2nCOjpuy5Zl6c8ffrg3n1dUFE3h555LbCmKrGbmAPADEODw/zYA7Zy26YikDthdKe03SwqTIUNogrjC0KGuBY2qY7J6ddLvW7VKKnwuX6Y6/fLL/LxgAX/33Xd05Lm5iVSpIn/4+QkAGRoQwGzVMmV0LkJ4OFV1i8W8KpsIzYyFC6nWAzRFzp9P3X3ILMTGsszlmTMkzW3dSlNkyRLWwf3xR/qSli0T+fNP3tN16+g/2LaNjtFTp3gfsnII22pl5O+TT1jus2BBXXgULcoaOfPn3/vzuHaN0b4uXXTTKEcOHnPx4iwpTEqBPYQPgP2ExyV879hrWAPwKYCTYAg5mb/EecmSwmTQIEp1V3jnHd5q57aYqhbKO+8k/X7mTH5//Lj+Xb9+jBbducOXvXp1OjgvXaIQq11bxNtbXi1YUADI0urV6Stp3VqnVN+5wxAgwGhSSoiIoFbk5cXozJAhj05k5n4hIoJdyiZM4LNybNBVuDBf7tmzDQs/pwnHjrHaXqNGernHIkU4Sa1enSQcneWESWYtWVKYPPEEw22uoISJKydYwYJULx1x6hS3/+gj/butW/ndnDn8rIhs33yjtybt00diAKmZM6fkBuR89+78/pVX9P3ExNDsAigoUjNAT5+mBuTuTgH1/PNJBV02Uo8LF2iW/O9/NFVVtEXTGDIePJia1pkz9y484uLY5HzkSFaoV0KqenU68PfuNTxGtjB5kKhSJWnvYUdMmMBb7SpKUbs2+R7OqFiRBXQU7HZGPxo10j8HBzN0HBvLoklBQSL16skxf3/xt1iksaZJvKq+9vnn+r6sVr1/8YABqff0nz1L7cTLizNbz5401R7WokaZjTt3GIX78ENGhFTDeICmRbNmFOgrVmRcqPzaNZqxPXvqVHyloc6cmaqwsEi2MHlwsFr5gr32muv1yoHqCl26JK9uL8LZP0+epDPH5Ml8ZMeO8fNvv/Hz3Lm0szWNNrWPj8yvXFkAyPg8eSisLJakdWrtdoYNAZo+acnzuHyZ16ro3yVKkEtx7lzq9/Ffw927eq+ap59O3qumSBFqhB99xFBuRoWO7XaRffuYQVy/vn7MAgU4hn75hUItNbh9m2OpQ4dsYfLAoFiuygRxxpgxnB1cYfBgvpTOUPwSR4bsxYvUCBQBzm4n3b5MGQq0gQNpiowZIwLIMyVKiBsgG1q0oMDKmTNpmUcRzmIeHtR60mq6REfTCaz8MJrGGfCTT2iq/RdhtzNysm6d3qumcuWkLScKFWKE6O23GaFKb2sRI9y9y4lk4ED6VdRx69ThMffsSb22GB5Oh+5jj3EcACLFimULkweGP//krXTMAnbEyJFUa11BhYcjI5N+v2+fJEZsHNGqFWc9BVX1/pdfqOLmzMm8i4YNJcLLS8oCUhiQGxMm0EFcokTywb1pE7Wg3LkZIUkPTp1iIqJqfQEw9+i116jqP2wV6u123qe1aykcBw1iKw/H0CzAouAdO4qMH89IUmbUfFHh4Rkz+GxVqDwggL66b76htpha3LpFDeSxx/Qs8CJFWMh6x46sm5uTGUuWEyaqMLQRc/Tpp/mwXEFFbm7eTPp9fDwFkKPz1HF7ZepYrRQQypfy4Ydc/+WXIp6estffXzwA6eTmJvbvvuM+69ZNLrxOnKD2oml0zt0Lv+LYMar7LVuKeHiQTq9tk9AyL5IXMWcO+SzpLNyTYYiPpw9h/Xq+XBMmkG/SqBEFq6PQyJWLmeEvvMBrW72agiOzwsxmOTXDhjEClJZkwBs3WPG+fXtdAylalPvaujWZJpMtTB4UunY1FhYiVD9btXK9Tr38ruqC1KzJfjiOUJGeDz/Uv1NJg7t2cYCVKUPVOyGKNMNiEQDycf78zGLWNCblOQuMyEj6XACe770QoRTCwyWk/A0BREJyHExaFczLi+ZV27YUMpMns9bLtm2ciS9c4H1Ji4M3Npaz9D//UOP69VcKr6lTSQTs25caRrFiSU0TZaYVKkRhMnAgBfeaNQy/ZzY3xWqlL2XiRAqtdObUJMHVq5xUWrfW91eyJLXFnTtdX9PlyyIzZmQLkweCqCjyPxx7CzvCbueAUIQzZ0ydKi7NHBGqodWrJ/++SpWktVPCw6ny9unDzz/+yH3Oni1Ss6bYAwLkMUA8AfmrY0e9PelLL7keUF9/TV5Jvnw0n+4RoaHMOQsNFQqGo0dJihs5kmp6cHBy88F58fOjU7FUKWpiRYrwc9681BoCApL3oXFePD05GzduTKHy5psUNKGh9Bc584AyGw7FiBI1Icecmk2b0u6ovXyZmdrNm+vCskwZ5hgZhYLDwpih3qpV4m/SI0zuY1LBfxRr17LwcZcurtdfucLciQoVXK+3WvnXVX5H4cIsGu2M9u1ZrDo2lnkpOXIAffsCX3/NqvTduwN16zIZ5ocfoDVqhG+LFkX18+fRa/ly7O3QAf6jRjEXp2BB4M03k+7/uef4+2eeYYX8Pn2ATz5hjkw60LIlF8INKFeOy1NPJd3w7l0WjT57lq1CXBWMjopKXhDa8f8cOZIWgHYsCu3n92ALQd+4obeaWLeOHQoAtvDo1Ik5VK1bA3nzpm2/Fy6wmPlPPwGbN1N0VqgAjBvHsVC1avLrjo0F/vwT+O475nHFxgKlSkFefx37a9bk79KKtEqfrLBkKc1EFQcymtXWrxeXlHkFxUFx5aOYNInrnJ2XP//M73fu1L/bk5CRq/gkGzbw8wcf0LEGyHo/P9EA6efmxspdzzwjif4VV4iLY3TA3Z1awIIFWZuKntUQFkan7KuvkpCmNCR/f5qv77+fWIwozTh2jOUcVKtXgBrrhAk08VzBamUU6vnndQ5KvnwiQ4fKhWXLZNrUqVKlShVBOjWTBy4Y0rNkGWFis/El69HDeJsPPuBtVl38nPHCC1TxXUExW535G+fP83vH0ox2O1Ph69bVv2vVioPl6lUmBxYpIm9yoMjCAgXo9O3Qgaq1YwV9Z+zbp5coaNCAvplsJMedO6zlMno0/WTKzPD2pjP63XfZ8jM9HBPFJ3nzzaQlFuvUYQmHI0eMf7d3L0szqlwrf3+RZ56Ru0uXyoK5c6V169aiaZoAkAY+PvJZtjB5ANiW0O7SsYm4M5o04Ututt654LSCqlLvLIjsdtLwnQtYq6iSivRs3szPH33EilqAxNepI40A8QfkeIcO5Cw0aULtY9ky4/O02WhX58/PfT71VHLOyqOG8+fp0B4yhM5yJTzc3enEfestaqbpDYvbbIy0jBhBxynAYzRtSuewWZ7UiRN05ioavYeHSOfOYl28WEL/+EOeeeYZ8UtIDC3h4yNvAnJMTRZDh2YLk/uO3r0p5Y2o0Ddu8OG/8YbxPvLnT56bo/DZZ3xErkhP7doxt8MRZ85IskhPy5Y8RnQ0Izg+PnKucGEJ1DQJBiR2xgw6cOvUYXRl1Srzaw4P58zr60uNpnt3zpj/dVitNEk+/ZSObsdwrZ8f7/Nbb/H+uSiDmGrExtIkHjSIWq9yHHfoQGfxtWvGv710ic51R9OnaVORL7+Uf7ZuldGjR0vhwoUFgOTw8JABHh6yCRBb0aIkT770Eh3b2ZrJfcbx4xQUo0YZb6MqzO/e7Xq9yhqeNs31+o8/5vobN5KvGzCAQsIZVasmjfSomidz5lDY+PiING0qvyaYO8Pd3Oh7uXmTkSMvr9QR165dIxNXhXqbNaNP5b9QSNpu5/NdvJgRp2bNkmbzFixI03bGjHtqdJWIsDAeq1cv/X76+fEYixebtxO5fZvRt5Ytdc2oRg2RadPk6t69MmPGDKlVq5YAEIubm3T09ZUfAIny9aVQHD2aAkdpPR06iCxZki1M7iteeIEvnhnzsFs32qlGDrYtW/gIjGq0Kg6KK81n/HiucyYujRtHToHKtbHbObgqVeL/inHbu7e8nCBQlgcFUWDduMFtvbzI6k0Nbt+mzZ5Q71Vy5uQMt27d/StdeC+wWmkW/vgjJ4aWLRlqduTC1K3La8rIRlfnzpFZ26qV3lQrXz46R5ctMxfKUVE838cf1xmspUuLvPmmRO/bJz/88IM89thjYkngF9Xy95cZgFxVQv+NN1jSQgnIMmX4DM+epTncuXO2MLlvuHCBNujgwcbbhIdzdhk0yHgb9eI7s18VXnmF+3A1eJU/5cyZpN+rKM7y5fp3SkNauZLmTokSIpUrS3TFilLNYpG8gFxs3pw2+s2bNJ88PJgqn1rYbBQgTz2l8z1y5BDp3l1CX1spITWjJHTNA4wEWa3kt/z6K8l8ffrompgjD6V2bT6zr76i+ZbRCXkTJtC/oo5ZvjyF2Nat5qzj+HhmF/ftS9MaoBn06qti37lTNm3cKC+88ILkzJlTAEhhX18Z7eEh/yhhMWoUuSbKh+LrS4Lipk0UIhMmkIOTINSyhcn9wrBhFAJmjETVrtOs0VWVKlQxjeDKL6KgugKq0owKERFUV8eP17+LjWV5gm7d+Fnl84wZI4ctFvF1c5PmgFhVVfzbt+lA1DQyL9OKiAi+tAMGiBQsKCHYIoBIiGUHi2GPG0dBtXMnI00ZMdNHR3O237KFDvF336X22KYNKejOpR+LF6dKP3IkHct79mR8r5qoKGp4gwfzeIqUFhJCsqJRBEbBZqMTffBgPj+l+T3/vMjatXL8yBF56623pGTJkgJA/Dw85Bk/P1kDiDUggNtNmcKkQ8WAbdiQJu+tW9SIO3bUzaM2bTg2YmOzhcl9waVLlOrOVeUdERVFldVVnRIFlWnsWADJGSVK0MnrCqp0o6tM36pVKYgc8eqr1DauX+fL27gxfS6vvy5fJ5g77yjtRYSM3I4deYyxY9Nfs8Rul9AvT0hIqcsS2v598i3UwFaLjw9nzAYNOKCfeIKz5tChXF56iYKhf3/e9x49aI7UqMHZ1Nc36f7UEhREx3L37qSQf/stw9qpTclPD86fpwB+7DG9HKKvL1/oOXNSziBWodyRI3VNwdub1/zrr3Lz0iX5/PPPpUGDBgJANE2T1rlyyXxAIpRA+PBDRoCUAzd/fmomR45Qk33zTT1MXKAAn6/K8r5wQWT8+Gxhkumw21l/xNtbD7+6wqxZvLUbNhhvo8K4J0+6Xh8VxVns7bddr1cNulxpR/36cbA44u+/uf3Mmfyswtpvvin2atWkt5eXWADZEhCg7zMuji8xwBfSFeU/PYiM5PksW0a/wYgR3H+rVizyVKkSX6RcuVieISiIg75IEQrY8uUpeB57jNc6fDhNly++oClw6NC9RVTSAquVJsrYsUmJaSVLMmS8YkXqQsOHDzMaVK6cJIaXO3YUWbhQYm/ckKVLl0q3bt3E09NTAEjlnDllqru7nAcoiN9+m0KkcWP+3mJhsa6lS3m/f/2VSX6axqV9e34XF8dxvW4dhbjFIkLOSdYSJgCKAlgP4DBYA/YVF9s0AxAOtg/dD+CtlPb7wISJavVpFH0R4cMpVoyqrJH6ruqQmPFPVqyQZL4PR3zxBde7KiysGoI5RxmqViWnRKFbNzrhNmyQcItFSvn5STFNk1vVq+svgN1O4p2m0ceQXa6R2cLOOTUWC03WadNYaT81ptuZMzRDatSQRBOoeXOR2bPFfuOG7Ny5U4YMGSJ58uQRAJLP11de9fOTvwCx58yp9yEeMEB3ppYty31eukRfyJtv6kWpCxXiZ+VnCw+nMK9Yketz56YGd+JElhQmBZHQOxhAAIBjACo5bdMMwB9p2e8DESbXrjGprE4d81DgN9/wtv7xh/E2K1dym6+/Nt5m0CA6X41mNVVwyZWD8PPPuX9nspty+KpQ88GD3O6NN0TGjZNdgLi7uUk3QOz9+iV9IZYvp5aQI0fK7TL+azBp8C39+tF/ldpKdZcuUTsMCdH3U68ew8yXLsmZM2fk3XfflfLlywsA8XJ3lyfz5JHlgMRpms43mTrVdfP4+Hg+q06dOKGo3yxdqo/bv//WxxfAaNXcuUkiSFlOmCQ7GPAbgNZO3z0cwqRXL768Bw8ab6METr165jNTkyZU2Y0cfjYbZxHlMHWFrl2TFklyxK+/ikvn744d/N6Rsdu9OwXElSsiFSvK+7lyCQD5DKC55ogzZ3RCVN++aSv1+DDBZmMZzA8+SLHBd6pw7RoFfLNmelnFatXoJD55UsLDw+Wbb76RZs2aSUJejDQJCpKv3N3lNkCzb8oUCq2ePfV6JPXqMeUiPJwUhXff1R29+fPT9HI0WX/4QTeDvL3pg3LkQMXHU/uuWzdrCxMAJQCcA5DD6ftmAG6C7TBWAKhs8PsH1x5UvZwTJ5pv9+STfNBGiVYiOsVd+S5cYdcubmPW36ZiRfIMXGHpUv7+r7+Sfm+1MqzoWA5h/35JNN02bxYbIG2LFxcvNzf522JJ7veJjWWkyN2d6vPixQ9/8p/dTj/LrFkU4I6FkcqXN23wbYhbt6iltmmjO5wrVOC9O3RI4uPjZcWKFdK7d2/xSWjvWjYwUCYGBMgpgFrgyy/zuG+/rQuJ3LnpTD94kEIvNJQTguKqtGzJXkJqorp2jf4kVd6xVCkKSUciZFgYv1Os3rJls64wAeAP9hru5mJdDgD+Cf93AHA8pf3dV83kwgWqs9Wrm3MOfvmFt3PSJPP9tW9P7cXMmdm/P2dDI/7JrVscoEY0fRX6/fvv5OuaNSOXwhEtWlBTiosTGThQrmqaFMidWyp6esrdPHlcO3n37tX5Eo0aJc1gzupwbnSloh4qZNy/P6NlRsmZRggP5+8c66mWLEl+x/79Ina77N+/X4YPHy4FChQQABLo6ysv5c8v2wGxu7nR6bp4MbWIDh30sG2rVtQaYmIoCD74QC+RmTs3TbCjR/Vz2buXETEVEm/ThoLJMSp36hQFk+KtNG/OFio2W9YUJgA8AKwCMDyV258BkNdsm/smTO7eJc/D39/1i6lw8yYHZI0a5gJHlQ6YOtV4m7NnOcsMHWq8zaJF3M/27a7XL1zI9a4S8UaO5EB3HFTLlnH7n3+moAoKkjUVK4qmaTLAw4MquatwqtVKcpfiQLRpQ00mq2kq4eGcwd9+23Wjq9696YdITxHsO3f4PLp00V9cVU911y4Ru10uXrwo77//vlSrVk0AiIe7uzxepIj84uEhMcqMef99cmRGj9aTKQsVop9LRfx27aJvRB2nYUMKL+VXU6ZMw4Zc7+dH7cZxHNjtjD498QQFlbs7TVYnLTbLCROwW998ADNMtikAJDZQr5tgCmlm+70vwsRmo8rr5mbuTBXhA7ZYkpsVjrh5k4MkJYEzdCgfsFlGaM+e3JcR9+Ptt2mbu9J+VNjaMQ0gPp4mi+r7k8CYfb1jRwEg32sa1xkxNO/coYBM6FMrFSty5rx0yfgaMgt37ohs3MgwaZ8+eqhVRUsyotFVRAS1h65d9Re7cGHO8tu2idhscvfuXVm4cKG0bdtW3NzcBIDUK1RIPs2VS24ADHsPHkwBsmgRNUYVFercmRpCfDydot98o5eA8Pcn78ZxcnNlykyfnjQNIz6egqZePW4TGEiNyVn7CgsTmTw5SwqTRgkOpb8dQr8dkLQ96JCEsPEBADsAhKS03/siTF5/nbcnpTaaKkRrlhkswkJE7u7mGbZXr9K86d/feJvISM6sAwYYb9OnD9V1VzByzo4apUd67HaRkBCJCwqS+l5eksPNjXa8c3FrV+f21Vd643NApHZtCe07V0Iq3ZbQ3yLMf59a2O0k3+3cSb7Na6/RJHBscqU0hMcf54u2cuW9NboKD2engG7ddIdswYLszLdli4jNJjabTdatWyfPPvus+Pv7CwApniePvFGkiBwBODG1b8+8mn37yKRW5SpLlqQDVVW3P3GCWqTy31SqxInAMelv377kpszvvycV+mFhFKzKH1KmDDOfnf0/Fy/yPiZobVlOmGTWkunCZO5c3poXXzSfudas4QvYoYN5XsXy5akTOH37UuCY0axVdfrNm423qVLFmH27di1/7+xYVU7fBQv4efduEUBON28uOQGp5+8vcSk5jh1x6BBfjgYNdDo9ttDO79CBs/LUqYxGLF5M7W/1atLPly2jyfXNN9RwXn+dBZ47dGA4VIU01eLpSY2jd2/6rDKqV83t2+wp07mz/sIWLEgy2qZNiZrhoUOH5PXXX5eiRYsKAAnw9ZXny5SRDV5eYlOO1ylTKCAWLNAjKu7udJ6uXs19Wa0UBu3aSaKW0qMHa6KocWizcZvmzSXRlBk8OHlaxaVLNJlUFnKTJnTMO4/To0c5MXl6Utj16iXy11/ZwiRDsGkTfQotW5qbI4cPM0+iShXzFPErV6h+VqpkXrA4oXiRvPmm8TaxsZxtGzc23kbVNHn/fdfrV63i+q1bk36vqsb16qV/l5C092ONGgJARhcrRlPhxx+Nj+8Cob/eoWbSfyFn9po19W6AqVnc3WlC1axJTePVV8nLWLqUL0NGtsy4fp38nw4ddCdqkSI8ZoIGIiJy7do1+fjjj6V27doCQCwWi7QvU0YWBwVJJMCXeOBA+rX++Ydanbrm0qUpSJXAu3WLz0tpVgULMurj2IMnMpLhZWW2FSnCCNzt20nP/9gxspaVcOjZ03UJjN276TfRNIaJX3opScO3bGFyr9i3j2pl+fLmHIobNzgg8uUzT/aLjqbK7+ND77oRwsIocCpXNhc4X3/NR7ZihfE2qgaKEd3fyMwR4SxYsqT++dQpvlD9+8tAX18BIKsqVuRANappmxbcucPkvH//JQdm0ya+fLt381mcOkVBndkO3fPnGdlxrOhevDgjJNu3JwqQ6OhoWbJkiXTq1Enc3d0FgNQoXlw+KltWLivB17IlHeC3btEXorQQDw9SB9au1X1dBw9S4CizqUkTCmrHSezSJTphlblTuzbNLeeJbtcuXTh4eVGrdsVW3rZN13xy5SIXxYUWly1M7gW7dvHmFi2atCWnM2Jj+dC9vPhgjGC3U+0GUmaMvvACB7FZeFUJnOBgc5p+/fp0gBrh/fd5Tq6E5bRpXOcYkh44UMTTUyK//loqA5LPx0cuV6zIF8As9yir48gRagfKIakcx+PGJWkJYbfbZcuWLfLiiy9KrgRCX8G8eeW16tXlbxVSVb2Wz5yhEB85kuF/pYVMm6ZXSLNaKdBbtOB6b29m9+7fn/T89u+nY9/DgwKia1cKW8dnb7dT01T7ypmTJqGrGjubNzO8DPDcpkxxHaE7flykf/9sYZJubNlCx1OpUuaaRlwcZ28gedtOZ7z9Nrd77z3z7VStEaOm5wpK4BhVbRPRa5k4M1cd8dxz1KhcQfFTHJ3Ep0/zuKNHyz9PPinegLSuXl1s5cuTxr1unfl5ZxVYrXzOo0ZR83RwEMt77yULo584cUImTJggpUuXFgDi6+srTwcHy6rixcWqhMDTT/P6Y2J479TLarHQnFO+EJHkpkzRonyhHcljdjv9cGo/fn6M7jlPboqpqng+BQty387mtt1Of4vyr+TLRx+UK/LdkSP02bm5iXh7ZwuTdGHdOj60cuXMSUoxMeQSAElrrLqCqjXyzDPmKrryz7RoYe6fWbOG+zMrESlCb36+fMZVuux2DuLOnV2vV5XfnM2o7t2ptV2/Ll/kzSsAZOqoUfQDeXnRd5EVER7Ol7x/f50L4+FBrsmsWcmq/t+6dUu+/PJLadiwoQBM729ZvbrMCw6WO4rFWq8eW4OEhSWm6ycm0hUtSpa0o6/j6FE6SFWZhCZNeE6Ofh6rlWHbWrV04TBlSnLtMS6OTunSpbld+fLkxzibxkooKROrYEH6mFxRBU6fpgbk5sZzHDFC5PLlbGGSZqxYwRmmcmXz8otRUbqdaTbrizCS4OVFVqiZ/+PECYYFy5VL2T9TvDi3Myvlt3o1z2/KFONt/vqL2xglGKqyBM41YJUwW7xY7Js2SXcwIXDHqlXM1XFzY0mFB01WUw2+P/oosc9xom+gd2/O5k7h4bi4OFm2bJl0795dvLy8BIBULF1aJjdrJucUeSwoiCS0f/7R0/W7d09M15f27RmBUpESpRF07sz1np4M4TqbMtHRpBYo4VCuHEPrzuMmNpZRL6XV1KzJaJczz8hu55hu0IDbFS5MX5CrZNErV6j1eHhwvI4YkaQdbLYwSQuWLuVDrlnTuOG4CAlKzZtzUMyZY77PX37hw6lVy3URaIXbt2mfBwaa10WJjWVau5cXHZRGuHOHAqd8efPaGYMH8/yMeggrgbRxY9LvbTbObt278/QHDZLigJQoUEDCLl6kPQ8wcuAcXchsGDX4rlSJmtzGjcmiPXa7XXbv3i1Dhw6VvAmaVt68eeV/bdrInuBgsStOSIcOfGljY12n648cmdQEiY1l6FeZH3nzskaJ80Sl6uYqYVW3Lo/jHLaNiWEER11XnToMCzsLbeU7UUmYxYrxd64ms9u36XT19aUwHDjQZRmLbGGSWsydy3BjvXrmWkF4OKnJbm46/8II333Hh1O/vvkLFRbG9HN3d85eRrDb6d8AzPvyiJAyrWk0U4xw+TKF0gsvGG+j2LGOarrCgAF08NlsInfvyraCBcUCSM/u3cVutfLlsFg4Gy5Zknlayp071JxGjaK/Q0VfAgIYNv78c0Na/Llz52Ty5MlSsWJFASBeXl7So3Vr+b1jR4lTHe5KlSLJTb1gR46QV6KcrbVrkyjnqCWGhVEjVNXLKlSgJuGsSV69Su6HovO3bUstx/leRUdTcClGa/361Dhc3dPt23X2bPHi1GxcZaNHRvIcVXi6d2/jiSybZ5IKxMeTL6BCeGb8kOvXKend3flymOHbb/kyN2liXhLw5k0ORg8PzkRmUA3NzXgnIrS1gZTZqS+9xBfPrLhR795U6V0NWlWnRRHqVq6U98B0+a+++orf7d6tVxurV4/XeC8FmVWT80WL+ELXqKELDw8P+gTGjzdt8H3nzh2ZO3eutGjRIrFrXaMGDWT200/LbVWUyMuLrOF163hMm43matu2+rH69k0ebbt0iUJNCYeWLUm6czY/Ll+mGaF6DfXq5ZoJHRVFxrXywTRqRBPT1fP4+2+aUQA1nE8+ca2JxMaS8aqSGTt2TG5uKZw6xfA1shmw5rh5k443gBRoM6LT33/TPvXyMu9yJ6LT6Vu1Ms8EvnqVL5qnJ1VVM6jmWz16mNde3bmTPp+GDc39MyrK8+qrxtvExFDzMKLyu3DO2rp3l1ZubuLj7S3//vsvv4yPT2rf580roa2nSEjpKxI66zBDpI7qfGwsX7b9+1ngeMYMqt6NGunsTYAvYosWNBtWrzYtB2C1WmXVqlXy1FNPJab3ly5dWia88IKc6NlTZ9BWqUJGrwqF37nDz2XKSKLjcuLE5DwMZ2LYk0+65hFdvEgh7+3N7fr2dc1uVvdMaSLNmrnWWESY9PfUUxRKOXOSZezqXlit1KZVC5LGjY1Z03fuMKTs5cX7PG5ctjAxxL//coB4eJhXNxPhbOrnR5XVjPdhtbLaGUDb2sxXcfEibW0fn5TJXqo27GOPme/z7FnOSCVLmnd5Cwujg69UKXNh9+WXPG5oqOv1qoasI/v17Fm55OUlQV5eUqVKFYlyVOvj4yk0+/SREMsOARLo9Eo4eHomLyytlsBACpOXXuLzOnAgVSzXv//+W0aOHCkFCxYUAJIrVy55sX9/2TpihNiVxuTrS17Hjh36y3r+PEPzytSpX58ROWdzYe9e+o0MWKOJOHeOpqeXF6+xf3/XGqHNRqewKiVQv75xqP3SJR7P3Z3jaPRo4xIVGzfqfpsaNagtuRJMVivvr/LdPPNMYkQzW5i4wrJlVEPz509OIXeEzaY3tqpf3zzjNTxcr9w+YIB5i4TDh/ky+/unTPJSDbKeeMJ8n2fOcJ85clBQGiE+nuFid3eaAkaIjKQmUaeOsa9DFXVybs41apSsSDB3XnrpJZc/DV1tY9+ct7dw5n/7bb4M48Zx5v/0U5qSu3ZRMKbB33L58mX58MMPpXr16gJA3N3dpXPnzvLT1KkS3b+/roXUqEF/iqNpu2cPzRt3d5167qqsw86dFO4ABY4Ba1SuXKHW6+nJfb7wguuC4XY776N64atUoVbm6rrDw6k1+Phwny+9ZDw2T52isFNh6u++M9ZsHQVOgwbJJs5sYeII1b1O08gadVV4WSEiQo9IPPusuUZw7Bi1DIuFDkuzgf/HH3zhg4KMa4+ocx07lsd/6inzWfjECXrrc+Uyj/DY7YzeAClHoYYN43Zmwk45Z52dm9evi3h7y8iEWh0/p+QLygBERkbKd999J+3bt09M769Tp4588uGHcm3mTApFgC/gc88l1hUREf5dvZr+DeW4HT7cNVnRkXqeOzfHkys/2+3bFIx+fhwXAwYkb46msGWLzv8oWZKmiKskUauVz82xUbxRJ4M7d6gle3pS85o0yZhG4CxwDCrlZQsThYiIREeS9Oljzs84dYpV293c6PwyEw6rVvElzpMn5UjM5MkUZDVrmtcmiYjQH+6AAebZx4cP0/zKk8c818du5wsCpMys/e03nqeBVpGIpk1pKrq6P88+K7H+/lKndm3JlSuXnDF6ke4BNptN1q9fL88995wEBAQIAClatKiMHTtWDq9enTxd/+OPk0bVrNbkrNFp01wLh+3bdRZqUBCd4a4c684RkiefTFrtzBEnT5IVC9AZ+tlnxtrnhg16xfqQEApDVzAxU5IhPDypwJk40djsPXUqW5iICNW1MmX4gkydai4c/vyTL2ZgoLkvw26nL8PNjYLHrCJXZKSek/Pkk+Z+itOn6ZR1cyMdOqVyB4GBZLiaFbW22XSNZMgQcwfuX39xNq1Tx/w8N23i/owykRNqzp74/nsJCAiQkJAQic+gTN4jR47IuHHjpFixYgJA/P39pX///rJu7VqxrV6tE8NcpeuLUMubN0/PtjVijYrwviqWc1CQMfXcaqXDVEVd2rc3LowVHs6Ij6cn7/WkScb32lHgFCtG4Wc0JhzNFDOBY7MlFTh9+xoLnLg4hvgfeTq91coHZbHwQTgTrxwRFcUXDaBwMAuX3ryp5+N060ZNwghnz5Kwpml8KGbCYf16CrKcOc2zgO12RjgsFjJ1zZIQo6KoiSmNxOz4e/ZwJi9a1Nw/FBHBF7B4ceNrP3aMx5w/X7777jsBIOPGjTPeZwq4fv26zJo1S+qySrq4ublJu3btZNGiRRJ58ybD1KrlRFAQ68Q4vyDx8Ump59Wrk9zmSriePMk8G02jWTppkvG1rlunh79DQox9Uaqkpao+16+fa/6OCAXO6NGpM1MuXNDHo4mZIiKscaJMKhd+kSTYvp3vAiDStesjLExOndLrXvbubU4a27eParAKlZr5R9atY7jO3Z3qrNks//33FAwBAebhZJuN5pS7O19SI7VYhLNn//481y5dzDksFy7opf3efddckGzcyJemRAljO1yEM1WnThRka9cab5dQSEllR/fv3180TZO1Zr9xQkxMjPz888/SpUsX8fDwEABSrVo1+eCDD+TSpUvJWaPVqpHf4/z8bDZyb5QmEhxs7NwMC6PQ9fBgdGbUKGPm8okTul+teHFGtYzu8fr1FF4Ax6VRcqbdznGjOCBmAsdqpd8qIIDnOnGiscCJiWEwwcOD2uycOcZjNzxcJz0WKZKYZ/XoCRO7nQ6sgAC+HGZMUZuNaqunJ9XTVauMt42J4SDTNL7wZs3Hw8OpOqookJnmcOmSToTq3Nm8jODRo7pweOstc0G2bRsHpL9/ykl3X37JQVauXLJEtySIjdVV7k8/Nd/nK69Q4CRoOHfv3pXy5ctLwYIF5ZpJ2Nput8u2bdvkpZdeksDAQAEgBQoUkBEjRsh+Ray6fp3OaUfWqBGRa+1aXfWvXJn3wigk+uWX1Go0jQLbSPW/e1fXGvz8KKiNXuLr1+m3SI2ZcuYMKQVK4BmZKSIMjatSCa1bm4+xDRv0rOg+fYwrztntpEEUKsR78L//JZmssqQwAdAOwFEAJwCMcbFeA/Bxwvq/kdAB0GwJDg7mTNWrFy+hUSPz0gEXLuje+8cfN8/FOXRId34NGmTeK2XrVnrk3dw4E5j5CZYuZa6Gjw+db2Y1SWbPprqbOzdrXxjB0bQrVcrclxIVxYI5ACMUZtrblSu6pjdjhvF2IoxYubmRaOaA/fv3i5eXl3To0EFsToLw1KlTMnHiRClTpowAEB8fH+nTp4+sXLlS97WEhZH96+fHwd6zp7Ff4vhxPlelNcyfb+zI3rFDf76NGplPFMuX6/1qnn3W2By02+mXyZOHGufYscYCJz6ek5qvL69t+nTjcRMVRaepuzvHzsKFxuPm5k3yZ1SUyDlZ0xHnzuns2erVXQqyLCdMAFgAnARQCoBnQtHoSk7bdACbb2kA6gPYmdJ+g8uXp+S3WJhHYRYB+eknvpS+vrRhzV7iTz+lCpk3r7mpEh9PbcHNjQ/OjL9y967+EtesmbxWpyOuX9dfilatzEsiOJp2ffqkbNqpBLXXXjO/X2vX0rTz8eHMaob58zlj16rl0scwa9YsASAfffSRhIWFyVdffSWNGzeWhCLj0rx5c/nmm28k3DGiYrNRa1DFhXr0MObSxMSITJigaw3vvWdstjqq84UL0xQyGguXL+vRwIoVzevtHj+uT1QNGpgL9D17dM3pscfMo3yrV+vs1f79jc0vu518knz5+D6MGmXs4LVayfPx9+fzff99Q0GWFYVJAwCrHD6/DuB1p22+BNDb4fNRAAXN9uuH8hJayEWuhCNu3ZLQttMkBFsktNxL5r6J8+cltN5YbltnjHk5gn//ldAKL3Pb1lNM83tCZ/4jId57JRQt+ZBNiGih72yTEI9dEurelvVSjMwau11CR66QEMsOCfXtZG7axcVJ6IDFEqJtk9DcPcxNu7t3JbT757yuIv3Mq+jfuCGhLd/jttWHGyZLxsbGSv369UXT3ASgH6R8+fLy7rvvug4fHz0qoZWGcL9VXjHVGkI/OZRwb1tQQzVxIoe+u4P3Fi2TqfNJYLdL6KiVvLfubemXMEpTsFol9IWEe+vbiYQ4o2cWHZ1wb7fyOZglQt6+LaGtJic8h2fMKQgXLkhonTHctvxg82d26JCElh+sj3GTiGTozH8EqHhXspgw6Q5gjsPnvgBmOW3zB4BGDp/XAqjtYl+J7UGBYAmpZzC7OkjqEGyl072+wUN2kNQh2jZu28DgIUdHU+328NDp4SEGT+P2bZFBg/Sq7JXDDDYUOty6d9e3rW5iVh05ItKsmb5tLRPn8c6dItWq6dvWMUm4+/NPkRIl9G2N7m18PF+aPHkM763dbpe9e/fKK6+8Ivny5UvQQtwEyCk+PuvFbvQSLVki4uur31uj55Bg2iUev6JJ1ndUlMjLL+vXVTWFJEzH51DThJt09qxI06b6tsEmeVF//y1Staq+bV0TU3jjRpFixVIetyKJGneI23bz++WgcYe47zTfNkHj5vGDRbKYMOnhQph84rTNchfCJNhsv35+wa5TSE6d0h2cdepI6BfHJSTEIN1k/36dKdmunYQuuGS87fr1enTg6acl9Kfbrre12/lSFCwo4uYmoU98JiH1rK73abPRd5Ijh4iXl4Q+t0hC6ttcbxsTQwq6p6dIrlwSOuwPCWlgd73t7dssepOgzodO2Gx8XUePMloDiFSoIKHT/3a9rc3G66pQgds2aSKhs08m2fb8+fMyZcoUqVy5sgAQT09PeeKJJ+S3336TF15YK4Cb1Kz5lGthMncuzzckREK/v258vrdvk9MBSGiLdyWkbrxhKpH8+2+iaRf6xGfG91ZEN+08PKhtGN1bET1q5+9PLcZoW5uN/iYvL5F8+agdGV1XbCwp85omUqaMhM46bLxtRIRemqJ2bQmde85426tX9bSPdu0k9MebxtueOJHo5A1tPSVLaiaZYuYkCw3HxZHN6ONDe3DmTGO/QGQk/QYWC+3M774zVjlv3NBDs6VKmRPbzp7V8zdq1jR37B08qFfDatnSnOeyapXu73jySWMTLD6ewilPHg7KIUOMTbBr1xiBcXdnlGTKFNfqvGKNKk5HpUosAJVwvyIiImTevHnSqlWrxPT+kJAQ+eKLL+SWk+nz9ttvCwCZO3du0mPs38/zaNXKPEx/6RIjNO7uzNQ2w7JlHAf585s/M5WPpaJ2Zqzi9EbtHnvMuBiVCLXN4GBu+/zz5jymHTvIm9E00vfNyjssX87x7eVFNrCZr/Cbb3i/cuVK9JNlRZ+JO4BTAEo6OGArO23T0ckBuyul/SYRJrt26TH9Ll3Mw50rV9JhCpC6bpR1qULOefNy8I4ZY+zUio+nR97Pj07eDz809s7fuqW/xHnyMAJg9JCPHNGFU6lSHBxGWL2aLxrA9HUj2zk8nI5jf386jwcMcC2cVI8WRfiqWJE1RaxWsVqtsmbNGunbt6/4JrS/KFmypIwfP16OmwhFq9UqTZs2FT8/PznimIbftSvvhdGzEKHwK1uW9zgl7soHH6Q9H+uZZ8wZwCofy82NDl+zqN2ff6YvameW0xQfT/+NxcLoklnSZoJpJ4qLY+YQvnmTSaUA0yUc3p0sJ0x4TugA4FhCVGdcwneO7UE1AJ8mrD/oyl/ivAQHB9OJptT5QoU4YxrhyhWd4l6+vDk79vBhve5J/frmDctXrtQJcB06GIen4+NJOMqdmwPyhReMw9O3bpFMp7SGadOMnYC7duk5JKVKJdEakiA8nPwI1Yqye3fXTc0vXSKbVG1Xty73abPJP//8I6NGjZJChQoJAMmZM6cMHDhQNm/ebOwHccKFCxckT548UqNGDYmJiaFm4OvLwW+EuDiGcL29zaMqdrue9d2jh3k+1unT1LbSmo9lVonfbqeGp2mc3Myidnfv6izWVq2MiWrqXFXU7qmnzLlJjlG74cPNa9w4mHYyZUoyTT5LCpPMWIJLlyZbT9M4EI1ucFwcVbzAQPobJkwwvsE3b9LT7+5OP8annxqbSocP64Sj0qXJBTEakKtW6QKneXPjKlcxMTTPlJkycKAx4ejQIZ1QljcvtSGjep8TJ+qJaB07umZj7tzJgap6tHTpIrJhg1y5fFmmT58utWrVEoBd6x577DH58ccfJdrMJDHBsmXLBID873//o1mjGLtGeOcdbrNokfmOJ0yQxDCqWej7yBG+RLlypZyP9eGHmZuPNXWqORlR5WPlyGF+/Sp3zMOD/jqz64qJYdU3NbEamHaPjjABOLOYpfX/+afuMGzZ0vVMLEKB88knutYwaJBxsaGbN3UzJUcOqtRGwunIEd35Vbq0ORvz2291clSLFsYC5/hxvS1BQABfIFd+Eacm1NKlS3IfTnQ0HZ+qCHFAgMgrr0jU33/L999/Lx07dhSLxSIAJDg4WGbOnClXzWz/NOCVV14RAPLbb7/xvj//vOsNr12juZBQyNoQCxfyGp591vzlPHiQfpSgILJKjWCz6blbKeVjnTuX+nysDRso/HPmNCeVpTUf6+mnea6PP25eyFzljilCponQe3SESZEixs6nQ4cSPf5StiydcUYPeOVKXS1s2dLYpHE2U1580dipdv48zRiLhQLn/fddCxy7neE9dfzatY1p4k4NkmT4cNdmkqsm1M7+kxMnKGiUKVOhgthmzJCNK1bIgAEDJEeOHAJAihQpImPGjNHLMWYgYmJipGbNmpI7d245364dzVRX2sSUKTxHo4lAhO0nvLzoKzIrKHXsGF/kQoXM92e18h4qU8FMOO3dS+GUUj6WCP0nac3Hevzx1OdjvfNOyj2agoI4JlPR5+jRESausoZv3qQPxWKh5P/wQ+PB5WimlCljnAhmt3OdMlNatDCe0a5f5+Dz8uLL/L//uRY4qi2BGgQVK9L55ur4//5L9VnTkjRISra/jRupfahygoMHJ03gi4vjMZQvyGIReeIJOTZvnrz5xhtSokQJASB+fn7Sr18/CQ0NFauZqZABOHr0qPj5+UmTSpXYIc/VAG/cmDOpEeLjGd4PCjIvXXnzJkP7efKYtxaJj9dn+TfeMH85t23jOCtePOVqd4oB3aGDuc/j0iX66YCU87G2b09fPpaZIHXAoylMHP0iqTFTHP0iH3xgLHDWr9fDt2XLGpspd+7Q3AgI4PH793ddZcu5QVKxYjRvXL20jk2o/fyYaOYsmOLiaEersGKePCTVOW53+jRDiKruRpEicnPMGPnsvfekfv36AjC9v02bNrJgwQK5a5aH5ApxcRRaW7ZQ6C5cSNNp0SLWf9271/TlmTdvngCQCblyuS4ZGRjI52mEb7/ldS1ebLxNfDwnAU/PlB24/fpJ4ixvho0b+RKXLm1OiXdMlhw1ytyXs2sXfTm+vqnrhuDpyeP/84/xdnFxem2blPKxnPDoCRNnv4iRmXL3Lp18OXOmLHD27GHdVIAPd/Zs1yZVdDSdXiqHpFs31x58u50vliLIGTVIcm5CnSsXBYGzOXP9Om1zVcm8QgXOPCp6obSQtm0pjDRNYtu1k1/HjpWujz+emN5fpUoVmTZtmlwwy/9RsFp5bQsWUDtq147XodpOpLCEFn5GQvKfkNCpe5LNtk8//bS4aZpsBJK/RG5u1BBcwWajyVCjhrkGoZqxf/ut+TWq+rvjx5tvt2YN/TgVKphHYaKjdZ/Z9Onm+1y1ivssXtzYXyZCwfjKK5IYBUopnN60KbdNKR/LBR4dYVKxol6b08wvEhND56qqgdG5s7HAOXxYj7nnyUMzyVXEwrlBUuvWrtPH7XZqM8rhVaKE6wZJ8fGcWVUCWKFC1JicbeW//qLW4+WlC8/ly/WX89gxajAJ12ovXFh2PP+8DO7bV3Lnzi0AJF++fDJs2DDZt2+feTg3IoIRgfHjeRzVgArg8WvUYHLhW2+R8LRyJaNER49SUzlyhPfkl19E3n1XQgIPCZBQnb5EiSQlC+/cuSNlypSRIu7ucqNo0aQhXR8fmo6uoFpvmDVHO3KE59utm7nAUQ3b+/Qx327rVu6valVzIlpkpG5SfvaZ8XYiHLuengwnm+3z7l2dCPfqq+ZclwMHKJi8vFJuHmeAR0eYAFSBjcwUq5Xqturd0qwZbVxXOHuWL6mbG1+a8eNdR0giI5M3SHJFoFJNqBWRrnRpvnDO2k1EBM0zlRlavjzL6zlqLHFxLMKjqmX5+lKrUqptVBQHi5qBLBY507q1vPP001KuXDkBIN7e3tKrVy/5888/jUspxsSQQ/HGG3q3QYD3pEYNqsrffstoSDrKMYaGMs8kdNw63SdQrlzi/duzZ494uLtLZ0DsjhXaypenL8gVRo2iH8CskVr79hwnZombhw9TaDVoYM7APX2avpnSpc1LWNy5w+ehaXzuZliyhPe6Th1zLePOHY4BNzdOSGb49Ve9VYtZjZQU8OgIkwIFXNt/djtnQ+UwDQ6mCulqtrl4UW9L4OXFCu2uTJ+7dxmRUeX3mjVLXmdUhAPRkTVarhxT9J1fvvPn+SLkyiWJxLgEYlgiLl2iH0a1myxZkpqSoqjv3Ut+TcI+wkuWlK+7dZOmDRqISu9v2rSpzJkzR8Jc+SxUg+8ZM/jC+fgkCiOpV495IitXmr+o6YXdTo1K9Yr53/9EYmNl+vTpAkA+cXPTWZt9+tDJ6Or5tWjBczXC3r3c/+TJxttYrRQigYEptzapXJn328yBGRdHTc5iYZqGGRYupHAICTF3yt6+zTFisaRcEmLBAu6zXj3z60kFHh1h4iqas2aN7peoUIGqq6tBeP580gZJzz/vmoJ/5w4HovKJtGrlmjmb2ibUu3cn7dHSo0dSbcluJw+hRw9dM2jXjo5Nq5WC5JNPEgv7xHt6yp/Nmkmv5s3FmwWApVy5cjJp0iQ57YqJGx1NB/DLL+samzIThwyhup0ZwsMIUVEUJAkC1X75snRs3Vo8AdlXuTKvWbUkdVUUqUgR0uCN8OSTdLKbvagzZnD/8+cbb2O1UuC6uxs3KBPh8xs4kPtzzj9yxpw51FyaNzfnsNy4QTPZw8O8SJYINRZNo5BNqyPdBR5NYbJjh+60LFaMA9CVKn7mDE2ElBok3b5Nb75qm9CuneviRxcvsr2CYzlBZ43FauUgaNJEEolhw4Ylpd2Hh5Ntq3JrAgPp5Dx+nNpKaCjDw15eYgdkX/nyMqx5c8mfkN6fO3duGTx4sOzYsSO5H+TyZQ6yTp1oIgHUQjp1Mm3wfV+xZEmi8/Ha9u1SKDBQygMSMXEiNUV3d95nZ+TMadxfOTKS+zSj6V+5wnvSoYO5n2T6dN63zz83vw7VifH11823+/lnfbyYUf6vXCEx09s7eeMzZ3zyCffZvr35PtOAR0uYbNumO6SCgjjLuCKHnTxJEpK7OwXJSy+5Dt06s0Yfe8x18aVDh6jNKGJY797JiWG3btE0UhpA8eIcbI4z/969FGiq41xwMAVhZCSFzfjxiazYizlyyLQGDaRqQplDDw8P6dq1q/z6668S6+gzstt5fpMnUzXWNP34L7/MQZlBgy1DsWcPn2FQkKybPVs0QJ51cyN/o0sXmpjOzzYgwLh38m+/8brXrDE+5pgxfH5mBLITJyiUHnvMXOAsW8Z73a2bOTdkz57U+WcuXKC/yNc35cRGFa3q0sU8FyeNeHSEiXrh8+ZlfoMrVfHYMdKrLRaaNEOGuM4iPXIkqXB48snkarXdzoGpmLWuiGEidIy++KKuBTRpQnNLaUp379LJ6txxbudOvuSLFiWWALwLyMIqVaRNtWqJXevq168vn332mdxwpEzb7fz96NF6zRUlnCZOpGc/De02HxiOHqXpkju3vPH00wJAFpUsSdPMVWi3ZEmaja4wYgSfkRGH6PZtCqMnnzQ+H5uN/rEcOczLZx4/Tsd9cLC5eXHhAn1gxYsb51yJcMKpUoXnlxIvZuJE3psnnzQvR5AOPDrCxN3duEHS4cNkMSrq+auvuuYD7NzJmcSsCXVMDAex6ieSPz8LODt6861WhoBVHVBvbwoIR23ln3/IzlWNsVXHuVu3eB4vvSSSM6dYAQnNn1/61agh/n5+AkBKlCghb775phx1nEHj4xl9GTqULyBAzat1a5pMZqn3WRknTogUKSLxefNKw5IlJQCQ4y+/zMhY+fJJ/VBNmhiXuuvYkb8xgjJJzMocKkKcWWtVm42Rmxw5zEtf3L3L0H9AgHlJgPh4TlgWS8r+GdVO9pln0swhSQ0eHWFSs2byq3dkjfr60s52RT1fsYIzjiKGuWpCff06hYZyqlatysHlqEbevJnUlClShOaFEjRRUaxX0qgR13t6cibdtInn9cEHiX6Sfz09ZUylSlIkKEgASI4cOWTAgAGyceNGvbJ7XBwjLAMG6E5hb2/mb8yfb1iHNcMRH09TLDyc9+DGDfOcmLTi6FGR/PnlbIECEujhIbUBiVUzsCNn4n//43N25R8rV46ObCM0amQubOLjGZULDjbX6j77LHUCp2tXTm5mNWlEdEJaSsWflB9n4EBzs+oe8OgIE8e+OY6s0Zw56QBzDvHGxzNUp7gfRsSww4dppnh7S6JDyzn5bvdu8lLUNo0a0YmoBvXBg9QYVOi3TBnatRcv0hnbubOIxSJXAZlZvLgEJ7S9tFgs0qFDB/n+++8lSvk1YmLYSuLZZ/UyAgEB9NP89FOGeO0TYbfzHDdt4ks7aRLNvzZtGFEoXjwpec158fSkkKtWjdf4v/+RPbxvX9q5Kfv2iQQEyC+FCwsAGZEzJ/dbooQu0Bct4nFdVbQrUoTaoStcu8YX+623jI+/eDH3bVaw6MwZ3o9WrcwFzqefcl8ffWS8jQgdvIqQZobff0+df+Ye8egIk1q1krNG338/eWgzLIz8jGLFxJAYZrNRW1H+EC8vOkYdk7ec0/X9/Ch0FPU5MpKai8q78fTkC79+PX0Ww4aJBAVJNCA/5solj5UuLe7u7gJAatasKdOnT5crSjuKiaED8emnqT4rIfnMM3T0pbOOSBLcukUhOX06NZ2QEF34OS4FCvCaO3Zk1vKrr7IW7dSpvK8zZnAfkybRZzNoEJ2VVasmETyhnu0lJNe/EjrkV/P+Ro5Yu1bEw0MG580rAORPlRj54Ydcf/UqX6oJE5L/tkCBZH18EvHdd9yPWZe9atVILzB7WTt25Dgwu57Tp7lN69Yp95G2WLhPM5PlwAHe11q1MnYicYFHR5h4ehoLBxE+xGHD9MhMkybUChwHx507DKkpp2WBAnxRHCnNp04lS9eXjz/WuQv799MRq176ChU42I8c4Xa1aokdkM0WiwwsWVJyJvhBChUqJKNGjZKDyn6OieGM07evvq/cuTm7/vnnvZkR4eH0r0ybxkZWinGrljx5eH9eeon3Y+VKnv+9Rn3sdjonv/tOQgqcFCCBTg/QzzB/vnkRIRGR+fMlCpCqXl4SBMil4GAKVqV5hoTwxXdG8eK8l67wxht8cY3u6fbtPMevvzY+r127uM2UKcbb2O0UIv7+rqOHCqdP85qqVjUvN3DlCifFQoXMHcIZhEdHmPj5JWeNinAg9OhBNdbdnT4K5xno2DGq4ErQ1K9PlVkNLpuNtm2HDpz5LBaqlGvXcoCEh9OmVTOllxe1iHXrqDk88YSIh4ecAGR8gQJSKk8eASC+vr7St29fWbNmDdP7Y2N5nH79dMdsYCAFyMqV6ffOnz3L2XfwYJp1jsl4JUqw0NDkycy9yaBiRykhNJTvfejc8yT4KZZwUBDPxYws99ZbcggQX4tFWtSrJ1ZFNBTR+RXOjtQGDWj6ukL37iTqGWHMGD5zMx9U167U5Mxe/q+/5rmZtVa122km+fubc36iojhOfX3Ni15nILKUMAHwPoAjYMvPXwHkMtjuTELt1/2pvYAkpLX4ePoslImRKxfp6o7edZuNL6iqYeLhQQHgyCO5coUDXTlUCxSgXX3+PB/61q30laiwb9Wq1D62bGEoMn9+uQXIF/7+ElKwoAAQTdOkVatWMm/ePImIiNCdqP3762ZFzpz0iaRHA7HbKRy/+IKFkFRkB+AAbd2aZsDKleYVuO437HYKX5WsGRjIkpWuBKjdLvLkkzIHTBN4V+UpbdtGB7CnJ8P+jujRg74qV6hZkyatESpWNBZEIjR/AXOfy+XLfK5Nm5qbSnPmcF9mhDi7XW+Da1bnOIOR1YRJGwDuCf9PBTDVYLszAPKmZd/BwcGczaZP11/+0qU5UzlyToxMGRXlUYO6Z0+dwt68OXMgYmOpTn/4oV4Nzd+f/pTVq7nf2rUlFpDf3NzkiUKFxDPBD1KpUiWZMmWKnD9/njbwunX0sShzKSCAPpA//kg70ejsWfpvnnkmqfAoVIh8g48/Jk8mHQl5DwS7d+sZthUruu5gFxkp9po1pZe7u1gsFtkaFERhHhfHScHfPyltfuJEapWuqPQVKhhHek6e5HnMnGl8vs88wwnFTDiPHJkyIe7CBZq0zZqZC5yvvuI5vfee8TaZgCwlTJIcBOgKYJHBurQLk/z5dd9C48b0hzg6rg4fZkTFyJS5eZPedSVkAgPpYzl8mA929WoOOA8Prm/QgJGJH38UeeIJsXt4yC5AhuTJI3kS/CBBQUHyyiuvyN69e8VutVKTGTqUAgzgAOzVi+eaFidqWBijCi++qJsHACMnPXtyVjt69OEgphlBVbRT/pyhQ5M7GM+dk7CgICnp7i7F8uSRW8pnoRL6PvhA33blSn7nij1apgyd466gShEYmRLR0XSoGjl3RTi2/PxYoNvseh97jKRFs/qu58+nTuBkArKyMPkdwNMG604D+AvAXgADTfaR2B40GOCAcEyxjo/nS6fCxJ6eSU0ZZar07avXBAkJoSMwKooz/oQJemHn3LkZvfj550Qz5iwg7/n6SoXAQAEgXl5e0rNnT/njjz8kLjaWYcqRI0WKFpVEf0q3bixJkFrvu9XK65o4kS0OLBZdm+nUiRGUv/++74PrvuDuXT35r2zZ5IWCNm+WnRaLuGuaPFGokNi9venkbd6cmpkS0jdvUjNxZYqUKmXMnFUFkoyS75Yv5/oVK4yvQVXJNyOnqbC2WdEku51meUoCJ5Nw34UJgFAA/7hYujhsMy7BZ6IZ7KNQwt98YJOuJikdN7hqVf2qL13ii6eKFRUvTqeeci6GhdEJplisAQF0Th44QBPj++/JpUioSiatW7Ny2UcfidSqJeGAfOPmJs3z5k3sWte4cWP56quv5Pbt28yFefNNPaXew4MhvgULUp+Fe+UKQ8u9eummEEDS1NixzFbOYLp0lsb69RQO3t7JCWGffSbTEvwnX3h7c9Zeu5b365NP9O1CQugfcUbdunzGrqBMRyMMHEiTysg0jYjgJNS5s/E+4uOZChAcbB4GXrAgZYGTichymgmAfgC2A/BN5fYTAIxMabvg4GC+YE8+qfs62rZlNMVqpVTfsYOREeUwrVmTpkpEBAXJ//6nZwYXK8aQ4ezZIt26Sby7u6wEpE9goPgklDksU6aMTJw4UU6dOsVw3uTJDEsCtI9btqR9a1bkRsFmo69gwgQ9T0f5dJ55hjOXWYHkRwFXr+opCi+/rPuA7Hax9esnbQHxtljkIEAHdOPGFEAq3Kwq2zunFnTtynQGV2jUiE5TV7DZWBjLrO2Gii6ZtWBRbTl++814m8uXaXo3aJApVPnUIEsJEwDtABwCEGSyjR+AAIf/twFol9K+gxX7NFculvVTbSlv3+YDVVqInx9JWTt3MtT3+ed6SNfTkz6Hzz6jYAkKkgOAjPDxkQIJbS8DAwNl0KBBsm3bNrFfukTHnIoaKV/MzJnmlbwU7tyhyfTcc7ofRdO4j3feYXjzfvk97Ha+rLt2kd8ydy4dzWPHknjWowdV7DZtaDaqPJh69agJdOpEU2HQIEbOZsxgftL+/WkqWpwirFaajWqyUA7VqCi5UrWq5Nc0qeTjI5H+/vRnOXI/jhzh5/ffT7rPIUPoh3B1r+vUYYTJFZRzdvZs4/Nt2pRJekaw2ZhCUaWKuZn6xBM0kVNZST4zkNWEyQkA5xNCvvsBfJHwfSEAfyb8XyrBtDkA4F8ktA9NaQn29WUcPzJS94X066dXDAsOpqkSHk6V+emndfp71apkbE6aJFK9ulwC5AM3N6meM6cAEHd3d+nSpYv8/PPPEnP5MtXsli11vka1atRKUlML5PRpvmitWunO3Jw5KcTmz89c7SMujibYkiU83xdf5ItSoYJ+n5wXNzc6dsuX5z2sX5+zdbNmvIY2bagB1KhBZ7BqjO28n9y5uf3rrzOcea8kq6++ogZavbouuE+fltX+/gJAXlBJjh078v6qSEv9+tRCHAWHyqdxxVytWZOC0hWUv2TLFtfrb9ygf8ux7KQzli7lPsy6823bxm1SqpCfychSwiQzl+DgYJoTM2fqRYX8/fnC7NlDjsmkSXr0I0cOhnSnThXp1EkiLRZZBEi7HDnELcEPUrduXZk1a5ZcP3uWDtMuXai9AIwAvPmmeX8UEQ7a3btpMikTCOALPHIkK6lltO8jPp7ntWQJzaYePXhPlPBSS548FBDdujFyNXMmB/fOnXTw3b6dPqeu3c4Xac8eRkM++ID3ukaNROdxKFpIiPdeCe3+OcPk6bkHq1bRZC1dWhfkf/4poxP8Jz8CfEZubtQ0RTihADR5Ffbs4Xc//JD8GFWr0gxyBZVpbFT/df58rndVA0eE96luXTqAzcL2nTrxWZlVYLsPeHSESe7cuqZRty61hxs3OEBUiweAM+qECSKDBoktMFDWAdLfx0cCPD0FgBQrVkzGjRsnRw4eJGmsb189p6RgQb50u3ebmx+xsQxFvvSS7gR2c+MM/sEHugmWEbBaKTjmzWP4NCQkqZahaXzZOncmk3PBAoY5H9TAjIoS2bZNQkpcZPBM26YL9/79OcunxbTbvp2+hEKFaMaISNwbb0g9QHK6uckpHx+aXxYL71NYGJ+nI7U+NpbalKuq98HBxmaO4gkZ4YknOGaMBLISYmaM2L//5jZvv228zX3CoyNM3NwYkdm3j8vQoboztWhRhnRHjxapXFkOAzLWYpFiCX6QgIAAee6552T92rVi27iRQkCl9OfKRR/L2rXmjq/bt6mq9uypc1l8fTmrzZ1rXr08tbDbqYovXMjradRIr8qm/EGNG3Pd/PkkqmXFKmriQKf/PYraUP/+utCuUIF5Q6ktoXDwIM2rAgXoU7Ba5VSjRpIDkHpubhJXty6fY5s2vIeq6ZpjTZvGjV1Hejp1cp3rI0IhU7u263V2O82rAQOMz/uNNzjJmI2NPn14X1LjxM9kPDrCpGpVOlpV1rCXF73sY8aItGkj1zVNPgGkTgKhzM3NTdq3by/fLVokkdu3U9CoTGIfH0aFfvvNnI16+TKjBm3a6BGk/Pk5gH7//d5f5Ph4zl4zZ1JIqcr06hxDQvhizJ/PWfcBefkzDBER9HuFhEhiyN5V+QhX+Pdf3vv8+ekXun5dfkzoDTQGYI0XgP6akyepsY0dq/9+6lSudy5oNGiQsfbRsiUFuitcuiTJQtPOqFqVjmwjnDhBYfPaa8bb3Ec8OsJEmTE1a9IU6dtXYgIC5CdAOnt7i3tCmcPq1avLhx9+KJd27CAhSbXAsFgYrVi40DxZ6/Rp2sqNGummU5kyfODbtt0bcezOnaSNrhy1jmLFSMr79FNGSB4Wanx6sX8/BaimUXAOG5by7Hz4MIVJoUIUGNu3ywsJ/q/V7u5kNxctSiLcE09QW1HRoEOHeJ+dG2RNmsTvXU0MTZrQbHaF9ev5u9WrXa8/dYrrVfkEVxg4kJPiPbaoyCg8OsIkTx6RQYPEXqqUbAXkRXd3CUzwgxQsWFBGjhwpB9av58uoZj6AQuHTT41nP7uds96kSbrWAzCK8PbbVLHTG76NiqK+P3as3gdF+Vdq1GDIcvFi8/J//3UcPsyonJsbzdZPPzUXpAcPcrsSJUTOn5fIqVOlEiD5NU2ulCzJ+ztqFE1AgORGET7DMmWSk9dUeNlVwaWQEEaoXEEVNjJ6dqqlhhGTNTaWmln//sbXep/xyAiTKoC8DUhpHx8BID4+PvLUU0/JqqVLxTp/PrUO9bJWqcLQqFERG7udfIsxY5IWZA4JoQPVVTuM1CA2lgWB336b/AMVGbJYyFUZO5YRivvZq+ZhwYEDpMir57dpk/G2u3fzRSxfXuTaNTnYsqV4A9IGEFu1ajRJ//6bTmlH7eSNN6gJOZLazpwRlxqLCDk2RszZV1+lz8xIU23d2pgoJ8JxANBcziJ4ZISJSu9v0aKFfDtnjtz54QeaBYrtWqwY/SJGfYVtNnJThg3T82gUV+Gzz8wbUhvBauXAnjKFfhV1LprGylgjRzJiZGZWZUOH3U6SX4kSvIevvGJcTGnjRkb36tYVuXBBPk+ozjZNOdXr1+ezAfTG5CdO8PO77yY9Zr581I6c0aGDa6etCP11FSoYX0v+/MZlJEUYTPD1zVIO9PQIE42/e7hQpHBh2fbxxygWGgosWQLcvAnkzg307An06QM0bAi4uSX9kc0GbN0K/PQT8PPPwKVLgKcn0K4d0K0b0KkT95EWXLoErF4NrFoFrFnD8wCAypWBFi2A5s2Bpk3Tvt/MgNUK3LrFc1R/b98G4uK4ztXi4QEEBCRfcuQAChfm38xGZCQwZgwwaxZQtizw7bd8vs747Tc+x7ZtIW+9hR4hIfhNBFtz5ULdsDDg00+BdeuAlSuBkyeB/Pn5bC5dAo4dAzSN++nUCThxAjh8OOn+Bw7kMa5eTX7sxx8HTp8GDhxIvi4sDAgMBKZNA157Lfl6EaBYMaBOHeCXX9J4czIPmqbtFZHaafpRWqVPVlgSyzb6+DA57vffXRcWio9nmPell/RK897eDOEuWpR2EyM6mvU6R47UKftA0pwas54omQG7nWHVfftY3mDGDKrdXbuSHl6ypF7JLaOXwEDO1l278pgzZtAJmRmV8teto5bi5kafliuTYvZsnle/fnJr2jQpBkhJQMIKFKAptHEjzUzV6W/uXG7v2KxLNbVy9n+obGBX46xjR2qfrrBjB39nlIuj+CcptRS9z8CjopnUzplT9nz6KWcEf/+kK+PjgfXrqYH8+itw4wbg6wt07Ah07w506JD8N0YQ4ay1ahVntA0bgOhoajSNGgFt23KpVk2f2TILd+4Ahw5x+fdf4Phx4MwZ4OxZrnOEnx9QvDhQtCiQNy+QJw+1ozx59CV3bs6YXl6Au3vyxWLhvYyISL6EhwMXL/L4jktUlH4OZcpwtlVL3bq8b/eCiAjgpZeARYuoUS5YwOtzxMSJwPjxwPjx2LZ5M5qsW4fuABZ7ekJr3Zr3ZM4c3seiRXmfatcGli/n748cASpWpCYzeLC+3zlzgBdeoAZSokTSY7ZrRw1kx47k5zxvHvDss8DRo0C5csnXjx8PvPMONR7na3mAeHQ0E+fG5aqeav/+eksIf39qLT//nHLhYkdERlLTefFFvbYJoDf4/uOPzK0MHhbGsPNXX9Gn06ZN0opqSruqWpUkq6FDGXL86SfOctevP5hCSSp5cM0aVgXr2jXxvEPRQkLctktoo/HklqQmMdLsOJ9/Tod20aJJqfJq/bPP8j59/rm8m1B7Zo5KL/j0U46NLl24vdI4VFKdivQ4M2HXreN2rmqZtGzJ2jOuMGYM/XFGKQTt2hn7Yh4g8Kg4YIODg/lwVqxIWk81Rw4m9S1dmjZn1pkzHGTt2+s0fX9/kp8ys8H39euk4r/zDl8+5Qx2JKvVqsVrmjyZqvKJEw8XYe3yZQkpf4MBMo9d+rXVrk3nZ3qTAPfupQnn5cWQuiNiY/mCu7uL9bPPpIWmiQ8gh/z9SUobN04SzZurV7mPF1/Ufz98OIWVo7M8PJyOYBVedkS7dqTiu8Kzz5rXSClRwrjy2wPEoyNM8ubV6fM5cjD3Ytmy1NdTjY9nuDGBcp84wMuUoe2/Zk2GNoEWEb1XzeTJJFE5aj1K8+ndm+uXLWNI+mESGiZIpNOvsZOg9u67jLAonk2HDtQg01pQ+/p1UuMBahiOGtnt2wzH5solF4cPlyBAqnl5SZSHBwV3qVJ89vHxTEz09tYJY5s2cZ/z5yc9XqVKLLfojOefZ16OKzz9NI/lClFRFFBZIBfHGY+OMHFz40NKiwC5cYOM1969dVPI3Z31Oj780Lz4b1oRF0fuyvTprnvVlCrF76dNo4M4I2uAPEw4cYJagkodCAoiySwtLNCYGIZyAea2OI6HU6eYd1WmjPxZr54AkMGKBa3qpHzyCQW3xaJ307Pb+YyaN096rGefZejY2Yx86y0KRVcEu169yF9yhQMHeA7ff5/6671PeHSEiZHn3BF2OxuGT55Me1bVI8mXj4NiyRLX1cvTg7AwmitvvMEBqDgmivPyxBM8jzVrskQSV5ZDfDx9Uaonr5cXuRdmzascYbfTTwPQvHE0T7ZsocnSpImMSKhZ80vJkjR3mjZlpOvKFY4JR+1E9Td2NHFVLRRnIqMqdeCqYXz37saENcW4da51mwXw6AgTI/s0Pp7hv+HDk1Zyr1WLs8fOnRlTiPnsWYaBVaMrlbfj5sZjDR3Kcgj3ofPafw7HjzN50sODmmP//uwNlBrMm0cNo06dpNm5CbVGYrt2ldqaJrnc3eWshwed2x4eNJOPH+dvhw3jb86e5XNVJDcRvXqbcxmBP/7g987OYBH63YwykZXASkuA4D7h0RQmd+8yO7RfP70Ys6cnnWKff37vL7TdzkH0+eeuG121akV7fc2abHZrRuLcOWZJe3tTqAwfnjpzcNky/qZChaRckdGjRQA53r69+APSMDBQ4gFGdQAWrurXj79VGkarVnSKO5ovZcuyZo4j/vmH+1iwIPn5dOrE3CtXUA3usyAeHWFSrRoJSh076mUDAwPpR1my5N5f6jNnRL75hvtzLAXwMDa6io3lS3j9OkOy588zT+n4cQrJf/6hmr1/P30YV6/SMZhV+vBcvkwHqabRpzJ7dsqO6Y0b6ZgvXlzPybJaOV4sFllUvrwAkDeDgvgyFylCU+ToUU5Eqv3or7/yuf/4o77vESOSR3qsVmZ9qwpvjjCL5rz0Eq8pCyI9wiTTSGuapk0A8AKA6wlfjRWRP11s1w7ATAAWAHNEZEpK+66tabIHIHmoSxcujRqR/p0eXL5Motu6dVxOn+b3QUGkxauldOnMJ6elhKgoEtXOnQOuXUtKj3f+/+ZN4O7d9B3HYklOoS9UiPfccSle/P7Q6vftA155Bdi8GahRg6SykBDj7ffuBVq35rlv2ACULEmyXf36wNWr6B8bi3lRUVjr6YnmVaty+/feA65fB2bOBA4eBMqX5xIUBGzfzv1u3Ag0a8aUjG7d9OM1bcrUBLWdwhtvAFOmALGxvKeOGDgQ+OMPUvqzGNJDWstsYXJXRD4w2cYC4BiA1gAuANgNoLeIHDLbd+3ChWXPihVA1arpe7lv3uQAU8LjyBF+nysXB4rKq6lc+f4Lj4gICgvFblXsUvX/9evJf6NpZLQ6s1wV09XPj4LWFdNVLXY7BY8rxqsz6zU6Ounxc+fmy1qzps54rVIl/cLdCCLMxRo5ErhwARg+HJg0CfDxcb39vn1Aq1a8/g0bgFKlyByuWxd3AwNR+/RpRLi7Y7/ViqDgYOCffygs2rThOPjtN+YEDR0KbNsGNGjAnKV8+cioXrBAP9ZrrwGffEI2siPT9/PPyaS9dAkoWDDp+T33HBAayokhiyFLMWCRih44ABoAWOXw+XUAr6e0b0MHrBFiYvScmho1dLPFz49EtfffJwnqfvE67HaaGytXMiz9/PPkXSjujOPi5cXQYps2VPfffZfO382b6Zi8efP+dvdTTNedO+lknjqV6nrLlnrIXbF069cXGTpUQkevlpDgGAkNzaBziIjgMQGWHjDrU7NvH+9rkSJ6PZEVK0Q0TfaVKyeegHT08xO7jw99YC1a6J39Nm3isXLlStqfeOBARuwcTR0Vmdm9O+nxf/vN9fcidPyWLJnu25CZQFbymSQIkzMA/gbwDYBAF9t0B00b9bkvgFkG+0tsD1qsWDHzO6GcpjNn6i0WAXrumzVjotjWrZnfJc9upxNwxQoKjeee4wum+iSrJSiI5zVoEF/O779nZODy5YerDajdTl/M4sV0mDZuLOLnJyHYQgas7z6Gz3fsyJjrWrOGoXc3N/JTjDhH+/fTOV+8uO6UnTxZBJCP8+UTADLd01OPAM6eTeFTsyYnmNGjkzYi376d2zl2G1S1UJybnqtEvqVLk59Xr1506GZB3HdhApP2oADyg34QNwDvAvjGxe97uBAmn6R0XJeaiWrwPXCg65ya33/P3CrtcXEctHPmcNasV08vNq2WfPnIQ3n5ZXIWNmz473fus1ol9MsTElLsvIRWeUXn+wQFMdN6yZK0NXJ3Rng4NTZF0TdKfdi7l0K8XDnySux2kR49xA5IZw8P8QBkD8Cxkzu3Xj1tzhxu7+PD8xXhbytUSJ6PozRIR1y7xv24avPZty+FYRZEltJMkhwEKAHgHxffp9/MsVqpartq8J3ZOTXx8SwZ+O23FAz16+s5PYri37QphdhnnzG6kBEV6/8LuHmTZlqfPrpZFBjIe7VvX/r3u3QpCWi5cjE87ApbttA8qVqV5xERIVK1qtzw85MigJRxd5c77u7UYHv04LjKl48Ca/hwCkLFeVFFqR2Z0yNG8LfOpS2KFHHdLH30aIa9s6D2maWECYCCDv8PA/C9i23cAZwCUBKAJ9jZr3JK+w4ODLx/Db6tVhYgnj+foT/nXjX+/iw2PHy4yHffcXBlwcGRJREfT1OlVy+9rGWtWiSFpSfF4ORJvXbvqFGuQ/dr1vBYdevS53HihEiuXLIxd25xA6Svh4fevvWDDxiSfu01mpze3noVtkuXOIGNHq3ve+NG/m7JkqTH7NaNeV/O+Phjbn+/a+CkAllNmCwAcDDBZ7JMCRc4tAdN+NwBjOicRGrbg7q7Z16D79u3WV5x3DhqF45V4319WZT61VdJUDp8OFtwZBRu3uTLVb26JDpwBwwwLsJshOhoZgAD9NlcvZp8m99+oyBo1Yo8nAQG64SEdqPzlElasKDIU09R2zhyhOxYi0VvrNa1K00ixWCNj6eWpcwhBdVE3TmV4uef+f3evWm7xvuALCVMMnNJczTHCI6NrgYNovqrqPEWC23wIUNI0/7nn/9MFm+Wht1Op6Vq/aCSOlNqzeqMhQspkEqUcP1bVWWtVy9OCBMmiBWQJm5u4gfIUZUe0bs3zadWraiN+PhQwIjo2cWff67v9+mnqTU7akVGtVB27uT3RmbZA0S2MEkJ8fGcBVw1ugoIIE164kQ+/MwsgJSN1OHSJfohfH0p5J94gszj1GLXLposOXIwDO8M5fd45ZVEhux5NzfJDUhNNzeJUQ70wYP59/vv2SgM4HnY7ZxwypfXNdSffuL6Vav049y547rUwIULyYVRFkF6hMnDWbaxdm3Zs2dPyhtGRLCU3tatwJYt/D8ykuuKFWNh4oYNyZ6tUiU5QzErwGpl4WdnZqvz57Awllm0Wlk826hItM3Gxdc3eaFof39z1mtqy11mNG7cICv1k09InuvSBfjgA5aGTAnnzrFI9L//Ah9/nLQUowgwYgQwfTpZqgMHAnXqYNnFi+gSE4NXAMzImZPXnTcvGcc7dpCcV7s2y3l+9x3w1FMs+9ihA5muhQqR+LZ4sX6sqlWBAgVYeFzBagW8vUl4mzw5w25XRiBLMWAzE4bCJCqKTEXFbN2zhy+OmxvrtCrB0bAh638+aMTGAufPu2a7XrqkCwkjWCw66zVXLjIvzViuFgv/urnxXhmxXZXAdUbevMnp9CVKkHJ+P1INwsMpUKZO5b0bNgwYNy5lOn9EBLsW/PEHKfkffaR3L7Dbgb59KRTmzQOqVwfq18crHh74OCICywB00jQKpN9/B159lWNn+HAKhqZNyf6tUIFsVoCM2a++YppGYCC/GzUKmDGDgtHxfGvUYKX8Vasy9l7dIx49YRIXB+zcqQuPHTv4nbs7Cxg3bw40acJ8jPuRP+KM+Hjg1KnklHj19/Jlzo4Kbm56kePChV1T4x3/z5EjeUuPjIDNRoFy5w5p686Fo9U1xMTov8mVi7O1otPXrg0UKZI5AubKFWDsWLa9yJ+fOTXPPmt+L2w2aiEzZ1KwzJ2r0/3j4qhVbNrEwuEXLiC2Xz/Uz5cP52/dwgGLBYVjY4H27dnaZPt2oEcPPofdu6kljR4N7NrFa//rLyA4GPjsMxbABphT1KQJ0wG6d9fPa+BAfnfr1oPP+3JAlqLTZ+YSXLhw8kZXwcEM4a1YkbnkNFeIjaWD9scfWf+iRw9moaoG52pxd9crePXvTxt67lyS106fznxGbkbCbmdIc/t2Fr8eOJBhWcdrzp9f5LHHJLTffAmpEi6hqzLYgb1rF7sjKnqAGa1enbMqotShQ9I6Irdvs4xjzpx02L74ohwBxM/TU5q6uYnV15c+tqAghpXnzeN+VL/q3Ln1ko52O2uY1Kmj7z8+ntv07Zv0nObM4X5SW7PlPgGPjAMWYNvI//2PZKXM6NPiCrGxJKv98AOLLXXvLlKxYtIXSNPIKejcmc66efOYR3P+/KMRDYqOJl3+k08YIq1YUafTu+9kJGTx4ox7ZnY7KQKFCzP68tprKTNqv/ySz6lhw6R8ljNn6LAtXpz/16kjc729BWA7WtE00gUA9giqXZvHvXuXRcEBvU/x9On8fPCgvn9XkZ6//9aFUhbCoyNMjCpXZSSuXCH/YMIERhEqVNBZtipsWLYsi+uMHcvBsG9flmrxmFUQujRCQspdl9BWk1mTVYXemzRhROWff+69fkp4OLUjgAJ+507z7X/8kfyRatWS1pzds4cab3CwyOHDYs+dW57KlUvcANmkaAO1apGs+Msv/Pzmmzx+YCAnEREynr28SDlwPKZKIFRQtVCGDr23689gPDrCJKN4JgrXrtE8mjSJwqFwYUmiaZQtS4r+2LGcBffvzxYa6YXVyr5A48bpBDWAJsZHH9172sGqVaSvu7nxeZkVHF+1ioKjXLmk/aV//52/79JFZPlyuQNIGX9/KeLmJjcDAigQfXy4vndv8llOn9bLMCoS2vPPc7sbN/g5PJwCbPjwpOfRpAlzubIQsoVJanDjBgfRu++S5lysmD6gVUp7nz4c2Js2ZZdizGycO0f6fN26vP8eHjQfV6xIv1kYFsYMbWUOm+X8bNlCLaNs2aQlPhXVfdQokbfekj2AeLi5SZfGjcXu7q4Lwi+/pMDo0YPHzZVL105UOUfH5uidO9OUcjR1xoxxndPzAJEtTJxx+zabtkyZwgFaokRSwVGmDBmQH3wgsn59xlWrz0b6cPAgUxVU3lXRojQh0puwuXw5KfFeXnzpjUyprVtJWixTRq//arfrNVO+/lqkVSv5yN1dAMisTp34fYkSFAxjxvDzhg26dqKcwW3a8BxUTyBVCvL33/Xjb97M7xzLQz5gPNrCxG5note8ebSdVX8UtaheNVOnPtq9ah4GxMTwxWrblmampvHZpaclxLVrfKEBRlKMmM3bt5MpW7q0XvMkLo40eg8Pkd9+E3uhQtLBx0e8vLxkf5UqjPxYLHQ0FyvGdIzbt5nX07gxx+TKlTy2augVF8eIULdu+rGtVgrQp59O+/VlEtIjTB5ensn27cD+/Tq7detW8g8AIGdO1gdt2JB8k+Bg8jIeJsTHk3vgXN81KsqY2WrEevX0NGa5ulrutcF4RuLcOZY+/PRTks86dSJRrV691O/DZgPefReYMAGoVIm8jooVk2+3cyeZq3nzsiZwsWIkDTZoQPbrF1/geu/eqO7hgRwFC2LvpUvwK1IEOHmSDcjffpuEOC8v4OWXSZLr0IHsai8v1pnVNPJdPvmExETVrPyZZ8iivXqVPKkHjEeHZxIQkLTRVYkSDDl+/jlDbVkxk9dq5Yy3eTOduDNnMrz88st04rVpwwhCiRLJK7GldvHwoP0eEMDIQlAQ1fDcubkutfvJm5dhz+7dWepy1ixGtv7558HlLN26RV6OKm3ZqhXNirREgdas4T3x82O5CFfYuZMaR5kyepTn+HEet0IFkXfekbWAaJomz6kWp/nzM5zcpg3v/Zkz/H2VKnzuX30lSfJ1Dh6UxPCywpIl/G7jxvTcnQwHHhnNxM9P9gwYoOfWFC78oE+JGoBiizozXVU1eas1+e9y5TJmubpivfr7J6fHK4p8SoiNNabQq+XOHb1wtFpiY5PuJyhIp9LXqEG2a+3a90f7i4gAvvgC+PBDzuING7KodPPmqfv9xYtAr17UZl99lexV55ys7dtZ2b5ECRaYzpOHf1u1oubi4YE3li3DuyL4rnJl9D55kveob1/ghx+Arl1Zub5nT7J0e/dmHlGxYjyuppEpGxdH7VrTeN/z5iXd//33M/impR2PjmaS0aHh1CIsjGHNr7+mY7BvX4b1VB1S5xm+YEEyNHv1opPuiy8Ypfj3X4ZAHwYSm83GwkDbt5NsNnky64W0bZu0ayLAz716sd7tpk2Zq8VERZEYp5qideuWekdtfDwdvQB76biK2K1dS8dt7dp6lEWVcnzlFYkvWVIaenpKgJ+fnPD357PWNIaDAWpBdevy/KKi9Naia9ZwX+qzI2u3TRvjvsT3GXikHbAZidu36eGfPZuDrnXrpNwTRVorWpTFkp5+mryJr77iYDl27N7qmj5MuH2b1zx5Ml/ookWT3qMqVSS07TQJKXtNQn8Jy/jjR0WRferry5d/7NjUp1N8/jkdqFWrsh2oM/74g+zmxo116v2gQby2SZPkjIeH5HJ3l9olS0oswLBwmTJ09pcvT7NGhYZjYihYGjWiaRYRQXOqZ0/9eErA3Ev5ygxCtjBJK27fJs9g9mzWtGjdOmmNE4A+iFq1qIVMnsxCNseP66G+bCTHlSsMfb71lkj79hLivpN0emwlhf2990QOHMjYroEXLlCoA3yGCxakzne2ejVf6vz5XfcK/uEHCsU2bSgQYmNJqff2Fhk9Wn4GBICMLFdOZ0h365YocKRrV/poLl6k7wmQxJ4fr73GfasG7TdvsqTkK69k0E1JP7KFiRFiY0mTnj2b+TytWlEtdRQavr5UaZ95huHj339nqDkrOnMfMoSutklI1TsS+vS3FMzqnhctypn+jz8yrnn3tm1MsANY6NtVvxpnHDrE/jXe3hQezvjmG11IWK0MNxcvzjHUtau8lCBQVgQG0vGtaey/4+VFweHpyXEVHU1Bp8LG585RAI0YoR+rRw86wB/wZJWlhAmAHwDsT1jOANhvsN0ZsFbs/tRegKkwiYtjFazZs2nbBwfrxYoBzhJ16rAw8LRpHMinT2cLjfuJixdpEj7+uF5j19ub9PSlS+89e9pmYzZ2gQKc+V9/3ZxWL0IB0bAhz2XatOTrZ8zgukGDKAgOHEgcS1FlykhVd3cJypFDLgGMxpUqRY2nWbPEpumyY4fOrF27lvvt3ZvbK7/Nn39y/c8/39s9uEdkKWGS5CDAhwDeMlh3BkDetOwvUZjEx/Ohfv012Yp16uiNzJUN27IlH+aSJdmaRlZETAx9C0OH6lXhCxTgM3NsI5EehIXpDtHKlfWMXrNzefJJSaTRO5thSiio8ouqROMTT8i/Xl7i4+YmLYsUEZsaf61b8++nn/Ka6tWjBla4MP+321lGwTFMbLVSe1HlDB4QsqQwAaABOA+grMH6tAuToCDXvWqaNycv4vvvWdU8I23ybGQ+4uKomXTurPsfGjemlnEvkaE//+QLarEwCmdmQlitupN1wICkETe7nRotwMiciF4Ttm9f+SrB3Hkvb15d46pWjRwVpZHMn8/JD6AwEqFTtmRJ/VhjxvBcL19O/zXfI7KqMGlidmIATgP4C8BeAANTs89gNze9V82iRf+tXjWxsVS5L15khOHkSV7fv/9SC/vrL85m27Yx/LpuHZ2If/5J5/DvvzPPaM8e/u7SJUYOHjbBeukSc6rKluUwDQhgmkR620LcuqULgmrVzCMmdjuFjvKTOEbm4uJYWMnNjSUIrFb2q3Z3F3u7dtITEIubm2xzc6NAKV6chME+fag5FyyoF2IqW5b7U6UMvv+exzh6VAzNrfuE9AiTeyKtaZoWCqCAi1XjROS3hG0+B3BCRD402EchEbmkaVo+AGsADBWRTS62Gwj2G0bxYsWCz5w9m+7zvi+w20nFNioAbfT57t3MOR9NM6bSOxeOLlGCpSN9fTPnXNICEZY8/Ppr0uCjo4G2bUmpb9w47fv7/XeWSrxxA3jrLZZ/NCokPnMmiW0tWgBLl/JeASxp2bIlCWerV7NYdJ06QEQEwgMCUOP0aYi/P/bfuYNcAFM6du1i4ephw4CRI1nCsXNnpgoMHMh9AMDBgyQgNmzIczx8OHNKc6aALFcDVtM0dwAXAQSLyIVUbD8BwF0R+cBsu1RXp88siJB96YrpeuYMc4Ru3+Z2rqBpZIsasV2di0M7Ml1TWuz2lFmurlivly6RkemIfPlcF5BWi49PJt1gA4SH8+X76CPg+nW+kG+8QWZqWuqn3roFDBnC6vGtWwOLFpHV6woLFgD9+7Mi/YoVei7NzZssTn75MouY2+2sNVymDHYeOYJGcXF4PDAQP0ZEQIuPJ/tV01iA+rvvgH37WCn/2DHgxAnm5fTqRQZtz57Jq97fZ2RFYdIO7B3c1GC9HwA3EYlI+H8NgIkistJsv5kuTOx2DhIjarxzMWWAgqB4cb5kBQsa0+GVsHgAs40p7HYKQVfFo9V1Owobi4UJbKqAdJ06/KyKNGcmoqJY/f3990mPr1OHQqVTp9QLFRFqO0OGUED88AO1AVf4/Xe+4GXLAmvX6oLn7FkKEC8vvaVK9+5A06aYtnEjRgP40tMTAy0Wan/XrgHPPw/88gvv1eTJPObEidSQqlbluPj7byYnlizJhETH9hj3CVmOTg9gLoBBTt8ltgcFUArsL3wAwL9IbXvQjCCtqTj/ihWkfz//PLkBpUu7TooLCqLN65j89vvvTNp6FAoo2Wz042zdSj/VuHEkcqnm4wAjafXq6V0QDx3KXF9WTAzrlJQsyeNXrUq/Q1rSFPbt4zO3WFjXxsi3tGYNHf5VqiRtOepY5vHu3USHrC0kRNoA4m2xyEGA3JPKlfl31CjdGdutG4szXb3K5ENA71U8eTI/HziQ7luUXiArOmAzY0mTMLHb6cj8808Olv79OeBVtzZHYVG/PnNLRo8m1XrFCr4Q2d39jGG3M3K2eDHJV02aJO3PHBAgodWHS0jRcxL6+bHMcQTHx/PFrFCBx6xRI23Zt2FhOmv18ceNa92EhpIRXblyUoGybBkdsp060YHeurWIh4dcKVZM8muaVPb3l0hVPzZfPubf1KnD/3ftImX/xRcpBCtUoMCy2ciI9fXlmL3PeLSFic1GWvLy5SLvvy/y7LNMtPL3Tyo08uenBjJkCAXGxo33Xnc0G0lhtbJcwbffigweLCF++xPo9FvIsXjhBTYPz2ghbbVSa1LJfz17us65cQW7nRXlVTsSo6jRunUUKJUqMW1AQVHlhwzheCpeXKRAAVnl5SUAZKCvL8mTSpPr358CaMgQ5n9pGrWkRYskSdh48GD+7j6HiR8dYVK1Kpmr06Yx3FenTtLZUBGfWrYkGeqLLxhGVYV9s3FfERoqElI7VkJHrmClf6UVenmJtGvH7N/0lmZ0hchI9i/y9uYyfnzq6frbtlEYeXvr5oYz1q+nxlCxYtKXfMQIXtdHH1EYeXuLVKwooxL4Jz+q5MeSJWlW9ezJz+vXk0LftCm1rPLlabLZbAwTaxrznO4jHh1h4ig0Chak0Pjf/yg0Nm+mepiNrIvYWEqYYcOo8qtnWakS/QmOvWbuBWfP6ozWIkVoiqXGzHKk1r/7ruvfbNzICax8eb2Iks1GYalppMN/+60IIHFVqkg9QHJ6eMhpR3Z2lSrUlIODyZJVGonynSxYwP126kRhcx87Ijw6wqRoUQqN+9V8KxuZi2PHaGK0bKk3NKtbl87VjKjYvmkTuw0CFBJ//ZXyb6Kj9SzkZ55xnduzaRPN6EqVdFM5Koq+N19faicJRalPFSwoOTRN6ufKJXFVquiadI8e/PvhhyTTlShB869mTZpKMTHUXBQt/z7h0REmD6o4UlaC3f5wFFdKK65do5lQuTKHp68vTdlNm+7NeavKJ+bLRxNj3LiUk//sdr3afKNGrn1rGzbQnKlVS+9ucOUKC2YVKsQk0vr1RXx85HsPDwEgr7/wAv0uBQpQeDZqROGyeDGP9c47ZDUrk8luFwkJYZZ1SuecQUiPMHk4yzY+aNKaGazWpMxXI7brrVtkcxoVgTZaVOFom43HcywWnZZFsV6LF7//5LPUQISs0a+/Br7/nuS6cuWA554D+vUDCrgiXqcCYWFkoc6dS17HvHkkpJnhhx/YGL1QIRaJdi5GvWIF0KUL+S6rVwN+fuSKNGwIlC/P82/YELBY8MLly/gawOohQ9Bq1iw+i6JFyeVp0YKclRUrSGbr359FqE+dYrHrtm1J2hs0KH3XngZkOdJaZuG+C5OwMJ20duECac5GdPiwMOP9WCxJyWu+vknZq2lhulosJDhFRqbMcL17l6Q0I+TPb8xyzQrCJjKSVPqvv2YNVYsF6NiR9VKbN08b+1Xhjz+AF17gs3zjDZLGzAh3O3eS/h4by3Np3Trp+p9/JrGteXPu29ubfzt3Zj3Yl14C2rRBVJEiqH3uHG7nzo0DVasi3/btJAO2bQusWgXMmsXq9U88wb/BwcDrr7O6fkgImcrHj2d6B4FsYZIeiFAQGLFdz5whjdsZOXOmXADa+bscOdI38O8VImSNplQ4Wl1zfHzS3zsLG5WLUqbM/WfyHj0KfPMNNYtr19iGYtw4Us7Tem9v3aJAWriQhbHnzQOqVTPe/uxZsmwPH+a2ffokXT9vHjWYTp0oXDw8SP0fMYLnGBAAjBmDg7lzo86tW2geEoLl//4LN3d3jrGyZZmG0bcv2b3r15Pp++uvFCAHDwLt2wNffsl8nkxElmPAZtaSZtLa1atsYfDDDwwnDx7MQsKVKyfnoagM1apV6UUfMoRktyVLWLXrypWkrR3/a7DZWAJxyxY2Y3/nHabit2rF+qaO7OCcOek0HTOG0YuzZ+9fdnJ0NGumqvauNWrwGaWHcfvrr/SleHjwes2eb3g4Cx5pGkPazlBRmV696Kex23n/FOP18cdFLBb5LKE74Pu9e3Odigy5uXH7kiXp2D18mOc1YAD3Va8erzmTK7HhkXbAxsYypPjDD4zJd+/Oh+Hjk1xYBAZy8D3+OOttTp/OAfXXX4wQPWzp+vcTqiDVnDlkbdaqpUdgFMOzY0eRCRNIILx2LXPPJy6OIVhVqqBiRb60aRX416/rYeTatcnqNUJ0NKvCAbxO5/EyZQrXKQEQG8taO56eLARVtqzYc+SQboC4u7nJrlat9O4GISH8q/YxdSppD25uLEOhKrHNnp3WO5UmPDrCpFIl5mC8+Sbj+hUq6MV0FDGobFkW2RkxgoVpli3jS5DdTzjjER1NzW/WLEZeKlXizJ3wPELz9ZaQfMcldHwmNoK3WhkNqVqVxy1ViqHltEY/fvyRHJCcOVlnxAjx8WRZK9ars0Y0bhzXjR3LzzdvUrPLn5+RGm9vuZU3rxQDpFSBAhJesCA1Ynd3RoHKl2e1NV9ftkXNmZO5UDYbSZolStx7eUsTPDrCxFFolCtHDWPcOFKR9++/r+SebBjgzh2GTd9/X0JyH9bp9B4eNJlmzGCV/4yGzUaqvioqXbgwC0Knxfw5dYraCcACXEYvrc3G9QBruTpuZ7czbQDgZCbCPK+AAGpzc+aIALIlIEAsgPRu2FDsKhReqhR/N3gwNetu3ag9A0wuXb5cMls7eXSEScmS1DIeld40DzlCQ0VC6tsk9IN9zLhWCXkAZ+Dhw1lgOSP9AHY7TYp69STRdNm2LfW/j4nhy6xMj/PnjY/z3nvcrn37pLT9+HiaQ5qmV73/4w9+7tUrkdD2jpubAJBvWrXS70v16tS2X35ZFyLly1Pjjokhd6VQoYyr6u+ER0eYZJPWHn6cPMkZu00bvXtAjhz0dX37bcalRNhspKWrfkhPPWUsGFxh8WI6R/Pm1XsFu8KXX1JINGyY1JSOiiIpzcND75ej/CETJ4rUqSNWLy9pDoivh4ccKl2a2oibm0iePGTFlivHMglLl/J3H3xAEh9AQZYJyBYm2Xg4ERFBB/iAAXo/Iy8vzt5r1mRMTZSICJrCXl40JSZNSr05fPiwXotk/Hhj5vGPP9LnUbdu0jIGt24xDycggE5+u501YTWN5k6uXHIxIEDyAlKtaFGJtlgoUFQXyeee49/x41l/NkcORhU7deL/mZDAmh5hks0zySjEx7smsoWHp47Rmpp16WW7enk9GH5LeiBC1ue8eeR/hIWR29K/P5eiRe9t/6dPA6+9Rh5I8eJsXP7EEynfn8hIllmcP59lIhcv1ks4OuK334AePYDq1cmGDQzk9xcvknQWG8uKbIUKsYbt0aPAtGnA4MH408cHHaOj8XLt2pilxne5cuT/tG7NimtLl5II9+yzrE9brRq5Mh99dG/3xQnZpLWMgN1OAZBSAWjn7yIiUn+M1LJdHbexWMiUdCSeOZeONDueEizOhaPVUqzYg2e6OiMmhoStr79muURNA9q0YenDzp0pJNOL9ev5Mv79N+uyzpxJAWAGERLmXn4ZKFKEDNcKFZJv9/vvLN9YtSoFgBIoR46QVh8YSIESHw/Urs1C348/Dnz4IYZrGqaL4NeKFfH4iRPcJkcOEtpOneI+g4OBGTModGfNotA9epTPMYOQLUxSg7t39TqujqxP9ff6dXPqeWCgOcvV+XPOnEmLQ7u5ZZyWEB/P63FmtqZUOFpdq9WadH8FChjT6h+0sDl9Gvj2Wy4XLlAr6NuXgqVy5fTt02oF5swhnT4sjCzVceNSpqrv2MFcnLg44KefWKneGcuXk0ZfpQoFSu7c+m9btOA5b9hAYdasGbUUux2xmzahoc2GU97eOODmhqJ2O4XNjRssOP3998DHHwOTJlGQLVzI/J/u3Vn8OoOQLUwAahXOdHjH/2/eTLq9pydfFJWHkj+/saAIDDRui/CwwWbTi2a7Ws6dS06rL1KEs2Lt2qTT167N+3K/z3v1amory5bxHOvXB0aPpraSHnr/rVvUUhYsoHYyb17KWsqZMzq1/rPPXNPbV6wAunZlYmBoqH6vfv+dmshjj7G49Pz5TGB86SVg6VKcCAtDzeho1ChWDOvPnYM7wDF6/To1k2PHKABHjmQV+/37Sb//6y+mBWQA7judHkAPsBC0HUBtp3WvAzgB4CiAtga/zw1WpD+e8DcwNccNrliRhKKPPiKDtUsXMlpz5ZJkbFdvb4Yi27YlY3PyZBaf2baNBZL/K827MhpWK6MemzczGjJpEiMh5csnvb+lSpE5+sEHLBgUEXH/zvHaNdYBUbyMypX5bNNbmmHpUpLK3N0ZaUmJFBYezkpxAAs9uTruihV0+lavnrSEgSrzOHQoHbIqBDx+vIibmyxMKFfwVoUKOgHQz48hbk9PVmmrXZtlDM6cIau7bdv0XbcL4H5HcwBUBFAewAZHYQKgElhx3gtASQAnAVhc/H4agDEJ/48BMDU1x01Sac3fn4OoY0c+kGnT6FXfuZM5OdnU+IxHWBh5IVOmkIGs8mMADvxKlciEnTWLzyGza3DExzOPqFIlnkOZMmzBmR7eyo0bJKABJJelVPUtPp4CAWB0xZUwXbWKk1rVqknTCxThbfp0Cq7GjRnFSfi+HyAaIOvz5qUAyZFDEvksAIl/msYJ9f33+d3q1Wm/Zhe478IkcSfJhcnrYL8c9XkVgAYufncUQMGE/wsCOJqa4wWXLs0qVjduZAuLrIKrV0nIGj+egj1fvkQBE2ppIyE5Dkroiz+KHDmSec/MZmPCYa1aPHbRohRo6WFE//QTOxZ4epLLkVKuz6xZJJlVr84WKs5QrTJq1NArBNpsZLeqMo9XrrC8ZPHiIh06SISmSTlACvn5yXUlqAsWpMApXZrbDRhATsqOHdTQKlfOkETU9AiTDPGZaJq2AcBIEdmT8HkWgB0isjDh89cAVojIT06/CxORXA6fb4tIoMExEtuDAqgC4J97PvGsh7wAbjzok8h4VKwA+PoBUZHA4SMP+mwyGP/RZ4byIhKQlh+4p7RBavoJu/qZi+/uSWqJyGwAsxPOaY+k1Tn0EOC/el3Af/fa/svXldbfpChMRKRVOs7lAgBHdlERAJdcbHdV07SCInJZ07SCAK6l41jZyEY2sgAyq0zWMgC9NE3z0jStJICyAHYZbNcv4f9+AIw0nWxkIxtZHPckTDRN66pp2gUADQAs1zRtFQCIyL8AfgRwCMBKAC+LiC3hN3M0TVNq4RQArTVNOw6gdcLn1GD2vZx3FsZ/9bqA/+61ZV9XAh5K0lo2spGNrIf7XA04G9nIxn8V2cIkG9nIRobgoREmmqb10DTtX03T7A4+F7XudU3TTmiadlTTtLYP6hwzApqmTdA07aKmafsTlg4P+pzuBZqmtUt4Lic0TRvzoM8nI6Fp2hlN0w4mPKcsVhMj9dA07RtN065pmvaPw3e5NU1bo2na8YS/LvlfjnhohAlIUusGYJPjl5qmVQLQC0BlAO0AfKZp2sOejTddRGokLH8+6JNJL/7fzv27RhGEYRz/PhhsxC4mhZVFeqs0abSJWMWUtlrY5A+ws0gXsApiIYSk0WATEAtR0tiKIIKoICoaIljYRwKPxczBEXQvp3N3O+v7afbHccsML/fezvK+m+NwB7hMarG4muPVJRdznGquNdkk/Xb63QR2bc8Bu/m4UTXJxPZb2+9/89ESsG37wPYnUnPh/HhHF/5gHvhg+6Ptn8A2KV6hRWw/B34cOb0EbOX9LeDKoOtUk0wanAW+9h3v5XM1W5H0Ot9+Dry9bLEuxqafgaeSXuZ2jy6Ztf0NIG9nBn1hYAXsOLWldH/UmuYJ3AVWSXNYBW4D18Y3uqKqi82QFmzvS5oBnkl6l//l/0utSiYjLt1vjePOU9I94PGIhzNK1cVmGLb38/a7pB3Ssq4ryWToVpcuLHOOW7pfhRy4nmXq7o5+AcxJOifpJOlB+aMJj6kISackne7tA4vUHaujhm51adWdSRNJy8A6cIZUuv/K9iXbbyT1SvcP6Svdr9SapPOk5cBn4MZER/MPbB9KWiG9z+YEsJFbLbpgFthRep/vFHDf9pPJDunvSHoAXACmc3vMLVJry0NJ14EvpLcqNl8nyulDCCV0YZkTQmiBSCYhhCIimYQQiohkEkIoIpJJCKGISCYhhCIimYQQivgFpb6B0ZP6ZZYAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from scipy.linalg import expm, eig\n", "\n", "# Example matrix\n", "A = np.array([\n", " [.6,.8],\n", " [.4,-.2]\n", "], dtype='float')\n", "\n", "# Plot parameters\n", "xlim = 10\n", "ylim = 10\n", "plt.axis('square')\n", "plt.xlim(-xlim, xlim)\n", "plt.ylim(-ylim, ylim)\n", "plt.title(\"Phase Diagram\")\n", "\n", "# Displays the contour with the initial condition\n", "def displayCountour(z0):\n", " z = np.zeros((2, t.size))\n", " for i in range(0, z.shape[1]):\n", " z[:,i] = expm(A * t[i]) @ z0\n", " plt.plot(z[0,:],z[1,:],'r') # show contour in red\n", " plt.plot(z0[0], z0[1], '.b', markersize=3)\n", "\n", "# Plot solutions to differential equation in phase diagram\n", "t = np.linspace(-4,4,num=51)\n", "for x0 in range(-xlim, xlim+1):\n", " displayCountour(np.array([x0,0]))\n", " displayCountour(np.array([0,x0]))\n", "\n", "# Plot eigenvectors (invariant subspaces)\n", "t = np.linspace(-20,20,num=301)\n", "l, T = np.linalg.eig(A)\n", "for i in range(0,l.size):\n", " plt.plot(T[0,i] * t, T[1,i] * t,'k')\n", "\n", "# Show the final plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An Engineering Application\n", "\n", "Note that in the above example, we are reliant on scipy.linalg.expm() to calculate the matrix exponential. Perhaps we would like to avoid such a computationally expensive function in a DSP application. Perhaps importing a library to do these computations would be too large on a small microprocessor. We can utilize invariant subspaces to simplify the computations.\n", "\n", "Assume we can diagonalize $A=T\\Lambda T^{-1}$. Then we can write\n", "$$\n", "e^{At} = I + At + \\frac{1}{2} A^2 t^2 + \\ldots \\\\\n", "e^{At} = TT^{-1} + T \\Lambda T^{-1} t + \\frac{1}{2} T^2 \\Lambda^2 T^{-2} t^2 + \\ldots \\\\\n", "e^{At} = TT^{-1} + T \\Lambda T^{-1} t + \\frac{1}{2} T \\Lambda^2 T^{-1} t^2 + \\ldots \\\\\n", "e^{At} = T ( I + At + \\frac{1}{2} A^2 t^2 + \\ldots ) T^{-1} \\\\\n", "e^{At} = T ( e^{\\Lambda t} ) T^{-1} \\\\\n", "$$\n", "\n", "Now we can substitute that expression into the solution to obtain\n", "$$\n", "x(t) = Te^{\\Lambda t}T^{-1}x_0 \\\\\n", "T^{-1}x(t) = e^{\\Lambda t}T^{-1}x_0 \\\\\n", "z(t) = e^{\\Lambda t} z_0\n", "$$\n", "where $z(t) = T^{-1}x(t)$ and $z_0 = T^{-1}x_0$.\n", "\n", "Now we can calculate everything in terms of $z$ easily because $e^{\\Lambda t}$ is a diagonal matrix. In scalar form, we obtain\n", "$$\n", "z_1(t) = e^{\\lambda_1 t} z_{10} \\\\\n", "\\vdots \\\\\n", "z_n(t) = e^{\\lambda_n t} z_{n0}\n", "$$\n", "\n", "Once we obtain our $z$, we can use $x(t) = Tz(t)$ to transform our coordinate system back into something in terms of $x$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/svg+xml": "\r\n\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2020-12-07T11:17:52.666787\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.3.1, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEICAYAAAB8uBDgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz3klEQVR4nO3dd3hUxfrA8e9LEkIXlN4REBSF0EU6KAgWFERCwkbxh4hXrCDWK1zFgqiIVAEpKSQgICgoinhBAREpCRZEyqX3TkINeX9/nAVCTCDZ7O7Z3cznec6T3T1zZt9NeXNmzswZUVUMwzByK5/dARiGERhMMjEMwy1MMjEMwy1MMjEMwy1MMjEMwy1MMjEMwy1MMvETIrJERPrYHce1iMg3IvKI3XG4SkSWi0j9a5SpKyIrvBWTvzDJxIeIyDYROS0iySKyX0SmiEgRu+O6SERURFKc8R0WkcUi0iN9GVXtpKrT7IoxN0TkPuCkqq67WjlVXQ8cc5Y3nEwy8T33qWoRoAHQGHjd5ngyqueMrxYwFRgtIoM9/aYiEuzp9wD6ATHZLBsHPOHBWPyOSSY+SlV3A98At6Z7uYrzNPykiHwnIiUv7hCRz0Vkn4gcF5EfRaROun2dReRP53G7RWRgun33ikiiiBwTkRUiUjeb8R1S1RjgSeAVEbnBWd+l5piIVBeRH5xnMYdEJE5Eiqd77wYiss4Z1+ciMkNEhjr3tRGRXSLykojsA6aISAkRmS8iB0XkqPNxxXT1LRGRoc7PkSwiX4nIDc73PSEiv4pI1cw+j4jkB9oBS9O9dsxZT7LzjEzTHb8EaC8iodn5fuUFJpn4KBGpBHQG0p9yRwC9gdJAfmBgun3fADWd+9Zi/ee86DPgCVUtipWcfnC+RwNgMtZ/2BuAT4Evc/gHMg8IBppk9jGAd4HywM1AJWCI873zA19gnd1cD8QDD2Y4vqxzXxWgL9bv6xTn88rAaWB0hmPCAQdQAagO/Ow85npgA5DVWVRNIE1Vd118QVWLq2oR55nYSOAnYLdz327gPNYZmgGgqmbzkQ3YBiQDx4DtwFigoHPfEuD1dGX/BSzMop7igALXOZ/vwEoYxTKUGwe8leG1jUDrLOpVoEYmr+8DItPF2SeL4x8A1jkft8L6w5R0+5cBQ52P2wDngAJX+X6FAUfTPV8CvJbu+YfAN+me3wckZlFXc2BfFvt6OH82pTK8vhtoZffvja9s5szE9zyg1n/EKqr6L1U9nW7fvnSPTwFFAEQkSETeE5EtInIC6xcf4GIzqBvWWc52EVkqIs2cr1cBBjhP54+JyDGss4fy2Q1WREKAUsCRTPaVFpEEZ9PqBBCbLqbywG51/lU67cxQxUFVPZOuvkIi8qmIbHfW9yNQXESC0h2zP93j05k8z6pD+yhQNJPPUB/r7OdBVT2YYXdRrMRvYJo5gSIC6ALcCVwHVHW+LgCq+quqdsFqAs0FZjr37wTediavi1shVY3PwXt3AVKBVZnsexfrbKauqhYDel2MCdgLVBARSVe+UobjM05pH4DVrGjqrK9V+s+ZS5sAEZEKF18QkVJYTbH+muEKj4iUx2pqbnTDewcEk0wCQ1HgLHAYKAS8c3GHiOQXkUgRuU5VzwMngAvO3ROBfiLSVCyFReQeEfnHf+iMROR6EYkExgDDVPVwFnElY11GrQC8mG7fz844+otIsIh0IfN+l4z1nXbWdz1Z93/kmPN78z3QGi5dPZoNxKnqjEwOaQP8oKpn3RWDvzPJJDBEY/Wx7Ab+BFZm2O8AtjmbBv2wzhBQ1dXA41in8UeBzcCj13ivJBFJdpbtAzyvqm9kUfY/WJe4jwMLgDkXd6jqOaAr8H9YTYVewHyspJiVj4GCwCHnZ1x4jVhz6lOs7xVARaAl8Fy6KzrJIlLZuT8SGO/m9/drcmWT1TDsIyK/AONVdYqNMSwDns7YrMlQ5jZggqo2y6pMXmSSiWEbEWmN1edwiMv/6W9U1b22Bma4xC3NHBGZLCIHROT3dK9dLyKLRGST82uJLI69W0Q2ishmEXnZHfEYfqMWkITVDBoAPGQSif9yy5mJiLTC6miLVtVbna+9DxxR1fecSaKEqr6U4bgg4G/gLmAX8CvQU1X/zHVQhmF4lVvOTFT1R/45zqALcHHC1zSsAUsZNQE2q+pWZ4dcgvM4wzD8jCcnT5W5eMqqqntFpHQmZSpw5UClXUDTzCoTkb5YQ6opXLhww9q1a7s53Kz9/fffpKSkcPPNN1OgQIF/FkhLg40b4fRpuOkmKOIzE30NwyVr1qw5pKqlcnKMN2ZiXk1mg40ybXep6gRgAkCjRo109erVnozrCnv27KFevXoEBQWxbNmyzBPKwYPQvDkcOABffAE33+y1+AzD3URke06P8eQ4k/0iUg7A+fVAJmV2ceWox4rAHg/G5JLy5cszdepU1q9fz4svvph5oVKl4NtvITQUOnaE3bu9G6Rh2MyTyeRL4OIdtx7Bml2a0a9ATRGp5pxFGu48zufcc889PP/884wePZq5c+dmXqhaNfj6azh2DO6+2/pqGHmFO2YLYk0f34s1JXsX1qjGG4DFWHMeFgPXO8uWB75Od2xnrCs6W0g34/NqW8OGDdUOZ86c0YYNG2qJEiV0x44dWRf8/nvVkBDVVq1UT5/2XoCG4SbAas1hHvDLQWve7jNJb/PmzdSvX5+wsDD++9//EhycRbdTQgL07An33w+zZkFIiHcDNYxcEJE1qtooJ8eYuTk5VKNGDcaNG8eyZct46623si4YHg5jxsCXX0Lv3tYVH8MIYHZfzfFLvXr14vvvv2fo0KG0bduWNm3aZF7wX/+C48fh1VehaFEYOxbEHbPlDcP3mDMTF40ePZoaNWoQGRnJoUOHsi74yivw0kswfrz12DAClEkmLipSpAgJCQkcOnSI3r17c9W+p3ffhX79YNgweO897wVpGF5kkkku1K9fn+HDhzN//nw++eSTrAuKWP0nERHW2cnYsd4L0jC8xCSTXHr66ae57777GDRoEOvWXWXtpnz5YOpUuO8+eOopiI31WoyG4Q0mmeSSiDB58mRKlSpFjx49OHnyZNaFQ0Jg5kxo2xYefdR6bBgBwiQTNyhZsiRxcXFs2bKF/v37X71wgQLW5eJmzaxmz+efeydIw/Awk0zcpHXr1vz73/8mOjqamJhrrDBZpIg17P72262BbbNmeSdIw/Agk0zc6PXXX6dly5Y8+eSTbNq06eqFixaFb76xEkp4OMye7Z0gDcNDTDJxo+DgYOLi4sifPz/h4eGcPXuNVRAuJpSmTa2EMmfO1csbhg8zycTNKlWqxJQpU1i7di2vZGeQ2sWE0rgx9Ohh3QvFMPyQSSYe0KVLF/r378+IESNYsGDBtQ8oVgwWLrQSysMPmzMUwy+ZZOIhw4cPp169ejz66KPszs6NktInlO7dYdq0ax9jGD7EJBMPKVCgADNmzODUqVP06tWLCxcuXPugYsXgu++gXTtrHMrVRtUaho8xycSDatWqxejRo1myZAnvvvtu9g4qUgTmz4cHH4Rnn4U33wQ/vOeMkfeYZOJhjz76KBEREQwZMoRly5Zl76DQUGt07COPwODB8MIL5n4ohs8zycTDRIRx48ZRtWpVIiIiOHIk4/JCWQgOhsmT4Zln4OOPoU8fSE31aKyGkRseTSYiUktEEtNtJ0TkuQxl2ojI8XRl3vBkTHYoVqwYCQkJ7Nu3jz59+lz9dgXp5ctnJZLBg2HKFOvS8enTHo3VMFyW05vGuroBQcA+oEqG19sA83NSl103lM6tDz74QAEdM2ZMzg8eMUIVVJs1Uz1wwO2xGUZ6uHBDaW82c9oDW1Q1x4v7BIrnn3+eTp068cILL7B+/fqcHfzcc9YcnnXrrCH4f/3lkRgNw1XeTCbhWEtiZKaZiCSJyDciUseLMXlVvnz5mDp1KiVKlKBHjx6kpKTkrIJu3WDJEkhOtmYdL1niiTANwyVeSSbOBbbuBzKbb78Wq+lTDxgFzM2ijr4islpEVh88eNBjsXpa6dKliY2NZePGjTz77LM5r6BpU1i5EsqVgw4dzOA2w2d468ykE7BWVfdn3KGqJ1Q12fn4ayBEREpmUm6CqjZS1UalSuVoPWWf0759e1555RU+++wzEhIScl5BtWqwYgW0bGkNbnvjDTMWxbCdt5JJT7Jo4ohIWRFr/QcRaeKM6bCX4rLNkCFDaNasGX379mXr1q05r6B4cWuCYO/e8NZb1iC348fdHqdhZJfHk4mIFALuAuake62fiPRzPn0I+F1EkoBPgHBnb3JACwkJYfr06QQFBdGzZ0/OnTuX80ry54fPPrMuHy9YAI0aQU47dg3DTTyeTFT1lKreoKrH0702XlXHOx+PVtU6qlpPVW9X1RWejslXVK1alUmTJrFq1Spef/111yoRsYbd//e/kJJiXekxN6s2bGBGwNqsW7duPPHEEwwfPpxvv/3W9YpatIC1a61Zxw6HdQd8V852DMNFJpn4gBEjRnDrrbcSFRXFvn37XK+obFn4/nsYMMBam6dVK9ieZ4f1GF5mkokPKFiwIAkJCZw8eRKHw0Fabib1hYTABx9Yd73/80+oW9e6fBz43VCGzUwy8RF16tRh5MiRfP/99wwfPjz3FT70ECQlQb161uXjrl3hwIHc12sYWTDJxIf06dOH7t2789prr7Fy5crcV1itmtUxO3y4tbTGbbfBvHm5r9cwMmGSiQ8RESZMmEDFihXp2bMnx44dy32lQUEwcCCsWQPly8MDD8Bjj5kxKYbbmWTiY4oXL05CQgI7d+6kb9++2b9dwbXceiv88gu8+qrVh3LzzRAXZ/pSDLcxycQH3X777QwdOpTPP/+cSZMmua/i/Pnh7bfh55+hQgXo1QtatzYD3YwrufgPxiQTHzVo0CDuvPNOnn32Wf744w/3Vt6kiTVZcMIE64pPgwbWwDd3NKsM/6MKW7bApEkQGQkVK7pUjUkmPipfvnzExMRQtGhRevTowWl332EtKAgefxz+/hv69oVRo6BWLWt4vrk9ZODbts26e19UFFSpAjVqWL8Pixdb45NckdO7KfnC5q93WnPFwoULFdB+/fp59o3WrLHu4gaq1aurTpmiev68Z9/T8I7z51XXrlUdPVq1Z0/VKlWsnzOoliql2r276tixqn/+qZqWpqqu3WnN9sTgypaXkomq6osvvqiAzpo1y7NvlJam+uWXqg0amKTiz44cUf36a9XXX1dt1061cOHLyaN8edWHHlIdNUr1t98uJY+MXEkmon7Ym9+oUSNdvXq13WF4zblz52jRogWbNm0iMTGRKlWqePYNVa21e4YMseb7VK8Or71mtafz5/fsexs5owqbNln3t1m+3Pr655/WvqAgCAuDO+64vFWqZE0OvQYRWaOqjXISikkmfmLr1q3Ur1+fOnXqsHTpUkJCQjz/phmTSunS1hiVxx+HG2/0/Psb/3TqFKxebSWNi9th5+1/ihe/MnE0bmwt6uYCV5KJ7U0WV7a81sy5KD4+XgF99dVXvfvGaWmqCxeqPvCAalCQdbrcsaPqnDmmCeRJJ06o/vijtTKBw6Fap45qvnyXmyy1aqn27q06caLV33HhgtveGtPMCXx9+vRh8uTJLFq0iPbt23s/gN27rSs+EyfCrl3WqNreva2bXYeFZesU2sjEkSPWygNr117e/v778v5y5axL+A0aWGcczZpByX/c3dRtTDMnD0hJSaFRo0YcO3aMpKQkSpcubU8gqanWbSM//dT6mpYGlStDly7WkP2WLa0ZzMaVVGHvXmugYPrE8b//XS5TufLlxNGwIdSvbyUTLzLJJI9Yv349TZo0oV27dsyfP598+WweLnTwoHXbyLlz4dtv4cwZq/1+771w333WjZvKl7c3Rm9LTYWtW631jTZssLaLj0+cuFyuevXLSaNBAytxePCMI7tMMslDxo4dy1NPPcWHH37ICy+8YHc4l6WkwKJF1uzkr7663DlYubJ1an777dbX+vUD48pQSgps3PjPpLFp05V3uitXDmrXtuZE1a5tzZWqX99Kuj7IJ5OJiGwDTgIXgNSMATrvTD8S6AycAh5V1bVXq9MkE6vjvFu3bsyfP58VK1bQqFHOOt69IjXVmq3888/WtnIl7Nhh7QsNtf4b161rjb6sXt36euONUKiQvXFfdP487Nljxbxzp7VdfHzxa/qF6PPlsz7HxaRxMXHUru2zSSMrvpxMGqnqoSz2dwaexkomTYGRqtr0anWaZGI5cuQIYWFh5M+fn7Vr11KsWDG7Q7q23butpHIxuWzYcOUfJFhNoovJpXx56w+xeHEoUeKfj6+7zhpPcfkaR+bb+fPW3KPMtqNHr3y+d6+VKPbutfqC0itRwhqrUamSdbZVqRLUrGkljJo1rSQZAPw1mXwKLFHVeOfzjUAbVd2bVZ0mmVy2bNkyWrduTXh4OLGxsYg/Xk05etSaaLZ58z+/HjgAFy7kqvr3eZHG/ArArzRmEJncya5IkcuJqkyZy4kifdKoVMnlcRv+xpVkEuypYNJR4DsRUeBTVZ2QYX8FYGe657ucr12RTESkL9AXoHLlyp6L1s+0aNGCIUOG8MYbb3DXXXfx6KOP2h1SzpUoYa35k1lTTdVaWznjGUT6x2lp1iXpLLbGW6vQ7rP3Afhh8FJo0eGfZzfB3vhTCHA5HZiS0w0o7/xaGkgCWmXYvwBoke75YqDh1erMq4PWspKamqpt2rTRQoUK6V9//WV3OD7nhx8ut3d++MHuaPwDLgxa88YiXHucXw8AXwBNMhTZBVRK97wisMfTcQWSoKAgYmNjKViwID169ODMmTN2h+RTfv0VfvjB2n791e5oApdHk4mIFBaRohcfAx2A3zMU+xKIEsvtwHG9Sn+JkbkKFSowbdo0kpKSGDRokN3h+JRBg6BtW2sz3xrP8fSZSRlgmXMd4VXAAlVdmGGt4a+BrcBmYCLwLw/HFLDuuecennvuOUaNGsU8cxd6w8vMoLUAc/bsWZo1a8b27dtJSkqioou34DPyNleu5pjbNgaY0NBQZsyYwblz54iIiCDV3ILR8BKTTAJQzZo1GTt2LD/99BNDhw61OxwjjzDJJEA5HA6ioqJ46623WLp0qd3hGHmASSYBbMyYMVSvXp3IyEgOX5xwZxgeYpJJACtSpAgJCQkcOHCA3r1744+d7Yb/MMkkwDVo0IDhw4fz1VdfMXr0aLvDMQKYSSZ5wDPPPMO9997LwIEDWbdund3hGAHKJJM8QESYMmUKJUuWJDw8nOTkZLtDMnyUqrr8D8ckkzyiZMmSxMXFsWnTJvr37293OIaP2b17N8OHD6du3bo0aNDApTpMMslD2rRpw+uvv860adOIi4uzOxzDZikpKcTGxtKhQwcqVarEoEGDKFq0KGPHjnWpPjOcPo9JTU2lbdu2JCYmsm7dOmrUqGF3SIYXXbhwgSVLlhAdHc3s2bNJSUmhatWqOBwOHA4HNWvWBHz35kiGDwkODmb69OnUq1eP8PBwVqxYQf5AuLGzcVV//PEHMTExxMbGsnv3booVK0bPnj2JioqiefPm7lnhIKc3QPGFzdwcKfe++OILBfSFF16wOxTDQ/bv368ff/yxNmjQQAENCgrSe+65R2fMmKGnTp266rG4cHMk2xODK5tJJu7x1FNPKaALFiywOxTDTU6fPq0zZszQe++9V4OCghTQBg0a6Mcff6z79+/Pdj2uJBPTZ5KHnTlzhqZNm7Jnzx6SkpIon9cWygoQqsqyZcuIiYlh5syZHD9+nAoVKtCrVy8cDgd16tTJcZ2mz8TIkQIFCjBjxgwaNmxIr169WLRoEUFBQXaHZWTT5s2biYmJISYmhv/9738ULlyYbt264XA4aNu2rfd/ljk9lfGFzTRz3Ouzzz5TQIcOHWp3KMY1HD58WMeNG6fNmjVTQEVE77rrLo2OjtaTJ0+67X0wfSaGK9LS0rRnz54aFBSky5YtszscI4OzZ8/q3LlztWvXrpo/f34FtE6dOjps2DDduXOnR97TlWTi0T4TEakERANlgTRggqqOzFCmDTAPuLgM/BxVffNq9Zo+E/c7ceIE9evXJzU1lcTEREqUKGF3SHmaqvLrr78SExNDfHw8hw8fpnTp0kRERBAVFUVYWJhHF1xzpc/Eo2cQQDmggfNxUeBv4JYMZdoA83NSrzkz8YxVq1ZpcHCwdu3aVdPS0uwOJ0/atm2bvv3221qrVi0FNDQ0VHv06KELFizQc+fOeS0OfG3dHFXdq85FyFX1JLABa7U+wwc1btyYd999lzlz5jB+/Hi7w8kzTpw4wZQpU2jbti1Vq1bltddeo0yZMkycOJF9+/aRkJBA586dCQkJsTvUq8tp9nF1A6oCO4BiGV5vAxzGWu3vG6BOFsf3BVYDqytXruz2TGxYLly4oB07dtTQ0FBdv3693eEErPPnz+s333yjPXv21IIFCyqgNWvW1DfffFO3bt1qd3i+2wELFAHWAF0z2VcMKOJ83BnYdK36TDPHs/bv369ly5bVm2++WZOTk+0OJ6AkJibqCy+8oGXLllVAS5QooU8++aT+/PPPPtW09MlkAoQA3wIvZLP8NqDk1cqYZOJ5ixYtUhHRPn362B2K39u9e7cOHz5c69atq4CGhIToAw88oHPmzNEzZ87YHV6mfC6ZAIJ1Nefjq5Qpy+XZy02cTSG5Wr0mmXjHK6+8ooAmJCTYHYrfSU5O1tjYWO3YsaPmy5dPAW3atKmOGTNGDx06ZHd41+RKMvH0peEWwE/Ab1iXhgFeBSoDqOp4EekPPAmkAqedZzArrlavuTTsHefPn6dVq1b8+eefJCYmUq1aNbtD8mlpaWksXbqU6OhoZs2aRXJyMlWqVMHhcNCrVy9q1apld4jZ5sqlYTM3x7iqbdu2ERYWRu3atfnpp598/4qCDTZs2HBpev/OnTspWrQoDz/8MA6Hg5YtW7pner+X+dw4E09tppnjXTNnzlRAX3rpJbtD8RkHDhzQTz75RBs1anRpen+nTp00Pj5eU1JS7A4v1/C1PhNPbSaZeF/fvn0V0G+//dbuUGxz+vRp/fzzz/W+++7T4OBgBTQsLEw/+ugj3bt3r93huZUrycQ0c4xsOXXqFE2aNOHgwYMkJSVRtmxZu0PyClVlxYoVxMTEMGPGDI4dO0a5cuUuTe+/7bbb7A7RI8wtCAyPKVSoEDNmzKBRo0ZERUWxcOFCv+wLyK4tW7YQGxtLTEwMW7ZsoVChQnTt2hWHw0H79u3NrRoyk9NTGV/YTDPHPuPHj1dAhw0bZncobnfkyBH99NNPtXnz5pem97dv316nTZumJ06csDs8r8L0mRielpaWpg899JAGBwfrypUr7Q4n186dO6dffvmlPvTQQxoaGqqA3nzzzfruu+/qjh077A7PNq4kE9NnYuTYsWPHLk2BT0xM5LrrrrM7pBxRVdasWUN0dDTx8fEcOnSIkiVLXpre36BBA49O7/cHps/E8IrixYsTHx9Py5Yt6du3LwkJCX7xx7dz507i4uKIjo5mw4YNhIaGcv/99xMVFUXHjh3NGJrcyumpjC9sppnjG9555x0FdOLEiXaHkqUTJ07o1KlTtV27dioiCmiLFi10woQJevToUbvD81mYZo7hTWlpaXTs2JHly5ezevVqbrnlFrtDAqxV6xYvXkx0dDRz5szh9OnTVK9e/dKw9urVq9sdos8zw+kNr9u7dy/16tWjTJkyrFq1ioIFC9oWy2+//UZ0dDRxcXHs3buX4sWL06NHD6KiomjWrJlfNMV8hekzMbyuXLlyREdH06lTJwYMGODyoteu2rdvH9OnTyc6OpqkpCSCg4Pp3LkzUVFR3HPPPRQoUMCr8eRpOW0X+cJm+kx8z8CBAxXQ2bNne/y9UlJSdPr06dqpU6dL0/sbN26so0aN0gMHDnj8/fMCTJ+JYZdz587RokULNm3aRGJiIlWqVHFr/Wlpafz444/ExMTw+eefc/LkSSpVqoTD4cDhcFC7dm23vl9eZ5o5hm3y589PfHw89evXJyIigqVLlxIcnPtfr40bN15atW7Hjh0UKVKE7t2743A4aN26dUAP6fc7OT2V8YXNNHN81/Tp0xXQ1157zeU6Dh48qKNHj9YmTZoooPny5dO7775b4+LiAmJ6vz/ADKc3fEHv3r1VRHTx4sXZPubMmTM6e/Zs7dKli4aEhCigdevW1Q8++ED37NnjwWiNzLiSTEyfieF2KSkpNGzYkBMnTpCUlESpUqUyLaeqrFy5kpiYGBISEjh69Chly5YlMjISh8NBvXr1vBy5cZFP3mkNuBvYCGwGXs5kvwCfOPevx7kC4NU2c2bi+xITEzU0NFQ7d+6sFy5cuGLf1q1b9c0339QaNWoooAULFtSIiAhduHChnj9/3qaIjfTwtWYOEARsAW4E8mMttJVxedDOWItvCXA78Mu16jXJxD+MHj1aAf3oo4/02LFjOnHiRG3ZsqUCCmjbtm118uTJevz4cbtDNTJwJZl4+u70zYAhqtrR+fwV59nQu+nKfAosUdV45/ONQBtV3ZtVvaaZ4x/OnTtH69atWbVqFUFBQZw/f55atWoRFRVFZGSk2y8fG+7ji5eGKwA70z3fBTTNRpkKwBXJRET6Yi0RSuXKld0eqOEeqsq6desuTe8/cOAA+fLlo1ChQsydO5fWrVubYe0BytMX6TP7rcl4KpSdMqjqBFVtpKqNsurQM+yza9cuhg0bxm233UbDhg0ZN24cLVu2ZN68eSxatIiTJ08yadIku8M0PMjTZya7gErpnlcE9rhQxvBBycnJzJkzh5iYGBYvXoyqcscddzB+/HgefvhhSpQocans4MGDGTx4MHfddRePPPKIjVEbHpPTTpacbFjJaitQjcsdsHUylLmHKztgV12rXtMBa5/U1FRdtGiROhwOLVSokAJarVo1HTx4sG7atOmqx7Vu3VoLFy6sf/31lxcjNlyBCx2wHj0zUdVU5/Kf32Jd2Zmsqn+ISD/n/vHA11hXdDYDp4DenozJcM0ff/xBdHQ0sbGx7Nmzh+uuu+7Scg/Nmze/Zj9IUFAQcXFx1KtXj/DwcFauXEloaKiXoje8IqfZxxc2c2biHfv27dMRI0ZogwYNLq1ad++99+rMmTP19OnTLtX55ZdfKqDPPPOMm6M13AlfG2fiqc0kE885deqUJiQk6D333KNBQUEKaMOGDXXkyJG6f/9+t7zHs88+q4DOmzfPLfUZ7udKMjHD6Q3S0tJYtmwZMTExzJw5kxMnTlCxYsVLzRh3347x7NmzNGvWjO3bt5OUlETFihXdWr+Re744zsTwYZs2bbo0vX/btm0ULlyYhx56CIfDQZs2bTy2al1oaCgJCQk0aNCAyMhIfvjhB7NCXiDI6amML2ymmeO6w4cP69ixY/X222+/NL2/Q4cOGhMTo8nJyV6NZdq0aQrokCFDvPq+xrVh+kyMzJw9e1a/+OILffDBBy9N77/11lv1/fff1127dtkaW69evTRfvny6dOlSW+MwruRKMjF9JgFKVVm1ahXR0dEkJCRw5MgRSpcuTWRkJFFRUdSrV88nhrWfPHmSBg0acObMGRITE7nhhhvsDsnA9JkYwPbt24mNjSU6Opq///6bAgUK8MADDxAVFcVdd93lllspulPRokVJSEigWbNmPPbYY8ydO9cnkpzhgpyeyvjCZpo5Vzp+/Lh+9tln2rp160vT+1u3bq2TJk3SY8eO2R1etowYMUIBHTVqlN2hGGqaOXlKamoqixYtIjo6mrlz53LmzBluuummS6vWVa1a1e4Qc0RVue+++1i0aBG//PILYWFhdoeUp5kV/QKcqpKUlER0dDTTp09n//79XH/99YSHhxMVFUWTJk38uolw8OBBwsLCKFq0KKtXr6ZIkSJ2h5RnmT6TALVnzx7i4uKIiYnht99+IyQkhHvvvZeoqCg6d+5M/vz57Q7RLUqVKkVsbCzt27fn6aefZsqUKXaHZORETttFvrDlhT6T5ORkjY2N1Q4dOlxate7222/XsWPH6qFDh+wOz6Nef/11BTQuLs7uUPIsTJ+Jf7tw4QJLliwhJiaG2bNnk5ycTNWqVS/1g9x00012h+gVqamptGnThvXr17N27Vpq1Khhd0h5jmnm+Kk///yTmJgYYmNj2bVrF8WKFSM8PByHw0GLFi3y3Kp1wcHBTJ8+nbCwMHr27Mny5csDpikX0HJ6KuMLWyA0c/bv368jR47Uhg0bXpre37lzZ01ISNBTp07ZHZ5PmDNnjgI6YMAAu0PJczDD6X3b6dOndebMmXrvvfdqcHCwAlq/fn0dMWKE7tu3z+7wfNK//vUvBfTrr7+2O5Q8xZVkYvpMPExVWb58OTExMcyYMYPjx49Tvnz5S9P7b731VrtD9GmnT5+madOm7Nu3j6SkJMqVK2d3SHmC6TPxIVu2bLk0vX/r1q0UKlSIbt26ERUVRdu2bc2U+2wqWLAgM2bMoFGjRvTq1YvvvvvOfO98VU5PZbK7AcOBv7CW/PwCKJ5FuW3Ab0Ai2Ty18tVmzpEjR3T8+PF6xx13KKAionfeeadOmzZNT548aXd4fm3SpEkK6Ntvv213KHlCdv8W02+eTCYdgGDn42HAsCzKbQNK5qRuX0omZ8+e1Xnz5mm3bt00f/78Cugtt9yi7733nu7cudPu8AJGWlqahoeHa1BQkC5fvtzucAKeK8nEK30mIvIg8JCqRmaybxvQSFUPZbc+u/tMVJXVq1dfWrXu8OHDlCpVioiICKKioqhfv75fD2v3VcePH6d+/fpcuHCBxMTEK9blMdzLlT4Tj52ZpN+Ar4BeWez7H7AWWAP0vUodfYHVwOrKlSu7NQtn1/bt2/Wdd97R2rVrK6ChoaH68MMP6/z58/XcuXO2xJTX/PLLLxocHKzdunXTtLQ0u8MJWHi7mQN8D/yeydYlXZnXsPpMJIs6yju/lsZapKvVtd7Xm82c48eP6+TJk7Vt27YqIgpoy5YtdeLEiXr06FGvxWFc9v777yug48ePtzuUgOX1ZHLNyuER4GegUDbLDwEGXqucp5PJ+fPndeHChRoREaEFCxZUQGvUqKFvvvmmbt261aPvbVzbhQsXtGPHjlqgQAH97bff7A4nIPlUMgHuBv4ESl2lTGGgaLrHK4C7r1W3p5JJUlKSDhgwQMuWLauAlihRQvv166crVqwwp9Q+Zt++fVqmTBm95ZZbNCUlxe5wAo6vJZPNwE6sS76JwHjn6+WBr52Pb3Q2bZKAP4DXslO3O5PJnj179IMPPtB69eopoMHBwdqlSxedPXu2njlzxm3vY7jfd999p4A+/vjjdocScHwqmXhyy20ySUlJ0bi4OL377rsvTe9v0qSJjh49Wg8ePJirug3veumllxTQmTNn2h1KQHElmeSZ4fRpaWksXbqUmJgYZs2axcmTJ6lcuTIOhwOHw0GtWrU8FK3hSefPn6dly5b89ddfrFu3jmrVqtkdUkAww+kz8ddff12a3r9jxw6KFi1K9+7dcTgctGrVKs9N7w80ISEhxMfHX7pdwU8//URISIjdYeVNOT2V8YXtWs2cgwcP6qhRo7Rx48aXVq3r1KmTTp8+3XTWBaiZM2cqoC+//LLdoQQE8nKfyZkzZ3TWrFl6//33X5reX69ePf3www91z549rn9XDb/x+OOPK6Dfffed3aH4PVeSiV/3magqP//8M9HR0cycOZOjR49Srlw5IiMjcTgc1K1b1+5QDS86deoUjRs35vDhwyQlJVGmTBm7Q/JbeWapi9tuu027d+9OdHQ0W7ZsoWDBgnTt2pWoqCjat29vpqjnYb///juNGzemVatWfPPNN6ZPzEV5JpmIiIoIbdu2xeFw0K1bN4oWLWp3WIaPGD9+PE8++STvv/8+L774ot3h+KU8k0wqVqyoK1asoHLlynaHYvggVaV79+7MmzeP5cuX06RJE7tD8jt5JpnYfQsCw/cdPXqUsLAwgoKCWLduHdddd53dIfkVV5KJaVAaAalEiRLEx8ezY8cOnnjiCfzxn6a/McnECFh33HEHb775JjNmzGDy5Ml2hxPwTDIxAtpLL71Eu3btePrpp9mwYYPd4QQ0k0yMgBYUFERMTAxFihQhPDyc06dP2x1SwDLJxAh45cuXZ9q0aaxfv56BAwfaHU7AMsnEyBM6derEgAEDGDt2LF988YXd4QQkk0yMPOOdd96hUaNGPPbYY+zYscPucAKOSSZGnpE/f37i4+NJTU0lIiKC1NRUu0MKKCaZGHlKjRo1+PTTT1m+fDlvvvmm3eEEFI8lExEZIiK7RSTRuXXOotzdIrJRRDaLyMueiscwLoqIiODRRx9l6NCh/Pe//7U7nIDh6TOTEaoa5ty+zrhTRIKAMUAn4Bagp4jc4uGYDINRo0Zx00030atXLw4ePGh3OAHB7mZOE2Czqm5V1XNAAtDF5piMPKBIkSIkJCRw6NAhevfubYbbu4Gnk0l/EVkvIpNFJLOFYStgLYdx0S7na/8gIn1FZLWIrDb/SQx3CAsL44MPPmDBggWMHDnS7nD8Xq6SiYh8LyK/Z7J1AcYB1YEwYC/wYWZVZPJapv8iVHWCqjZS1UalSpXKTdiGcUn//v25//77GTRoEGvWrLE7HL+Wq7vTq+qd2SknIhOB+Zns2gVUSve8IrAnNzEZRk6ICJMnTyYsLIzw8HDWrl1rbrTlIk9ezSmX7umDWAuaZ/QrUFNEqolIfiAc+NJTMRlGZm644Qbi4uLYunUrTz31lN3h+C1P9pm8LyK/ich6oC3wPICIlBeRrwFUNRXoD3wLbABmquofHozJMDLVqlUr3njjDWJiYoiOjrY7HL9k7rRmGE4XLlygXbt2rFmzhrVr13LTTTfZHZJtzJ3WDCMXgoKCiIuLIzQ0lPDwcM6ePWt3SH7FJBPDSKdixYpMmTKFdevW8dJLL9kdjl8xycQwMrj//vt55plnGDlyJF999ZXd4fgNk0wMIxPvv/8+YWFh9O7dm927d9sdjl8wycQwMhEaGkpCQgJnzpwhMjKSCxcu2B2SzzPJxDCyUKtWLcaMGcPSpUt5++237Q7H55lkYhhXERUVRWRkJP/5z3/46aef7A7Hp5lkYhhXISKMGzeOG2+8kYiICI4cOWJ3SD7LJBPDuIaiRYuSkJDA/v37eeyxx8ztCrJgkolhZEPDhg0ZNmwY8+bNY+zYsXaH45NMMjGMbHruuefo3LkzAwYMICkpye5wfI5JJoaRTSLC1KlTuf766+nRowcpKSl2h+RTTDIxjBwoVaoUsbGx/P333zzzzDN2h+NTTDIxjBxq164dr776KpMnTyY+Pt7ucHyGSSaG4YIhQ4bQvHlznnjiCbZs2WJ3OD7BJBPDcEFwcDBxcXEEBQURHh7OuXPn7A7JdiaZGIaLqlSpwmeffcbq1at57bXX7A7HdiaZGEYudO3alSeffJIPPviAhQsX2h2OrTx220YRmQHUcj4tDhxT1bBMym0DTgIXgNTs3CrO3LbR8CWnT5+madOm7Nu3j6SkJMqVK3ftg3ycT922UVV7XFwaFJgNzLlK8bbOsjkK3jB8QcGCBUlISCA5ORmHw0FaWprdIdnC480cERHgYcBcQzMC1i233MInn3zC4sWLGTZsmN3h2MIbfSYtgf2quimL/Qp8JyJrRKSvF+IxDI/4v//7Px5++GH+/e9/8/PPP9sdjtflqs9ERL4Hymay6zVVnecsMw5rcfLMlgdFRMqr6h4RKQ0sAp5W1R8zKdcX6AtQuXLlhtu3b3c5bsPwlOPHjxMWFoaqkpiYSPHixe0OySWu9Jl4dN0cEQkGdgMNVXVXNsoPAZJV9YOrlTMdsIYv++WXX2jRogUPPPAAM2fOxGrp+xef6oB1uhP4K6tEIiKFRaToxcdABzJfRtQw/EbTpk15++23mTVrFhMnTrQ7HK/xdDIJJ0PHa/rlQYEywDIRSQJWAQtUNW9frDcCwsCBA+nQoQPPPvssv/+eN/4/muVBDcND9u/fT7169ShZsiSrVq2iUKFCdoeUbb7YzDGMPKtMmTJER0fzxx9/8Pzzz9sdjseZZGIYHtShQwcGDRrEhAkT+Pzzz+0Ox6NMMjEMDxs6dChNmzbl8ccfZ9u2bXaH4zEmmRiGh4WEhBAfH4+q0rNnT86fP293SB5hkolheEG1atWYMGECK1euZPDgwXaH4xEmmRiGl/To0YM+ffrw3nvv8f3339sdjtuZZGIYXjRy5Ehq166Nw+HgwIEDdofjViaZGIYXFSpUiBkzZnD06FEeeeSRgLpdgUkmhuFlt912GyNGjGDhwoV89NFHdofjNiaZGIYN+vXrR9euXXnllVf49ddf7Q7HLUwyMQwbiAiTJk2ifPnyhIeHc+LECbtDyjWTTAzDJiVKlGD69Ols376dfv364Y/z5NIzycQwbNS8eXP+85//EB8fz9SpU+0OJ1dMMjEMm7388su0bduW/v37s2HDBrvDcZlJJoZhs6CgIGJjYylUqBDh4eGcOXPG7pBcYpKJYfiA8uXLM23aNNavX8/AgQPtDsclJpkYho/o3Lkzzz//PGPGjGHu3Ll2h5NjJpkYhg959913adiwIY899hg7d+60O5wcMcnEMHxIaGgoCQkJnD9/noiICFJTU+0OKdtylUxEpLuI/CEiaSLSKMO+V0Rks4hsFJGOWRx/vYgsEpFNzq8lchOPYQSCGjVqMH78eJYtW8Zbb71ldzjZltszk9+BrsAVi2aJyC1Yd6avA9wNjBWRoEyOfxlYrKo1gcXO54aR50VGRvLII4/w1ltvsWTJErvDyZZcJRNV3aCqGzPZ1QVIUNWzqvo/YDPQJIty05yPpwEP5CYewwgko0ePpmbNmkRGRnLo0CG7w7mmYA/VWwFYme75LudrGZVR1b0AqrrXuURoptIvDwqcFZFAXIykJOD7vzWuCdTP5pXPVapUKU+/RUa1cnrANZNJdtYTzuywTF7L1cQDVZ0ATHDGtDqna3r4g0D9XBC4ny2QP1dOj7lmMlHVO12IZRdQKd3zisCeTMrtF5FyzrOSckBg3XrKMPIQT10a/hIIF5FQEakG1MRa/jOzco84Hz8CZHWmYxiGj8vtpeEHRWQX0AxYICLfAqjqH8BM4E9gIfCUql5wHjMp3WXk94C7RGQTcJfzeXZMyE3cPixQPxcE7mczn8vJL9caNgzD95gRsIZhuIVJJoZhuIXfJJPcDt33FyIyRER2i0iic+tsd0y5ISJ3O38um0UkoEY4i8g2EfnN+XPK8aVUXyEik0XkQPqxW65MdfGbZELuh+77kxGqGubcvrY7GFc5fw5jgE7ALUBP588rkLR1/pz8eazJVKy/nfRyPNXFb5KJG4buG97XBNisqltV9RyQgPXzMnyIqv4IHMnwco6nuvhNMrmKCkD6Gz9kNXTfn/QXkfXO009/nkkdiD+b9BT4TkTWOKd7BJIrproAWU51uchTc3Nc4itD9z3tap8TGAe8hfUZ3gI+BB7zXnRu5Xc/mxxqrqp7nHPKFonIX87/8nmSTyUTDw/d9xnZ/ZwiMhGY7+FwPMnvfjY5oap7nF8PiMgXWM26QEkmOZ7qEgjNnOwO3fcLzh/cRQ9idTz7q1+BmiJSTUTyY3WUf2lzTG4hIoVFpOjFx0AH/PtnlVGOp7r41JnJ1YjIg8AooBTW0P1EVe2oqn+IyMWh+6mkG7rvp94XkTCs5sA24Albo8kFVU0Vkf7At0AQMNk51SIQlAG+EBGw/o6mq+pCe0NyjYjEA22Aks7pMYOxprbMFJH/A3YA3a9ZjxlObxiGOwRCM8cwDB9gkolhGG5hkolhGG5hkolhGG5hkolhGG5hkolhGG5hkolhGG7x/0vWtxmPTSpGAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } } ], "source": [ "# initial condition\n", "x0 = np.array([\n", " [3],\n", " [5]\n", "],dtype='float')\n", "\n", "# transform x initial condition to z coordinates\n", "z0 = np.linalg.inv(T) @ x0 \n", "\n", "# solve in terms of z\n", "# note that we do not need to call scipy.linalg.expm() here. Because we know the eigenvalues and eigenvectors, we can avoid that messy computation.\n", "t = np.linspace(-4,4,num=51)\n", "z = np.zeros((2, t.size))\n", "for i in range(0, t.size):\n", " z[0,i] = np.exp(l[0] * t[i]) * z0[0]\n", " z[1,i] = np.exp(l[1] * t[i]) * z0[1]\n", "\n", "# transform z solution to x coordinates\n", "x = np.zeros((2, t.size))\n", "for i in range(0, t.size):\n", " x[:,i] = T @ z[:,i]\n", "\n", "# plot solution with initial condition\n", "plt.plot(x[0,:], x[1,:], 'r')\n", "plt.plot(x0[0,0], x0[1,0], color='b', marker='x', markersize='3')\n", "# plot eigenvectors for reference\n", "t = np.linspace(-20,20,num=301)\n", "l, T = np.linalg.eig(A)\n", "for i in range(0,l.size):\n", " plt.plot(T[0,i] * t, T[1,i] * t,'k')\n", "\n", "# plot paramters, show plot\n", "xlim = 10\n", "ylim = 10\n", "plt.axis('square')\n", "plt.xlim(-xlim, xlim)\n", "plt.ylim(-ylim, ylim)\n", "plt.title(\"Phase Diagram (z)\")\n", "plt.show()" ] }, { "source": [ "This matches the solution we obtained by using the np.linalg.expm() function call." ], "cell_type": "markdown", "metadata": {} }, { "source": [ "## Challenge Problem: Markov Chains\n", "\n", "An engineer is analyzing data from their popular video-based social media company TikTak to better understand how users interact with their platform. A user can do one of three things: watch a video, view a profile, or create a video. The engineer denotes each of these activies as states which are denoted W, P, and C, respectively. The engineer wants to know the distribution of types https requests (Q, V, M) to model the amount of data that must be processed by the servers to fulfill the users' requests.\n", "\n", "The engineer observes that if the user just made a W request (the user watched a video), there is a 90% the next request will be a W (to watch another), a 5% chance that the next will be a P (view a profile), and a 5% chance the next will be a C (create a video). If the most recent request was a P, there is a 70%, 0%, and 30% chance for a subsequent W, P, or C request, respectively. Likewise, a recent C request means that there is a 60%, 20%, and 20% chance for a subsequent W, P, or C request, respectively. These transition probabilities can be summarized by the transition matrix\n", "$$\n", "M = \\begin{bmatrix}\n", ".90 & .70 & .60 \\\\\n", ".05 & .00 & .20 \\\\\n", ".05 & .30 & .20\n", "\\end{bmatrix}\n", "$$\n", "Such a matrix is a Markov matrix because there are no negative entries and the sum of each column is equal to one.\n", "\n", "The engineer would like to determine the proportions of https requests from the users as a whole.\n", "\n", "1. Write a difference equation that represents the probabilities of the next request in terms of the probabilities of the previous request.\n", "2. In the light of invariant subspaces, determine the steady-state probabilities associated with each of the three requests.\n", "3. What percentage of https requests will be a user wanting to view a profile (P)?\n" ], "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3.8.5 32-bit", "metadata": { "interpreter": { "hash": "8276f2cb448e58cabb0cf7fb45a799f61ea300db3f13689733add06a41a529f5" } } }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "3.8.5-final" } }, "nbformat": 4, "nbformat_minor": 2 }