{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "----\n", "\n", "This Tutorial will give you an example of how to compute Least Squares and Recursive Least Squares.\n", "\n", "We will be using numpy for matrix computation and matplotlib for visualization." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy.linalg import inv\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model and data\n", "----\n", "\n", "The example of resistance measurement in the slides is an estimation of a single constant. In this code, we will extend the example to estimate a vector $\\mathbf{x}$ with length of 2.\n", "\n", "You will be fitting a linear model between voltage and current: $V = RI + b$. You are given the measurements of voltage $V$ and current $I$, while your goal is to estimate the value of resistance $R$ and offset $b$. If Ohm's law ($V = RI$) holds, you expect this offset $b$ to be near zero.\n", "\n", "The measurements you get is in the following table:\n", "\n", "| Current (A) | Voltage (V) |\n", "|-------------|-------------|\n", "| 0.2 | 1.23 |\n", "| 0.3 | 1.38 |\n", "| 0.4 | 2.06 |\n", "| 0.5 | 2.47 |\n", "| 0.6 | 3.17 |\n", "\n", "How good are these measurements? You may assume that the current values are known exactly, and that the voltage measurements are corrupted by additive, independent and identitically distributed zero-mean Gaussian noise with a standard deviation of $0.15~V$ (i.e., a variance of $0.0225 ~ V^2$). You may also assume that your initial estimates for $\\hat{R}$ and $\\hat{b}$ are uncorelated (i.e., the off-diagonal elements of the $2 \\times 2$ covariance matrix are zero). " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "I = np.array([[0.2, 0.3, 0.4, 0.5, 0.6]]).T\n", "V = np.array([[1.23, 1.38, 2.06, 2.47, 3.17]]).T" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEJCAYAAACOr7BbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeSklEQVR4nO3df5QcZZ3v8feHMCxjBolLYCBDYnAXooiEMBF/BCXjRRJYXYKLR2I2nnVl5+CVHK9KVsLxgq7HFTfqvSpirossugvM6iUJGAMhyoSIiJJJQn4QgzEGzIQrgoIZdpRM+N4/qiLNUN3T1Zmansx8Xuf0SddTz9P1mc6kv6muqqcUEZiZmfV3WL0DmJnZ8OQCYWZmmVwgzMwskwuEmZllcoEwM7NMLhBmZpapsAIh6UhJP5X0kKStkj6V0WeepE3p435JU0vW7ZK0WdJGSeuKymlmZtkOL/C1/wi8LSJ6JDUA90m6MyIeKOnzS+CciPidpPOBrwNvKFnfFhFPFpjRzMzKKKxARHIFXk+62JA+ol+f+0sWHwBOPJhtjh8/PiZPnlzT2GeffZaxY8cezOYL4Vz5OFc+zpXPSMzV1dX1ZEQcm7kyIgp7AGOAjSSF4nMD9L0CuKFk+ZfAeqALaK9me62trVGrzs7OmscWybnyca58nCufkZgLWBdlPlMVQzDVhqRxwDJgQURsyVjfBlwPnB0RT6VtEyJij6TjgNXp2LUZY9uBdoDm5ubWjo6OmjL29PTQ1NRU09giOVc+zpWPc+UzEnO1tbV1RcT0zJXlKsdgP4BrgCsy2k8HfgGcUmHsJ7PG9n94D2LoOFc+zpWPc+VT1B5EkWcxHZvuOSCpETgX+Fm/PpOApcD8iHikpH2spKMOPAfOA16y52FmZsUp8iymE4BvShpDcjrttyNihaTLACJiCXA1cAxwvSSAvkh2dZqBZWnb4cAtEXFXgVnNzKyfIs9i2gRMy2hfUvL8UuDSjD47gan9283MbOj4Smozs0PU8g3dzLj2HjZ3P8OMa+9h+YbuQX39Ir9iMjOzgizf0M2ipZvp3bcfJkL3070sWroZgDnTWgZlG96DMDM7BC1etT0pDiV69+1n8artg7YNFwgzs0PQnqd7c7XXwgXCzOwQNGFcY672WrhAmJkdghbOmkJjw5gXtTU2jGHhrCmDtg0fpDYzOwQdOBCdHHPYS8u4RhbOmjJoB6jBBcLM7JA1Z1oLc6a1sGbNGhbMmznor++vmMzMLJMLhJmZZXKBMDOzTC4QZmaWyQXCzMwyuUCYmVkmFwgzM8vkAmFmZplcIMzMLFOR96Q+UtJPJT0kaaukT2X0kaQvS9ohaZOkM0vWzZa0PV13ZVE5zcwsW5F7EH8E3hYRU4EzgNmS3tivz/nAyemjHfgaQHof66+m608F5ko6tcCsZmbWT2EFIhI96WJD+oh+3S4EvpX2fQAYJ+kE4CxgR0TsjIjngI60r5mZDZFCj0FIGiNpI/AEsDoiftKvSwvwq5Ll3WlbuXYzMxsiiuj/n/oCNiKNA5YBCyJiS0n794DPRsR96fIPgH8EXgXMiohL0/b5wFkRsSDjtdtJvp6iubm5taOjo6aMPT09NDU11TS2SM6Vj3Pl41z5jMRcbW1tXRExPXNlRAzJA7gGuKJf2/8B5pYsbwdOAN4ErCppXwQsGmgbra2tUavOzs6axxbJufJxrnycK5+RmAtYF2U+U4s8i+nYdM8BSY3AucDP+nW7A3hfejbTG4FnIuJx4EHgZEknSToCuCTta2ZmQ6TIGwadAHwzPSPpMODbEbFC0mUAEbEEWAlcAOwA/gt4f7quT9LlwCpgDHBjRGwtMKuZmfVTWIGIiE3AtIz2JSXPA/hQmfErSQqImZnVga+kNjOzTC4QZmaWyQXCzMwyuUCYmVkmFwgzM8vkAmFmZplcIMzMLJMLhJmZZXKBMDOzTC4QZmaWyQXCzMwyuUCYmVkmFwgzM8vkAmFmZplcIMzMLJMLhJmZZSrshkGSJgLfAo4Hnge+HhFf6tdnITCvJMtrgGMj4reSdgF7gf1AX5S7qbaZmRWiyFuO9gEfi4j1ko4CuiStjoiHD3SIiMXAYgBJ7wQ+EhG/LXmNtoh4ssCMZmZWRmFfMUXE4xGxPn2+F9gGtFQYMhe4tag8ZmaWz5Acg5A0meT+1D8ps/5lwGzgtpLmAO6W1CWpvfCQZmb2IoqIYjcgNQH3Ap+JiKVl+rwH+NuIeGdJ24SI2CPpOGA1sCAi1maMbQfaAZqbm1s7OjpqytnT00NTU1NNY4vkXPk4Vz7Olc9IzNXW1tZV9hhvRBT2ABqAVcBHB+i3DHhvhfWfBK4YaHutra1Rq87OzprHFsm58nGufJwrn5GYC1gXZT5TC/uKSZKAbwDbIuKLFfodDZwD3F7SNjY9sI2kscB5wJaispqZ2UsVeRbTDGA+sFnSxrTtKmASQEQsSdsuAu6OiGdLxjYDy5Iaw+HALRFxV4FZzcysn8IKRETcB6iKfjcBN/Vr2wlMLSSYmZlVxVdSm5lZJhcIMzPL5AJhZmaZXCDMzCyTC4SZmWVygTAzs0wuEGZmlskFwszMMrlAmJlZJhcIMzPL5AJhZmaZXCDMzCyTC4SZmWVygTAzs0wuEGZmlskFwszMMrlAmJlZpiLvST1RUqekbZK2SvpwRp+Zkp6RtDF9XF2ybrak7ZJ2SLqyqJxmZpatyHtS9wEfi4j1ko4CuiStjoiH+/X7YUS8o7RB0hjgq8Dbgd3Ag5LuyBhrZmYFKWwPIiIej4j16fO9wDagpcrhZwE7ImJnRDwHdAAXFpPUzMyyDMkxCEmTgWnATzJWv0nSQ5LulPTatK0F+FVJn91UX1zMzGwQKCKK3YDUBNwLfCYilvZb93Lg+YjokXQB8KWIOFnSu4FZEXFp2m8+cFZELMh4/XagHaC5ubm1o6Ojppw9PT00NTXVNLZIzpWPc+Uz3HI93buPXz/zB15xxPP87rnDaD76SMY1NtQ71p8Mt/frgIPJ1dbW1hUR0zNXRkRhD6ABWAV8tMr+u4DxwJuAVSXti4BFA41vbW2NWnV2dtY8tkjOlY9z5TOcci1bvzte/Yk745UfXxFf/o/l8cqPr4hXf+LOWLZ+d72j/clwer9KHUwuYF2U+Uwt8iwmAd8AtkXEF8v0OT7th6SzSL7yegp4EDhZ0kmSjgAuAe4oKquZ1d/iVdvp3bf/RW29+/azeNX2OiWyIs9imgHMBzZL2pi2XQVMAoiIJcDFwAcl9QG9wCVpReuTdDnJ3scY4MaI2FpgVjOrsz1P9+Zqt+INWCAkHUfyYT+B5EN8C8kuyfOVxkXEfYAG6HMdcF2ZdSuBlQPlM7ORYcK4RrozisGEcY11SGNQ4SwmSW2SVgHfA84HTgBOBT5BslfwqfQgs5nZQVs4awqNDWNe1NbYMIaFs6bUKZFV2oO4APiHiHis/wpJhwPvILmQ7baCspnZKDJnWnIme3LMYS8t4xpZOGvKn9pt6FUqEJ+PiF9nrYiIPmB5IYnMbNSaM62FOdNaWLNmDQvmzax3nFGv0llMD0laLenvJR09ZInMzGxYqFQgWoDPA28BHpG0XNJ7JPmIkZnZKFC2QETE/ohYFRHvByYC/wbMAX4p6eYhymdmZnVS1YVykUyY9zDJhHu/JzmbyczMRrCKBULSJEkLJa0HVpBctHZhREwbknRmZlY3Zc9iknQ/yXGI7wDtEbFuyFKZmVndVTrNdRGwNp36wszMRplKXzG9BRhXbqWkt0l6R7n1ZmZ2aKu0B7EZWCHpD8B64DfAkcDJwBnA94F/LjqgmZnVR9kCERG3A7dLOplksr4TSM5g+g+SYxKeYtHMbAQbcDbXiPg58PMhyGJmZsPIkNyT2szMDj0uEGZmlqnIW45OlNQpaZukrZI+nNFnnqRN6eN+SVNL1u2StFnSRkm+BsPMbIgNWCAknSLpB5K2pMunS/pEFa/dB3wsIl4DvBH4kKT+U3T8EjgnIk4HPg18vd/6tog4IyKmV7E9MzMbRNXsQfwryUVz+wAiYhNwyUCDIuLxiFifPt9LMo9TS78+90fE79LFB4ATq49uZmZFqqZAvCwiftqvrS/PRiRNBqYBP6nQ7QPAnSXLAdwtqUtSe57tmZnZwdNAM2lIuhO4HPhORJwp6WLgAxFxflUbkJqAe4HPRMTSMn3agOuBsyPiqbRtQkTskXQcsBpYEBFrM8a2A+0Azc3NrR0dHdXEeomenh6amppqGlsk58rHufJxrnxGYq62trausl/jR0TFB/Aqkqum/wvoBu4DJg80Lh3bAKwCPlqhz+nAL4BTKvT5JHDFQNtrbW2NWnV2dtY8tkjOlY9z5eNc+YzEXMC6KPOZWs2FcjuBcyWNBQ6L5HjCgCQJ+AawLSK+WKbPJGApMD8iHilp/9O20ufnAf9UzXbNzGxwDFggJH203zLAM0BXRGysMHQGMB/YLOlAv6uASQARsQS4GjgGuD593b5IdnWagWVp2+HALRFxV7U/lJmZHbwBCwQwPX18N13+K+BB4DJJ34mIf8kaFBH3Aar0whFxKXBpRvtOYOpLR5iZ2VCppkAcA5wZET0Akq4B/i/wVqALyCwQZmZ2aKvmNNdJwHMly/uAV0Yym+sfC0llZmZ1V80exC3AA5JuT5ffCdyaHjx+uLBkZmZWV9WcxfTp9FqIGSTHFC6LF+5PPa/IcGZmVj/V7EEQEeskPUZyRzkkTYqIxwpNZmZmdVXNZH1/LennJBPr3Zv+eWflUWZmdqir5iD1p0lmY30kIk4CzgV+VGgqMzOru2oKxL5I5kc6TNJhEdEJnFFsLDMzq7dqjkE8nU64txa4WdIT5JzN1czMDj3V7EFcSDJR30eAu0gm1ntHkaHMzKz+qikQV0fE8xHRFxHfjIgvAx8vOpiZmdVXNQXi7RltVd0Lwmy0W76hmxnX3sPm7meYce09LN/QXe9IZlUrewxC0geB/w68StKmklVH4bOYzAa0fEM3i5ZupnfffpgI3U/3smjpZgDmTGsZYLRZ/VU6SH0LyfUOnwWuLGnfGxG/LTSV2QiweNX2pDiU6N23n8WrtrtA2CGhUoEYA/we+FD/FZL+3EXCrLI9T/fmajcbbioViC7gwA2r+9/XIUhuRWpmZUwY10h3RjGYMK6xDmnM8itbINKrps2sRgtnTXnhGESqsWEMC2dNqWMqs+pVcxbTgfmYPp8+qroGQtJESZ2StknaKunDGX0k6cuSdkjaJOnMknWzJW1P113Zf6zZcDdnWguffdfraEn3GFrGNfLZd73Oxx/skFHNPamvBV4P3Jw2fVjSjIhYNMDQPuBjEbFe0lFAl6TVEVF6D4nzgZPTxxuArwFvkDQG+CrJKba7gQcl3dFvrNmwN2daC3OmtbBmzRoWzJtZ7zhmuVQz1cYFwBkR8TyApG8CG4CKBSIiHgceT5/vlbQNaOHFNxm6EPhWRATJTYnGSToBmAzsSO9NjaSOtK8LhJnZEFHy2VyhQ3INxMwDZy1J+nNgTUScXvVGpMkkczmdFhG/L2lfAVwbEfelyz8guUp7MjA7Ii5N2+cDb4iIyzNeux1oB2hubm7t6OioNtaL9PT00NTUVNPYIjlXPs6Vj3PlMxJztbW1dUXE9MyVEVHxAcwFHgVuAr5Jcj+ISwYaVzK+ieSMqHdlrPsecHbJ8g+AVuDdwA0l7fOBrwy0rdbW1qhVZ2dnzWOL5Fz5OFc+zpXPSMwFrIsyn6mVrqS+DrglIm6VtIbkOISAj0fE/6umMklqAG4Dbo6IpRlddgMTS5ZPBPYAR5RpNzOzIVLpLKafA1+QtAv4H8BjEXF7juIg4BvAtoj4YpludwDvS89meiPwTCTHLh4ETpZ0kqQjgEvSvmZmNkQqXQfxJeBLkl5J8gH9b5KOBG4FOiLikQFeewbJV0ObJW1M264CJqWvvwRYSXIQfAfJlOLvT9f1SbocWEVyRfeNEbG1pp/QzMxqMuBZTBHxKPA54HOSpgE3AteQfHBXGncfL70Cu3+fIGMqj3TdSpICYmZmdTDghXKSGiS9U9LNJJP3PQL8TeHJzMysriodpH47yRlMfwX8FOgA2iPi2SHKZmZmdVTpK6arSKb8viI8c6uZ2ahT6SB121AGMTOz4aWqyfrMzGz0cYEwM7NMLhBmZpbJBcLMzDK5QJiZWSYXCDMzy+QCYWZmmVwgzMwskwuEmZllcoEwM7NMLhBmZpbJBcLMzDINeMOgWkm6EXgH8EREnJaxfiEwryTHa4BjI+K36W1O9wL7gb6ImF5UTjMzy1bkHsRNwOxyKyNicUScERFnAIuAe/tNK96WrndxMDOrg8IKRESsBaq9j8Rckntdm5nZMFH3YxCSXkayp3FbSXMAd0vqktRen2RmZqObIqK4F5cmAyuyjkGU9HkP8LcR8c6StgkRsUfSccBqYEG6R5I1vh1oB2hubm7t6OioKWtPTw9NTU01jS2Sc+XjXPk4Vz4jMVdbW1tX2a/yI6KwBzAZ2DJAn2XAeyus/yTJbU8H3F5ra2vUqrOzs+axRXKufJwrH+fKZyTmAtZFmc/Uun7FJOlo4Bzg9pK2sZKOOvAcOA/YUp+EZmajV5Gnud4KzATGS9oNXAM0AETEkrTbRcDdEfFsydBmYJmkA/luiYi7isppZmbZCisQETG3ij43kZwOW9q2E5haTCozM6tW3c9iMjOz4ckFwszMMrlAmJlZJhcIMzPL5AJhZmaZXCDMzCyTC4SZmWVygTAzs0wuEGZmlskFwszMMrlAmJlZJhcIMzPL5AJhZmaZXCDMzCyTC4SZmWVygTAzs0wuEGZmlqmwAiHpRklPSMq8n7SkmZKekbQxfVxdsm62pO2Sdki6sqiMZmZWXpF7EDcBswfo88OIOCN9/BOApDHAV4HzgVOBuZJOLTCnmZllKKxARMRa4Lc1DD0L2BEROyPiOaADuHBQw5mZ2YAUEcW9uDQZWBERp2WsmwncBuwG9gBXRMRWSRcDsyPi0rTffOANEXF5mW20A+0Azc3NrR0dHTVl7enpoampqaaxRXKufJwrH+fKZyTmamtr64qI6ZkrI6KwBzAZ2FJm3cuBpvT5BcDP0+fvBm4o6Tcf+Eo122ttbY1adXZ21jy2SM6Vj3Pl41z5jMRcwLoo85lat7OYIuL3EdGTPl8JNEgaT7JHMbGk64kkexhmZjaE6lYgJB0vSenzs9IsTwEPAidLOknSEcAlwB31ymlmNlodXtQLS7oVmAmMl7QbuAZoAIiIJcDFwAcl9QG9wCXp7k6fpMuBVcAY4MaI2FpUTjMzy1ZYgYiIuQOsvw64rsy6lcDKInKZmVl1fCW1mZllcoEwM7NMLhBmZpbJBcLMzDK5QFguyzd0M+Pae9jc/Qwzrr2H5Ru66x3JzApS2FlMNvIs39DNoqWb6d23HyZC99O9LFq6GYA501rqnM7MBpv3IKxqi1dtT4pDid59+1m8anudEplZkVwgrGp7nu7N1W5mhzYXCKvahHGNudrN7NDmAmFVWzhrCo0NY17U1tgwhoWzptQpkZkVyQeprWoHDkQnxxz20jKukYWzpvgAtdkI5QJhucyZ1sKcaS2sWbOGBfNm1juOmRXIXzGZmVkmFwgzM8vkAmFmZplcIMzMLFNhBULSjZKekLSlzPp5kjalj/slTS1Zt0vSZkkbJa0rKqOZmZVX5B7ETcDsCut/CZwTEacDnwa+3m99W0ScERHTC8pnZmYVFHnL0bWSJldYf3/J4gPAiUVlMTOz/IbLMYgPAHeWLAdwt6QuSe1FbtjTV5uZZVNEFPfiyR7Eiog4rUKfNuB64OyIeCptmxAReyQdB6wGFkTE2jLj24F2gObm5taOjo6q8z3du4/u3/XyfATNjfDrXjhMouUVjYxrbKj6dYrU09NDU1NTvWO8hHPl41z5OFc+B5Orra2tq9xX+XW9klrS6cANwPkHigNAROxJ/3xC0jLgLCCzQETE10mPX0yfPj1mzpxZ9fZnXHsP3U8ncwt97HV9fGFz8na0jBvDj66s/nWKtGbNGvL8TEPFufJxrnycK5+ictXtKyZJk4ClwPyIeKSkfaykow48B84DMs+EOlievtrMrLzC9iAk3QrMBMZL2g1cAzQARMQS4GrgGOB6SQB96W5OM7AsbTscuCUi7ioi44RxjXRnFANPX21mVuxZTHMHWH8pcGlG+05g6ktHDL6Fs6a8cAvNlKevNjNLjOrZXD19tZlZeaO6QICnrzYzK2e4XAdhZmbDjAuEmZllcoEwM7NMLhBmZpbJBcLMzDIVOhfTUJP0G+DRGoePB54cxDiDxbnyca58nCufkZjrlRFxbNaKEVUgDoakdcPx3hPOlY9z5eNc+Yy2XP6KyczMMrlAmJlZJheIF/S/5elw4Vz5OFc+zpXPqMrlYxBmZpbJexBmZpZpxBcISbMlbZe0Q9KVGevnSdqUPu6XNLXasXXMtUvSZkkbJa0b4lwXppk2Slon6exqx9YxV93er5J+r5e0X9LFecfWIVdh71c12STNlPRMuv2Nkq7O+3PVIVddf8fSbBslbZV0b56xFUXEiH0AY4BfAK8CjgAeAk7t1+fNwCvS5+cDP6l2bD1ypcu7gPF1er+aeOGrydOBnw2T9yszV73fr5J+9wArgYuHw/tVLleR71eOv8uZJPeyr+nnGupc9f4dA8YBDwOT0uXjBuv9Gul7EGcBOyJiZ0Q8B3QAF5Z2iIj7I+J36eIDwInVjq1TriJVk6sn0t8+YCwQ1Y6tU64iVfszLwBuA56oYexQ5yrawfzcw+E9G2rV5HovsDQiHgOIiCdyjK1opBeIFuBXJcu707ZyPgDcWePYocoFyYff3ZK6JLUPUqaqc0m6SNLPgO8Bf59nbB1yQR3fL0ktwEXAkrxj65QLinu/qsqWepOkhyTdKem1OccOdS6o77/JU4BXSFqTbv99OcZWNNJvGKSMtsz/WUpqI/kgPvDdddVjhzgXwIyI2CPpOGC1pJ9FxNqhyhURy0juG/5W4NPAudWOrUMuqO/79b+Bj0fEfulF3ev9fpXLBcW9X9VmW08y/UOPpAuA5cDJVY6tRy6o7+/Y4UAr8N+ARuDHkh6ocmxFI30PYjcwsWT5RGBP/06STgduAC6MiKfyjK1DLiJiT/rnE8Aykl3JIctVkmMt8BeSxucdO4S56v1+TQc6JO0CLgaulzSnyrH1yFXk+1VVtoj4fUT0pM9XAg3D4XesQq56/47tBu6KiGcj4klgLTC1yrGVDfZBleH0IKmsO4GTeOEgzWv79ZkE7ADenHdsnXKNBY4qeX4/MHsIc/0lLxwMPhPoJvmfSr3fr3K56vp+9et/Ey8cpK7r+1UhV2HvV46/y+NL/i7PAh4bJr9j5XLV+9/ka4AfpH1fBmwBThuM92tEf8UUEX2SLgdWkRzRvzEitkq6LF2/BLgaOIbkf1AAfRExvdzYeucCmkm+RoHkF+CWiLhrCHP9DfA+SfuAXuA9kfyW1vv9yswlqd7vV66x9c5Fgb9fObJdDHxQUh/J3+Ulw+R3LDNXvX/HImKbpLuATcDzwA0RsQXgYN8vX0ltZmaZRvoxCDMzq5ELhJmZZXKBMDOzTC4QZmaWyQXCzMwyuUDYqCXpeEkdkn4h6WFJKyWdMsQZZkp6c4X1c0pnDU3bHpJ0a7+2z0t6W1E5bXRygbBRSclJ68uANRHxFxFxKnAVyXUA1b7GmErLVZpJMnNvOf8IXF+yjdeQ/Lt9q6SxJf2+Agzq9NdmLhA2WrUB+0ovGouIjRHxw/R/9SsOtEu6TtLfpc93Sbpa0n3AuzOWz5P0Y0nrJX1HUlPJuE+l7ZslvVrSZOAy4CNK5vJ/S2nAdG/mj5FMn3DAe4F/B+4G/rok+6PAMZKOH9R3yUY1FwgbrU4Dumoc+4eIODsiOkqXge8DnwDOjYgzgXXAR0vGPZm2fw24IiJ2kcym+r8i4oyI+GG/7cwgmSCu1HuA/wRuBeb2W7c+HWM2KEb0VBtmBfnPMstvBE4FfpROu3AE8OOSfkvTP7uAd1WxnROA3xxYkPR64DcR8aik3cCNkl4RL9w35AlgQp4fxKwSFwgbrbaSzK2TpY8X710f2W/9s2WWBayOiP7/sz/gj+mf+6nu314vcHTJ8lzg1ekMrAAvJ5mD6oaSnL1VvK5ZVfwVk41W9wB/JukfDjQouT/zOcCjwKmS/kzS0STz7FfjAWCGpL9MX+9lVZwVtRc4qsy6bSSz1CLpMODdwOkRMTkiJpPcHay0GJ1CMpOn2aBwgbBRKZ0d9CLg7elprluBTwJ7IuJXwLdJZse8GdhQ5Wv+Bvg74FZJm0gKxqsHGPZd4KKsg9Qk8/pPS8+4eivQHRHd/dafKukESQ0kxWRdNVnNquHZXM2GMUlfAr4bEd8foN9FwJkR8T+HJpmNBt6DMBve/pnkJjADORz4QsFZbJTxHoSZmWXyHoSZmWVygTAzs0wuEGZmlskFwszMMrlAmJlZJhcIMzPL9P8BLFKmhZXhnaAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(I, V)\n", "plt.xlabel('Current (A)')\n", "plt.ylabel('Voltage (V)')\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 1: Least Squares\n", "----\n", "\n", "Let's first suppose you got all the measurements at once. In this case, you can use least squares to estimate the results.\n", "\n", "As you may find in the slides, the equation to be used for this estimation is:\n", "\n", "\\begin{align}\n", "\\hat{\\mathbf{x}}_{LS} = \\left(\\mathbf{H}^T\\mathbf{H}\\right)^{-1}\\mathbf{H}^T\\mathbf{y}\n", "\\end{align}\n", "\n", "All you need to do is to appropriately define $\\mathbf{x}$, $\\mathbf{H}$ and $\\mathbf{y}$." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slope and offset parameters of the best-fit line (i.e., the resistance and offset) are [R, b]:\n", "Estimates value for R is: 4.970000000000006\n", "Estimates value for b is: 0.07399999999999807\n" ] } ], "source": [ "# Least Squares Estimation\n", "H = np.ones((5, 2))\n", "H[:, 0] = I.ravel()\n", "x_ls = inv(H.T.dot(H)).dot(H.T.dot(V))\n", "print('The slope and offset parameters of the best-fit line (i.e., the resistance and offset) are [R, b]:')\n", "print('Estimates value for R is:', x_ls[0, 0])\n", "print('Estimates value for b is:', x_ls[1, 0])" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuJ0lEQVR4nO3deXxU9fX/8dcBAgQChDVAWMK+oxAWFRdCsSyigEuLtdqqLdXq19+3rRDEvbYVpdqv1oVa69Za05YgKCKokCCiFAEhG4R9S9iXkEASksz5/TFjTWOWIcnNzM2c5+MxD+be+7l33o7JnNy5954rqooxxpjQ1SDQAYwxxgSWFQJjjAlxVgiMMSbEWSEwxpgQZ4XAGGNCXKNAB7hQ7dq105iYmGqte/bsWZo3b167gRzkprxuygruyuumrOCuvG7KCjXLu3HjxuOq2r7charqqkdsbKxWV1JSUrXXDQQ35XVTVlV35XVTVlV35XVTVtWa5QU2aAWfq/bVkDHGhDgrBMYYE+KsEBhjTIizQmCMMSHOCoExxoQ4KwTGGBPirBAYY0yIs0JgjDFBrqjEw0vJO9mdU+LI9h0rBCLSVETWi8gWEUkXkcfLGTNWRHJEZLPv8YhTeYwxxo3SsnKY9uJanl6eyYbDzhQCJ1tMFALjVDVPRMKAz0TkQ1VdV2bcGlWd4mAOY4xxnYKiEv64agcLVu+mdbPGvHzLcMJPZDryWo4VAt8lzXm+yTDfw26HZowxVdiw9ySzE1PYfewsN8V24aFrBtKqWRjJyc4UAlEHb1UpIg2BjUBv4EVVjS+zfCyQCBwEsoH7VTW9nO3MBGYCREVFxSYkJFQrT15eHhEREdVaNxDclNdNWcFded2UFdyVN9iy5hcrC7efZ9X+Yto0FW4f3JjB7b75e70meePi4jaq6ohyF1bUhKg2H0AkkAQMLjO/JRDhez4Z2FHVtqzpXHByU1ZVd+V1U1ZVd+UNpqzJmUf1sidXasycpfrokjTNKyj61hinms7VSRtqVT0tIsnARCCt1PwzpZ4vE5GXRKSdqh6vi1zGGBNop8+d54mlW0ncdJBe7Zvzr59dyoiYNnWawbFCICLtgSJfEQgHxgNPlRnTETiiqioio/CexXTCqUzGGBNMPkw9xMNL0jl17jz3xvXm3nG9aRrWsM5zOLlH0Al403ecoAHwT1VdKiJ3AajqAuBG4G4RKQbygRm+XRhjjKm3jp4p4JEl6SxPP8zg6Ja8ecdIBnVuFbA8Tp41lAIMK2f+glLPXwBecCqDMcYEE1Vl4caDPLE0g4JiD/ET+/PTK3rQqGFgr+113a0qjTHGjQ6cPMfcd1NZs+M4o2LaMO+GIfRsHxxnLFkhMMYYB5V4lLe+2Mv8FZkI8MTUQdwyujsNGkigo/2HFQJjjHHIzqO5xCemsnHfKa7q257fXT+E6MjwQMf6FisExhhTy4pKPPxp9S6eX7mTZk0a8uz3LmL6sGhEgmcvoDQrBMYYU4vSsnKYtTCFrYfOcM3QTjx27SDat2gS6FiVskJgjDG1oKCohP/7ZAd/XrObts0b86dbY5kwqGOgY/nFCoExxtTQ+j0nmZOYwu7jZ/n+iK7MvWYArcLDAh3Lb1YIjDGmmnILinh6eSZ/XbePrm3CefsnoxnTu12gY10wKwTGGFMNSZlHeXBRKofOFHDHmB7cP6EvzRq78yPVnamNMSZATp09zxNLM1j0VRZ9OkSQePdlDO/WOtCxasTuWWyMCTmLv8pizLxVpGblMGbeKhZ/lVXlOqrK0pRsxj+7mve2ZHPfuN4sve9y1xcBsD0CY0yIWfxVFg8sSiW/qAS6QtbpfB5YlArAtGHR5a5z5EwBDy9O46OMIwzt0oq//WQ0Azq1rMvYjrJCYIwJKfNXZHqLQCn5RSXMX5H5rUKgqvxzwwF+88FWzhd7mDu5P3eMCXyTuNpmhcAYE1KyT+f7NX//iXPMWZTC57tOMLpHG566YSgx7ZrXRcQ6Z4XAGBNSOkeGk1VOMejs6wFU4lHe+Hwvv1+RScMGwm+nD+bmkd2CqklcbbNCYIwJKbMm9PvmGIFPeFhDZk3ox/YjucxemMLmA6cZ178Dv50+mE6tgq9JXG2zQmCMCSlfHweYvyITyCU6MpxfjO/D/pPnmLVwCxFNGvHcjIu57qLOQdskrrZZITDGhJxpw6KZNiya5ORkrrzyYuITU9h2OJfrLurMo9cOpG1EcDeJq21WCIwxIamgqIR/ZJ5nxYq1dGjRlFdvG8H4gVGBjhUQVgiMMSHni10neGBRCntPFHHzqG48MLk/LZu6p0lcbXPsZFgRaSoi60Vki4iki8jj5YwREXleRHaKSIqIDHcqjzHGnCkoYu67qdz853V4FGaPbMqT1w8J6SIAzu4RFALjVDVPRMKAz0TkQ1VdV2rMJKCP7zEaeNn3rzHG1KpV244wd1EaR3ML+OkVPfjl1f349+drAh0rKDhWCFRVgTzfZJjvoWWGTQXe8o1dJyKRItJJVQ85lcsYE1pO5BXy66UZLNmcTb+oFiy4NZaLu0YGOlZQEe9nsEMbF2kIbAR6Ay+qanyZ5UuBear6mW96JRCvqhvKjJsJzASIioqKTUhIqFaevLw8IiIiqrVuILgpr5uygrvyuikrBE9eVeXfh0t4O6OQc8Vwba8wpvQMo1GpC8OCJau/apI3Li5uo6qOKHehqjr+ACKBJGBwmfkfAJeXml4JxFa2rdjYWK2upKSkaq8bCG7K66asqu7K66asqsGRN/v0Ob3zjfXaPX6pXvfCZ7rt0JlyxwVD1gtRk7zABq3gc7VOzhpS1dMikgxMBNJKLToIdC013QXIrotMxpj6x+NREr48wJPLtlLk8fDQNQO4fUwPGtbj9hC1wbFCICLtgSJfEQgHxgNPlRn2HnCviCTgPUico3Z8wBhTDXuPn2XOohTW7T7JpT3bMu+GIXRvWz+bxNU2J/cIOgFv+o4TNAD+qapLReQuAFVdACwDJgM7gXPA7Q7mMcbUQyUe5bXP9vDMx5mENWjAvOuH8P2RXUOmPURtcPKsoRRgWDnzF5R6rsA9TmUwxtRvmYdzmb1wC1sO5jB+QAd+M20IHVs1DXQs17Eri40xrnO+2MOLSTt5KXknLZuG8cebhzFlaCfbC6gmKwTGGFf5av8p4hNT2H4kj+nDonl4ykDaNG8c6FiuZoXAGOMK584X88xH23lt7R46tmzKaz8ewbj+odkkrrZZITDGBL3Pdx5nzqJU9p88xw8v6Ub8xP60CPH+QLXJCoExJmjl5Bfx5LKtJHx5gJi2zUiYeQmX9Gwb6Fj1jhUCY0xQ+jjjCA8tTuVYbiE/u6onvxjfl6ZhDQMdq16yQmCMCSrH8wp57L10lqYcon/HFvz5thEM7RIZ6Fj1mhUCY0xQUFUWb87i8fczOFdYwq+u7stdY3sR1tCx26YYHysExpiAyz6dz4PvppKUeYxh3SJ5+oah9IlqEehYIcMKgTEmYDwe5e31+3nqw22UeJRHpgzkR5fFWJO4OmaFwBgTEHuOnyU+MYX1e05yee92PHn9ELq2aRboWCHJCoExpk4Vl3h49bM9/OHj7TRp1ICnbxzKTbFdrD1EAFkhMMbUmYzsM8QnppCalcOEQVE8MXUwHVpak7hAs0JgjHFcYXEJL6zaycvJu4hsFsZLtwxn0uCOthcQJKwQGGMctXGft0nczqN5XD88moevGUhraxIXVKwQGGMccbawmN9/lMkbn++lc6tw3rh9JGP7dQh0LFMOKwTGmFq3ZscxHliUysFT+fzo0u7MmtifiCb2cROs7P+MMabW5Jwr4i+phaxZvp6e7Zvzr7suZWRMm0DHMlWwQmCMqRXL0w7z8JI0TuQV8/OxvbjvO32sSZxLONbEQ0S6ikiSiGwVkXQR+X/ljBkrIjkistn3eMSpPMYYZxzNLeDnb2/krr9tpH1EEx65pCmzJ/a3IuAiTu4RFAO/UtVNItIC2CgiH6tqRplxa1R1ioM5jDEOUFUWbcri10szyC8qYdaEfsy8sidr13wa6GjmAjm2R6Cqh1R1k+95LrAViHbq9YwxdefgqXP86PUv+dW/ttC7QwTL7ruC6Mhwxs5PJjUrhzHzVrH4q6xAxzR+ElV1/kVEYoBPgcGqeqbU/LFAInAQyAbuV9X0ctafCcwEiIqKik1ISKhWjry8PCIiIqq1biC4Ka+bsoK78gZTVo8qq/YXs3D7eRS4qW9jxnVrxJmCYrJO5eNRJSocjuRDAxGiW4cTGR68t5QMpvfWHzXJGxcXt1FVR5S3zPFCICIRwGrgt6q6qMyyloBHVfNEZDLwnKr2qWx7I0aM0A0bNlQrS3JyMmPHjq3WuoHgprxuygruyhssWXcdy2NOYgpf7j3FlX3b87vpg+nS2tskbsy8VWSdzgfgV0OKeSbV+61zdGQ4a+eMC1jmqgTLe+uvmuQVkQoLgaNnDYlIGN6/+N8uWwQASu8dqOoyEXlJRNqp6nEncxlj/FdU4uGVT3fz3ModhIc15Pc3XcQNw6P/qz1Etq8IlFXRfBNcHCsE4v0p+QuwVVWfrWBMR+CIqqqIjMJ7zOKEU5mMMRcmLSuH+MQU0rPPMHlIRx67bhAdWny7SVznyPD/7BGUnW+Cn5N7BGOAW4FUEdnsmzcX6AagqguAG4G7RaQYyAdmaF0ctDDGVKqgqITnV+7gT5/upnWzxiz44XAmDu5U4fhZE/rxwKJU8otK/jMvPKwhsyb0q4u4poYcKwSq+hlQaWtBVX0BeMGpDMaYC7dh70lmJ6aw+9hZbortwkPXDKRVs8oP+E4b5j0hcP6KTCCX6MhwZk3o95/5JrjZlcXGGADyCouZv3wbb63bR3RkOH+9cxRX9Gnv9/rThkUzbVg0ycnJ/M8tY50LamqdFQJjDKu3H2PuolSyc/L50aUxzJrQj+bWJC5k2P9pY0LY6XPn+fXSDBZtyqJX++YsvOtSYrtbk7hQY4XAmBC1LPUQjyxJ4/S5Iu6N682943pbf6AQZYXAmBBz9EwBjyxJZ3n6YQZHt+TNO0YxqHOrQMcyAWSFwJgQoar8a+NBfrM0g8JiD3Mm9ecnl/egUUPHWo4Zl7BCYEwIOHDyHHPfTWXNjuOMimnDvBuG0LO9e3rsGGdZITCmHivxKG99sZenl2fSQOCJaYO5ZVQ3GjSo9BIfE2KqLAQi0gHvVcKd8V79mwZsUFWPw9mMMTWw82gusxemsGn/acb2a89vpw8h2lo+mHJUWAhEJA6YA7QBvgKOAk2BaUAvEVkIPFO6cZwxJvCKSjz8afUunl+5k2ZNGvKH71/EtIv/u0mcMaVVtkcwGfipqu4vu0BEGgFTgKvxdhc1xgSB1IM5zFq4hW2Hc5kytBOPXTeIdhFNAh3LBLnKCsHvVfVIeQtUtRhY7EgiY8wFKygq4Q+fbOfVNXto27wxr9way3cHdQx0LOMSlRWCLSKSCrwDJKpqTh1lMsZcgH/vPsGcRansOX6WGSO78sDkAbQK4ruCmeBTWSGIBsYDM4AnReQLvEXhPVW1u00YE2C5BUU8tXwbf1u3n65twnn7J6MZ07tdoGMZF6qwEKhqCbACWCEijYFJeIvCcyKyUlVvqaOMxpgykrYdZe67qRw+U8Cdl/fgV9/tS7PGdja4qR6/fnJU9byIZABbgVhgoKOpjDHlOnn2PL9+P53Fm7Pp0yGCxLsvY3i31oGOZVyu0kIgIt2A7wM3A82BBGCqqm6tg2zGhLTFX2Uxf0UmM7rmMvfJlVw9MIqlKYfIyS/ivu/04Z64XjRpZE3iTM1Vdh3B53iPE/wLmKmqG+oslTEhbvFXWf+59WNeR8jOKeDNL/bRtXU4b/90NP07tgx0RFOPVLZH8ADwqd1D2Ji6N39F5n/u//vG9m/+6i/xqBUBU+sqazt4BRBZ0UIRGSciUypZ3lVEkkRkq4iki8j/K2eMiMjzIrJTRFJEZPgFpTemnso6/c2Jee1LdYU4lFMQgDSmvqtsjyAVWCoiBcAm4BjeFhN9gIuBT4DfVbJ+MfArVd0kIi2AjSLysapmlBozybe9PsBo4GXfv8aEpBKP8vraPQjw9a7493qU8Gya91e1s/UKMg6o7PTRJcASEemDt+lcJ+AM8De8xwwqvZZAVQ8Bh3zPc0VkK95jDqULwVTgLd/XT+tEJFJEOvnWNSakZB7OZXZiClsOnGZQ55bsOppHQbGHr1sEhYc1ZNaEfoENaeolqYtDACISA3wKDC7dpE5ElgLzVPUz3/RKIL7sgWkRmQnMBIiKiopNSEioVo68vDwiItzTg91Ned2UFYIrb7FHWbq7iPd3FdGsEdwyoAmjOzUkp6CYIzkFtG7s4dT5BkS1akqkC64YDqb3tipuygo1yxsXF7dRVUeUu1BVHX0AEcBG4Ppyln0AXF5qeiUQW9n2YmNjtbqSkpKqvW4guCmvm7KqBk/ezftP6XefXa3d45fqfe9s0uO5Bd8aEyxZ/eWmvG7KqlqzvHhvH1Du56qjlyKKSBje7qRvq+qicoYcBLqWmu4CZDuZyZhgkH++hGc/zuQvn+2hQ4umvHrbCMYPjAp0LBOiHCsE4m1+/hdgq6o+W8Gw94B7RSQB70HiHLXjA6ae+2LXCeYsSmHfiXP8YHQ35kzqT8umwf+Vj6m//LlDWV+8Z/NEqepgERkKXKeqv6li1THArUCqiGz2zZsLdANQ1QXAMrz3PdgJnANur85/hDFucKagiCeXbeOd9fvp3rYZf//paC7rZU3iTOD5s0fwZ2AW8CcAVU0Rkb8DlRYC9R4ArvSWSL7vre7xL6ox7vVJxhEeXJzKsdxCZl7Zk1+M70t4Y2sPYYKDP4WgmaquL3Obu2KH8hhTr5zIK+Tx9zN4b0s2/Tu24JVbR3BR18hAxzLmv/hTCI6LSC9817eIyI34rg8wxpRPVXlvSzaPvZdOXmExvxjfl7vH9qJxo8ou5jcmMPwpBPcArwD9RSQL2AP80NFUxrjYoZx8Hno3jZXbjnJx10ievnEofaNaBDqWMRWqshCo6m5gvIg0Bxqoaq7zsYxxH49HeefL/Ty5bBvFHg8PXTOA28f0oGGDSg+VGRNw/pw19Msy0wA5wEZV3exMLGPcZc/xs8xJTOHfe05yWa+2zLt+KN3aNgt0LGP84s9XQyN8j/d909cAXwJ3ici/VPVpp8IZE+yKSzy8tnYPz3y0ncaNGvDUDUP43oiulDm5wpig5k8haAsMV9U8ABF5FFgIXIm3dYQVAhOSth46Q3xiCikHc7h6YBS/mTaYqJZNAx3LmAvmTyHoBpwvNV0EdFfVfBEpdCaWMcGrsLiEF5N28VLSTlqFh/HCD4ZxzZBOthdgXMufQvB3vC2il/imrwXe8R08zqh4NWPqn037TxG/MIUdR/OYPiyaR6YMpHXzxoGOZUyN+HPW0BMi8iHelhEC3KXftIm+xclwxgSLc+eLeeaj7by2dg8dWzbl9R+PJK5/h0DHMqZW+NV0TlU3iMh+vHcoQ0S6qep+R5MZEyTW7jzOnEUpHDiZz62XdGf2xH60sCZxph7x5/TR64BngM7AUbzHDLYBg5yNZkxg5eQX8bsPtvKPDQfo0a45/5h5CaN7tg10LGNqnT97BE8AlwCfqOowEYkDbnY2ljGB9VH6YR5anMaJs+e566pe/O/4PjQNsyZxpn7ypxAUqeoJEWkgIg1UNUlEnnI8mTEBcCy3kMfeT+eDlEMM6NSSv/xoJEO6tAp0LGMc5U8hOC0iEXjvOfy2iBzFuo+aekZVWbw5i8ffz+BcYQn3f7cvP7uqF2ENrUmcqf/8KQRTgXzgF3jPEmoFPO5kKGPqUtbpfB58N5XkzGMM7+ZtEte7gzWJM6HDn0LwiKrGAx7gTQDfV0PxTgYzxmkej/L2v/cx78NtKPDYtQO59dIYaxJnQo4/heBqvv2hP6mceca4xuGzHma8so71e09yRZ92/G76ELq2sSZxJjRVWAhE5G7g50BPEUkptagFsNbpYMY4objEw5/X7OGZtfk0a1zE/BuHcmNsF2sPYUJaZXsEfwc+BJ4E5pSan6uqJx1NZYwDMrLPMDtxC2lZZ4iNasjLd15FB2sSZ0ylhaAhcIZybi4vIm2qKgYi8howBTiqqoPLWT4WWIL3jmcAi1T11/7FNsZ/BUUlvLBqJwtW7yKyWWNevmU44ScyrQgY41NZIdiI7z7FeHsMlaZAzyq2/QbwAvBWJWPWqOqUKrZjTLVt3HeS2QtT2HXsLDcM78LDUwYQ2awxycmZgY5mTNCosBCoao+abFhVPxWRmJpsw5jqOltYzPwVmbz5xV46twrnzTtGcVXf9oGOZUxQElWtepC339CVvslkVV3q18a9hWBpJV8NJQIHgWzgflVNr2A7M4GZAFFRUbEJCQn+vPy35OXlERERUa11A8FNeYMpa9rxYl5PO8/JAmVct0bc2Lcx4Y3+e6c2mPJWxU1ZwV153ZQVapY3Li5uo6qOKHehqlb6AOYBK4E7fI+PgSerWs+3bgyQVsGylkCE7/lkYIc/24yNjdXqSkpKqva6geCmvMGQ9fTZ8/qrf27W7vFLNe73Sbp+z4kKxwZDXn+5Kauqu/K6KatqzfICG7SCz1V/riOYDFysqh4AEXkT+Ap4oFpl6ZsCdKbU82Ui8pKItFPV4zXZrglNy9MO8fCSdE6ePc/Px/bivu9Ykzhj/OXX/QiASODrs4RqpQOXiHQEjqiqisgooAFwoja2bULH0dwCHl2SzodphxnUuSWv/3gkg6OtSZwxF8KfQvAk8JWIJOE9e+hK/NgbEJF3gLFAOxE5CDwKhAGo6gLgRuBuESnG28tohm/3xZgqqSqJm7J4YmkG+UUlzJ7Yj59e0bPKJnGLv8pi/opMZnTN5cF5q5g1oR/ThkXXUWpjglNlVxa/APxdVd8RkWRgJN5CEK+qh6vasKpWes8CVX0B7+mlxlyQAyfPMffdVNbsOM7ImNbMu2EovdpXfQBt8VdZPLAolfyiEujqbTb3wKJUACsGJqRVtkewA3hGRDoB/wDeUdXNdZLKmHJ4PMpbX+zl6RWZCPDrqYP44ejuNPCzSdz8FZneIlBKflEJ81dkWiEwIa2y6wieA54Tke7ADOB1EWkKvAMkqOr2OspoDDuP5jEnMYUN+05xZd/2/G76YLq0vrAmcdmn8y9ovjGhospjBKq6D3gKeEpEhgGv4f2+307JMI4rKvHwyqe7ee6THYQ3bsgzN13E9cOjq9UkrnNkOFnlfOh3jgyvjajGuFaVt18SkTARuVZE3sbbhG47cIPjyUzIS8vKYeoLa5m/IpOrB0bxyS+v4oYadAqdNaEf4WVOKQ0Pa8isCf1qI64xrlXZweKr8d6k/hpgPZAAzFTVs3WUzYSogqISnlu5g1c+3U2b5o1Z8MNYJg7uWOPtfn0cYP6KTCCX6MhwO2vIGCr/amgu3lbU96u1nTZ15Mu9J4lfmMLu42f53oguPDh5IK2ahdXa9qcNi2basGiSk5P5n1vG1tp2jXGzyg4Wx9VlEONOtXVefl5hMU8v38ZbX+yjS+tw/nbnaC7v086BxMaYsvy9stiYb6mt8/KTM4/y4LtpZOfkc/uYGO7/bj+aN7EfTWPqiv22mWqr6Xn5p86e54kPMli0KYveHSJYeNdlxHZv7VRcY0wFrBCYaqvuefmqyrLUwzz6XhqnzxVx37je3DOuN00a2RnJxgSCFQJTbdU5L//omQIeWpzGRxlHGBLdirfuGM3Azi2djGmMqUKV1xEYU5ELOS9fVfnnlwf4zrOrWb39GA9M6s+7P7/MioAxQcD2CEy1+Xte/oGT53hgUSqf7TzOqB5tmHf9EHr60STOGFM3rBCYGqnsvPwSj/Lm53uZvyKThg2E30wbzA9GdfO7SZwxpm5YITCO2HEkl/jEFDbtP01cv/b8dvoQ6+ljTJCyQmBq1fliD39avYs/rtpJ8yYN+b/vX8zUiztXuz+QMcZ5VghMrUk5eJrZC1PYdjiXay/qzKPXDqRdRJNAxzLGVMEKgamxgqIS/pF5nhUr1tK+RRP+fNsIrh4YFehYxhg/WSEwNbJu9wnmJKaw90QRN4/qypxJA2gVXntN4owxznPsOgIReU1EjopIWgXLRUSeF5GdIpIiIsOdymJqX25BEQ++m8qMV9bhUZg9silPXj/UioAxLuTkBWVvABMrWT4J6ON7zARedjCLqUWrth3hu3/4lHfW7+cnl/dgxf9eycC21h7CGLdy7KshVf1URGIqGTIVeEtVFVgnIpEi0klVDzmVydTMybPn+fX76SzenE3fqAheuuUyhnWzJnHGuJ14P4cd2ri3ECxV1cHlLFsKzFPVz3zTK4F4Vd1QztiZePcaiIqKik1ISKhWnry8PCIi3HNFa7DkVVX+fbiEtzMKOVcMU3qGcW2vMBqVujAsWLL6y0153ZQV3JXXTVmhZnnj4uI2quqIcheqqmMPIAZIq2DZB8DlpaZXArFVbTM2NlarKykpqdrrBkIw5D10Ol/vfONL7R6/VK/74xrdeiin3HHBkPVCuCmvm7Kquiuvm7Kq1iwvsEEr+FwN5FlDB4Gupaa7ANkBymLKUFUSvjzA7z7YSpHHw4OTB3DH5T1oaO0hjKl3AlkI3gPuFZEEYDSQo3Z8ICjsO3GWOYmpfLH7BJf0bMO864cS0655oGMZYxziWCEQkXeAsUA7ETkIPAqEAajqAmAZMBnYCZwDbncqi/FPiUd5fe0efv9RJmENGvDk9UOYMbKrtYcwpp5z8qyhm6tYrsA9Tr2+uTCZh3OZnZjClgOnGT+gA7+ZNoSOrZoGOpYxpg7YlcUh7nyxh5eSd/Ji0k5aNA3j+ZuHce3QTrYXYEwIsUIQwjYfOE38whQyj+Qy9eLOPHrtINo0bxzoWMaYOmaFIATlny/hmY8yeW3tHjq0aMpffjSC7wywJnHGhCorBCHm813HmZOYyv6T57hldDfiJ/WnZVPrD2RMKLNCECLOFBTx5LKtvLP+ADFtm5Ew8xIu6dk20LGMMUHACkEI+CTjCA8uTuVYbiE/u7In/zu+L+GNrUmcMcbLCkE9djyvkMffz+D9Ldn079iCP982gqFdIgMdyxgTZKwQ1EOqypLN2Tz+fjp5hcX88uq+3HVVLxo3crLruDHGrawQ1DPZp/N5aHEaq7YdZVi3SJ66YSh9o1oEOpYxJohZIagnPB7l7+v3M+/DbZR4lEemDORHl8VYkzhjTJWsENQDe46fJT4xhfV7TjKmd1uenD6Ubm2bBTqWMcYlrBC4WHGJh1c/28MfPt5O40YNePqGodw0oou1hzDGXBArBC6VkX2G+MQUUrNyuHpgFL+ZNpioltYkzhhz4awQuExhcQkvrNrJy8m7iGwWxos/GM7kIR1tL8AYU21WCFxk475TxCemsPNoHtcPi+bhKQNpbU3ijDE1ZIXABc4WFvP7jzJ54/O9dGrZlNdvH0lcvw6BjmWMqSesEAS5NTuO8cCiVA6eyue2S7sze2J/IprY/zZjTO2xT5QglXOuiL+kFrJm+Xp6tGvOP392KaN6tAl0LGNMPWSFIAitSD/Mw4vTOJ5XzN1je/H/vtOHpmHWJM4Y4wwrBEHkWG4hj72XzgephxjQqSU/Hyz8eGL/QMcyxtRzjnYhE5GJIpIpIjtFZE45y8eKSI6IbPY9HnEyT7BSVRI3HmT8s6v5OOMIsyb04717xxDTyvYCjDHOc2yPQEQaAi8CVwMHgS9F5D1VzSgzdI2qTnEqR7A7eOocD76bxurtx4jt3pqnbhhK7w4RgY5ljAkhTn41NArYqaq7AUQkAZgKlC0EIcnjUf7273089eE2FHjs2oHcdmkMDaxJnDGmjomqOrNhkRuBiar6E9/0rcBoVb231JixQCLePYZs4H5VTS9nWzOBmQBRUVGxCQkJ1cqUl5dHRETg/9o+lOfh9fRCtp/yMKhtA348qAntm337W7pgyesPN2UFd+V1U1ZwV143ZYWa5Y2Li9uoqiPKXaiqjjyAm4BXS03fCvyxzJiWQITv+WRgR1XbjY2N1epKSkqq9rq14Xxxib6YtEP7PLhMhzy6XP/55X71eDwVjg903gvhpqyq7srrpqyq7srrpqyqNcsLbNAKPled/GroINC11HQXvH/1ly5CZ0o9XyYiL4lIO1U97mCugEjLyiE+MYX07DNMGtyRx6cOokMLaxJnjAk8JwvBl0AfEekBZAEzgB+UHiAiHYEjqqoiMgrvWUwnHMxU5wqKSvjjqh0sWL2b1s0a8/Itw5k0pFOgYxljzH84VghUtVhE7gVWAA2B11Q1XUTu8i1fANwI3C0ixUA+MMO3C1MvbNh7ktmJKew+dpYbY7vw0DUDiGxmTeKMMcHF0QvKVHUZsKzMvAWlnr8AvOBkhkDIKyxm/vJtvLVuH51bhfPWHaO4sm/7QMcyxphy2ZXFtWz19mPMXZRKdk4+P7o0hlkT+tHcmsQZY4KYfULVktPnzvPE0q0kbjpIz/bN+dfPLmVEjDWJM8YEPysEteDD1EM8vCSdU+fOc09cL/5nnDWJM8a4hxWCGjh6poBHlqSzPP0wgzq35M07RjKoc6tAxzLGmAtihaAaVJWFGw/yxNIMCoo9xE/sz0+v6EGjho728DPGGEdYIbhAB06eY+67qazZcZyRMa2Zd8NQerV3zyXqxhhTlhUCP5V4lLe+2Mv8FZkI8MTUQdwyurs1iTPGuJ4VAj/sPJpLfGIqG/ed4qq+7fnt9MF0ad0s0LGMMaZWWCGoRFGJhz+t3sXzK3fSrElDnv3eRUwfFo2I7QUYY+oPKwQVSMvKYdbCFLYeOsM1Qzvx2LWDaN+iSaBjGWNMrbNCUEZBUQn/98kO/rxmN22aN+ZPt8YyYVDHQMcyxhjHWCEoZf2ek8xJTGH38bN8f0RX5k4eQKtmYYGOZYwxjrJCAOQWFPH08kz+um4fXVqH87c7R3N5n3aBjmWMMXUi5AtBUuZRHlyUyqEzBdwxpgf3T+hLs8Yh/7YYY0JIyH7inTp7nieWZrDoqyx6d4hg4V2XEdu9daBjGWNMnQu5QqCqfJB6iEeXpJOTX8R943pzz7jeNGlkTeKMMaEppArBqQIPP/vrRj7KOMKQ6Fb87SejGdCpZaBjGWNMQIVMIUjadpS5n+XjoZAHJvXnzsutSZwxxkAIFYIe7ZrTO7Ihz//4Cnq0ax7oOMYYEzQc/ZNYRCaKSKaI7BSROeUsFxF53rc8RUSGO5Ulpl1zfjWiqRUBY4wpw7FCICINgReBScBA4GYRGVhm2CSgj+8xE3jZqTzGGGPK5+QewShgp6ruVtXzQAIwtcyYqcBb6rUOiBSRTg5mMsYYU4aoqjMbFrkRmKiqP/FN3wqMVtV7S41ZCsxT1c980yuBeFXdUGZbM/HuMRAVFRWbkJBQrUx5eXlERLjnJjJuyuumrOCuvG7KCu7K66asULO8cXFxG1V1RHnLnDxYXF6v5rJVx58xqOorwCsAI0aM0LFjx1YrUHJyMtVdNxDclNdNWcFded2UFdyV101Zwbm8Tn41dBDoWmq6C5BdjTHGGGMc5GQh+BLoIyI9RKQxMAN4r8yY94DbfGcPXQLkqOohBzMZY4wpw7GvhlS1WETuBVYADYHXVDVdRO7yLV8ALAMmAzuBc8DtTuUxxhhTPkcvKFPVZXg/7EvPW1DquQL3OJnBGGNM5Rw7a8gpInIM2FfN1dsBx2sxjtPclNdNWcFded2UFdyV101ZoWZ5u6tq+/IWuK4Q1ISIbKjo9Klg5Ka8bsoK7srrpqzgrrxuygrO5bWua8YYE+KsEBhjTIgLtULwSqADXCA35XVTVnBXXjdlBXfldVNWcChvSB0jMMYY822htkdgjDGmDCsExhgT4uplIQimG+L4w4+8/UXkCxEpFJH7A5GxVJaqst7ie09TRORzEbkoEDl9WarKOtWXc7OIbBCRywORs1SeSvOWGjdSREp8HX4Dwo/3dqyI5Pje280i8kggcpbKU+V768u8WUTSRWR1XWcslaOq93ZWqfc1zfez0KZGL6qq9eqBt53FLqAn0BjYAgwsM2Yy8CHe7qeXAP8O8rwdgJHAb4H7gzzrZUBr3/NJgXpv/cwawTfHyYYC24L5vS01bhXeK/ZvDNaswFhgaaDez2rkjQQygG6+6Q7BmrXM+GuBVTV93fq4R+C2G+JUmVdVj6rql0BRIAKW4k/Wz1X1lG9yHd6OsoHgT9Y89f02Ac0ppwV6HfLn5xbgf4BE4GhdhivD36zBwp+8PwAWqep+8P7O1XHGr13oe3sz8E5NX7Q+FoJo4ECp6YO+eRc6pq4EU5aqXGjWO/HueQWCX1lFZLqIbAM+AO6oo2zlqTKviEQD04EFBJa/PweXisgWEflQRAbVTbRy+ZO3L9BaRJJFZKOI3FZn6f6b379jItIMmIj3D4MacbTpXIDU2g1x6kgwZamK31lFJA5vIQjU9+7+3vToXeBdEbkSeAIY73SwCviT9//w3sGvRKS84XXGn6yb8Pa2yRORycBivPcmDwR/8jYCYoHvAOHAFyKyTlW3Ox2ujAv5PLgWWKuqJ2v6ovWxELjthjjBlKUqfmUVkaHAq8AkVT1RR9nKuqD3VVU/FZFeItJOVQPRhMyfvCOABF8RaAdMFpFiVV1cJwm/UWVWVT1T6vkyEXkpyN/bg8BxVT0LnBWRT4GLgLouBBfyczuDWvhaCKiXB4sbAbuBHnxzsGVQmTHX8N8Hi9cHc95SYx8jsAeL/Xlvu+G9v8RlLvg56M03B4uHA1lfTwdj3jLj3yBwB4v9eW87lnpvRwH7g/m9BQYAK31jmwFpwOBgzOob1wo4CTSvjdetd3sE6rIb4viTV0Q6AhuAloBHRP4X75kEZyrabqCyAo8AbYGXfH+5FmsAujv6mfUGvHfIKwLyge+r77csSPMGBT+z3gjcLSLFeN/bGcH83qrqVhFZDqQAHuBVVU0Lxqy+odOBj9S7B1Nj1mLCGGNCXH08a8gYY8wFsEJgjDEhzgqBMcaEOCsExhgT4qwQGGNMiLNCYOo9EekoIgkisktEMkRkmYj0reMMY0XkskqWTyvbodPXnuGdMvN+LyLjnMppQpMVAlOvifdihneBZFXtpaoDgblA1AVso2Fl034ai7cza0VmAy+Veo0BeH8/rxSR5qXG/RGosEW1MdVhhcDUd3FAUekLslR1s6qu8f2VvvTr+SLygoj82Pd8r4g8IiKfATeVM/1d8d4jYpOI/EtEIkqt97hvfqp47yURA9wF/MLXQ/6K0gF9eyeF+t/tF34A/BX4CLiuVPZ9QFvfRYbG1AorBKa+GwxsrOa6Bap6uaomlJ4GPgEeAsar6nC8V33/stR6x33zX8bbEmQv3o6hf1DVi1V1TZnXGYO3SVtp3wf+gbeXzM1llm3yrWNMrah3LSaMqUX/qGD6EmAgsNbXRqMx8EWpcYt8/24ErvfjdToBx76eEJGRwDFV3SciB4HXRKS1fnOfh6NA5wv5DzGmMlYITH2XjrfvTXmK+e+94qZllpft4/L1tAAfq2rZv9S/Vuj7twT/fsfy8TYR+9rNQH8R2eubbom3L9KrpXLm+7FdY/xiXw2Z+m4V0EREfvr1DPHe8/cqYB8wUESaiEgrvL3o/bEOGCMivX3ba+bHWUi5QIsKlm3F2wkVEWkA3AQMVdUYVY3Be4eq0kWnL97umMbUCisEpl7zdbycDlztO300HW8772xVPQD8E2/HybeBr/zc5jHgx8A7IpKCtzD0r2K194Hp5R0sBj4FhvnOcLoSyFLVrDLLB4pIJxEJw1s0NviT1Rh/WPdRY4KAiDwHvK+qn1QxbjowXFUfrptkJhTYHoExweF3eG+IUpVGwDMOZzEhxvYIjDEmxNkegTHGhDgrBMYYE+KsEBhjTIizQmCMMSHOCoExxoS4/w/j8WfL6/3NmwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Visualize results\n", "I_line = np.arange(0, 0.8, 0.1).reshape(8, 1)\n", "V_line = x_ls[0]*I_line + x_ls[1]\n", "\n", "plt.scatter(I, V)\n", "plt.plot(I_line, V_line)\n", "plt.xlabel('Current (A)')\n", "plt.ylabel('Voltage (V)')\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the offset parameter $\\hat{b}$ is near zero, while $\\hat{R}$ closely approximates the true resistance value of $R = 5~\\Omega$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2: Recursive Least Squares\n", "----\n", "\n", "Next, let's assume you don't get all measurements at once, but you get one set of $[V, I]$ measurements at a time, and you need to update your estimation of $[R, b]$ every time you get a new measurements.\n", "\n", "In this case, you will need to use Recursive Least Square.\n", "\n", "### Initialize\n", "----\n", "\n", "To use the recursive least squares formulation, you must have a prior estimate of the resistance and its associated uncertainty (otherwise, you won't know how to weigh the information you receive from a new measurement). You choose to set the initial parameters under the assumption that your prior estimate of the resistance, $R = 4$, is not very good. Also, since you are fairly sure that Ohm's law ($V = RI$) does, in fact, hold, you feel that it is safe to assume with high confidence that the offset term $b$ is close to zero. After some thought, you choose to intialize the recursive estimator as follows:\n", "\n", "$$\\hat{R} \\sim \\mathcal{N}(4, 9.0),~~\\hat{b} \\sim \\mathcal{N}(0, 0.2)$$\n", "\n", "Your initial guess is that $\\hat{R}$ follows a Gaussian or normal distribution (recall that you do not know the exact value of $R$, so it must be considered as a random variable) with a mean of $4~\\Omega$ and a standard deviation of $3~ \\Omega$ (i.e., a variance of $9~\\Omega^{2}$). Similarly, your intial guess is that $\\hat{b}$ should also follow a normal distribution with a mean of $0~V$ and a variance of $0.2~V^{2}$.\n", "\n", "With this initial guess, you can initiate the recursive least square.\n", "\n", "**Initialize the parameter and covariance estimates**:\n", "\n", "----\n", "$$\n", "\\hat{\\mathbf{x}}_0 = E\\left[\\mathbf{x}\\right],\\quad \\mathbf{P}_0 = E\\left[(\\mathbf{x} - \\hat{\\mathbf{x}}_0)(\\mathbf{x} - \\hat{\\mathbf{x}}_0)^T\\right]\n", "$$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# Initialize the 2x1 parameter vector x (i.e., x_0).\n", "# In this specific case, our x is [R, b]\n", "x_k = np.array([[4.0, 0.0]]).T\n", "\n", "#Initialize the 2x2 covaraince matrix (i.e. P_0). Off-diangonal elements should be zero.\n", "P_k = np.array([[9.0, 0], [0, 0.2]])\n", "\n", "# Our voltage measurement variance (denoted by R, don't confuse with resistance).\n", "R_k = np.array([[0.0225]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Recursive Updates\n", "\n", "After initialization, you can code the loop that updates your estimation with each new measurement.\n", "\n", "**For every measurement k**:\n", "\n", "----\n", " * Calculate the gain term: $$\\mathbf{K}_k = \\mathbf{P}_{k-1}\\mathbf{H}_k^T\\left(\\mathbf{H}_k\\mathbf{P}_{k-1}\\mathbf{H}_k^T + \\mathbf{R}_k\\right)^{-1}$$\n", " * Update the parameter estimate: $$\\hat{\\mathbf{x}}_k = \\hat{\\mathbf{x}}_{k-1} + \\mathbf{K}_k\\left(\\mathbf{y}_k - \\mathbf{H}_k\\hat{\\mathbf{x}}_{k-1}\\right)$$\n", " * Update the covariance estimate: $$\\mathbf{P}_k = \\left(\\mathbf{I} - \\mathbf{K}_k\\mathbf{H}_k\\right)\\mathbf{P}_{k-1}$$\n", " \n", "In this case, the initial parameter vector $\\hat{\\mathbf{x}}_0$ should contain $\\hat{R}$ and $\\hat{b}$. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slope and offset parameters of the best-fit line (i.e., the resistance and offset) are [R, b]:\n", "Estimates value for R is: 4.976925034352892\n", "Estimates value for b is: 0.06966257824825756\n" ] } ], "source": [ "# Pre allocate space to save our estimates at every step.\n", "num_meas = I.shape[0]\n", "x_hist = np.zeros((num_meas + 1, 2))\n", "P_hist = np.zeros((num_meas + 1, 2, 2))\n", "\n", "x_hist[0] = x_k.T[0]\n", "P_hist[0] = P_k\n", "\n", "# Iterate over all the available measurements.\n", "for k in range(num_meas):\n", "\n", " # Construct H_k (Jacobian).\n", " H_k = np.array([[I[k, 0], 1]])\n", "\n", " # Construct K_k (gain matrix).\n", " K_k = P_k.dot(H_k.T).dot(inv(H_k.dot(P_k).dot(H_k.T)+R_k))\n", " \n", " # Get the new data\n", " y_k = V[k, 0]\n", " \n", " # Update our estimate.\n", " x_k = x_k + K_k.dot(y_k - H_k.dot(x_k))\n", " \n", " # Update our uncertainty (covariance)\n", " P_k = (np.identity(2) - K_k.dot(H_k)).dot(P_k)\n", "\n", " # Keep track of our history.\n", " P_hist[k + 1] = P_k\n", " x_hist[k + 1] = x_k.T[0]\n", " \n", "print('The slope and offset parameters of the best-fit line (i.e., the resistance and offset) are [R, b]:')\n", "print('Estimates value for R is:', x_k[0, 0])\n", "print('Estimates value for b is:', x_k[1, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the Results\n", "----\n", "Let's plot out the solution at each step. Does the resistance value converge towards the batch least squares solution?" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABzgUlEQVR4nO2ddXhUx/rHPxN3JSQhAQIEDxAgQHB3p7TFS/XWfm1v5VK5FWh7a9Rue1tKDQja4lqkNLgGCRIkeILEXXd3fn+ckAaILMluBObzPDzk7Jkz53tONvvumfed7wgpJQqFQqG4f7GoagEKhUKhqFpUIFAoFIr7HBUIFAqF4j5HBQKFQqG4z1GBQKFQKO5zVCBQKBSK+xwrc59ACGEJHARipZTDbtsngK+BIUAWMFVKeai0/mrVqiUDAgLKpSUzMxNHR8dyHVsV1CS9NUkr1Cy9NUkr1Cy9NUkrVExvREREgpTSq9idUkqz/gNeBhYCa4vZNwTYAAggFNhXVn/t27eX5eWvv/4q97FVQU3SW5O0Slmz9NYkrVLWLL01SauUFdMLHJQlfK6adWhICOEPDAV+KqHJSGBegc69gJsQwtecmhQKhUJxK+bOEXwF/AswlLDfD7hSZDum4DWFQqFQVBJCmsliQggxDBgipXxWCNELeFXemSNYB3wkpdxZsP0n8C8pZcRt7Z4CngLw9vZuv3jx4nJpysjIwMnJqVzHVgU1SW9N0go1S29N0go1S29N0goV09u7d+8IKWVIcfvMmSzuCowQQgwB7AAXIcR8KeWkIm1igLpFtv2Bq7d3JKWcDcwGCAkJkb169bplf35+PjExMeTk5JQqyNXVFTs7u3JcStVQk/RWVKudnR3+/v5YW1ubUFXJhIeHc/v7qLpSk7RCzdJbk7SC+fSaLRBIKd8A3gAo8kQw6bZmq4HnhRCLgU5AqpTy2t2eKyYmBmdnZwICAtAKkYonPT0dZ2fnu+2+yqhJeiuiVUpJYmIiMTExNGjQwMTKFApFWVT6PAIhxNNCiKcLNtcD54Fo4Efg2fL0mZOTg6enZ6lBQFF9EULg6elZ5hOdQqEwD2afRwAgpQwHwgt+nlXkdQk8Z4pzqCBQs1G/P4Wi6qiUQKBQKBSK8pNvyGfuiblY5ZrnI1tZTJgIS0tLgoODadOmDe3atWP37t2ltk9JSeG7774rs99evXpx8ODBUtsYDAZeeOEFgoKCaNWqFR06dODChQulHhMQEEBCQkKpbebMmcPVq3/n7p944glOnjxZpmaFQmE6jsQd4aE1D/H1oa+JzIo0yznUE4GJsLe358iRIwBs3LiRN954g23btpXY/mYgePbZcqVFbmHJkiVcvXqVyMhILCwsiImJMcm0+Tlz5hAUFESdOnUA+OmnkuYFKhQKU5Oam8rXh77m9zO/4+Pow397/xdx3jxDqOqJwAykpaXh7u4OaHW/ffv2pV27drRq1YpVq1YB8Prrr3Pu3DmCg4N57bXXAPj0009p1aoVbdq04fXXXy/s7/fff6djx440adKEHTt23HG+a9eu4evri4WF9uv09/cvPP+iRYto1aoVQUFBTJs27Y5jL168SFBQUOH2zJkzee+991i6dCkHDx5k4sSJBAcHk52dfcvTSUn9Ojk58dZbb9GmTRtCQ0O5ceNGhe6lQnG/IaVkw4UNjFw5kmVnlzGlxRRWjVxF73q9zXbOe+6JYPqaE5y8mlbsPr1ej6Wl5V332aKOC+8Ob1lqm+zsbIKDg8nJyeHatWts3boV0OrjV6xYgYuLCwkJCYSGhjJixAg+/vhjjh8/XvgUsWHDBlauXMm+fftwcHAgKSmpsG+dTsf+/ftZv34906dPZ8uWLbec+6GHHqJbt27s2LGDvn37MmnSJNq2bcvVq1eZNm0aERERuLu7M2DAAFauXMmoUaPKvOaxY8fy7bffMnPmTEJCbp2DUly/a9euZfz48WRmZhIaGsqHH37Iv/71L3788Uf+/e9/G3GXFQrFlfQrfLj3Q3Zd3UWQZxDf9/ue5p7NzX5e9URgIm4ODZ06dYo//viDKVOmFBo6vfnmm7Ru3Zp+/foRGxtb7LfkLVu28Oijj+Lg4ACAh4dH4b4xY8YA0L59ey5evHjHsf7+/pw+fZqPPvoICwsL+vbty59//smBAwfo1asXXl5eWFlZMXHiRLZv317hay2u3127dgFgY2PDsGHDStWrUChuJd+Qz0/HfmL0qtEciT/C6x1fZ/6Q+ZUSBOAefCIo7Zt7ZU3Q6ty5MwkJCcTHx7N+/Xri4+OJiIjA2tqagICAYuvlpZQlllDa2toCWkJap9OV2Gbw4MEMHjwYb29vVq5cSd++fcvUamVlhcHwtxWUMbX8pdmSWFtbF15HaXoVCoXG4bjDzNgzg+iUaPrV68e0jtPwcfSpVA3qicAMnDp1Cr1ej6enJ6mpqdSuXRtra2v++usvLl26BICzszPp6emFxwwYMIBffvmFrKwsgFuGhsri0KFDhdU9BoOByMhI6tevT6dOndi2bRsJCQno9XoWLVpEz549bznW29ubuLg4EhMTyc3NZe3atYX7btd4k+L67datm/E3SKFQkJqbyvQ905myYQqZ+Zl80+cbvuz9ZaUHAbgHnwiqips5AtC+Mc+dOxdLS0smTpzI8OHDCQkJITg4mGbNmgHg6elJ165dCQoKYvDgwXz22WccOXKEkJAQbGxsGDJkCG+88YZR546Li+PJJ58kNzcXgI4dO/L8889jZ2fHRx99RO/evZFSMmTIEEaOHHnLsdbW1rzzzjt06tSJBg0aFOoDmDp1Kk8//TT29vbs2bOn8HVfX987+h06dGhFbp9Ccd9wMxn8yYFPSM1N5ZEWj/Bs8LM4WDtUraia9K+4hWlOnjxpxLIMUqalpRnVrrpQk/SaQquxv0dTUJMWJKlJWqWsWXorW+vl1MvyqU1PyaA5QXLcmnHyZMLdvefNtTCNeiJQKBQKM5Ovz2fuybnMOjoLKwsr3uj4Bg83fRhLi7uvYjQHKhAoFAqFGSmaDO5fvz/TOkzD29G7qmXdggoECoVCYQZSc1P5MuJLlp1dhq+jL9/0+YZedXtVtaxiUYFAoVAoTIiUkvUX1vPpgU+rTzK4DFQgUCgUChNxJe0K7+99nz3X9hDkGcSsfrMqbVJYRVCBQKFQKCpIvj6fOSfm8EPkD9UyGVwWakKZiRBCMHny5MJtnU6Hl5dXod3Cvc6OHTtKtN6WUvLCCy8QGBhI69atOXToUCWrUyjMx6Ebh3hwzYP89/B/6eHfg1UjVzGh+QSTB4GcfD35hpJn9VcEsz0RCCHsgO2AbcF5lkop372tTS9gFXDTPH+5lHKGuTSZE0dHR44fP052djb29vZs3rwZPz+/KtGi0+mwsqrch70dO3bg6elJly5d7ti3YcMGzp49y9mzZ9m3bx/PPPMM+/btq1R9CoWpuT0Z/G2fb+lZt2fZB94lUko2HL/Oh+ui6Oylo7/Jz2DeJ4JcoI+Usg0QDAwSQoQW026HlDK44F+NDAI3GTx4MOvWrQM0m+bx48cX7svMzOSxxx6jQ4cOtG3bttCO+uLFi3Tv3p127drdsqDNtWvXGDRoEMHBwQQFBRXaTzs5ORX2uXTpUqZOnQpos4BffvllevfuzbRp0zh37hyDBg2iffv2dO/enVOnThW2e+aZZ+jduzcNGzZk27ZtPPbYYzRv3rywL4BNmzbRuXNn2rVrx4MPPkhGRgagLWjz7rvvFtpqnzp1iosXL/LLL7/w5ZdfEhwcfIdV9qpVq5gyZQpCCEJDQ0lJSeHatWsmvPMKReUhpWTd+XWMWDmCldErmdpyKitHrjRLEDh1PY0JP+7j2QWHcLazItDdPENNZvvaWDCTLaNg07rgn3mea4qy4XW4fqzYXfZ6HViW45J9WsHgj8tsNm7cOGbMmMGwYcOIjIzkscceK/xQ/PDDD+nTpw+//PILKSkpdOzYkX79+lG7dm02b96MnZ0dZ8+eZfz48Rw8eJCFCxfSt29fZsyYgV6vL/QgKo0zZ86wZcsWLC0t6du3L7NmzaJx48bs27ePZ599ttAaOzk5ma1bt7J69WqGDx/Orl27+Omnn+jQoQNHjhzB39+fDz74gC1btuDo6Mgnn3zCF198wTvvvANArVq1OHToEN999x0zZ87kp59+4rHHHsPT05NXX331Dl2xsbHUrVu3cNvf35/Y2Fh8fX2Nuv0KRXWhaDK4Va1W/ND/B5p5NCv7wLskJSuPLzafYf7eS7jYW/P+yJaM71iPnTsq7h5cHGYdPxBCWAIRQCDwPyllceMBnYUQR4GrwKtSyhPm1GROWrduzcWLF1m0aBFDhgy5Zd+mTZtYvXo1M2fOBDSXz8uXL1OnTh2ef/55jhw5gqWlJWfOnAGgQ4cOTJ06FQsLC0aNGlXoY1QaDz74IJaWlmRkZLB7924efPDBwn03fYgAhg8fjhCCVq1a4e3tTatWrQBo2bIlFy9eJCYmhpMnT9K1a1cA8vLy6Ny5c+HxRW2xly9fXqYuWYxbqVqsXlGTuD0Z/GanN3moyUMmzwPo9AYW7b/M55vPkJadz6TQ+vyzXxPcHW1Mep7bMWsgkFLqgWAhhBuwQggRJKU8XqTJIaC+lDJDCDEEWAk0vr0fIcRTwFOguWWGh4ffst/V1fVvl8xub5Wop7wL0wBQjAvnnU3SGThwIK+88grr168nKSkJnU5Heno6er2eefPm0bjxrZf3n//8B3d3d3bu3InBYMDLy4v09HTatm3LunXr2LJlCxMnTuSFF15gwoQJCCEKrzU5OZn8/HzS09PJz8/HwsKC9PR00tLScHV1vWOI5mY7g8FAeno6WVlZWFtbF/an1+tJT0/H1taWXr168euvv95xvJSy8Jw5OTnk5uYWvn7z59vx9vbmzJkztGnTBoDLly/j4uJyR9ucnJw7frfmIiMjo9LOVVFqklaoWXqN0Xou5xyLkxZzPf86wQ7BjHUfi+t1V3Zcv3O1wIoQlahnQVQuMRmSZh4WvNLWnrrOCRw98Pfa4ua6t5WSUZRSpgghwoFBwPEir6cV+Xm9EOI7IUQtKWXCbcfPBmYDhISEyF69et3Sf1RUlFHrDJh7PQJnZ2eeeeYZateuTWhoKOHh4VhZWeHs7MzgwYP55Zdf+OabbxBCcPjwYdq2bUtOTg7169fH1dWVX3/9Fb1ej7OzM5cuXcLHx4f/+7//Q6/XF16jt7c3MTExNG3alD/++ANnZ2ecnZ2xtrbG3t6+cLthw4b88ccfPPjgg0gpiYyMpE2bNre0c3JywsLCovCe3NzXs2dPXn31VW7cuEFgYCBZWVnExMTQpEkThBA4OTnh7OyMo6MjlpaWODs74+LiQm5ubrH394EHHuDbb7/l0UcfZd++fbi7u98REEFbza1t27Zm+/0UJTw8nNvfR9WVmqQVapbe0rQWJoMvFSSDu5snGRyTnMV/1kex/th1/Nzs+W5icwYH+RT71Gyue2u2ZLEQwqvgSQAhhD3QDzh1WxsfUXC1QoiOBXoSzaWpMvD39+fFF1+84/W3336b/Px8WrduTVBQEG+//TYAzz77LHPnziU0NJQzZ84ULjofHh5O165dadu2LcuWLSvs8+OPP2bYsGH06dOn1DH2BQsW8PPPP9OmTRtatmxZmJw2Bi8vL+bMmcP48eNp3bo1oaGhhcnmkhg0aBArVqwoNlk8ZMgQGjZsSGBgIE8++STfffed0VoUispGSsna82vNngzOztPz5eYz9P18G1tPxfHPfk3485WeDGnlW/lDpyXZklb0H9AaOAxEoj0FvFPw+tPA0wU/Pw+cAI4Ce4EuZfWrbKirJ8qG2nzUJK1S1iy9t2u9lHpJPrHxCRk0J0iOXzteRiVGmfycBoNBrjkaKzv/Z4usP22tfG5BhIxJziqX3ruBqrChllJGAnc850spZxX5+VvgW3NpUCgUCmPI1+fz64lf+eHoD9hY2vBWp7d4sMmDJk8Gn7yaxvQ1J9h3IYnmvi588XAwoQ09TXqO8qAsJhQKxX1NxI0IZuyZwfnU8wyoP4BpHadR26G2Sc+RlJnH55tOs2j/ZVztrflgVBDjO9bD0sK4ISApJZeOHSEnJdmkum6iAoFCobgvSc1NZWHiQvb8sYc6jnX4X9//0cO/h0nPodMbWLDvMl9sPkNGro4pnQN4qV9j3ByMLwe9Fn2anYvmcvl4JLVatIFRo02qEVQgUCgU9xmyIBk88+BMUnJSeLTlozzd5mmT20Tvjk5g+pqTnL6RTpdGnrw7vCVNfYyvWkyMucLOxfOIPrAHe2cXGrR/ACs/f5NqvIkKBAqF4r7hUtol3t/7Pvuu7aN1rdY86f4kk0ImmfQcV5Ky+HBdFH+cuI6/uz2zJrVnYEtvoyuB0hLi2P37Qk5u24qVrS0BwUNJvB7ItQuWeDmYp5pIBQKFQnHPk6/P55fjvzA7cvYtyeAd2003KSwrT8es8HPM2n4eSyF4pX8TnuzREDtr4xLOWWmp7F/5G0c2rgMhqBvUm9SkFly/ZE1A61p0GtGQ49EHTaa3KCoQmAghBJMmTSIsLAzQHEB9fX3p1KkTa9eurWJ15mfHjh24ubkV6z566tQpHn30UQ4dOsSHH35YrB+RQmEuzJ0MllKyJvIaH62P4lpqDiPa1OGNIc3wdbU36vi87Cwi1q3i4Nrl5Ofk4te8M5kZbYiLscO/mTudRjTEp6Gr1jjaZLJvQQUCE6FsqEu2ofbw8OC///0vK1eurFRNivub1NxUvoj4guVnl5stGXw8NpUZa06y/2ISLXxd+HpcWzo28DDqWF1+PpFbNrB3+RKy01LxadyO/LwOJFxzxLuBCwMeb4h/M+P6qihqYRoTomyoi7ehrl27Nh06dMDa2tqEd1uhKB4pJWvOrWHEyhGsil7Foy0fZcXIFSYNAokZubyx/BjDv91JdHwG/xndijX/182oIGAw6Dmx7U9+/ec/+GvObJzc6+DV8HFSEnph51ybIc+25oF/ta+0IAD34BPBJ/s/4VRS8XYI5TWda+bRjGkdp5XZTtlQF29DrVBUFrcng2f3n01Tj6Z3tFt5OJbPNp5mXN103vp4K68NbMqotmU/wefrDczfe4kvN58hM0/P1C4BvNS3Ca4OZX/JkVJy7uA+di6eR2LMZdzrBFA7cDJpiV64etkz4PGGBLavjTByboEpuecCQVWibKgViqohT5/Hr8d/LUwG/7vTvxnbZGyxM4NXHo7ljeXHyM7XQ12ITcnmjeXaGialBYOdZxOYvuYEZ+My6N64Fu8Ma0Fjb+PKQa+ciGTHorlcO3sa51q+eDd+mJT4Ojjb2dF7UgOadvbB0rLqBmjuuUBQ2jd3c7uPAowYMYJXX32V8PBwEhP/9s+TUrJs2TKaNr3128l7772Ht7c3R48exWAwYGdnB0CPHj34448/2LZtG5MnT+a1114rXOXrJjk5Obf0ddOwzmAw4ObmxpEjR4rVaGtrC4CFhUXhzze3dTodlpaW9O/fn0WLFpV6vKWlJTqdzpjbolCYjYPXD/L+3vc5n3qegQED+VeHf5WaDP5s42ktCBQhO1/PZxtPFxsILidm8cG6k2w6eYN6Hg7Mntye/i2MKwe9cT6anYvncfHoIRxcPfBpMprkuPpY5drS/aEAWvaog5WRVUXm5J4LBFXNY489hqurK61atbrFN3zgwIF88803d9hQp6am4u/vj4WFBXPnzkWv196gly5dwsvLiyeffJLMzEwOHTrElClT8Pb2JioqiqZNm7JixYpiA5uLiwsNGjTg999/v8OG2hhCQ0N57rnniI6OvsOGuiScnZ2LXYtAoTAXKTkpfBHxBSuiV9xVMvhqSrZRr2fm6vguPJofd1zAUgheG9iUx7s1MKocNPlaLLuWzOf0nh3YOjjh02QoyfENycqwJ3RUPVr39sfGrvp8/FYfJfcIpdlQv/TSS7Ru3RopJQEBAaxdu5Znn32WBx54gN9//53evXvfYkP9ySefYGtri5OTE/PmzQP+tqGuW7cuQUFBhUnc21mwYAHPPPMMH3zwAfn5+YwbN87oQFDUhvrmkNIHH3xQaiAYNGgQU6dOZdWqVXzzzTd07969cN/169cJCQkhLS0NCwsLvvrqK06ePImLi4tRehSKotycGfzZgc9Iy0vj0aBHebq18TOD67jZE1tMMKjjZl/Y/6ojV/l4wymup+UwKrgOrw9ujo+rXZl9ZyQlsmfZIo5t3YSllTU+jfuSktCcjFQ72g+qS9v+9bBzrH5FE0IWs4xgdSYkJEQePHjrpIqoqCiaN29e5rGVMTRkSmqSXlNoNfb3aArulcVTqiPm1Ht7Mvidzu8UmwwujaI5glda6fj8mBX21pZ8NKYVgbWdeG/1CQ5eSqaVnyvvjWhB+/plV+/kZGSwf/VSDq9fjcFgoFZAKOkpQQgLR4K6+9F+cAAOLhVfbrIi91YIESGlDClun3oiUCgU1Z48fR6/HP+FHyN/LDMZXBY38wCfbTwNpOPnZs/TPRuy93wi//ztCB4ONnw8phUPhtQt0x00PyeHQxtWc2D1MnKzs/Cq357MjLakJzvTrIsvHYY2wNmj7CeJqkYFAoVCUa05eP0gM/bO4ELqBQYGDGRah2l4OXhVqM9Rbf0Y1daPLVv/4tF69fl042my8/Q81rUBL/RtjKt96cM3ep2OY1s3sXfZIjJTkvGsG4R1bgjpqR4EhtSm47AGuPs4VkhjZWK2QCCEsAO2A7YF51kqpXz3tjYC+BoYAmQBU6WUh8ylSaFQ1ByKJoP9nPz4ru93dPfvXvaBRrL9TDzv7MrmamYUPZp48c6w5gTWLn14UxoMnNqzg91L5pNy4xpuPo1w9h5CZoY3Aa086TiiIV51a8ZwblHM+USQC/SRUmYIIayBnUKIDVLKvUXaDAYaF/zrBHxf8L9CobhPKS4Z/EybZ7C3Ms67pywuJWby/tootkTdoLaD4KcpIfRtXrvUclApJRePRLBj0VziL13AuZY/Lj4Pkp3jj39Td0JHNfrbD6gGYs6lKiVws6TFuuDf7ZnpkcC8grZ7hRBuQghfKeU1c+lSKBTVl4upF/lg7wfsu17+ZHBJZObq+PavaH7ecQErS8G0Qc0INFymXwvvUo+LPXWSHYvmEnvqBA6uXrj6jiInuwHe/i6EjmqEfzP3yl9s3sSYNUcghLAEIoBA4H9Syn23NfEDrhTZjil4TQUCheI+omgy2NbSlrdD32Zsk7FYiIrPtpVSsvJILB+tP0Vcei5j2voxbXAzvF3sCA+/UuJx8ZcvsnPxPM5H7MfO0RXXOoPJyWqCm4cLfUY0pEGbWjU+ANykUspHhRBuwArg/6SUx4u8vg74SEq5s2D7T+BfUsqI245/CngKwNvbu/3ixYtv6d/V1ZXAwMAydZTXa8gYXFxcePjhh/nxxx8BzQG0cePGhISE8Pvvv5erT3PqNTXbtm3Dzs6OTp3uHNlbsmQJX331FaDNfv7yyy8LbS2KEh0dTWpqqrmlApCRkXGLgV91piZphbvXezbnLIsTFxOni6OdQzvGuI/B1co0wywXUvUsiMojOsVAAxcLJja3IdD977+p4rTmpqVw9cBuks6cxMLaFhvnDkhDW2ydrfEKErjWo0r8gErSayy9e/eu2vJRKWWKECIcGAQcL7IrBqhbZNsfuFrM8bOB2aDNI7i9jjYqKsqoGnZz1uU7Ojpy+vRprKyssLe3Z8OGDfj7+2NlZVXuc5ZXb1XYUO/atQtPT0/69et3x74WLVqwY8cO3N3d2bBhA//85z/Zt+/2h0Ows7Ojbdu2lSG3RtXm1yStYLzelJwUPo/4nJWXVmrJ4E6mSwbHp+fy2cZT/B4Rg6ejDZ+ObcbYdv5Y3PYBXlRrZkoye5cv4eSWPxBC4FanG9mZrXB0daPD0ACadfGtUj8gQ14e2//6yyzvBbNdlRDCq+BJACGEPdAPuN0WdDUwRWiEAqk1OT+gbKiLt6Hu0qUL7u7ugGZfERMTY6pbrqiBSClZfW41I1aOYM25NTwW9BgrRq4wSRDI0xn4cft5+swMZ/mhWJ7o1oCtr/bioZC6dwSBm+RmZbJrSRg/v/AkRzevx6V2W6wcp4JlV7o91IpJ74fSsrtflQUBXXw88d98S3Sfvjj++adZzmHOr42+wNyCPIEF8JuUcq0Q4mkAKeUsYD1a6Wg0WvnooxU96fX//IfcqOJtqHV6PUnlGGqxbd4MnzffLLOdsqEu24b6559/ZvDgwWVei+Le5JZksFdr3gk1XTI4/HQcM9ae5Hx8Jr2aevH2sBY08ip5GCU/L5cbRw7wU9gP5GSk4+rTmuzs9hiopfkB9albpX5A2cdPkBw2j9T1GyA/H8eePYg3Ygi8PJizaigSuOM5vyAA3PxZAs+ZS0Nlo2yoS+evv/7i559/ZufOnUYfo7g3yNPn8fPxn/kx8kfsLO1Mmgy+kJDJB2tP8uepOBrUcuSXqSH0aVZyJZBBry0Ms3vpQjISE3Dxaoq07IDe4EPI4Lq0HVB1fkBSpyN982aS5oWRffgwFg4OuD/0EO6TJmLboAHnixhZmpJ7bmZxad/clQ21RlXYUEdGRvLEE0+wYcMGPD09jTpGcW9w4PoBZuyZwcW0iwwKGMS/OvyrwjODATJydXyz9Sy/7LyAjaUFrw9uxqNdA7C1Kv6pX0rJ2f272bk4jOSrMTh61MfGpRc6WZ/WffxoP7g+jq62xR5rbnTJyaT8vpTkhQvRXb+Odd26eL/xOq5jxmBZCX5j91wgqGqUDfWdXL58mTFjxhAWFlZqH4p7i8JkcPRKk84MNhgkKw7H8vEfp4hPz+WBdv5MG9SU2i4le/pcijzCjkVzuXH+LPYu3ti5jUIvG+ARKBj+eCgutUwzWe1uyTlzhuSwMFJXr0Hm5uLQORSfd97GqWdPRCVWDKpAYGKUDfWdNtQzZswgMTGRZ599FgArKytud5BV3DvcTAbPPDCT9Lx0Hgt6jKfbPG2SmcFHrqTw3uoTHLmSQpu6bsye3J629dxLbH89+gw7Fs3l8vGj2Dq6Y+82GANNaRziQ8fhDTh66kClBwGp15OxbRtJ88LI2rsXYWuL64gRuE+ehF0VfVFSNtTVmJqkV9lQm4+apPVi6kVe3vAyZ3PP0sarDe90focm7hX/cItLz+HTP06zNCKGWk62vD64GWPa+pVYCZQYe4VdS8I4u2831nZOWNl3wiBbEtDKm04jGuJVT3uvVua91aenk7JsGckLFpJ/5QpWPj64T5iA24NjsXIvOZgVRdlQKxSKakvRZLCVtDJZMjhPZ2DO7gv8989ocnV6/tGjIc/3CcTZrvhkblpCPHuWLuJE+BYsrKyxd++GQQbj27g2oSMb4hvoViE95SH3wgWS5y8gdcUKDFlZ2LdrR+1XXsa5Xz+EsfN9MuJh3yzcU52BXibXqAKBQqGoEEWTwYMDBtM1vysjm46scL9/nYrj/bUnOZ+QSZ9mtfn30OY0LKEcNCstlf0rf+fIpnVIg8TeLQS9oR1e9b3pNLIhdZt7VKodhJSSzJ27SAqbR+b2HWBtjeuQwbhPnoJ9UEvjO0q+BLu/gcNhoMvFtf7DZtGrAoFCoSgXyTnJfH7wc1adW4Wfkx/f9/uebn7dbimSKA/n4zN4f+1J/jodT8Najvw6tQO9mxW/GH1eTjYR61ZycM1y8nJycHBrjV7fAVdfX0JHNKRBcOX6ARmyskhdtYqksPnknT+PZa1a1Hr+edwffggrr7uolLpxAnZ+BceXgbCANuOg64tcPB5LgBl0q0CgUCjuisJk8MGZZORl8HjQ4/yjzT8qnAxOz8nn263R/LLrArZWlrw5pBlTuzTAxurO4SVdfj6RW/5g34olZKWm4ODWHKw74lTLj47DG9K4g3eJ+QNzkBcTS/KCBaQsW4YhLQ27oCDqfPoJzoMGYWFzF0tUXtoDO7+EsxvB2hFCn4HOz4FLnYIGsWbRrwKBQqEwmgupF3h/7/scuH7AZMlgg0Gy7FAMn/xxmoSMXB5s789rg5pS2/nOclCDQc+pndvY9dsC0uJvYO/SABvnQdi7BRAytAHNu/hiWUzgMAdSSrIOHCA5LIz0P7eCEDgP6I/H5CnYtw02/klESjizUQsAV/aCgyf0/jd0eBwcyl4v2RSoQKBQKMokT5/Hz8d+5sdjpp0ZfPhyMu+tOcnRKym0refGz4+E0Kau2x3tpJSci9jPrsXzSLhyCTunOlg7jcHGOZD2gwNo1dMPK5vKqbs35OaStnYtSWHzyT11Cks3NzyfeAL38eOw9vU1viO9Dk4s1wJA3ElwrQuDP4W2k8HGwXwXUAwqEJgIIQSTJk0iLCwM0BxAfX196dSpE2vXrq1ideZnx44duLm50aVLlzv2rVq1irfffhsLCwusrKz46quv6NatWxWoVJSH25PB/+r4L2rZ16pQn3FpOXz8xymWH4rFy9mWzx9sw+gSykFjTh5n+6I5XDtzChuHWlg7DsXaqTkdB9SnTZ+62NhXzsdY/o04khctJGXJb+iTk7Ft3Bif92fgOnw4FnZ3sUB9XhYcWQC7/wspl8GrOYz+AYIeAMuqsbZQgcBEODo6cvz4cbKzs7G3t2fz5s34+flViZaqsKHesWMHnp6exQaCvn37MmLECIQQREZG8tBDDxW6oSqqLyUlgytCrk7Pr7su8s2fZ8nXS57u2Yjn+wTiZHvn+zXu4nl2LprLhSMRWNu5YOPYH2uHIFr3qU+7AfWxc6qcD83sI0dICptP2saNoNfj1KcPHpMn49Cp490lorOT4cBPsHcWZCWAf0ftCaDxQLCoOntrUIHApNy0oR47dmyhDfVN99HMzEz+7//+j2PHjqHT6XjvvfcYOXIkFy9eZPLkyWRmZgLw7bff0qVLF65du8bYsWPJzMxEp9Px/fff0717d5ycnApnEy9dupS1a9cyZ84cpk6dioeHB4cPH6Zdu3Y8++yzPPfcc8THx+Pg4MCPP/5Is2bNmDp1Kvb29pw6dYpLly7x66+/MnfuXPbs2UOnTp2YM2cOoJnkvfvuu+Tm5tKoUSN+/fVXnJycCAgI4JFHHmHNmjXk5+fz+++/Y2dnxy+//IKVlRXz58+/Y2ZxUevszMzMe2ZVp3uV25PBT7R6gqdaP1WhZLCUkq0F5aAXE7Po17w2/x7agoBajne0Tb5+lV1L5nN693YsrR2wduyBlV0wQd3r035IQKX4Acm8PNI2biIpLIycyEgsnJzwmDgR90kTsalbt+wOipJ2Dfb+Dw7+CnkZ0HgAdPsn1OsM1eRv4Z4LBDt+O0PCleJtF8q74letuk50f6jshJiyoS7ZhnrFihW88cYbxMXFFa7ZoKh+FE0GB3sF807nd2js3rhCfZ6Lz2DGmpNsOxNPQy9H5jzagV5N7ywHzUhKZO/yxRzbugmwxMYxFAvrdjTrEkCHoQ0qxQpCl5hI8pIlpCxajC4+HpuAALzf/jduo0Zh4Xhn0CqVhGjY/TUcXQwGnTb00/VF8Llzdb6q5p4LBFWJsqEumdGjRzN69Gi2b9/O22+/zZYtW4w6TlE53J4MfqfzOzzQ+IEKJYPTcvL575azzNl9EXtrS/49tDlTOgfcUQ6ak5HBgdVLObRhNfp8HdYObRBWHQgMaUDH4Q3x8L3LD+BykBMVRdK8MNLWrUPm5eHYvTu+H36AY7duiLsdtrl6WEsAn1wNVrbQbgp0fh48GphHvAm45wJBad/clQ21RlXYUN+kR48enDt3joSEBGrVqljCUWEabkkGNxjMvzpULBlskJLfDlzh042nSMzM46H2dXl1YFO8nG8d0snPzeHwH2vZv/J3crMysXFogYVdKAGtGxE68m8/IHMhdTpsDx3i4k8/kX0wAuHggNvYB3CfNAnbhg3vsjMJF7ZpAeB8ONi6QveXodPT4FT8ZLjqxD0XCKoaZUN9J9HR0TRq1AghBIcOHSIvL0+tSVANSM5JZubBmaw+txo/Jz9m9ZtFV7+uFeoz4lIy7+/J4UJaJO3qufHL1A609ne7pY1ep+P4X5vYs3QxmSlJ2DgEYuPcGb9mgYSOakQdM/sB6VNSSFm2jKQFC3C7eg2dnx+1p03D7YExWLq43F1nBj2cWqsFgKuHwckb+k2HkEfBztU8F2AGzBYIhBB1gXmAD2AAZkspv76tTS9gFXCh4KXlUsoZ5tJUGSgb6jttqJctW8a8efOwtrbG3t6eJUuWqIRxFSKlZNW5VXx+8HOTJYNvpOXw8YZTrDgci5ut4KuHgxkZXOeW37M0GDi9Zwe7fptPyvVrWNv7Y+PUD+9GzQkd2ZC6LczrB5QbHU1S2HxSV69GZmfj0LEjCSNGEPp//3f33v+6XIhcAru+hsRocG8Aw76CNuPB+i5KSasJZrOhFkL4Ar5SykNCCGcgAhglpTxZpE0v4FUp5TBj+1U21NUTZUNtPkyp1dTJ4Fydnp93XuDbrdHo9JInujegtdU1BvXrXdhGSsnFo4fYsWgu8RfPY21XGyy74Fm3JaGjGtEw2MtsAUAaDGRs20Zy2Hwyd+9G2NjgMmI4HpMnY9e06d3f29x0iJgLe/4H6VfBp7VWAdRiJFiYf0JbjbOhllJeA64V/JwuhIgC/ICTpR6oUChMTq4+l5+P/cxPx37CzqriyWApJVui4vhg3UkuJWbRv4U3nRp48Ouui9jXTef9j7fy2sCmdHRMY8eiucScPI6VjRvWDoNwqxNMx+GNaNLRx2x+QPqMDFKXryBpwXzyL13Gytsbr5dewu3hh4z2/r+FzATY9wPsnw05KRDQHUZ+C436VJsS0IpQKQvTCCECgO1AkJQyrcjrvYBlQAxwFe3p4EQxxz8FPAXg7e3dfvHixbfsd3V1JTAwsEwd5S0frSpqkl5TaI2OjiY1NdVEikonIyPjlvkN1ZmKaj2Tc4YliUuI08XR3qE9YzzG4GJ5l2PhRbiaYWBhVB7HE/XUcRRMaG6Dv6MkNjkbg5R428O1G0kYTh3EcP0SwtIBS5tQbF1bUTvICreGYGFpng9Py7g47MPDsd+9B4ucHPIaNiSrT29y27aFYt6fZd1b25w46l5Zie+1zVga8oivFcrlemNId2la4jHmpCLvhd69e1fdwjRCCCe0D/uXigaBAg4B9aWUGUKIIcBK4I7nVCnlbGA2aENDtz8aRUVFGTUsUZOGWqBm6TWFVjs7O9q2bWsiRaVzPwwNFSaDb6zG38mfH0J/oIvfnTO/jSU1O5+vt5xl3p6L2NtY8vawFkzpXB9rSwu6fryV2BRLnPPTeEzsR3/lLFLYYGPXFQePDoQMbkJQLz+szeAHJKUka88ekuaFkbFtG1hZ4TJoEB5TJmPfqvSa/RLv7Y2T2vj/sd+1b/ytx0HXF/DyaspdmEmbHHO9b80aCIQQ1mhBYIGU8o6C86KBQUq5XgjxnRCilpQywZy6FIp7mduTwU+2epKnWj+FnVX5kph6g+T3g1f4bONpkrLyGNehLq8MaEotp7/LQZMTE+mecohWaScwCIGlbXsMdh3Y52DNnP/0MosfkCE7m9RVq0maH0Ze9DksPT2p9cwzuI17GOva5SzZvLxPqwA6swGsHbTyz87Pgqu/acVXM8r87QghagNdgTpANnAcOCilNJRxnAB+BqKklF+U0MYHuCGllEKIjoAFkFhcW4VCUTbnU8/z/p73OXjjIMFewbzb+V0C3cseNi2JgxeTeG/NCY7HphFS3525IzoS5Pd3WWRuVhYH167gkZilWBp0WNi0xNohlIN29uy31eHhYW3yIJB/9SrJCxeS/PtSDKmp2LZoju9HH+EydMjdef/fREo4u1kLAJd3g7079HoDOj5VaTbQVU2JvyEhRG/gdcADOAzEAXbAKKCREGIp8Hkxwz036QpMBo4JIY4UvPYmUA9ASjkLGAs8I4TQoQWZcbIykhYKxT1Grj6Xn479xM/HfsbOyo53O7/LmMZjyp0Mvp6aw0cbolh15Co+LnZ8PS6YEW3+LgfV5eVxdPN69i5fQk5GOjY2TbCw68wxe1faBOew7YwOe2tLXhtomrF0KSXZEREkzQsjfcsWzfu/Xz9t+Kddu/JVHel11L6xDWa9BTeOg4s/DPpYmwlsY/7ZzNWJ0kL1EOBJKeXl23cIIayAYUB/tKGfO5BS7gRK/e1IKb8FvjVabTVG2VCXbEN9kwMHDhAaGsqSJUsYO3ZsJaq7t9l3bR/v732fS2mXGNJgCK91eK3cM4Nz8rVy0P/9FY3OIHm+dyDP9GqEY4E7qEGv5+T2rez6bQEZSQlY2tTHxmU4zbsGk9LAnt/2XKCNrcTPzZ7XBjZlVNuKOfAacnNJW7eepPlh5J6MwsLVFc/HH8N9/His69Qpu4PiyM+Gw/Nh9ze0SLkEtZrCqO8haCxYleOJ4h6gtEAwU0p5o7gdUkodWmJXUYCyoS7Zhhq0qqJp06YxcODAStV1L5OUk8TnBz9n9bmCZHC/8ieDpZRsOnmDD9dFcTkpi4EtvXlrSAvqeToU7o8+sIedi+aRdDUGS2tfrJ3G0qRje80PqI72DXpszwaEh4fzfxN7Veja8uPiSFm8mOTFS9AnJWHbOBCf6dNxHTEcC/tyTnzLToGDP8Pe7yEzHvxCOO43gaAH/lXlNtBVTWmfFkeFEMeARcAyKWXl1PXVYJQNdfE21ADffPMNDzzwAAcOHKi8X8g9ipSSldEr+TziczLzMiucDD57I53pa06yMzqBxrWdmP94J7o1/vuJ4vLxo2xfOJcb585gae2JteNwAtp0JHRUI2rXL38ZanFkR0Zq3v9//AE6HU69euExeRIOnTuXf9JZ+nXY+x0c+AXy0iGwnzYJrH5XErZtu++DAJQeCPyAfsA44CMhxB60oLBaSpldGeLKw19zZhN36Xyx+/Q6PZZWd1++Vrt+Q3pPfarMdsqGungb6tjYWFasWMHWrVtVIKggRZPBbWu35Z3Qd8qdDE7NzuerLWeYt+cSjjaWvDu8BZNCtXJQgBvno9m+cA6Xjx3BwsoFK4cB+LcIpfOoxtRpXI5JWSUg8/NJ27SJ5HlhZB89ioWjI+7jx+ExcSI29euXv+PEc9oqYEcWajbQLUdD15fAt7XJtN8rlBgIpJR6YCOwUQhhAwxGCwpfCyH+lFJOrCSNNQZlQ108L730Ep988kmNmRxXHbmZDP7p2E/YW9lXKBmsN0iWHLjCzE2nSc7KY3zHerzSvwmeBeWgSVdj2LUkjDN7dyEs7bGy74l3o650HtOUeib0A9IlJZHy228kL1yELi4O6/r18H7zTVzHjMayIpP9rh6BXV/ByVVgYQ1tJ0GX/wOPu3QUvY8waiBZSpknhDgJRAHtgRZmVVUBSvvmrmyoNSrbhvrgwYOMGzcOgISEBNavX4+VlRWjRo0q81gFnM4+zczVM7mUdomhDYfyasir5U4GH7iYxHurT3DiahodAtx5d/jf5aDpiQnsWbqI439tBmGJpV0otep1J3RUcxq19UKYyA4i5/RpkubNI23NWs37v0sXfGZMx6lHj7v3/r+JlHBxh1YCem4r2Lpoi8B0egacvU2i+16m1EAghKgHPAyMBxyBxcBIKWVUJWirkSgb6ju5cOFC4c9Tp05l2LBhKggYQWEyOG41dZ3r8kP/H+hSp3zJ4Gup2Xy0/hSrj17F19WO/45vy/DWvgghyE5PY/+qpRzesAa93oClTWvcfHsSOiqIJh29sbCs+Bi61OtJ37qV5HlhZB04gLCzw3X0aDwmT8LWCHuYEjEY4PQ6LQDERoBjbej7LnR4vEbZQFc1pc0j2I2WJ/gdeEpKebCktoq/UTbUd9pQK+6O25PBA1wG8OHwD8uVDM7J1/Pj9vN8F34OvZS80CeQp3s1wsHGivycHCLWr2L/qqXk5+RgYdMcV+8edBoRTItudbC0qngA0KelkbJ0GckLFpAfG4tVHV9qv/Yqbg88gKWbW/k71uXBsd80G4iEM+AeAEO/gOAJYG3+JS3vNUo0nRNC9AS2V7cJXsqGunqibKhNw/mU88zYO4OIGxGFyeCYozF3rVVKycYT1/lgXRQxydkMDvLhzSHNqevhgF6XT+SfG9nz+yKy01OxsG6Eg0cPOgwNoVVv/wr7AYWHh9O5Xj2SwsJIXblK8/4PCcF9ymSc+/RBVKS0OTcDDhXYQKfFgncr6PYStBgFlnffb3V9H5REVdhQdwcigeQSOu0DOEgp7/3ZUgqFmcnV5/Jj5I/8fPxn7K3sea/ze4xuPBoLYUEMMXfV1+nr6Uxfc4Ld5xJp6u3Mwic60SWwFtJgIGrHX+xYFEZ6YhwWVn44eAyl3ZBQgvvVw7aCVhDSYCBz507c/vsN50+eRFhb4zJsGB6TJ2HXooJpxcxE2P+DZgWdkwL1u8Hw/0Jg33vCBrqqKe03fwxYK4TIQXMJjUezmGgMBANbgP+YW6BCca+z99pePtj7QYWTwalZ+Xy55Qxhey/hZGvF9BEtmdipHpYWgvOHDrB9wRwSYy4hLL2wdX2A4P5daT8oAHvnis2m1WdkkrpyJcnz55N38SJWrq54vfgCbg89hFVFlyRNuaJ9+z80F/KzoOlQ7QmgbseK9au4hdLKR1cBq4QQjdF8g3yBNGA+Ws6g2s4lUChqAkk5Scw8MJM159dUKBmsN0gW7b/M55tOk5qdz4RO9Xi5f1M8HG2IOXWC7fPncO1sFMLSDRunoQT17kmHoQ1xcq/Ykop5V66QPH8BKcuWYcjIwK51a+p89hmH7e0I6tevQn0Td6rABvo3bbvVQ1oVUO1mFetXUSxlPgtKKc8CZytBi0JxX3BLMji/YjOD951P5L01J4m6lkbHBh68N7wlLeq4EH/pAsu+nsvFowcRFo5YOfSjebc+dBoRiKuXQ4W0Z+3bp3n///UXWFriMnAgHpMnYX9zrkuRarm75soBrQLo9DrNBrrDk9D5OXCrW/4+FWVSuYY0CsV9TtFkcLva7Xg79O1yzQyOTcnmP+ujWBd5jTqudnw7oS1DW/mSGneDtV//wOnd2xHCBiv7bgR26E/oqKZ4+pV/kpYhJ4fUNWtIDptP7pkzWLq74/mPpzTzN+8K1ulLCdF/agHg0k6wc4Oer2s20I4VHFpSGIUKBApFJVA0Gexg5cD0LtMZFTjqrmcG5+Tr+WHbeb7fFo2U8GLfxjzdsxGGrDT+/Pl7Iv/ciJQCS7sQ6rcZQJcxLfEOKL8fUP61ayQvXETKb7+hT03FtlkzfD/8EJdhQ7EoMhmxXOh1cHIl7PwKbhwD5zow8D/Q7hGwrRnLiN4rqEBgIiwtLWnVqhX5+flYWVnxyCOP8NJLL2FRykzJixcvsnv3biZMmFCJShWVze3J4NdCXsPTvuxvuisPx/LZxtOMq5vOmx/9yYCWPmw+eYPYlGyGtvLljSHN8LKV7F+2gIi1K9HrdFjatsK3aX+6jW2DX9Py+QFJKck+fJiksDDSN20GKXHu2wf3yZNx6NCh4hYT+TlwZIHmA5R8ETwbw8j/aXmA+9QGuqoxZoWyJsD3gLeUMkgI0RoYIaX8wOzqahD29vaFlg5xcXFMmDCB1NRUpk+fXuIxFy9eZOHChSoQ3KNUJBm88nAsbyw/Rna+nngvuJqaw5zdF/F1tWPRk6GE1HXiyB9rWb7sN/JzMrGwbopP0/50fTCE+kGe5fqwNuTlkb5hA0nzwsg5cQILFxc8HnkE9wkTsPE3gaV6TiocuGkDHQd+7WHAB1olkHIArVKMeSL4EXgN+AFAShkphFgIlBoIhBB1gXmAD2AAZkspv76tjQC+RlsEJwuYKqU8dLcXcbfc/KZ1NSWbOiZaQKMotWvXZvbs2XTo0IH33nuPS5cuFWs1/frrrxMVFUVwcDCPPPIIo0ePvqXdp59+Sr+KVl8oKh2DNLAyeiVfRHxBZn4mT7V+iidbPXlXyeDPNp4mO1+zGwk7+/cELwtpwOHCAX78aAHZ6clYWAXgGTCerg92JrBd7XL5Aeni40le8hvJixejT0jApmFDfN59B9eRI7FwKH9iuZD0G5oN9MFfIDcNGvXRbKADuqs5ANUEYwKBg5Ry/23fMMp2GtPavCKlPCSEcAYihBCbpZQni7QZjDYvoTHQCe3Jo5Nx0stH0W9aoCXd3lh+DMCkwaBhw4YYDAbi4uJKtJr++OOPmTlzZuEKZllZWbe0e/jhhzl0yOxxUWFCzqecZ/qe6RyKO0S72u14p/M7NHJrdFd96PQGYlP+rs5u7SE5mgiBWefpdmU/W46mICx9ca0zka4P9qRpJ59y+QFlHz9Bctg80tZvQObn49izBx6Tp+DYpXP5zd+KknQedn8DhxeAIR9ajNRsoOsEV7xvhUkxJhAkCCEaARJACDEWuFbWQVLKazfbSSnThRBRaN5FRQPBSGBegY3FXiGEmxDCt+BYs1D0m9ZNsvP1fLbxtEkDAWhjrQD5+fnFWk3fjrHtFNUPUyWD95xLZPqaE3+/ICV9rC8TFHsQj/x4hIUnDl4P0HlMP4K6+2FpfXf9S52O9C1bSJoXRvahQ1g4OOD28MO4T5yAbYMGd9VXSTiln4el8+DECrCw0vx/urwAnncXEBWVhzGB4DlgNtBMCBELXAAm3c1JhBABQFtg3227/IArRbZjCl4zWyC4mlL8PLiSXi8v58+fx9LSktq1azN9+vRiraZv58svvzSqnaJ6sefqHj7Y+wGX0y8zrOEwXg151ahkcFFikrP4z/oo1h+7jp+bPVO7BPDn9gN0ituD7mIsnhYuSMeB1OranfGTWmFte3d+QLrkZFJ+X0rywoXorl/Hum5dvN94HdcxY7A0hZ+VlHBpF+z8kpDoLWDjrK0BEPosOPtUvH+FWTFmQtl5oJ8QwhGwkFIW7zVcAkIIJ7QF7l+SUqbdvru4UxbTx1PAUwDe3t632DsDuLq6lmiBXBS9Xo+Piy3X0nLv2OfjYmtUH6Vx8/iEhASeeOIJnnzySTIyMoiPj8fPz4/MzEzmz5+PXq8nPT0dCwsLUlJSCo8rqV1NwBRac3Jy7vjdmouMjIwKnytdn86K5BUcyDyAl5UXz9V+jmb6Zhzbd8zoPnL1kvXn81l/IR8BjA60prdbEjd2rGbUlWgQ9lg59sKqUXN8W9ng4ZLCrj07jO7fMjYWh7/Csd+3D5GfT26zpmSPHkVuq1bEWFhAREQ5rrwI0oBn4gHqXV6Ga9pp8qxdOe/3EAkBI9FZO0HEKeBUxc5hRkzxPqhMzKXXmKqhl2/bBkgFIqSUR8o41hotCCyQUha3lFUMUHTKoD9w9fZGUsrZaE8lhISEyNvd96KiooxyvkxPT2fa4Oa35AgA7K0tmTa4eYXcM7Ozs+nevXth+ejkyZN5+eWXsbCw4KWXXuKBBx5g9erVhVbTzs7OdO7cGVtbW7p168bUqVNLbFcTMIX7qJ2dHW3btjWRotKpiIvjzWTw5wc/J0uXVa5ksJSSdceu8Z91UVxNzWdYa19eCK3FmTVLOb35LxDWWDt0pU3/YeR7JNBvUG/j+9brydi2jaR5YWTt3YuwtcV11CjcJ0/CrhQr8btClwfHl2pzABJOg1t9GPo5NsETub5rX41x9Lyf3EdLw5ihoZCCf2sKtocCB4CnhRC/Syk/Le6ggoqgn4EoKeUXJfS9GnheCLEYLUmcas78APydEDZ11dDNBWWKo3HjxkRGRhZuf/TRRwBYW1vz559/3tK2aLs333yzQpoUpudcyjlm7JlRmAx+t/O7NHS7uyUQT15NY/qaE+y7kERzXxc+GdqQ3AObWPvWeqQBrOza0aLnCDqPbomTu53R3wD16emkLl9O0vwF5F+5gpWvL16vvIzb2LFYuZtojeG8TDg0D3Z/C2kxULsljPlJWw+4HDbQiuqBMb85T6CdlDIDQAjxLrAU6AFEAMUGAjSjusnAMSHEkYLX3gTqAUgpZwHr0UpHo9HKRx8t11XcJaPa+pk8May4t8nR5fDjsR/55fgv5U4GJ2Xm8cXm0yzcdxlXe2tmDAkkIPYAEZ9+jT4/F0vblgR2Gk63B9viVtv4ss3cCxdInr+A1BUrMGRlYd+uHbVfeRnnfv0q5v1flKwk2D9bs4HOToJ6XWDYl9C4vyoBvQcw5l1SD8grsp0P1JdSZgsh7hxsL0BKuZPicwBF20i0ZLRCUW2paDJYpzewYN9lvth8hoxcHVM6+tNXRnP453c5kJOBhXVjGnYYSo/xoUb7AUkpydy1m6SweWRu2655/w8ZgvvkydgHtSzvpd5JaoxmAx0xR7OBbjJYs4GuF2q6cyiqHGMCwUK00s5VBdvDgUUFyeOTJR9WuUgpKz71XVFlVLOF8ABIzE7ks4Ofse78Ouo512N2/9l0rtP5rvrYHZ3A9DUnOX0jna6N3JnqlciZtV+wPz0RC6t6+LeaSM+JPfBpYNz6uoasLFJXrSIpbD55589jWasWtZ5/HvdxD2NVq3wL2hdL/GnNBjpyiVYR1OpBzQbau4ILzCiqJcZUDb0vhNiANtQjgKeLrF880ZzijMXOzo7ExEQ8Pcs3tV5RtUgpSUxMrDblsgZpYMXZFXwR8QVZuiz+0fofPNn6SWwtjTdZu5KUxYfrovjjxHX83ez4NMSC5K1ziEy8irD0pnbgZHpP7o9/Mw+j+suLiSV54UJSli7FkJaGXVAQdT79BJdBgxA2JvTniTmouYCeWgdWdhDyOHR5Htzqme4cimqHUQOIUsqDQojLaCuUIYSoJ6W8bFZld4G/vz8xMTHEx8eX2i4nJ6fafNgYQ03SW1GtdnZ2+Pv7m1BR+ahoMjgrT8es8HPM2n4eSyH4Z0uB06FVXDl8HmHhjkfdsfSaPJSA1rXK/NIipSTrwAFcZ/3AuchIEALnAf3xmDwF+7bBpvvSIyWc26oFgIs7NBvoHq9Bp3+AowmfMhTVFmPKR0cAnwN1gDi0nMEpwIQDkRXD2tqaBkbMigwPD6+08kRTUJP01iStxZGjy2F25Gx+PfErDlYOzOgyg5GBI41OBkspWRN5jY/WR3EtNYex9SHo0nYS154gRTjh4j2UHpNG0iTEt0w/IENuLmlr15EUFkbuqVPYODri+cQTuE8Yj7WPCSdnGfRwcpUWAK5HgrMvDPgQ2j8CtjWjbFlhGox5IngfCAW2SCnbCiF6A+PNK0uhqDx2X93NB3s/4Er6FYY3HM4rIa/cVTL4eGwqM9acZP/FJELcdTxhdYiU8AgShR0OHn3oNm4MLbvVK9MPKP9GHMmLFpKy5Df0ycnYNm6Mz/szOOrqSssBAyp6mUVOlANHF2k5gOQL4BkII76F1g+BVQXXGFDUSIwJBPlSykQhhIUQwkJK+ZcQ4hOzK1MozMztyeAfB/xIqK/x1TBJmXnM3HSaRfsvU8c6j9esT5JzaB8pWGLn0oXQMWNp068RVtal20FkHzlCUth80jZuBL0epz598Jg8GYdOHbXhH1PNJM1J0xxA934HGTegTlvoHwbNhoLF3VlWKO4tjAkEKQU2EduBBUKIOIxzH1UoqiUGaWB3+m7eWvlWuZLB+XoD8/de4svNZ9BlZ/IP62hszu4hx6DHxjGYkOFjaT+4BTZ2Jf95ybw80jZuIiksjJzISCycnfGYNAn3iROwqWvi9Xkz4rQ1AA78DLmp0LAXjJkNDXqqOQAKwLhAMBLIBv6JViXkCpS82opCUY0pTAYnlS8ZvPNsAtPXnODC9WTGWETjG7MXgy4HS7sWtOk/ltDRbbFztC7xeF1iIslLlpCyaDG6+HhsAgLwfvvfuI0ahYWjoyku8W+SLhTYQM8HfR60GKHZQPu1M+15FDUeYwLBO1LKaWiLy8wFKBgammZOYQqFKSmaDHa0dmSCxwSmDZpmdDL4cmIWH6w7yZYTV+mVf44h8fsw5GUgbBoS1G8M3R/ugoNLyWWcOVFRJM0LI23dOmReHo7du+P74Qc4dutmGu//olw/pnkAnVgOwhKCx0OXF6FWoGnPo7hnMCYQ9OfOD/3BxbymUFRLiiaDRzQawSshrxC5N9KoIJCVp+O7v84xe/s5WmZE83zKAWROMlj50aTLFHpN7oOzR/Fls1KnI/3PrSSFzSP7YATCwQG3sQ/gPmkStg3vzp+oTKSEy3u0CqCzm8DGCTo/p9lAu9Qx7bkU9xwlBgIhxDPAs0BDIURkkV3OwC5zC1MoKkpidiKfHviU9RfWU9+lPj8N+IlOvsYtgCelZPXRq3y0LgqH62d4Iv0gltlxYOlFQNtH6TN1EO4+xQ/l6FNTSVm6lKQFC9BdvYa1nx+1p03D7YExWLq4mPISwWCAsxu1AHBlHzjUgj7/hg5PgL2JjOYU9zylPREsBDYAHwGvF3k9XUqZZFZVCkUFuH1m8NNtnuaJVk8YnQw+HpvKe6tPcO1UFMPSDmCfFYOwcMW3+cP0e2wEXvWKt4PIjY4mKWw+qatXI7OzcejYEZ8338Spd2+EpYmrcvT5cHyZNgQUHwWu9WDITAieCDYmWGdYcV9RWiCwBNIoxhROCOGhgoGiOhKdHM2MvTM4HHeY9t7teSf0HaOTwYkZuczcdJrNO48wOPkgnbPOg3CkVsAw+j3+IH5N7pxbIA0GMrZvJ3leGJm7dyNsbXEZPgyPyZOxa9rU1JcHeVlwOExLAqdegdotYMyPBTbQJSepFYrSKC0QRPD3amG315hJwMSDnApF+SlMBh//FUcbR2Z0mcGowFFG2TDk6w3M23OJn9YdoOv1fYzPOg3CFrc6fenz6DgatPa94xh9Rgapy1eQtGA++ZcuY+XtjddLL+H28EOm8/4vSlYSHPgJ9s2CrESo1xmGfg6NB6gSUEWFKTEQSClNs5K1QmFmdl/dzft73icmI6YwGexhZ5yZ2/Yz8Xy8fD+Nzu7kgYwTCAROnp3pNWUCTToF3BFI8i5dImnBAlKXLceQmYl9cDC1X3wR5/79EdZm+EaeGkuj6F9g1xbIz4TGA6HbP6H+3bmgKhSlYZTpXIHfUI+CzXAp5VrzSVIojCMhO4HPDnxWrmRwXJaBf/y8E7lnC30yjiKkDnvXYLqOm0DrXs1u8QOSUpK1Zw9J88LI2LYNrKxwGTwIj8mTsW/VykwXdxZ2fQVHl+Bv0EOrsZoNtE+Qec6nuK8xxnTuY6ADsKDgpReFEF2llG+YVZlCUQIGaWD52eV8EfEF2brsu0oGZ+bq+N/mKE6u2UvL1ENYyBxsHJsTOnoC7Ya0wbKIH5AhO5vU1WtInh9G7tloLD09qfXMM7iNexjr2rXNc3GxEVoFUNRazfen/VT2WXQgdPDD5jmfQoFxTwRDgGAppQFACDEXOAyUGgiEEL8Aw4A4KeUdX2OEEL2AVcCFgpeWSylnGK1ccV9yRzK48zs0dC07XSWlZMWhy6wNW0HTG3tpZcjAyrYB7YaOo/OY0Fv8gPKvXiV54UKSf1+KITUV2xbN8f3oI1yGDsHClN7/f4uD8+Gw8wu4sB3sXKH7K9DpaXDyIsdUXkMKRQkYu6CpG3CzSsi4pZRgDvAtMK+UNjuklMOM7E9xH3N7Mvj9ru8zstFIo5LBR68kM/uHZfif30ZLfTIW1r7UChrMw69MLPQDklKSHRFB0rww0rds0bz/+/XDY8pk7Nu1M8+CRwY9RK3RngCuHQEnH+j/PrSfCnYmnm+gUJSCMYHgI+CwEOIvtOqhHpTxNAAgpdwuhAiomDyFAnbH7ub9vXefDE7IyOW/P67B9uB6GuriEJYeNOn2GH0fHca+g7uxsbPCkJdH2rr1JIXNI/dkFBaurng+/hju48djXcdMM3J1uXB0sWYDnXQOPBrB8P9Cm3HKBlpRJYiS1ooVQnwLLJRS7hZC+KLlCQSwT0p53ajOtUCwtpShoWVADHAVeFVKeaKEfp4CngLw9vZuv3jxYmNOfwcZGRk4ORm3OHh1oCbpNYfWNH0aK5JWcDDrIF5WXjzs8TBN7cuuzdcZJNsOxWB9dDeOeVfAwhnn+l0I6N4CG0dtCCj76lVqHYzAfscOLNPT0dXxJat3b7I7dQJzDP8Alros6lzdiH/MKmzzkkl3asTleg8Q7xWqeQKVQE16H0DN0luTtELF9Pbu3TtCShlS3L7SAsGLwDjAF1gCLJJSHrmbE5cRCFwAg5QyQwgxBPhaStm4rD5DQkLkwYMHy2pWLOHh4fTq1atcx1YFNUmvKbUapIFlZ5fxZcSX5OhyeKLVEzze6nGjksHr/jxExIIwHDPPgrDHM7API5+fiLuPNtSSfewYSfPCSF2/HmEw4NSrFx5TJuMQGmq+9a4z4rX6/wM/Qk6qZv/c7Z+aHbQR56xJ7wOoWXprklaomF4hRImBoLR5BF8DXwsh6qMFhF+FEHbAImCxlPJMudT83X9akZ/XCyG+E0LUklImVKRfRc3mbPJZZuyZwZH4I4R4h/B257eNSgYfO3GB1d/8iF3yMRyxwtGvF2NemErtgFrI/HxS160jOWw+2UeOYOHoSFbPHrSeNg2b+vXNdzHJF2H3t9pMYF0uNB8O3V4Cv/bmO6dCUQ7KzBFIKS8BnwCfCCHaAr8A76JZUJQbIYQPcENKKYUQHQELILEifSpqLjm6HH6I/IE5x+fgZONkdDL4xvUE5n/6AzJ2P3aApUcHRj73KA2C6qJLSiJh1g8kL1qE7sYNrOvXw/utt3AdPYodBw+aLwjcOKF5AB1fBsIC2jysrQNQq8wHXoWiSjBmHoE1MAjtqaAvsA0jFqYRQiwCegG1hBAxaMHDGkBKOQsYCzwjhNChLXwzTpY0TqW4p9kVu4sP9n5QmAx+NeRV3O1Kt2nISs9g3syfyTi9DSHzwaklvaZOoX33FuScPs3Vt94ibc1azfu/a1d8pr+HU48erDp6jc++3c+4uum89fFWXhvYlFFt/UxzIZdu2kBvBGtHCH1Gs4F2NVH/CoWZKM2Guj/aIvVDgf3AYuApKWWmMR1LKUtd4F5K+S1aeaniPiUhO4FPD3zKhgsbCHAJ4OcBP9PRt2Opx+Tn5rH0uwXE7t+AMGRhsAskaNQ4Bg0PITM8nEtTPiFr/36EvT2uY0bjMWkStoHagiwrD8fyxvJjZOfroS7EpmTzxvJjAOUPBlLCmZs20HvB3gN6v6XZQDsYZ3OhUFQ1pT0RvIlmRf2qchpVmJLbk8HPtHmmzGSwXq9j07xVnNy8HPSpGKz9qdX7cSY92In0FSs4P/Df5MfGYlXHl9qvvYrb2LFYut465eWzjae1IFCE7Hw9n208ffeBQK/TVgDb+SXEnQTXujD4U2g7CWxMvOSkQmFmSksW965MIYr7g6LJ4A4+HXg79G0auJbsbyilZM/KLexdvhCZF4/BsjYET+HR4a3RLVvC+T5vaN7/ISHUnvYvnPv0QVgV/7a+mpJ9V68XS14WHFkAu/8LKZfBqzmM/gGCHlA20Ioai7EzixWKCpGty2Z25OzCZPAHXT9gRKMRpSaDI//cx1/zf0WXFYO0cCO1wUgmhtbFYe0yEsa+i7C2xmXYMDwmT8KuRYsyNdRxsye2mA/9Om72RlxAsmYDvXcWZCWAf0ftCaDxQDD1msMKRSWjAoHC7BRNBo9sNJJXQl4pNRl8Zv9xNv/8MzkpZ5HCkXjPPowIsMdn63LyVl4k18sLrxdfwO2hh7DyvHOxmJJ4bWDTv3MEBdhbW/LawFImqaVdg73/g4O/Ql6G5v/f7Z/aegBqHQDFPYIKBAqzkZCdwKf7P2XDReOSwZdPnGfjDz+RdiMShC2Jjh3pYZdDg/2LMGzNwKJ1a+p89hkuAwcgyjH792Ye4LONp4F0/NzsS64aSoiG3V9rVhAGnTb00/VF8DGT7bRCUYWoQKCoECsPx/LZxtO3lGOOCPZl6ZmlfBXxFTn6HJ5t8yyPt3ocG8viP7xvXLjKH9//SsKlfYAFGTZBtM24zpCji8HSEqeBAzXztzZtKqx3VFs/RrX1Izw8nP+b2OvOBrGHtHUATq7WfH/aTYHOz4OHWqdJce+iAoGi3BRbjrl2Mz9Gb+By1skyk8HJ15PY8N1crp3eBhjIs2pAq6sXaHRtFZbu7rj94ynN/M3b27wXIiVc2KZVAJ0PB1tX6P5ygQ20mdYdUCiqESoQKMpN0XJMnczDxmsLlp7buZJhzwfdS04GZ6RksvH7+Vw8uglkLtKiLsEXo/FL3oRt06Z4PP8hLsOGYmFrZidOgx5OrdUCwNXD4OQN/aZDyKPamgAKxX2CCgSKcnOz7NLS8Qy/61dgWyuZ/JT25MYNYeSjI+9on5OZy+affufs3jVIQyYWwps2ly7ik7odh9698Zr6IQ4dOpjP/O0mulx8rm2G/70CidHg3gCGfQVtxoO1nXnPrVBUQ1QgUJQbH498ku2WYu16FAtqkXXpSfRZjfC7rRxTl6tj67w1nAhfjkGXjBUetIhNplbGEZxGP0D9xx/Bxr8SbBhy0yFiLuz5H83Sr4JPaxj7K7QYCRYVss5SKGo0KhAo7hqDNLD0zFJ0vl9gpc8hN74fY3268VWW3S3lmHqdnl2/beXQhkXo8+KwxIVm1yX2uVdwGjeBNk9NxMLBwfyCMxNg3w+wfzbkpEBAd442eJI2o/+pSkAVClQgUNwlZ5LPMGPPDI7GH6WjT0dCXZ5kzvUsLMXf5Zgj29Rh/+rd7F0WRn7OFSxwpFGCJemGHLLHP0r/x0dhY10Jb72Uy5oN9KF5oMuGZsO0OQD+ISSHh6sgoFAUoAKBwiiyddn8cPQH5p6Yi5ONEx92+5DhDYcjhODJztqCGc9P6EnklsP877FPyM0+j8CWOmn2XBZOnHjwUf4xqTe1nCphKcYbJ7VlII/9rn3Ytx4HXV8Ar7JXN1Mo7kdUIFCUyc7YnXyw9wNiM2IZFTiKl9u/fMfM4JSo68ya+wpZWWcAazxz3Dhq7cOBAQN588GOtPKvhCqcy/u0CqAzGzQb6E5PQ+dnwdXf/OdWKGowKhAoSiQhO4FP9n/CHxf/IMAlgF8G/kIHnw63tDn7x162Ll5GRvZpQOAsvdlm14hzrUJ4fWhLRgbXMW8VkJRwdrMWAC7v1myge70JHZ9UNtAKhZGoQKC4g5vJ4MKZwcHP8njQ3zODpV7PhWWb2LpmK6l50YAOOwtvNnu046RbA57s0YDZvQJxtDXj20uvgxMrCmygT4CLPwz6BNpNVjbQCsVdYra/VCHEL8AwIK6ExesF8DUwBMgCpkopD5lLj8I4bk8G/zv034Uzg/Xp6Vyat4xtOyNJkpdB5uDoUJetPl3YLz3o38KbzUObU9/TjB/E+dlweD7s/gZSLkGtpjDqewgaC1Z37z+kUCjM+0QwB20Fsnkl7B8MNC741wn4vuB/RRWQrctm1tFZzDsx745kcO6FC8TM+Y1dx6+SYH0DKTNwcKrLsSZ9WZ9sT6CXE6/Wy+f5sSFmFJgCB3+Gvd9DZjz4hcCgj6DJYGUDrVBUELMFAinldiFEQClNRgLzCtYp3iuEcBNC+Eopr5lLk6J4iksGu9m6kblzF1fnLeLgDbjumIa0Ssbe0Zfktg/z/VVb7LMteWdYEyZ3rs+uHdvNIy79Ouz9Dg78AnnpENhPKwGt31WVfyoUJqIqcwR+wJUi2zEFr6lAUEnEZ8Xz6YFPb0kGt3dpSeqKVZye/xuR+BDrkovBPg5bu1rYdX2M7686kRSbz7gO/rw6oCme5ioHTTynrQJ2ZKFmA91yNHR9CXxbm+d8CsV9jNC+kJupc+2JYG0JOYJ1wEdSyp0F238C/5JSRhTT9ingKQBvb+/2ixcvLpeejIwMnJycynVsVWAuvQZpYHfGblYnryZf5jPAdQADdW1xCd+Fzd79nKvVhkseFuj1MVhYO2PTPJTfCeRiOjR2s2BicxsCXG+1ZDCVVqf0c9S7vAyv+D1IYck1375cqTuKHHvfCvddlJr0XqhJWqFm6a1JWqFienv37h0hpSx2/LYqnwhigLpFtv2Bq8U1lFLOBmYDhISEyF69epXrhOHh4ZT32KrAHHrPJJ9h+p7pRCZF0sm7I2/YjsR2+WZSt/6H834dON8sFF3+OSyEA80GjGOTXRNWHYvHx8WOr8c1Y0Sb4stBK6RVSri4Q6sAOrcVbF2g24uITs/g5+yNOVyIatJ7oSZphZqltyZpBfPprcpAsBp4XgixGC1JnKryA+YjKz+LWZFaMthDOPNN1mjqzYok59TrnGvYjTMhg8jNiULok2necyTnGoby6q5Y9DKR53sH8mzvRjjYmPjtYjDA6XVaAIiNAMfa0O89CHlM2UArFJWIOctHFwG9gFpCiBjgXcAaQEo5C1iPVjoajVY++qi5tNzv7IjZwYf7PiTrWgxvXWhC693X0Cf/zpVWAzjZtTVZGUch9zqBHfth2W0QH22L5Ur4ZQa29OatIS2o52liYzhdHhz7DXZ+BYlnwT0Ahn0JbSYoG2iFogowZ9XQ+DL2S+A5c51foSWDPznwCed3/cGUo/YEHwcMp0jp9RBHsSU96QBknKNuUFeajHmImXsS2LkimibeTix4ohNdA2uZVlBuBhyaqxnBpV8F71Yw9hdoPhIs1dxGhaKqUH999yAGaWDp8UXsXvgFfffl8PhVAxZOkPvAU0SkQ9K1nSCz8G4UTJcpjzDvTB5vLDiNo40l7w1vwaTQ+lhZmrA2PzMR9v+gWUHnpED9bjDiGwjsq0pAFYpqgAoE9xinovfx17dv0HrnNf6RAaKeH5YvPs6+mGzizv6JNKTiXqcxfZ94gl0Zzoz9/TTJWXmM71iPV/o3MW05aMoV2PM/7SkgPwuaDoVuL0HdjqY7h0KhqDAqENwjpBw7zIFv3qP2rjP00kNmuyY4PvwSe4/HE7NjPVKfgJOHP30ef4kkj4Y8u/YkJ65epGOAB++OaEHLOiZMzsadKrCB/k3bbvUQdH0Rajcz3TkUCoXJUIGgBiN1OtL/3MrFn77F+thZalnDxR6NaPPIdM7svs6O3xZh0F3FzqkWPSa9gkebjny88Qxrju7F19WOb8a3ZVhrX9O5g145oFUAnV4H1g7Q4Uno/By41S37WIVCUWWoQFAD0aemkrJ0KQlhYRiu3yDZFfYPrUWXSTOw2Qcrvp2DIf8C1nYudBv/D1r0HcAvuy/z3Zc7MEjJC30CebqXicpBpYToPwk+/C6EHwd7d+j5OnR8Chw9K96/QqEwOyoQ1CByo6NJCptP6urVyOxsoupb8cdYG9oNfZbg00Ec/HIlutwoLK3t6TR6Eh1HjWRrdAoD/7uLmORsBgf58OaQ5tT1MEE5qF4HJ1dqJaA3jmFv4wkDP4J2U8C25szUVCgUKhBUe6TBQMb27STPCyNz926wseZwsAsLg/Ko07Iro5InED17E5FZ6xAWlrTpP4pu48ZxKUMyNewou88l0tTbmYVPdKKLKcpB83PgyALNByj5ItRqAiO/Y29ybXp27l/x/hUKRaWjAkE1RZ+Rif3WrZz76CPyL13GsnZtTj8Uwud1IrFytuEx/QekrjhBVPoXIAw07dKHXpMno7Nz5qPNZ5i/7zJOtlZMH9GSiZ3qVbwcNCcVDty0gY4Dv/Yw4ENoOgQsLJDh4Sa5boVCUfmoQFDNyLt8maT580ldthyXzEysgoNJnDSA6babiM08xjj9U7juyCYuJQxkLgHBXejz6FRcavuycP9lvtgUQWp2PhM61eOV/k1xd6zgYi3pNzQb6IO/QG4aNOqr2UAHdFNzABSKewQVCKoBUkqy9uwhaV4YGdu2gZUVLoMGcbxpPTYFXmDzxbl0vjGQYaf8yU7cTrbMxLdxa/o+/jjeDRqx73wi732zk6hraXRq4MG7w1vSoo5LxUQlnddWATu8AAz50GKUNgfAt40pLlmhUFQjVCCoQgzZ2aSuXkPy/DByz0Zj6elJrWeeweXhB1mRHM7n+z6n7uHmPHH+CXQJ+8k2ROHhH0jfxx6jXsvWxKZk89zCQ6yLvEYdVzu+ndCWoa0qWA56LRJ2faWtB2xhBcETocv/gWcjk123QqGoXqhAUAXkX71K8sKFJP++FENqKrYtmuP70Ue4DB3C2YwLvLT7FZKic3ggejw2ScfJ12/CuZYffR55nkYdOpGrM/DVljPM2nYOKeHFvo15umcj7G0syz55cUgJl3ZpcwCit4CNM3R5AUKfAWcf0168QqGodqhAUElIKcmOiCApbD7pW7aAlDj374/HlMnYt2tHti6br45+y8a9O+ga3RXn+HMYdFuxc/akx8SXaNGzN0JYsOH4dT5cF0VsSjZDW/nyxpBm+LuXsxzUYIAzG7QAEHMAHL2g77uaDbS9m0mvX6FQVF9UIDAzhrw80tatJylsHrkno7BwdcXz0am4T5iAdZ06AGyP2c43m34k8HgLBsXVw5AfjqWdM74dezP22Rewsrbm1PU03lt9gr3nk2jm48yiJ0Pp3KicE7Z0eXB8qTYHIOE0uNWHoZ9rw0DW9qa7eIVCUSNQgcBM5MfFkbJ4CclLlqBPTMS2cSA+06fjOmI4Fvbah21cVhxfbPkW/U5bul6rgz5vN8LaltAHJtBh+Ch279tPep7ki3XHWbDvEi721rw/KojxHeqWrxw0LxMOzdNsoNNiwDsIHvhZSwQrG2iF4r5F/fWbmOxjx0iaF0baH3+ATodTr154TJ6EQ+fOhUlcvUHP4gPLOLz6LH5XJIbc/RgsBMEDh9N57MM4uLii0xv483I+L20PJy07n0mh9Xm5fxPcHMpRDpqVBPtnw75ZkJ0M9bvC8K8gsJ8qAVUoFCoQmAKZn0/apk0kh80n+8gRLBwdcR8/Do+JE7GpX/+WtpGXTvLb/I14ns3AN/sQBnQ069aH7uMn4lKrNgB7ziUyfc0JTl3Po3NDT94d0YJmPuUoB02N0WygI+YU2EAPga4vQb1OFb9ohUJxz2DWQCCEGAR8DVgCP0kpP75tfy9gFXCh4KXlUsoZ5tRkSnTJyaQs+Y3kRYvQ3biBdf16eL/1Fq6jR2HpdKvfTnJKGnPnr0ZGxOCeHYFBZtOgbSd6TpqKp7/mzhmTnMV/1kex/th1/NzseS7Yllcf7nT35aDxpzUb6Mgl2narBwtsoJub4rIVCsU9hjnXLLYE/gf0B2KAA0KI1VLKk7c13SGlHGYuHeYg5/RpkubNI23NWmReHo5duuAz/T2cevRAWNw6dp+Xo2Pl8u3E/HkEmXUADOnUDmxGv0efxDewKQDZeXq+33aOH7adQwj4Z78m/KNnQ/bu2nF3QSDmoFYBdGodWNlBhycKbKDrmfLyFQrFPYY5nwg6AtFSyvMAQojFwEjg9kBQI5B6PRl//UXSvDCy9u9H2NvjOmY0HpMmYRsYeEd7XZ6ePZtPEbE6HEPafqQhEQcfX4Y8Po36rYO1PqVk3bFr/GddFFdTcxjW2pc3hjTHz+0uKnekhHNbtQBwcQfYuUHPf0HHfygbaIVCYRRCW0PeDB0LMRYYJKV8omB7MtBJSvl8kTa9gGVoTwxXgVellCeK6esp4CkAb2/v9osXLy6XpoyMDJyc7s4iWWRlYb9rNw7h4VgmJqL38CCrV0+yu3ZFOjre0V4aJEnnJFcOnceQuhepv47e0Z6GXfvg2bBZ4Tf8y2l6FkTlcTrZQF1nCyY1t6Gpx60TwkrVK/V4xe+h3uVlOGecJ9fGkyt1R3LNdwB6q8ovAS3Pva1KapLemqQVapbemqQVKqa3d+/eEVLKkOL2mfOJoLgxjdujziGgvpQyQwgxBFgJNL7jIClnA7MBQkJCZK9evcolKDw8HGOPzT1/nuT580lZuQqZlYVDSAju776Dc58+CKs7b5vBIDl74AbbfttG5vXtGHSX0NlZEzpxIt0HPYSFpfYhn5SZx+ebTrNo/2Vc7a35YFQLxnesh6XFnberWL35OXB0kZYDSL4Ano2h7/+wbfUQgVY23PlsUjnczb2tDtQkvTVJK9QsvTVJK5hPrzkDQQxQdI1Cf7Rv/YVIKdOK/LxeCPGdEKKWlDLBjLpKRBoMZO7cSdK8MDJ37kTY2OAybBgekydh17z4RKuUkgtHE9jx2z6SLm/FkH8GnZXAd3Bnxo1/BRtbOwB0egML9l3mi81nyMjVMaVzAP/s1wRXB2vjxOWkaQ6ge7+DjBtQpx0MeF9bEN6ighbTCoXivsacgeAA0FgI0QCIBcYBE4o2EEL4ADeklFII0RGwABLNqKlYDJmZpKxYSfL8+eRdvIiVlxdeL76A20MPYeVZ/Di7lJKYqGR2/n6YG9Fb0OUdx2Ah0XXw5cnHZ+Dl7lvYdnd0AtPXnOT0jXS6BnryzrCWNPVxNk5cRpy2BsCBnyE3FRr1gQd+goDuag6AQqEwCWYLBFJKnRDieWAjWvnoL1LKE0KIpwv2zwLGAs8IIXRANjBOmitpUQx5V66QPH8BKcuWYcjIwK5Na+rMnInLgP4Im5Inbl07l8ruZce4cnwL+tzDGISe64EWPDj1NUIb9yhsdyUpiw/XRfHHiev4u9sza1J7Brb0Nq4SKOkCjc/Mgh1bQZ8HLUdpcwDqBFf4uhUKhaIoZp1HIKVcD6y/7bVZRX7+FvjWnBqK0UTWvn0khc0nY+tWsLTEZeBAzfytTele+/FX0tmz4hQXIragyz2AlHlc8sshaNQwXur6HNaW2jBPVp6OWeHnmLX9PJZC8OqAJjzRvSF21ka4g14/pnkAnViOLxbQbpLmBKpsoBUKhZm4b2YWG3JysN+5kwtffEnumTNYurvj+fQ/cB83Dmtv71KPTb6eyd5VZzmzZyv63L1IQxYxXjnkd6vLv4ZMp56LVqcvpWRN5DU+Wh/FtdQcRrSpwxtDmuHrWkYVj5RwaXeBDfRmsHGCzs+zVwbTZeADproFCoVCUSz3TSBIW7cel/kLoFkzfD/8EJdhQ7GwtS39mMRsDqw9z4lt4ehydiP1qcS75XO6TT5PDnqVwQ0GFw7zHI9NZcaak+y/mETLOi78d3xbOgR4lC7KYIAzfxTYQO8Hh1rQ523o8DjYu5On1gFWKBSVwH0TCFyGDuFEUiJdnniizDH6zNRcDm64yLE/d5KXuROpjyfDBfY0uUFolyHMCXkZV1tXQCsHnVlQDuruYMNHY1rxUEjdYstBC9Hnw7Gl2kpg8ae0mb9DZkLbScoGWqFQVDr3TSCwsLMjv3HjUoNATmY+hzdd5vDGPeSm78Cgi0XnZMOuwHgsmvnwny4/0M67HQD5egPz917iy81nyMzTM7VLAC/1LaMcNC8LDodpawGnXoHaLWHMT9BytLKBVigUVYb69EHzA4rcGsPB9QfJTtmOIf88Fg72RLTIILpuBk+1fZpHWjxSmAzeeTaB6WtOcDYug+6Na/HOsBY09i6lHDQrCQ78pNlAZyVCvS4w9Ato3F+VgCoUiirnvg4Eunw9J7ZfZf/aI2QkbMeQdwpLW1tigm0Ir32ajvU6s6zTv6nros2Lu5yYxQfrTrLp5A3qeTgwe3J7+rcopRw0NVabAHbwV8jPhCaDodtLUC+08i5SoVAoyuC+DAR6vYFTu6+xf/UJUm9sR593DAtLC2QHfxa67cXe2YUPO3xcmAzOzNXxffg5Zu84j5WF4LWBTXm8W4OSy0Hjz8Dur+HoEpCGv22gvVtU7oUqFAqFEdxXgUBKyZn919m7Moqk2J0Y8g4Berw6teL3Wge4aNjJ2CZjeandS7jauiKlZNWRWD5af4rraTmMCq7D64Ob4+NqV/wJYiO0CqCotWBlCyGPQufnwb1+8e0VCoWiGnDfBIJr0SlEr88jMn4e+rwDSH029Tt0YH9gPL+mribQJZB5nefRtnZbQCsHfW/1CQ5eSqaVnyv/m9iW9vWLKQeVEs7/pQWAC9vBzhV6vKrZQDt5VfJVKhQKxd1z3wSCG+dPkB77K1KfQf3WbckIrc3MGwvQZeh4sd2LhcnghIxcZm48zZKDV/BwsOGTB1rxYPu6WNxeDmrQQ9RqLQBcOwrOvjDgA2g/FWyN9BFSKBSKasB9Ewj8m/vjUMuFoPGT+F/iAo7HrKRLnS78uyAZnK838PPOC3y15QzZeXoe69qAF/o2xtX+tnJQXe7fNtBJ58EzEEZ8A60f1oaDFAqFooZx3wQCBx8vonrbM+vsv3GzdePTHp8yKGAQQgi2n4lnxtqTRMdl0KOJF+8Ma05g7du+1eekQcSvsOc7yLgOvsHw0DxoNgwsjPAQUigUimrKfRMINl3axNa0rTzY5EFebPcirrauXErM5IN1UWw+eYP6ng78NCWEvs1r31oOmhEP+76H/T9pNtANe8GYH6BBTzUHQKFQ3BPcN4FgRKMRpJ9LZ3LnyWTm6vj0j1P8tOMC1paCaYOa8Vi3AGytinyzT76ozQA+PF8bDmo+HLr9E/zaVdk1KBQKhTm4bwKBhbDA38afFYdj+HjDKW6k5TKmnR/TBjXD26VIOej145oH0PHlICygzThtDkCtO1bQVCgUinuC+yYQnLyaxof7cohOOUprf1e+n9SedvXc/25w0wb67CbNBjr0Gej8HLjUqTrRCoVCUQmYNRAIIQYBX6OtUPaTlPLj2/aLgv1DgCxgqpTykDm0pGbnE5cl+XRsa8a289fKQQ0G7YN/55dwZS84eELvf2s20A5lWEgrFArFPYLZAoEQwhL4H9AfbSH7A0KI1VLKk0WaDQYaF/zrBHxf8L/J6dzIk5k97RkQUlezgT66XBsCijsJrvVg8GeaDbSNgzlOr1AoFNUWcz4RdASipZTnAYQQi4GRQNFAMBKYV7BO8V4hhJsQwldKec0cguzIg32zC2ygL4NXcxg9G4LGgGUp9tEKhUJxD2POQOAHXCmyHcOd3/aLa+MHmD4QnNlI6N4nID8N6naCIZ9B4wFgYWHyUykUCkVNwpyBoLgie1mONgghngKeAvD29ia8HEs42mclEODQkKsNHibVrYUWaq5tv+t+KpOMjIxyXWtVUJO0Qs3SW5O0Qs3SW5O0gvn0mjMQxAB1i2z7A1fL0QYp5WxgNkBISIjs1atXuQSFO/hR3mOrgvDw8BqjtyZphZqltyZphZqltyZpBfPpNee4yAGgsRCigRDCBhgHrL6tzWpgitAIBVLNlR9QKBQKRfGY7YlASqkTQjwPbEQrH/1FSnlCCPF0wf5ZwHq00tFotPLRR82lR6FQKBTFY9Z5BFLK9Wgf9kVfm1XkZwk8Z04NCoVCoSgdVTKjUCgU9zkqECgUCsV9jgoECoVCcZ+jAoFCoVDc56hAoFAoFPc5QivcqTkIIeKBS+U8vBaQYEI55qYm6a1JWqFm6a1JWqFm6a1JWqFieutLKb2K21HjAkFFEEIclFKGVLUOY6lJemuSVqhZemuSVqhZemuSVjCfXjU0pFAoFPc5KhAoFArFfc79FghmV7WAu6Qm6a1JWqFm6a1JWqFm6a1JWsFMeu+rHIFCoVAo7uR+eyJQKBQKxW3ck4FACDFICHFaCBEthHi9mP1CCPHfgv2RQoh2VaGziJ6y9DYTQuwRQuQKIV6tCo1FtJSldWLBPY0UQuwWQrSpCp0FWsrSOrJA5xEhxEEhRLeq0FlET6l6i7TrIITQCyHGVqa+2zSUdW97CSFSC+7tESHEO1Whs4ieMu9tgeYjQogTQohtla2xiI6y7u1rRe7r8YL3gkeFTiqlvKf+oVlenwMaAjbAUaDFbW2GABvQVkgLBfZVc721gQ7Ah8Cr1VxrF8C94OfBVXVvjdTqxN/Do62BU9X53hZptxXN1XdsddUK9ALWVtX9LIdeN7T11OsVbNeurlpvaz8c2FrR896LTwQdgWgp5XkpZR6wGBh5W5uRwDypsRdwE0L4VrbQAsrUK6WMk1IeAPKrQmARjNG6W0qZXLC5F23VuarAGK0ZsuCvCXCkmGVSKxFj3rcA/wcsA+IqU9xtGKu1umCM3gnAcinlZdD+5ipZ403u9t6OBxZV9KT3YiDwA64U2Y4peO1u21QW1UlLWdyt1sfRnryqAqO0CiFGCyFOAeuAxypJW3GUqVcI4QeMBmZRtRj7PugshDgqhNgghGhZOdKKxRi9TQB3IUS4ECJCCDGl0tTditF/Y0IIB2AQ2heDCmHWhWmqCFHMa7d/0zOmTWVRnbSUhdFahRC90QJBVY27G6VVSrkCWCGE6AG8D/Qzt7ASMEbvV8A0KaVeiOKaVxrGaD2EZmmQIYQYAqwEGptbWAkYo9cKaA/0BeyBPUKIvVLKM+YWdxt383kwHNglpUyq6EnvxUAQA9Qtsu0PXC1Hm8qiOmkpC6O0CiFaAz8Bg6WUiZWk7Xbu6r5KKbcLIRoJIWpJKavCe8YYvSHA4oIgUAsYIoTQSSlXVorCvylTq5QyrcjP64UQ31XzexsDJEgpM4FMIcR2oA1Q2YHgbt634zDBsBBwTyaLrYDzQAP+Tra0vK3NUG5NFu+vznqLtH2Pqk0WG3Nv66GtQd2lBrwPAvk7WdwOiL25XR313tZ+DlWXLDbm3voUubcdgcvV+d4CzYE/C9o6AMeBoOqotaCdK5AEOJrivPfcE4GUUieEeB7YiJaB/0VKeUII8XTB/lloFRdD0D6wsoBHq7NeIYQPcBBwAQxCiJfQKgnSSuq3qrQC7wCewHcF31x1sgpMvYzU+gAwRQiRD2QDD8uCv7JqqrdaYKTWscAzQggd2r0dV53vrZQySgjxBxAJGICfpJTHq6PWgqajgU1Se4KpMGpmsUKhUNzn3ItVQwqFQqG4C1QgUCgUivscFQgUCoXiPkcFAoVCobjPUYFAoVAo7nNUIFDc8wghfIQQi4UQ54QQJ4UQ64UQTSpZQy8hRJdS9o+63aGzwJ5h0W2vzRRC9DGXTsX9iQoEinsaoU1mWAGESykbSSlbAG8C3nfRh2Vp20bSC82ZtST+BXxX5BzN0f4+ewghHIu0+wYo0aJaoSgPKhAo7nV6A/lFJ2RJKY9IKXcUfEtfe/N1IcS3QoipBT9fFEK8I4TYCTxYzPYAoa0RcUgI8bsQwqnIcdMLXj8mtLUkAoCngX8WeMh3Lyqw4OkkV95qvzABCAM2ASOKaL8EeBZMMlQoTIIKBIp7nSAgopzH5kgpu0kpFxfdBrYA/wb6SSnboc36frnIcQkFr3+PZglyEc0x9EspZbCUcsdt5+mKZtJWlIeBJWheMuNv23eo4BiFwiTccxYTCoUJWVLCdijQAthVYKNhA+wp0m55wf8RwBgjzuMLxN/cEEJ0AOKllJeEEDHAL0IId/n3Og9xQJ27uRCFojRUIFDc65xA870pDh23PhXb3bb/dh+Xm9sC2CylvP2b+k1yC/7XY9zfWDaaidhNxgPNhBAXC7Zd0HyRfiqiM9uIfhUKo1BDQ4p7na2ArRDiyZsvCG3N357AJaCFEMJWCOGK5kVvDHuBrkKIwIL+HIyoQkoHnEvYF4XmhIoQwgJ4EGgtpQyQUgagrVBVNOg0QXPHVChMggoEinuaAsfL0UD/gvLRE2h23lellFeA39AcJxcAh43sMx6YCiwSQkSiBYZmZRy2BhhdXLIY2A60Lahw6gHESiljb9vfQgjhK4SwRgsaB43RqlAYg3IfVSiqAUKIr4E1UsotZbQbDbSTUr5dOcoU9wPqiUChqB78B21BlLKwAj43sxbFfYZ6IlAoFIr7HPVEoFAoFPc5KhAoFArFfY4KBAqFQnGfowKBQqFQ3OeoQKBQKBT3OSoQKBQKxX3O/wOiTIRwjcnI7AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(I, V, label='Data')\n", "plt.plot(I_line, V_line, label='Batch Solution')\n", "plt.xlabel('Current (A)')\n", "plt.ylabel('Voltage (V)')\n", "plt.grid(True)\n", "\n", "I_line = np.arange(0, 0.8, 0.1).reshape(8, 1)\n", "\n", "for k in range(num_meas):\n", " V_line = x_hist[k, 0]*I_line + x_hist[k, 1]\n", " plt.plot(I_line, V_line, label='Measurement {}'.format(k))\n", "\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resistance estimate $\\hat{R}$ should approach the true resistance value of $R = 5~\\Omega$ very closely (i.e., to within a few hundredths of ohms). As expected, the offset term $\\hat{b}$ should be small as well (less than 0.1 ohms). Try modifying the initialization (e.g., the intial uncertainty of the prior guess) - can you get a better final esimate?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 2 }