{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Extended Least Squares\n", "\n", "Example created by Wilson Rocha Lacerda Junior" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we import the NARMAX model, the metric for model evaluation and the methods to generate sample data for tests. Also, we import pandas for specific usage." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pip install sysidentpy" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from sysidentpy.model_structure_selection import FROLS\n", "from sysidentpy.basis_function._basis_function import Polynomial\n", "from sysidentpy.metrics import root_relative_squared_error\n", "from sysidentpy.utils.generate_data import get_siso_data\n", "from sysidentpy.utils.display_results import results\n", "from sysidentpy.utils.plotting import plot_residues_correlation, plot_results\n", "from sysidentpy.residues.residues_correlation import compute_residues_autocorrelation, compute_cross_correlation\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating 1 input 1 output sample data \n", "\n", "The data is generated by simulating the following model:\n", "$y_k = 0.2y_{k-1} + 0.1y_{k-1}x_{k-1} + 0.9x_{k-2} + e_{k}$\n", "\n", "If *colored_noise* is set to True:\n", "\n", "$e_{k} = 0.8\\nu_{k-1} + \\nu_{k}$\n", "\n", "where $x$ is a uniformly distributed random variable and $\\nu$ is a gaussian distributed variable with $\\mu=0$ and $\\sigma$ is defined by the user.\n", "\n", "In the next example we will generate a data with 3000 samples with white noise and selecting 90% of the data to train the model. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "x_train, x_valid, y_train, y_valid = get_siso_data(n=1000,\n", " colored_noise=True,\n", " sigma=0.2,\n", " train_percentage=90)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the model\n", "\n", "First we will train a model without the Extended Least Squares Algorithm for comparison purpose." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [], "source": [ "basis_function = Polynomial(degree=2)\n", "\n", "model = FROLS(\n", " order_selection=False,\n", " n_terms=3,\n", " extended_least_squares=False,\n", " ylag=2, xlag=2,\n", " info_criteria='aic',\n", " estimator='least_squares',\n", " basis_function=basis_function\n", ")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.5649873632814074\n" ] } ], "source": [ "model.fit(X=x_train, y=y_train)\n", "yhat = model.predict(X=x_valid, y=y_valid)\n", "rrse = root_relative_squared_error(y_valid, yhat)\n", "print(rrse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clearly we have something wrong with the obtained model. See the *basic_steps* notebook to compare the results obtained using the same data but without colored noise. But let take a look in whats is wrong." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Regressors Parameters ERR\n", "0 x1(k-2) 8.8975E-01 7.42588820E-01\n", "1 y(k-1) 3.0445E-01 8.42278843E-02\n", "2 y(k-2) -4.7924E-02 2.09392505E-03\n" ] } ], "source": [ "r = pd.DataFrame(\n", " results(\n", " model.final_model, model.theta, model.err,\n", " model.n_terms, err_precision=8, dtype='sci'\n", " ),\n", " columns=['Regressors', 'Parameters', 'ERR'])\n", "print(r)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Biased parameter estimation\n", "\n", "As we can observe above, the model structure is exact the same the one that generate the data. You can se that the ERR ordered the terms in the correct way. And this is an important note regarding the Error Reduction Ratio algorithm used here: __it is very robust to colored noise!!__ \n", "\n", "That is a great feature! However, although the structure is correct, the model *parameters* are not ok! Here we have a biased estimation! The real parameter for $y_{k-1}$ is $0.2$, not $0.3$.\n", "\n", "In this case, we are actually modeling using a NARX model, not a NARMAX. The MA part exists to allow a unbiased estimation of the parameters. To achieve a unbiased estimation of the parameters we have the Extend Least Squares algorithm. Remember, if the data have only white noise, NARX is fine. \n", "\n", "Before applying the Extended Least Squares Algorithm we will run several NARX models to check how different the estimated parameters are from the real ones." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\wilso\\AppData\\Local\\Temp/ipykernel_11404/3592060221.py:12: DeprecationWarning: setting an array element with a sequence. This was supported in some cases where the elements are arrays with a single element. For example np.array([1, np.array([2])], dtype=int). In the future this will raise the same ValueError as np.array([1, [2]], dtype=int).\n", " parameters[:, i] = list(model.theta)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD7CAYAAACBiVhwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApYUlEQVR4nO3de3BU9fk/8Pc5e0mIxFvcFYqUaa0lLSrWG6U64avVAEqwBGYKtKJTL1AtKF/FC6SkVSlUsdGOZaY4VMfbt0WKBRkMWqmMGH/a8ocURKVKUkFIliDCJtns7jmf3x+75+zZ3YQsSc45n+S8X/+Q3cTdx5PsPvt8ns9FEUIIEBERWahuB0BERPJhciAiojxMDkRElIfJgYiI8jA5EBFRHiYHIiLKw+RARER5/G4H0F++/LINut73JRv33/+/AIDf/vZ3fX6s/nL//f+LQMCHRx55zO1QpFdWNhStrVG3w5Aer1NhZLtO/fn+pKoKzjjjlG6/P2iSg66LfkkOLS0t5uPJoqWlBcGgX6qYZMbrVBhep8LIdJ2cfH/isBIREeVhciAiojxMDkRElIfJgYiI8jA5EBFRHiYHIiLKw+QwgAghEN/1d4iYPPOuiWhwYnIYQPTIPnQ2vIDY//s/t0MhokHO1uQQjUYxZcoU7N+/P+v+F154ATfeeKN5e8+ePaiursbEiROxZMkSJJNJO8MasES8PfVv9IjLkRDRYGdbcvjggw8wa9YsNDY2Zt3/n//8B6tXr866b9GiRVi6dCm2bNkCIQTWrl1rV1gDm5pe0K5r7sZBRIOebclh7dq1qK2tRTgcNu+Lx+NYunQpFixYYN534MABxGIxXHTRRQCA6upq1NfX2xXWwCZ0tyMgIo+wbW+lZcuW5d33+OOPY/r06TjnnHPM+1paWhAKhczboVAIzc3NJ/18ZWVDexdojmDQn46jtF8erz8YMZ1WGkQHgEDAJ1V8suG1KQyvU2Fkuk5Ovj85tvHeO++8g4MHD+LBBx/Ee++9Z96v6zoURTFvCyGybheqtTXaL5tRxeOpfkckcrzPj9Vf4vEkgkE/vjqamqWUSOpSxSeTUKiU16YAvE6Fke069ef7k6oqJ/xQ7Vhy2LRpE/bu3YsbbrgB7e3tOHz4MO6++24sWrQIkUjE/LnDhw9nDUVRhmCvgYgc4lhyWL58ufn1e++9h6eeegpPPPEEAKCoqAg7duzAJZdcgg0bNqCiosKpsAYWMzmcfGVFRAOfEEAvBlZ6RYrzHFauXImamhpEo1GMGTMGc+bMcTskOelsSBN52SefH8UpQwKOPJftyWHr1q15940bNw7jxo0zb5eXl2PdunV2hzLwGZWDUx8diEg6bR0JR56HK6QHEKGnFwcyORCRzZgcBhL2HIg8y+njSpkcBhL2HIg8K550drYik8NAwqmsRJ4VTzr74ZDJYQARbEgTeVaSyYG6ZTak+Wsj8hpdsOdA3TE23mPlQOQ5DucGJocBhT0HIs8SluzgRBXB5DCQMDkQeZY1HyQc6D8wOQwgZkPa6fqSiFxnrRY6E/Z/UGRyGEjM5MD1DkReY/1MGGdyoCysHIg8S2RVDhxWIiuz58DkQOQ11lc9KwfKIgQrByKv4rASdY89ByLP4rASdY89ByLPYuVA3TN2ZWVyIPIcTmWlbhmH/Qg2pIk8h5UDdY/DSkSeZe05dMSZHMiKyYHIs6wv+7aY/edI25ocotEopkyZgv379wMA/vKXv2DKlCmoqqrCgw8+iHg8DgDYs2cPqqurMXHiRCxZsgTJZNLOsAYuc5YSkwOR11h7DvH4AJ6t9MEHH2DWrFlobGwEAOzbtw9r1qzBn//8Z2zcuBG6ruOll14CACxatAhLly7Fli1bIITA2rVr7QprQBOcykrkWdZhJc2BI4NtSw5r165FbW0twuEwACAYDKK2thZDhw6Foij49re/jS+++AIHDhxALBbDRRddBACorq5GfX29XWH16MvOr7D36GfQZXwD5rASkWdZX/ZJzf73AL9dD7xs2bKs2yNGjMCIESMAAEeOHMGLL76I5cuXo6WlBaFQyPy5UCiE5ubmk36+srKhfQs47UjsCABgyKkqTi0u7ZfH7KtgMPVr8qtAHIDfpyAUkiM2GfHaFIbXqTCyXKcvjsbMr/0Bn+1x2ZYcutPc3Ixbb70V06dPx7hx47Bjxw4olpPNhBBZtwvV2hqFrvc9mxrZ+YuWI+gc0ueH6xfxeBLBoB/JdI8mmdAQiRx3OSo5hUKlvDYF4HUqjEzX6cuj7ebX0fZ4n+NSVeWEH6odna306aefYubMmZg2bRruvPNOAMCwYcMQiUTMnzl8+LA5FOUGYw1BUpevKW72HCDhkBcR2Sqr56AN4J5Drmg0iltuuQV33XUXfvazn5n3jxgxAkVFRdixYwcAYMOGDaioqHAqrG4lhYSnrrHnQORZ1pe91g+jJD1xbFhp3bp1OHz4MJ555hk888wzAICrr74ad911F1auXImamhpEo1GMGTMGc+bMcSqsbslYOTA5EHmXWTkoQNKBysH25LB161YAwM0334ybb765y58pLy/HunXr7A7lpCRlPK85HRO3zyDyHj2TG/qlv9oTrpDOoSDVDGflQEQyMSoHRVHgQG5gcuiOjD0HweRA5FmWUaWs5rRdmBy6oclYOQiukCbyKmvPQWdycI8u47i+4HkORF5l5gZFgQO7ZzA5dEfO7TOMpMDkQOQ1RrXAYSWXSZkcwMqByKvMlz2HldwlZXIw/iBkjI2IbCUslQOHlVwkW3LIKiNZORB5jrXn4MRaJyaHbsiWHIjI2/SsyoHJwTXyJYfMH4OQLjYislt2z8H+52Ny6IZ0ycH6x8BhJSLPyfQcFM5WcpMTswF6TebYiMgWZuGgcFjJVbp022eIbr4mIi/IXudg//MxOXRDyhXSAKD4OJWVyIOEpXTgOgcX6U5MJD4Zxh+DqnJYiciDstY5MDm4R5f1KE7Vx+RA5EGZdQ4cVnJJ6qprsg7dqBxWIvIi3borKxvSzjO3tpPtDdg86EMFG9JE3pPJDZzK6jjr2ga5Kwe3gyAip2VOguMiOMdZs7ETmflkmHuppCsH2eIjIntZX/IDviEdjUYxZcoU7N+/HwDQ0NCAqqoqVFZWoq6uzvy5PXv2oLq6GhMnTsSSJUuQTLpzCpuWVTnIts4hTfWlv2ByIPIS6xnSYiD3HD744APMmjULjY2NAIBYLIbFixdj1apV2Lx5M3bt2oVt27YBABYtWoSlS5diy5YtEEJg7dq1doV1QtaFb9J9MjfCMZKDbPERka30wTKVde3ataitrUU4HAYA7Ny5E6NGjcLIkSPh9/tRVVWF+vp6HDhwALFYDBdddBEAoLq6GvX19XaFdULWnoN8U1nTfxhq+lfG5EDkKU5PZfXb9cDLli3Lut3S0oJQKGTeDofDaG5uzrs/FAqhubn5pJ+vrGxo74NNOxbLfB0s8iEUKu3zY/aHYNAPPZGqavyBIOIAzjrrFKj+oLuBSUqW35vseJ0KI8t1KilJvd59PhVCCNvjsi055NJ1HYqimLeFEOmDsru+/2S1tkb7PPf3q85j5tftHZ2IRI736fH6SzyehD9dOSTTBc3hyDEo/iIXo5JTKFQqze9NZrxOhZHpOkWjnQBSaxx0gT7HparKCT9UOzZbadiwYYhEIubtSCSCcDicd//hw4fNoSinDYyprBxWIvKi3D6D3X1Rx5LD2LFjsW/fPjQ1NUHTNGzatAkVFRUYMWIEioqKsGPHDgDAhg0bUFFR4VRYWXSJp7IaFIUNaSIvMwZW7H4LcGxYqaioCCtWrMD8+fPR2dmJCRMmYNKkSQCAlStXoqamBtFoFGPGjMGcOXOcCiuLLvNUVuvGe6k7XAuFiJynC6SmKpm3BVSc/BB8oWxPDlu3bjW/Hj9+PDZu3Jj3M+Xl5Vi3bp3dofTIOkNJ1soBavpXJmt8RGQLIQQUZPLDoBlWGgiE1D2HnMqByYHIUzIv+VR6sHsdHJODhbXnIN0Z0oZ0z0G6jQGJyFZmpaDk3LYJk4NF1iI42d58za2V2JAm8iLzILj0v3afR8bkYGHtOci2QlrkDiuxIU3kKblrwITN7wFMDhbWakHahjSnshJ5Uu5L3u63ACYHC6l7DrlTWZkciDzFurcSYP9pcEwOFlL3HNIUbtlN5EmZ0Qwl57Y9mBwshNTJwfjYYAwryRYfEdkpZw0cp7I6KbNfuuLIfum9wmElIk+yHhNqvW0XJgcLo1pQFEW62Uo87IfI2wSQt32GnZgcLIyEkDppSbLkYE5lZXIg8qKcNXCcreQks3KAImFySDEa0kK2yoaI7DVYt+weCMyeg6JIuM7BXCKdvilbfERkJ12kqwaHtuwuKDnMnz8fDQ0N9kYiASFz5WDOYlOzbxORR2S/6KXoOVx77bVYtWoVJk6ciDVr1uDo0aO2BuUWHZnKQbqGtMFc5yBpfERki9R5DgoUmXZlnTp1Kl544QWsWrUKra2tmDFjBhYtWoSdO3faG53DdD11wI+C1NnWMhFsSBN5W972GRJUDgCg6zqamprQ2NgITdNQVlaGX/3qV/j9739vZ3yOyq4cJHvz5a6sRJ5mHvYj0zGhdXV1WL9+PUaOHInZs2fjySefRCAQQHt7O6666iosWLDA3igdIvf2GawciLwsdxjJ7r2VCkoOR44cwdNPP43y8vKs+0tKSvD444/bEpgbrCukpT1Mx2xISxofEdkkZyqrDFt2a5qWlxiMauHKK6/s/6hcIiwrpKU7JjR3hbRsw15EZCunF8GdsHKora1Fc3MzduzYgSNHjpj3J5NJfP75571+0g0bNmD16tUAgIqKCtx///1oaGjA8uXL0dnZicmTJ2PhwoW9fvzeyqyQlnedA3sORN6UnqxkNh3snsp6wuQwY8YM7N27Fx9//DEmTpxo3u/z+XDRRRf16gk7OjqwbNky1NfX49RTT8WsWbOwdetWPPTQQ3j++ecxfPhwzJ07F9u2bcOECRN69Ry9pcm8t5IhPawkX/IiIjvlvubtHtw4YXK44IILcMEFF+CKK67A2Wef3S9PqGkadF1HR0cHSkpKkEwmMXToUIwaNQojR44EAFRVVaG+vt7x5CCydmWVLTnkNKQ5rETkKbnDSq5WDnfddReefPJJ3HrrrV1+/9VXXz3pJxw6dCjuuusuTJ48GUOGDMFll12GlpYWhEIh82fC4TCam5tP6nHLyoaedCy5So4EAAA+nwodAqFQaZ8fsz8Eg34kY6mvTz39FMQAnH5aMYZIEp9sZPm9yY7XqTCyXKdg0A9FUeD3p0YPTjttiK2xnTA53HbbbQCAX/7yl/32hB999BH++te/4h//+AdKS0tx7733orGxMfvg7JyDtAvR2hrt89SuY8c7Us+vp9Z1RCLH+/R4/SUeT5ozB45F4wCAo1+2ITpEjvhkEgqVSvN7kxmvU2Fkuk6xzgQggKSmIwDgyJftfYpNVZUTfqg+4Wyl888/HwBw+eWXY/jw4bj88svR3t6Of/7zn/jOd77Tq4C2b9+O8ePHo6ysDMFgENXV1XjvvfcQiUTMn4lEIgiHw716/L7IDCvJu85BUdiQJvIikT4KLjNbSYKprEuXLsXTTz+NTz/9FDU1Ndi/fz8WL17cqycsLy9HQ0MD2tvbIYTA1q1bMXbsWOzbtw9NTU3QNA2bNm1CRUVFrx6/L6wNaQEhV9OXh/0QeVpeQ1qGFdK7du3CunXrsHr1akybNg333HMPqqure/WEV155JT788ENUV1cjEAjgggsuwPz583HFFVdg/vz56OzsxIQJEzBp0qRePX5fWKeyAqnqwWd8UncdG9JEXmYmg/SQu90fXgtKDkIIqKqKd955B/PmzQMAxGKxXj/p7bffjttvvz3rvvHjx2Pjxo29fsz+kFkVnU4OEJAsNXCFNJHHOTVbqaBhpa9//eu47bbbsH//flx++eW45557MHr0aFsDc4Oet7GVhG/AHFYi8iTj/Slz297nK6hyWL58Od544w1ccsklCAQCuPTSS/GjH/3I3shckGpCK4BlWEkaRrNcNfI5kwORF0nVkC4pKcGll16KY8eOYffu3bjwwgvx2Wef2RqYG3ToltQgWXIwcLYSkSfpQlg/u8rRkH7yySfxpz/9CWVlZeZ9iqLgzTfftC0wN+hCNzYvSd+W6Q1YpGJzqBlFRHLJfclLsWX3hg0b8Prrr/fbFhqy0i3rHADItb+SAADVuZM+iEguRk/UuGnz0xU0rDR8+PBBnxiAVANaAcyrL9ewkgBUxZIcZIqNiOyWKhQy6UGKqazjx4/Ho48+ih/+8IcoLi427x8zZoxtgbnBaEgrMjakAaRyORvSRJ5mfniVIDmsX78eAFBfX2/eNzh7DrlTxSR7A1Y5rETkVSJn2NvVLbsNW7dutTcKSehINaSlrBzyN1ZxMxoicpgQmc+GgCSL4Nra2vDQQw/hpptuwtGjR7F06VK0tbXZGpgbdKPnkCbdIjhFgcJhJSJPcnpvpYKSwyOPPILS0lK0traiqKgI0WgUS5cutTcyF5g9h3R6luscaQFF4bASkVflbK0kxyK4PXv2YOHChfD7/RgyZAhWrlyJPXv22BqYG0ROz0FI9OlcGAfIMjkQeVLeOgcZkoOqZv+Ypml59w0GRo9B2hXSCqeyEnmVgPHh1ZjKau/zFdSQvuyyy/DYY48hFovh7bffxgsvvIBx48bZG5kLNKGnh5QkbEhDpHZkNXdlZeVA5CXmnBSZhpXuvfdelJSUoLS0FE888QTKy8tx33332RqYG3ShIdVzMG5L9AYskE4M6U8NEg15EZH98oeV7H2+HiuHN954A2vWrMHHH3+M4uJijB49GhdffDGKiorsjcwFmjlbSdbKgT0HIq8SRumQ5uoiuNdeew11dXVYsGABysvLoSgK/v3vf2PZsmXo7OxEZWWlrcE5TdM1c6YSIFtyQHpYiT0HIi8yZysZt92sHJ577jk8++yz+NrXvmbed+6552Ls2LFYvHjx4EsOQktvnmEM3Uj2BmzphxCRt5izKWXoObS1tWUlBsM3vvENdHZ22haUW3IrBbl6DoKVA5GHSTWV1efr/gTlwXiegCZSw0pyLoJDOjbOViLyoswOOs5MZXVlscLWrVtRXV2NyZMn45FHHgEANDQ0oKqqCpWVlairq3MjrHRDOjNsI9f2GbmVA5MDkZfkzlB0dcvujz/+GBdffHHe/UIIxOPxXj3h559/jtraWrz88ssoKyvDTTfdhG3btqG2thbPP/88hg8fjrlz52Lbtm2YMGFCr56jtzRdAyDpIjiukCbytvQpoZktu+19uhMmhzfeeKPfn/CNN97Addddh2HDhgEA6urq0NTUhFGjRmHkyJEAgKqqKtTX1zufHISWvQhOtrUEWQ1piRIXEdnOelKlACDcPCZ0xIgR/f6ETU1NCAQCmDdvHg4ePIj/+Z//wXnnnYdQKGT+TDgcRnNz80k9blnZ0L4HpwI+VUUgkOq1DC0NIhQq7fvj9lEw6EccAv5AAGeFTkUUwCklQZwuQWwykuF3NhDwOhVGluvk86lQVQWBgA9xAENK7H1/Kmj7jP6kaRr+9a9/4fnnn0dJSQl+/vOfo7i4OGt9gRAi63YhWlujfT5wO5FMQtcFksnUp/KjX7UjUny8T4/ZH+LxJASApCZwuDUKAIhGY0hE3I9NNqFQKSK8Lj3idSqMTNcpkdShC4FEQoMCIBrt7FNsqqqc8EO148nhrLPOwvjx43HmmWcCAK655hrU19dnzYyKRCIIh8NOh5ZaBGdpScvVcxA5w0qSDXkRkb2EMN+dFEWxfQsdx2crXXXVVdi+fTuOHTsGTdPw9ttvY9KkSdi3bx+ampqgaRo2bdqEiooKp0NL9xwApw7wPmnceI/Is6wveUUBdBmOCe1PY8eOxa233orZs2cjkUjgiiuuwKxZs/DNb34T8+fPR2dnJyZMmIBJkyY5HVrmsJ/0benWOYCzlYi8ylopqKri7lRWu8yYMQMzZszIum/8+PHYuHGjG+GYjO0zDNKtc1C5zoHIq3Irh0G5CE5WqfMcYDbDdZmmiwogNdzF7TOIvEiITMdRURQ5ToLzAiEE9JwV0lLtrZSuHDLrMGSKjYjsJpDJDqoiyTGhXqAJLf1VZr6SVLOVAJh/GU7UlEQkFetLXlUUDis5xTw/Wsm/TwrG9hkAACYHIq+xNqAVxf6GNJNDmlE5KJZxfamSg7HxHpCuHGSKjYjsljeVlZWDMzRz0rBimRAk2adzRcn+l4g8w3pIqMrKwTlm5WDZL12qdQ5CZM5yUBTJptkSkd2EuUuC0XZkcnBEUk8CQPZsJZmmsgKWYSWVPQcij8keVlI4rOSUhJ4AgJwNAOVJDqlpbNaGtDyxEZH9rJWCqti/ZTeTQ1oiXTmo6cpBVVTJ1jkgUzmoqv0bqxCRdLIXwdn7XEwOaUZyMMb1U8lBojdgy1RWRfUBunbinyeiQUXPXecw2HZllVXSGFYyKgcociUH61RW1Q/B5EDkKdnrHLi3kmPiWnbPQVVUCRvSxtp5H5CudIjIGyyTlVLDSuw5OCN3tpIiW8/B+pfBYSUiz7G+HzmxZTeTQ5rZkE6/Aftk6zkg0w9hz4HIe1KVgnWdg73Px+SQlsirHGTrOSCrcmDPgchbsoaVwC27HZO7zkGFKtU6h+yGtA8QTA5EXqLnrnNg5eCMRO5sJel6DshODqwciDzD6C+Y6xxUVg6OSZo9B8s6B6lmK1nGG5kciDzFTASWw35YOTgkoSez9lVSZew5qFznQORFuRsiDPpjQn/729/igQceAAA0NDSgqqoKlZWVqKurczyWhJZAQPWbt2VbIS04rETkWZlhJQ/syvruu+/ilVdeAQDEYjEsXrwYq1atwubNm7Fr1y5s27bN0XgSehIBNWDelm6dg2XjPUX1ARoXwRF5Rf6w0iA9JvTo0aOoq6vDvHnzAAA7d+7EqFGjMHLkSPj9flRVVaG+vt7RmJJ6An5L5eBTZJutlFnnwMqByFtycoMjx4T6e/6R/rd06VIsXLgQBw8eBAC0tLQgFAqZ3w+Hw2hubj6pxywrG9qnmNRPFRQHixAMpi5J0O+HL6giFCrt0+P2h2DQjw4hUHJKMc4MlaK5pBidR3UpYpMRr0theJ0KI8N1irbHAQA+n4pg0I+ioB+6ELbG5nhyePnllzF8+HCMHz8e69evBwDoup5zjoLIul2I1tZon/Yaiba3QxUq4vHUcI2mCcQ644hEjvf6MfuLEVN7RwJa5Dg64wJaIiFFbLIJhUp5XQrA61QYWa7T8XRy0DQd8XgSyaSGhKb3KTZVVU74odrx5LB582ZEIhHccMMN+Oqrr9De3o4DBw7A5/OZPxOJRBAOhx2NK9VzyG5Iy3WGdGYRnOLzs+dA5CHmsFLWVNZBNqz0zDPPmF+vX78e77//Pn7961+jsrISTU1NOOecc7Bp0yZMnz7d0bgSehJ+S0NattlKADKzlXx+CO7KSuQZmckxme197P7s6krPIVdRURFWrFiB+fPno7OzExMmTMCkSZMcjSGpJxBUg+ZtudY55Hxs8AWA9BbjRDT4ddWQtnvLbleTQ3V1NaqrqwEA48ePx8aNG12LJaElUOIvMW+rUKHJtn+RuSsrh5WIvMRMBFwh7byueg7SrHMwK0pL5SB0CJ4jTeQJue9FTkxlZXJI67LnIM3eSunVkZaeAwBA59ASkRfkbbynZJ8pbQcmh7SknkDQZ60cZOo5pFlnKwEcWiLyCD27H51eIc3KwRHxnMpBkWm2Uu48Nl8qTsGmNJEnGD0Hc28lbtntnKSev/GeXOsckN2QBlg5EHlE7nsRG9IOEULkN6Qh4bCSmtNzYHIg8oTc/gIb0g5JpqesBqRdBJcz4GgMK7EhTeQJWnpmonmGtJI/g6m/MTkgNaQEQN7zHIy/AZUNaSIvSmrGjEXjGONBumW3bOLpN9n8qayy9Bxyp7IaDWkmByIv0LR05ZC+PagP+5FJV5WDIuVUVmMem1E5cFiJyAuSxmwlc1hJ4ToHJyTSm9gFfJnKwSfTsJK5txKHlYi8yKwczGEl9hwcYSaHrMpBouRg9qNzhpXYkCbyBE3LXiHNnoNDEuk32ayeA2Rc52DMVmLlQOQlyZyN9ziV1SFdz1ZSpNtbKTOslE5iTA5EnpDMGVZSFNi+ZTeTA7oeVlIVFZosw0oGY1gpHSe3zyDyBg4ruSSTHKwNaR90XZLzHPL2VuKwEpGXJPXcykGBsHmqPZMDUgf9ANmVg1/1Iyk0qfoOSu6wEhvSRJ5gVg5ZK6TtfU4mB2QqB2tD2p9OFEmZToPLqRy4CI7IG3IXwamKAsGeg/2MhnTQZ00OvvT3JHgDFtkNaXBXViJPySyCM7bsZuXgiEzlkD2sBMiRHMy/AWNYSVFSCYINaSJPyCyCS90etIf9PPXUU7j++utx/fXX49FHHwUANDQ0oKqqCpWVlairq3M0nkQXU1n9ikSVA3Ia0gDg80PI0jAnIlsZG+8ZlMF4nkNDQwO2b9+OV155BX/729+we/dubNq0CYsXL8aqVauwefNm7Nq1C9u2bXMspoSehAIFvnRCAKyVg0RvwIpq+ZKVA5FXJHUdfl/mw6GCQXgSXCgUwgMPPIBgMIhAIIBzzz0XjY2NGDVqFEaOHAm/34+qqirU19c7FlMifQqcYvlknmlIS1A55GzZDSDVlJaiqiEiu2magM/y+lfV1HuVnQnC3/OP9K/zzjvP/LqxsRGvvfYafvrTnyIUCpn3h8NhNDc3n9TjlpUN7XVMgf+qCPqDCIVKEQymLknZ6aUAgFNPK0LojNJeP3Z/8PsVxAGcccYpKA6lYukIBFEUUBEKuRubjHhNCsPrVBgZrlMg6EcwoALp96fS0iIAwJlnDkXAb89nfMeTg2Hv3r2YO3cu7rvvPvh8PjQ2NprfE0JkfYovRGtrtNfLyY+1tcMHHyKR44jHU5/G247HAQAtrcdwSvJ4rx63vySTqWbU0a9i8BWlYtGgItbWjkjE3dhkEwqV8poUgNepMLJcp+PRGFRVMd+f4rHUv4eav0JxsHdv46qqnPBDtSsN6R07duDmm2/GPffcg2nTpmHYsGGIRCLm9yORCMLhsGPxGMNKVjLNVsrdWyn1pR+QqR9CRLZJJAUCvszr35fuP+Q2qvuT48nh4MGDuPPOO7Fy5Upcf/31AICxY8di3759aGpqgqZp2LRpEyoqKhyLKaEns85yACTtOSjZPQfurUTkDUlNh9+SHIyvjSmudnB8WGnNmjXo7OzEihUrzPtmzpyJFStWYP78+ejs7MSECRMwadIkx2JKdlE5BGSsHNTMbKpUQ5qVA5EX5CUH1f7KwfHkUFNTg5qami6/t3HjRoejSUnoyaytMwDJprIaMxLU3GElGRIXEdktoekI+C2zKdOJwtiQzw5cIY2uew7GDq1xLe5GSF1SrDFyWInIM5JJHb7B3nOQUUJPZm3XDQBD/MUAgJjW6UZI2XL3VgJSQ0yaBFUNEdkuqWU3pJ3oOTA5wEgO2ZVDkS81jziWjLkRUo78noPiC3BYicgjEnkNaVYOjkhqibyeQ0D1Q1VUOSqHrhrSqp9bdhN5hKZlb59hDDElWTnYK6EnEfRlVw6KoqDYV4ROGZKDOapkna3kY+VA5BEJTWSthDZmK3FYyWYJPb9yAIBifzFiSQmSQ5eL4ALceI/II5LJrtc5JDisZK+ueg4AUOwrkmNYSXS1zsHHLbuJPMKNRXCeTw5CCCS7Sw7+Ikka0mlZySHAk+CIPCKp6V1vn2HjcXCeTw6a0CAg8qayAqkZS1IMK3UxlTV1ngOTA5EXJDQd/q4WwbFysE9Xp8AZSvxD0KF1OB1SFwQAJXunWtUHCA1C2PfHQURySCZFN9tnMDnYJq4Z50fnVw4lgRK0J9xPDkKI7CNCgdSwEsD9lYgGOV0X0EXOIji/0XPgsJJtjO0xinzBvO+dEihBW6IdugyfznOSg+JL9x84tEQ0qCXS1YHf38VspSQrB9sY6xiK/EV53zvFPwQCwv2+gxDIO/ooPQwmuNaBaFAzho6sw0pFgdSHw86EfSMHTA4nqBxKAiUAgPZku6Mx5RPZ+yoBmWElVg5Eg5pxEqR1hbTfp8CnKkwOdjIrB18XlUM6ObQlXE4OQs8fVlI5rETkBcbQkXWFtKIoKAr4EIszOdjmhJWDP105uN2UFgJKd5UDh5WIBrVYujowhpIMRUEfKwc7nbhyGAIAaEu0ORpTHqF3O6wkkvKcN0FE/a+zu+QQ8KGTlYN9YmZy6Gq20ikAgLak+5VD3rBSMJW4RNz9qbZEZJ94nJWDK6LxVFVQ4h+S9z3jvnbXew75DWmlKDXkJeJuN8uJyE6diVTPoSiYnRyK2XOwVzTRhlMCJfCpvrzv+VQfin1FOC7BsJKi5lYOqeQAJgeiQS2WSPUVcyuH4qAPsU77eo5SJYdXX30V1113HSorK/Hiiy868pzReBRDA0O7/X7ZkDPR2tHqSCzdEbpmrmswGMlBdDI5EA1m0fbUFj9DS7J3cThtaBGORu1bg5W/oZBLmpubUVdXh/Xr1yMYDGLmzJkYN24cvvWtb9n6vJGOVpxZfHq33w+XhHDg+Be2xnAiQtcBXcusiDYEhwCKD6LjK3cCIyJHHGuPQ1UUDB2SnRzOLC3CsfYEEkk9a5prf5EmOTQ0NOD73/8+Tj/9dADAxIkTUV9fj1/84hcF/feqmreGuEfH4seREHF856zzzP8+HA5nPd53zvwW9kf3IxI7jLNLQif9HH2lH4sgPPxrCJxxdvb/o+pH8OtjgPYjvfp/H8x4PQrD61QYt69TtCOJ0V8/HX6fmvX+9PXhpQifMQQtRzswMtz96Ed3evr/UoQQ9u3cdBL++Mc/or29HQsXLgQAvPzyy9i5cycefvhhlyMjIvIeaXoOuq5nbUkthMjeopqIiBwjTXIYNmwYIpGIeTsSiZglFBEROUua5PCDH/wA7777Lo4cOYKOjg68/vrrqKiocDssIiJPkqYhffbZZ2PhwoWYM2cOEokEZsyYgQsvvNDtsIiIPEmahjQREclDmmElIiKSB5MDERHlYXIgIqI8TA5ERJSHyUEyPW0++Pe//x033HADpk6dijvuuANffeXNvZUK3aTxrbfewtVXX+1gZHLp6Tp99tlnuPHGGzF16lTccsst/Hvq5jrt3r0b06dPx9SpUzF37lwcO3bMhSgdJkgahw4dEldddZX48ssvRVtbm6iqqhJ79+41v3/8+HFxxRVXiEOHDgkhhHjiiSfEww8/7Fa4runpOhkikYiYNGmSuOqqq1yI0n09XSdd10VlZaXYtm2bEEKIxx57TDz66KNuheuaQv6eZs2aJd566y0hhBDLly8Xv/vd79wI1VGsHCRi3XywpKTE3HzQkEgkUFtbi7PPPhsAMHr0aBw8eNCtcF3T03Uy1NTUFLxx42DU03XavXs3SkpKzMWm8+bNw09+8hO3wnVNIX9Puq6jrS11rktHRweKi4vdCNVRTA4SaWlpQSiU2fk1HA6jubnZvH3GGWfg2muvBQDEYjGsXr0a11xzjeNxuq2n6wQAzz33HL773e9i7NixTocnjZ6u03//+1+cddZZWLx4MaZNm4ba2lqUlJS4EaqrCvl7euCBB1BTU4Mrr7wSDQ0NmDlzptNhOo7JQSKFbj54/Phx3H777SgvL8e0adOcDFEKPV2nTz75BK+//jruuOMON8KTRk/XKZlM4v3338esWbPwyiuvYOTIkVixYoUbobqqp+sUi8WwZMkSPPvss9i+fTtmz56N+++/341QHcXkIJFCNh9saWnB7NmzMXr0aCxbtszpEKXQ03Wqr69HJBLB9OnTcfvtt5vXzGt6uk6hUAijRo3CBRdcAACYMmUKdu7c6XicbuvpOn3yyScoKioyt/P58Y9/jPfff9/xOJ3G5CCRnjYf1DQN8+bNw+TJk7FkyRLPbmne03VasGABtmzZgg0bNmD16tUIh8N46aWXXIzYHT1dp+9973s4cuQIPvroIwDA1q1bMWbMGLfCdU1P12nUqFE4dOgQPvvsMwDAm2++aSbUwUyajfeo+80Hb7vtNixYsACHDh3Chx9+CE3TsGXLFgDA+eef77kKoqfr5IUXbiEKuU5/+MMfUFNTg46ODgwbNgyPPvqo22E7rpDrtHz5ctx9990QQqCsrAy/+c1v3A7bdtx4j4iI8nBYiYiI8jA5EBFRHiYHIiLKw+RARER5mByIiCgPkwMREeVhciAiojxMDkRElOf/A41pv39A+9PEAAAAAElFTkSuQmCC", "text/plain": [ "