{ "cells": [ { "cell_type": "markdown", "id": "54086665-879b-438f-af50-2ee8ccb39d49", "metadata": {}, "source": [ "## MT09 - TP5 - Automne 2024\n", "### Problèmes aux moindres carrés\n" ] }, { "cell_type": "markdown", "id": "e5dd0ce4-cb08-4aba-a9ac-bd51b794385e", "metadata": {}, "source": [ "#### Exercice 1 - Approximation affine par morceaux\n", "\n", "On considère le nuage de $m$ points $(x_i,y_i)_{i=1,...,m}$ suivants donnés par les tableaux $x$ et $y$ du script code python ci-dessous :" ] }, { "cell_type": "code", "execution_count": 4, "id": "1d688136-b21c-44fb-a51f-a99f4f62a3f6", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGfCAYAAACX9jKsAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJ+RJREFUeJzt3X9wVfWd//HXTQi50CbXDZjcSwmY0i9KiLASNhB38UfdxIRuRlvnO7oulPZrnQG1Xc2XcRfdnRBHl51px2FbS6i7ulSzfvUPqjUjjWSmBewAjSFkJQapdlNx8YYUUu8NMAmXm/v9g72Ry70JN+Gec/K59/mYyYz35HPPfedtnLw85/P5HFckEokIAADAEFlOFwAAADARhBcAAGAUwgsAADAK4QUAABiF8AIAAIxCeAEAAEYhvAAAAKMQXgAAgFEILwAAwCiEFwAAYJRpVp583759+v73v69Dhw7J7/fr9ddf19133z3m+D179uj222+PO3706FHdcMMNSX3myMiIPv30U+Xl5cnlck22dAAAYKNIJKLBwUHNmTNHWVnjX1uxNLycPXtWS5cu1be//W3dc889Sb/v2LFjys/PH3197bXXJv3eTz/9VMXFxROqEwAATA2ffPKJ5s6dO+4YS8NLbW2tamtrJ/y+wsJCXXPNNZP6zLy8PEkXf/hLA9DVCoVC2r17t6qrq5WTk5Oy8yIWfbYPvbYHfbYPvbaHVX0OBoMqLi4e/Ts+HkvDy2TddNNNGhoaUmlpqf7hH/4h4a2kqOHhYQ0PD4++HhwclCTNmDFDM2bMSFlN06ZN08yZMzVjxgz+o7AQfbYPvbYHfbYPvbaHVX0OhUKSlNSUD1ckEomk7JPH+yCX64pzXo4dO6Z9+/apvLxcw8PDevnll7V9+3bt2bNHt9xyS8L3bN68WY2NjXHHX3nlFc2cOTNV5QMAAAudO3dO999/vwKBwBXvnEyp8JJIXV2dXC6X3nzzzYTfv/zKS/Sy06lTp1J+26itrU1VVVUkegvRZ/vQa3vQZ/vQa3tY1edgMKjZs2cnFV6m5G2jS61cuVLNzc1jfj83N1e5ublxx3Nyciz55bXqvIhFn+1Dr+1Bn+1Dr+2R6j5P5FxTfp+Xw4cPy+fzOV0GAACYIiy98nLmzBl99NFHo697e3vV1dWlgoICzZs3T5s2bdKJEyf00ksvSZK2bt2q6667TosXL9b58+fV3NysnTt3aufOnVaWCQAADGJpeOno6IhZKVRfXy9JWrdunXbs2CG/36/jx4+Pfv/8+fPauHGjTpw4oRkzZmjx4sV66623tHr1aivLBAAABrE0vNx2220abz7wjh07Yl4//vjjevzxx60sCQAAGG7Kz3kBAAC41JRfbQQg/YVHImrvHVD/4JAK89yqKClQdhbPJgOQGOEFgKNau/1qbOmRPzA0eszncauhrlQ1Zaw0BBCP20YAHNPa7deG5s6Y4CJJfYEhbWjuVGu336HKAExlhBcAjgiPRNTY0qNEU/qjxxpbehQesWUTcAAGIbwAcER770DcFZdLRST5A0Nq7x2wrygARiC8AHBE/+DYwWUy4wBkDsILAEcU5rlTOg5A5iC8AHBERUmBfB63xloQ7dLFVUcVJQV2lgXAAIQXAI7IznKpoa5UkuICTPR1Q10p+70AiEN4AeCYmjKfmtYsk9cTe2vI63Grac0y9nkBkBCb1AFwVE2ZT1WlXnbYBZA0wguQAmxvf3Wys1yqXDDL6TIAGILwAlwltrcHAHsx5wW4CmxvDwD2I7wAk8T29gDgDMILMElsbw8AziC8AJPE9vYA4AzCCzBJbG8PAM4gvACTxPb2AOAMwgswSWxvDwDOILwAV4Ht7QHAfmxSl4HYDTa12N4eAOxFeMkw7AZrDba3BwD7cNsog7AbLAAgHRBeMgS7wQIA0gXhJUOwG2xywiMRHfjdaf2864QO/O40YQ4ApiDmvGQIdoO9MuYDAYAZuPKSIdgNdnzMBwIAcxBeMgS7wY6N+UAAYBbCS4ZgN9ixMR8IAMxCeMkg7AabGPOBAMAsTNjNMOwGG4/5QABgFsJLBmI32FjR+UB9gaGE815cunh1KhPnAwHAVMRtI0xYuu2FwnwgADALV14wIem6F0p0PtDlP5s3DX42AEg3hBckLboXyuXXWaJ7oZg+6Zf5QABgBsILknKlvVBcurgXSlWp1+g/9swHAoCpjzkvSAp7oQAApgrCC5LCXigAgKmC8IKksBcKAGCqILwgKTwbCQAwVRBekBT2QgEATBWEFySNZyMBAKYClkpjQtgLBQDgNMILJoy9UAAATuK2EQAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhccDABlsJCL9pndAp89d4DlVAIxBeAEy1Nvvn1RjZ7Y+O9gxeszncauhrjTtnxAeHonwcFHAYIQXIAO1dvv13Vf/U5HLjvcFhrShuVNNa5albYBp7farsaVH/sDQ6LFMCW1AumDOC5BhwiMRNbb0/E9wib3aEA0zjS09Co9cHm3M19rt14bmzpjgIn0e2lq7/Q5VBmAiCC9AhmnvHYj7432piCR/YEjtvQP2FWWD2NAWK91DG5BuCC9AhukfHDu4TGacKTI1tAHpiPACZJjCPHdKx5kiU0MbkI6YsAtHXb7q46a5eU6XlPYqSgrk87jVFxhKeAvFJcnrubgCJ51kamgD0hHhBY5JtOrDm5+r1V6XVjtYV7rLznKpoa5UG5o7dfFmyeeTdqP/1FBXmnZLhzM1tAHpiNtGcMRYqz5OBof14m+z9Pb7Jx2qLDPUlPn0o/uW6prpsce9HnfaLpOOhjbp8jVW6R3agHTElRfYLplVH8/84gPVLvkSf0gsdOfiIoV+H9a1pSszZofdmjKfmtYsi7/ixz4vgFEIL7DdlVZ9SC75A8Nq7x1Q5YJZttWVibJc0oqSAuXk5Dhdim1qynyqKvWywy5gMMILbMeqDzgtO8tFMAYMZumcl3379qmurk5z5syRy+XSG2+8ccX37N27V+Xl5XK73fryl7+s7du3W1kiHMCqDwDA1bA0vJw9e1ZLly7Vc889l9T43t5erV69WqtWrdLhw4f1xBNP6Hvf+5527txpZZmwWXTVx9gX6SPyeXJZ9QEASMjS20a1tbWqra1Nevz27ds1b948bd26VZK0aNEidXR06Ac/+IHuuecei6qE3S5dquuSYibuRl8/WXsDcxAAAAlNqTkvBw4cUHV1dcyxO++8Uy+88IJCoVDCSYXDw8MaHh4efR0MBiVJoVBIoVAoZbVFz5XKc2ayO66frR/dt1RP7/pAfcHP//0V5edqtfecvrqwgF5bjN9pe9Bn+9Bre1jV54mcb0qFl76+PhUVFcUcKyoq0oULF3Tq1Cn5fPHLGLds2aLGxsa447t379bMmTNTXmNbW1vKz5nJ/q5U+l3QpWBIys+RFuSfVZaLPtuJXtuDPtuHXtsj1X0+d+5c0mOnVHiRJJcr9lZBJBJJeDxq06ZNqq+vH30dDAZVXFys6upq5efnp6yuUCiktrY2VVVVZdSyUrvRZ/vQa3vQZ/vQa3tY1efonZNkTKnw4vV61dfXF3Osv79f06ZN06xZiZc15ubmKjc3N+54Tk6OJb+8Vp0Xseizfei1Peizfei1PVLd54mca0o9HqCysjLuMtTu3bu1fPlyfhEBAIAki8PLmTNn1NXVpa6uLkkXl0J3dXXp+PHjki7e8vnmN785On79+vX6+OOPVV9fr6NHj+rFF1/UCy+8oI0bN1pZJgAAMIilt406Ojp0++23j76Ozk1Zt26dduzYIb/fPxpkJKmkpES7du3SY489ph//+MeaM2eOfvjDH7JMGgAAjLI0vNx2222jE24T2bFjR9yxW2+9VZ2dnRZWBQAATDal5rwAAABcyZRabQRgbOGRCE9CBgARXgAjtHb71djSI3/g8ydt+zxuNdSVqqYsfvNGAEhn3DYCprjWbr82NHfGBBdJ6gsMaUNzp1q7/Q5VBgDOILwAU1h4JKLGlh4lmvYePdbY0qPwyNgT4wEg3RBegCmsvXcg7orLpSKS/IEhtfcO2FcUADiM8AJMYf2DYweXyYwDgHRAeAGmsMI8d0rHAUA6ILwAU1hFSYF8HrfGWhDt0sVVRxUlBXaWBQCOIrwAU1h2lksNdaWSFBdgoq8b6krZ7wVARiG8AFNcTZlPTWuWyeuJvTXk9bjVtGYZ+7wAyDhsUgcYoKbMp6pSLzvsAoAIL4AxsrNcqlwwy+kyAMBx3DYCAABGIbwAAACjEF4AAIBRCC8AAMAohBcAAGAUwgsAADAK4QUAABiF8AIAAIxCeAEAAEYhvAAAAKMQXgAAgFEILwAAwCiEFwAAYBTCCwAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjTHO6gEwUHomovXdA/YNDKsxzq6KkQNlZLqfLAgDACIQXm7V2+9XY0iN/YGj0mM/jVkNdqWrKfA5WBgCAGbhtZKPWbr82NHfGBBdJ6gsMaUNzp1q7/Q5VBgCAOQgvNgmPRNTY0qNIgu9FjzW29Cg8kmgEAACIIrzYpL13IO6Ky6UikvyBIbX3DthXFAAABiK82KR/cOzgMplxAABkKsKLTQrz3CkdBwBApiK82KSipEA+j1tjLYh26eKqo4qSAjvLAgDAOIQXm2RnudRQVypJcQEm+rqhrpT9XgAAuALCi41qynxqWrNMXk/srSGvx62mNcvY5wUAgCSwSZ3Nasp8qir1ssMuAACTRHhxQHaWS5ULZjldBgAARiK8AECKJHpuGYDUI7wAQAqM9dyyJ2uvd7AqID0xYRcArtJ4zy377qv/qf88zZw2IJUILwBwFZJ5btnPfp/Fc8uAFCK8AMBVSOa5ZZ+dd6nj4z/aVxSQ5ggvAHAVkn9u2bDFlQCZg/ACAFch+eeW5VpcCZA5CC9JCo9E9GHApZb3/Drwu9PcvwYgKbnnll0zPaLl8//EzrKAtMZS6SS0dvu1+c331RfMlnqOSLq4BLKhrpQt/YEMF31u2YbmTrmkmIm70UDzjetG2EUbSCGuvFxBdAlkXzD2fnVfYEgbmjvV2u13qDIAU8V4zy370X1LtXQWV2qBVOLKyziutATSJamxpUdVpV7+rwrIcGM9t2wkfEG7Pna6OiC9EF7GkcwSSH9gSO29AzyrCEDC55aNhB0qBkhj3DYaR/JLIJMbBwAArh7hZRzJL4FMbhwAALh6hJdxJLME0ufhybEAANiJ8DKO6BJISXEBJvq6oa6UyboAANiI8HIF0SWQRfmxu2N6PW41rVnGPi8AANjMlvCybds2lZSUyO12q7y8XO+8886YY/fs2SOXyxX39cEHH9hRakI1ZT7t+b+36JHSsJ793zfq/z24Ur/+u68SXAAAcIDlS6Vfe+01Pfroo9q2bZv+/M//XD/5yU9UW1urnp4ezZs3b8z3HTt2TPn5+aOvr732WqtLHVd2lkv/yxPR6iU+5eTkOFoLAACZzPIrL88++6weeOABfec739GiRYu0detWFRcXq6mpadz3FRYWyuv1jn5lZ2dbXSoAADCApeHl/PnzOnTokKqrq2OOV1dXa//+/eO+96abbpLP59Mdd9yhX/3qV1aWCQAADGLpbaNTp04pHA6rqKgo5nhRUZH6+voSvsfn8+n5559XeXm5hoeH9fLLL+uOO+7Qnj17dMstt8SNHx4e1vDw588dCgaDkqRQKKRQKJSynyV6rlSeE/Hos33otT3os33otT2s6vNEzueKRCKWPTHs008/1Ze+9CXt379flZWVo8efeeYZvfzyy0lPwq2rq5PL5dKbb74Z973NmzersbEx7vgrr7yimTNnTr54AABgm3Pnzun+++9XIBCImfOaiKVXXmbPnq3s7Oy4qyz9/f1xV2PGs3LlSjU3Nyf83qZNm1RfXz/6OhgMqri4WNXV1Vf84SciFAqpra1NVVVVTNi1EH22D722h4l9Do9E1PHxH9U/OKzCvFwtn/8nRuxnZWKvTWRVn6N3TpJhaXiZPn26ysvL1dbWpq9//eujx9va2nTXXXclfZ7Dhw/L50u8LDk3N1e5ublxx3Nyciz55bXqvIhFn+1Dr+1hSp9bu/1qbOmJeSitz+NWQ12pMdtDmNJr06W6zxM5l+VLpevr67V27VotX75clZWVev7553X8+HGtX79e0sUrJydOnNBLL70kSdq6dauuu+46LV68WOfPn1dzc7N27typnTt3Wl0qAGS01m6/NjR36vK5BH2BIW1o7mRjTkwZloeXe++9V6dPn9ZTTz0lv9+vsrIy7dq1S/Pnz5ck+f1+HT9+fHT8+fPntXHjRp04cUIzZszQ4sWL9dZbb2n16tVWlwoAGSs8ElFjS09ccJGkiC4+EqWxpUdVpV4jbiEhvVkeXiTpoYce0kMPPZTwezt27Ih5/fjjj+vxxx+3oSoAQFR770DMraLLRST5A0Nq7x1Q5YJZ9hUGJMCzjQAA6h8cO7hMZhxgJcILAECFee6UjgOsRHgBAKiipEA+j1tjzWZx6eKqo4qSAjvLAhIivAAAlJ3lUkNdqSTFBZjo64a6UibrYkogvAAAJEk1ZT41rVkmryf21pDX42aZNKYUW1YbAQDMUFPmU1WpV+29A+ofHFJh3sVbRVxxwVRCeAEAxMjOcrEcGlMat40AAIBRuPICAA4Lj0S4TQNMAOEFAByUDg9CBOzGbSMAcEj0QYiXb8sffRBia7ffocqAqY3wAgAOuNKDEKWLD0IMjyQaAWQ2wgsAOGAiD0IEEIvwAgAO4EGIwOQRXgDAATwIEZg8wgsAOIAHIQKTR3gBAAfwIERg8ggvAOAQHoQITA6b1AGAg3gQIjBxhBcAcBgPQgQmhttGAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjEF4AAIBRCC8AAMAohBcAAGAUwgsAADAK4QUAABiF8AIAAIxCeAEAAEYhvAAAAKMQXgAAgFEILwAAwCiEFwAAYBTCCwAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjEF4AAIBRCC8AAMAohBcAAGAUwgsAADAK4QUAABhlmtMFAAAwGeGRiNp7B9Q/OKTCPLdumpvndEmwCeEFAGCc1m6/Glt65A8MjR7z5udqtdel1Q7WBXtw2wgAYJTWbr82NHfGBBdJOhkc1ou/zdLb7590qDLYhfACADBGeCSixpYeRRJ8L3rsmV98oPBIohFIF4QXAIAx2nsH4q64xHLJHxhWe++AbTXBfsx5AQBMyuUTZitKCpSd5bL0M/sHxwsuEx8HMxFeAAATlmjCrM/jVkNdqWrKfJZ9bmGeO6XjYCZuGwEAJmSsCbN9gSFtaO5Ua7ffss+uKCmQz+PW2Nd3IvJ5clVRUmBZDXAe4QUAkLRkJsw2tvRYNmE2O8ulhrpSSYoLMNHXT9beYPntKziL8AIASNqVJsxGJPkDQ5ZOmK0p86lpzTJ5PbG3hryeXP2fhSO6c3GRZZ+NqYE5LwCApE2VCbM1ZT5VlXrjdth9u/UXln4upgbCCwAgaVNpwmx2lkuVC2aNvg6FQpZ/JqYGbhsBAJJ2pQmzLl1cdcSEWVjJlvCybds2lZSUyO12q7y8XO+888644/fu3avy8nK53W59+ctf1vbt2+0oEwBwBclMmG2oK2XCLCxleXh57bXX9Oijj+rJJ5/U4cOHtWrVKtXW1ur48eMJx/f29mr16tVatWqVDh8+rCeeeELf+973tHPnTqtLBQAkYewJs241rVlm6T4vgGTDnJdnn31WDzzwgL7zne9IkrZu3aq3335bTU1N2rJlS9z47du3a968edq6daskadGiRero6NAPfvAD3XPPPVaXCwBIQqIJs3bssAtIFoeX8+fP69ChQ/r7v//7mOPV1dXav39/wvccOHBA1dXVMcfuvPNOvfDCCwqFQsrJyYn53vDwsIaHh0dfB4NBSRcnbqVy8lb0XEwIsxZ9tg+9tke693n5vHxJ+ZKkkfAFjYSdqyXdez1VWNXniZzP0vBy6tQphcNhFRXFrrkvKipSX19fwvf09fUlHH/hwgWdOnVKPl/s5cgtW7aosbEx7jy7d+/WzJkzr/IniNfW1pbycyIefbYPvbYHfbYPvbZHqvt87ty5pMfaslTa5Yq9jBiJROKOXWl8ouOStGnTJtXX14++DgaDKi4uVnV1tfLz86+m7BihUEhtbW2qqqqKu/qD1KHP9qHX9qDP9qHX9rCqz9E7J8mwNLzMnj1b2dnZcVdZ+vv7466uRHm93oTjp02bplmzZsWNz83NVW5ubtzxnJwcS355rTovYtFn+9Bre9Bn+9Bre6S6zxM5l6WrjaZPn67y8vK4S0ttbW26+eabE76nsrIybvzu3bu1fPlyfhkBAID1S6Xr6+v1b//2b3rxxRd19OhRPfbYYzp+/LjWr18v6eJtn29+85uj49evX6+PP/5Y9fX1Onr0qF588UW98MIL2rhxo9WlAgAAA1g+5+Xee+/V6dOn9dRTT8nv96usrEy7du3S/PnzJUl+vz9mz5eSkhLt2rVLjz32mH784x9rzpw5+uEPf8gyaQAAIMmmCbsPPfSQHnrooYTf27FjR9yxW2+9VZ2dnRZXBQAATMSzjQAAgFEILwAAwCiEFwAAYBTCCwAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjEF4AAIBRCC8AAMAohBcAAGAUwgsAADAK4QUAABiF8AIAAIxCeAEAAEYhvAAAAKMQXgAAgFGmOV0AAABTWXgkovbeAfUPDqkwz62KkgJlZ7mcLiujEV4AABhDa7dfjS098geGRo/5PG411JWqpsznYGWZjdtGAAAk0Nrt14bmzpjgIkl9gSFtaO5Ua7ffocpAeAEA4DLhkYgaW3oUSfC96LHGlh6FRxKNgNUILwAAXKa9dyDuisulIpL8gSG19w7YVxRGEV4AALhM/+DYwWUy45BahBcAAC5TmOdO6TikFuEFAIDLVJQUyOdxa6wF0S5dXHVUUVJgZ1n4H4QXAAAuk53lUkNdqSTFBZjo64a6UvZ7cQjhBQCABGrKfGpas0xeT+ytIa/HraY1y9jnxUFsUgcAwBhqynyqKvWyw+4UQ3gBAGAc2VkuVS6Y5XQZuAS3jQAAgFEILwAAwCiEFwAAYBTCCwAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjEF4AAIBRCC8AAMAohBcAAGAUwgsAADAK4QUAABiF8AIAAIxCeAEAAEYhvAAAAKMQXgAAgFEILwAAwCjTnC4AAIBMEx6JqL13QP2DQyrMc6uipEDZWS6nyzIG4QUAABu1dvvV2NIjf2Bo9JjP41ZDXalqynwOVmYObhsBAGCT1m6/NjR3xgQXSeoLDGlDc6dau/0OVWYWwgsAADYIj0TU2NKjSILvRY81tvQoPJJoBC5FeAEAwAbtvQNxV1wuFZHkDwypvXfAvqIMRXgBAMAG/YNjB5fJjMtkloaXP/7xj1q7dq08Ho88Ho/Wrl2rzz77bNz3fOtb35LL5Yr5WrlypZVlAgBgucI8d0rHZTJLw8v999+vrq4utba2qrW1VV1dXVq7du0V31dTUyO/3z/6tWvXLivLBADAchUlBfJ53BprQbRLF1cdVZQU2FmWkSxbKn306FG1trbq4MGDWrFihSTpX//1X1VZWaljx47p+uuvH/O9ubm58nq9VpUGAIDtsrNcaqgr1YbmTrmkmIm70UDTUFfKfi9JsCy8HDhwQB6PZzS4SNLKlSvl8Xi0f//+ccPLnj17VFhYqGuuuUa33nqrnnnmGRUWFiYcOzw8rOHh4dHXwWBQkhQKhRQKhVL002j0XKk8J+LRZ/vQa3vQZ/uY0Os7rp+tH923VE/v+kB9wc//dnk9uXqy9gbdcf3sKV2/ZF2fJ3I+VyQSsWRN1j/90z9px44d+u1vfxtzfOHChfr2t7+tTZs2JXzfa6+9pi9+8YuaP3++ent79Y//+I+6cOGCDh06pNzc3LjxmzdvVmNjY9zxV155RTNnzkzNDwMAQAqNRKTfBV0KhqT8HGlBfkSZfsHl3Llzuv/++xUIBJSfnz/u2AlfeRkrLFzq3XfflSS5XPH/JiKRSMLjUffee+/oP5eVlWn58uWaP3++3nrrLX3jG9+IG79p0ybV19ePvg4GgyouLlZ1dfUVf/iJCIVCamtrU1VVlXJyclJ2XsSiz/ah1/agz/ah1/awqs/ROyfJmHB4eeSRR3TfffeNO+a6667Te++9p5MnT8Z97w9/+IOKioqS/jyfz6f58+frww8/TPj93NzchFdkcnJyLPnlteq8iEWf7UOv7UGf7UOv7ZHqPk/kXBMOL7Nnz9bs2bOvOK6yslKBQEDt7e2qqKiQJP3mN79RIBDQzTffnPTnnT59Wp988ol8Pp73AAAALFwqvWjRItXU1OjBBx/UwYMHdfDgQT344IP6q7/6q5jJujfccINef/11SdKZM2e0ceNGHThwQL///e+1Z88e1dXVafbs2fr6179uVakAAMAglu7z8h//8R+68cYbVV1drerqai1ZskQvv/xyzJhjx44pEAhIkrKzs3XkyBHdddddWrhwodatW6eFCxfqwIEDysvLs7JUAABgCMuWSktSQUGBmpubxx1z6WKnGTNm6O2337ayJAAAYDiebQQAAIxCeAEAAEYhvAAAAKMQXgAAgFEILwAAwCiEFwAAYBTCCwAAMArhBQAAGIXwAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMMs3pAgAAgDXCIxG19w6of3BIhXluVZQUKDvL5XRZV43wAgBAGmrt9quxpUf+wNDoMZ/HrYa6UtWU+SZ1zvBIRL/pHdChUy7N6h1Q5VcKHQlDhBcAANJMa7dfG5o7FbnseF9gSBuaO9W0ZtmEA0xsGMrWSx92XHUYmizmvAAAkEbCIxE1tvTEBRdJo8caW3oUHkk0IrFoGLr0Ko70eRhq7fZPvuBJILwAAJBG2nsH4kLGpSKS/IEhtfcOJHU+K8LQ1SK8AACQRvoHxw4ukxmX6jCUCoQXAADSSGGeO6XjUh2GUoHwAgBAGqkoKZDP49ZYa4BcurjqqKKkIKnzpToMpQLhBQCANJKd5VJDXakkxQWY6OuGutKklzinOgylAuEFAIA0U1PmU9OaZfJ6Yq+GeD3uCS+TTnUYSgX2eQEAIA3VlPlUVepNyQ670TB0+aZ3Xof2eSG8AACQprKzXKpcMCsl54qGoQMf9Wv3O79R9aoV7LALAACmtuwsl1aUFOj00YhWOPicJOa8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjpN0Ou5FIRJIUDAZTet5QKKRz584pGAwqJycnpefG5+izfei1Peizfei1Pazqc/TvdvTv+HjSLrwMDg5KkoqLix2uBAAATNTg4KA8Hs+4Y1yRZCKOQUZGRvTpp58qLy9PLlfqnrkQDAZVXFysTz75RPn5+Sk7L2LRZ/vQa3vQZ/vQa3tY1edIJKLBwUHNmTNHWVnjz2pJuysvWVlZmjt3rmXnz8/P5z8KG9Bn+9Bre9Bn+9Bre1jR5ytdcYliwi4AADAK4QUAABiF8JKk3NxcNTQ0KDc31+lS0hp9tg+9tgd9tg+9tsdU6HPaTdgFAADpjSsvAADAKIQXAABgFMILAAAwCuEFAAAYhfCSpG3btqmkpERut1vl5eV65513nC4p7ezbt091dXWaM2eOXC6X3njjDadLSjtbtmzRn/3ZnykvL0+FhYW6++67dezYMafLSktNTU1asmTJ6EZelZWV+sUvfuF0WWlvy5YtcrlcevTRR50uJa1s3rxZLpcr5svr9TpWD+ElCa+99poeffRRPfnkkzp8+LBWrVql2tpaHT9+3OnS0srZs2e1dOlSPffcc06Xkrb27t2rhx9+WAcPHlRbW5suXLig6upqnT171unS0s7cuXP1z//8z+ro6FBHR4e++tWv6q677tL777/vdGlp691339Xzzz+vJUuWOF1KWlq8eLH8fv/o15EjRxyrhaXSSVixYoWWLVumpqam0WOLFi3S3XffrS1btjhYWfpyuVx6/fXXdffddztdSlr7wx/+oMLCQu3du1e33HKL0+WkvYKCAn3/+9/XAw884HQpaefMmTNatmyZtm3bpqefflp/+qd/qq1btzpdVtrYvHmz3njjDXV1dTldiiSuvFzR+fPndejQIVVXV8ccr66u1v79+x2qCkiNQCAg6eIfVVgnHA7r1Vdf1dmzZ1VZWel0OWnp4Ycf1te+9jX95V/+pdOlpK0PP/xQc+bMUUlJie677z7913/9l2O1pN2DGVPt1KlTCofDKioqijleVFSkvr4+h6oCrl4kElF9fb3+4i/+QmVlZU6Xk5aOHDmiyspKDQ0N6Ytf/KJef/11lZaWOl1W2nn11VfV2dmpd9991+lS0taKFSv00ksvaeHChTp58qSefvpp3XzzzXr//fc1a9Ys2+shvCTJ5XLFvI5EInHHAJM88sgjeu+99/TrX//a6VLS1vXXX6+uri599tln2rlzp9atW6e9e/cSYFLok08+0d/+7d9q9+7dcrvdTpeTtmpra0f/+cYbb1RlZaUWLFign/70p6qvr7e9HsLLFcyePVvZ2dlxV1n6+/vjrsYApvjud7+rN998U/v27dPcuXOdLidtTZ8+XV/5ylckScuXL9e7776rf/mXf9FPfvIThytLH4cOHVJ/f7/Ky8tHj4XDYe3bt0/PPfechoeHlZ2d7WCF6ekLX/iCbrzxRn344YeOfD5zXq5g+vTpKi8vV1tbW8zxtrY23XzzzQ5VBUxOJBLRI488op/97Gf65S9/qZKSEqdLyiiRSETDw8NOl5FW7rjjDh05ckRdXV2jX8uXL9ff/M3fqKuri+BikeHhYR09elQ+n8+Rz+fKSxLq6+u1du1aLV++XJWVlXr++ed1/PhxrV+/3unS0sqZM2f00Ucfjb7u7e1VV1eXCgoKNG/ePAcrSx8PP/ywXnnlFf385z9XXl7e6BVFj8ejGTNmOFxdenniiSdUW1ur4uJiDQ4O6tVXX9WePXvU2trqdGlpJS8vL27O1he+8AXNmjWLuVwptHHjRtXV1WnevHnq7+/X008/rWAwqHXr1jlSD+ElCffee69Onz6tp556Sn6/X2VlZdq1a5fmz5/vdGlppaOjQ7fffvvo6+h91HXr1mnHjh0OVZVeosv9b7vttpjj//7v/65vfetb9heUxk6ePKm1a9fK7/fL4/FoyZIlam1tVVVVldOlARP23//93/rrv/5rnTp1Stdee61WrlypgwcPOvZ3kH1eAACAUZjzAgAAjEJ4AQAARiG8AAAAoxBeAACAUQgvAADAKIQXAABgFMILAAAwCuEFAAAYhfACAACMQngBAABGIbwAAACjEF4AAIBR/j/LddYQm2FvGAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "m = 20\n", "# Data : \n", "xi = np.linspace(0.0, 5.0, m) + 0.05*np.random.randn(m)\n", "yi = 0.8 + 0.5*xi-xi**2/6 + 0.2*np.random.randn(m)\n", "plt.plot(xi, yi, 'o'); plt.grid()" ] }, { "cell_type": "markdown", "id": "54b3f1d5-16ba-4ade-8fb9-787fc1146231", "metadata": {}, "source": [ "On souhaite représenter ce nuage de points par une ligne brisée en minimisant les écarts aux carrés entre les points et la ligne brisée.\n", "Dans un premier temps, coder la fonction dite chapeau ```y = phi(x)``` définie par\n", "\n", "$$\n", "\\phi(x) = \\max(0, 1-|x|).\n", "$$\n", "\n", "Représenter graphiquement la fonction $\\phi$ sur l'intervalle $[-3,3]$." ] }, { "cell_type": "code", "execution_count": null, "id": "1d356719-f874-4991-ab99-db793ff66ead", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "95453e70-1bf1-48a9-8b93-8d9ce34e40d8", "metadata": {}, "source": [ "On cherche ensuite une approximation $f(x)$ de la forme\n", "\n", "$$\n", "f(x) = \\sum_{j=0}^5 \\, u_j\\, \\phi(x - j),\n", "$$\n", "où les $(u_j)_{j=0,...,5}$ sont les inconnues. Réécrire le problème aux moindres carrés\n", "\n", "$$\n", "\\min_{(u_0,...,u_5)} \\quad \\frac{1}{2}\\sum_{i=1}^m\\, [f(x_i)-y_i]^2\n", "$$\n", "\n", "sous la forme\n", "\n", "$$\n", "\\min_{u\\in\\mathbb{R}^6} \\quad \\frac{1}{2} \\|Au - y\\|^2\n", "$$\n", "en déterminant la matrice $A$. Coder en python la matrice $A$, résoudre les équations normales et tracez sur le même graphique le nuage de points et la fonction d'approximation $f(x)$ ainsi trouvée. Afficher le vecteur $\\mathbf{u}$. Comment interprétez-vous les valeurs $u_i$ de $\\mathbf{u}$ ?\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2a4618a0-88c6-4cd5-be76-909561690c97", "metadata": { "tags": [] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "833c28d6-c8ec-4589-af4a-12c87100048e", "metadata": {}, "source": [ "#### Exercice 2 - Recherche d'une ellipse passant par un nuage de points\n", "\n", "On considère le nuage de points $(x_i,y_i)$, $i=1,...,m$ générés par le script python ci-dessous :" ] }, { "cell_type": "code", "execution_count": 3, "id": "1d84359f-80e0-4bcc-9168-7ff9880dee37", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAGwCAYAAACjPMHLAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOyxJREFUeJzt3X90VPWB9/HPNcSBVH5IIkkggcRqBaVVG6yGA8K0JQJbFoiyRc5SaNFTRKqYh9IHbWvoHmVrKY1WBd1S44+CbCW0npZqss8mgCtUsKRaWznaBkNDskjUBEGTYZjnj3GmN8wkmcncmTt35v06J2eYm3vJN9/MJJ/7/Wn4fD6fAAAAIEk6z+4CAAAAJBPCEQAAgAnhCAAAwIRwBAAAYEI4AgAAMCEcAQAAmBCOAAAATAbZXYBkd/bsWR07dkxDhw6VYRh2FwcAAETA5/Pp5MmTGj16tM47L7q2IMJRP44dO6bCwkK7iwEAAAbg6NGjKigoiOoawlE/hg4dKklqamrSyJEjbS6Ns3k8HtXW1qqsrEyZmZl2F8exqEfrUJfWoS6tQT1a57333lNxcXHw73g0CEf9CHSlDR06VMOGDbO5NM7m8XiUlZWlYcOG8aaPAfVoHerSOtSlNahH63g8Hkka0JAYBmQDAACYEI4AAABMCEcAAAAmhCMAAAATwhEAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgIljwtH69et1zTXXaOjQoRo1apTmzZunw4cP93vd7t27VVJSosGDB+viiy/W5s2bE1BaAADgVI4JR7t379btt9+u/fv3q66uTmfOnFFZWZlOnTrV6zVNTU2aPXu2pk6dqkOHDunuu+/WHXfcoR07diSw5AAAwEkG2V2ASL3wwgs9nj/xxBMaNWqUXn31VV1//fVhr9m8ebPGjh2rqqoqSdKECRN08OBBbdiwQTfeeGPYa7q6utTV1RV83tnZKUnyeDzyeDwWfCfpK1B/1GNsqEfrUJfWoS6tQT1aJ5Y6dEw4OldHR4ckaeTIkb2es2/fPpWVlfU4dsMNN2jLli3yeDzKzMwMuWb9+vVat25dyPH6+nplZWXFWGpIUl1dnd1FSAnUo3WoS+tQl9agHmN3+vTpAV/ryHDk8/lUUVGhKVOmaOLEib2e19bWptzc3B7HcnNzdebMGZ04cUL5+fkh16xdu1YVFRXB552dnSosLJTb7VZ2drZ130Qa8ng8qqur04wZM8IGU0SGerQOdWkd6tIa1KN12tvbB3ytI8PRypUr9dprr+mll17q91zDMHo89/l8YY8HuFwuuVyukOOZmZm8UC1CXVqDerQOdWkd6tIa1GPsYqk/x4Wjb33rW3r++ee1Z88eFRQU9HluXl6e2traehw7fvy4Bg0aRCsQAAAIyzGz1Xw+n1auXKmamhr993//t4qLi/u9prS0NKTftra2VpMmTSKRAwCAsBwTjm6//XY988wz2rp1q4YOHaq2tja1tbXpo48+Cp6zdu1afe1rXws+X758ud555x1VVFToL3/5i37+859ry5YtWr16tR3fAgAAcADHhKNNmzapo6ND06dPV35+fvBj+/btwXNaW1vV3NwcfF5cXKxdu3apoaFBV111lf7t3/5NDz30UK/T+AEAABwz5igwkLov1dXVIcemTZumP/zhD3EoEQAASEWOaTkCAABIBMIRAACACeEIAADAhHAEAABgQjgCAAAwIRwBAACYEI4AAABMCEcAAAAmhCMAAAATwhEAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgAnhCAAAwIRwBAAAYEI4AgAAMCEcAQAAmBCOAAAATAhHAAAAJoQjAAAAE8IRAACACeEIAADAhHAEAABgQjgCAAAwIRwBAACYOCoc7dmzR3PmzNHo0aNlGIZ+9atf9Xl+Q0ODDMMI+XjzzTcTU2AAAOA4g+wuQDROnTqlK6+8Ul//+td14403Rnzd4cOHNWzYsODziy66KB7FAwAAKcBR4WjWrFmaNWtW1NeNGjVKI0aMsL5AAAAg5TgqHA3U1VdfrY8//liXX365vvvd78rtdvd6bldXl7q6uoLPOzs7JUkej0cejyfuZU1lgfqjHmNDPVqHurQOdWkN6tE6sdSh4fP5fBaWJWEMw9DOnTs1b968Xs85fPiw9uzZo5KSEnV1denpp5/W5s2b1dDQoOuvvz7sNZWVlVq3bl3I8a1btyorK8uq4gMAgDg6ffq0Fi1apI6Ojh5DayKR0uEonDlz5sgwDD3//PNhPx+u5aiwsFCtra3Kzs6Opchpz+PxqK6uTjNmzFBmZqbdxXEs6tE61KV1qEtrUI/WaW9vV35+/oDCUVp0q5ldd911euaZZ3r9vMvlksvlCjmemZnJC9Ui1KU1qEfrUJfWoS6tQT3GLpb6c9RUfiscOnRI+fn5dhcDAAAkKUe1HH344Yd6++23g8+bmprU2NiokSNHauzYsVq7dq1aWlr01FNPSZKqqqpUVFSkK664Qt3d3XrmmWe0Y8cO7dixw65vAQAAJDlHhaODBw/2mGlWUVEhSVqyZImqq6vV2tqq5ubm4Oe7u7u1evVqtbS0aMiQIbriiiv029/+VrNnz0542QEAgDM4KhxNnz5dfY0fr66u7vF8zZo1WrNmTZxLBQAAUknajTkCAADoC+EIAADAhHAEAABgQjgCAAAwIRwBAACYEI4AAABMCEcAAAAmhCMAAAATRy0CCQCwgdcr7d0rtbZK+fnS1KlSRobdpQLihnAEAOhdTY10553S3//+j2MFBdKDD0rl5faVC4gjutUAOI/XKzU0SNu2+R+9XrtLlJpqaqSbbuoZjCSppcV/vKbGnnIBcUY4AuAsNTVSUZHkdkuLFvkfi4pS4w91MoU+r9ffYhRuP8vAsVWrCKZISYQjAM6Ryi0ZyRb69u4NrWczn086etR/HpBiCEcAnCGVWzKSMfS1tlp7HuAghCMAzpCqLRnJGvry8609D3AQwhEAZ0jVloxkDX1Tp/pnpRlG+M8bhlRY6D8PSDGEIwDOkKotGcka+jIy/NP1pdCAFHheVcV6R0hJhCMAzpCqLRnJHPrKy6XnnpPGjOl5vKDAf5x1jpCiWAQSgDMEWjJuuskfhMxjdJzckhEIfS0t4ccdGYb/83aFvvJyae5cVshGWqHlCIBzpGJLhhO6rzIypOnTpZtv9j8SjJDiaDkC4CyxtmSY9gkzLrooOab+B0JfuG06qqqcGfoAByMcAXCeQEtGtM7ZJ2yQpLLsbBmPPir9y79YWsSo0X0FJA3CEYD4S4Zd3QMLLZ4zrmdwe7u0cKE0aJD9LTQDDX0ALMWYIwDxlQzbYvSx0GJwlI9TV9cGYDnCEYD4SZZtMfpZaNFw6uraAOKCcATAel6v9P/+n3TrrcmxLUayLrQIICkRjgBYK9CN9uUvS++91/t5iWytSeaFFgEkHQZkA7BOL4Oe+5SI1pp+Flr0GYYMOxdaBJBUaDkCYI2+dpfvSyJaa/pYaDFYWrsXWgSQNAhHgNN5vVJDg7Rtm//RrhlX/e0uf65E74XWy+raH+XkyPvss/ZP4weQNBwVjvbs2aM5c+Zo9OjRMgxDv/rVr/q9Zvfu3SopKdHgwYN18cUXa/PmzfEvKJAoyTBNPiCa7jG7tsUoL5eOHJHq66WtW3Wmrk51jz0m3/z5iSsDgKTnqHB06tQpXXnllXr44YcjOr+pqUmzZ8/W1KlTdejQId1999264447tGPHjjiXFIg/Y+fO5JgmHxBN95ide6GZ9gnzTZtGVxqAEI4akD1r1izNmjUr4vM3b96ssWPHqqqqSpI0YcIEHTx4UBs2bNCNN94Yp1ICCeD1KqOiovdp8obhnyY/d27i/vj3t7u8JGVnS9u3s3kpgKTmqHAUrX379qmsrKzHsRtuuEFbtmyRx+NRZmZmyDVdXV3q6uoKPu/s7JQkeTweeTye+BY4xQXqj3qMjcfjUfaf/yyjpaX3kz6ZJn+mvt7fOpIgxo9/rIyFCyXD8C+sGCjOJ91o3kcfle/666WzZ/0fNuM1aR3q0hrUo3ViqcOUDkdtbW3Kzc3tcSw3N1dnzpzRiRMnlB+mG2D9+vVat25dyPH6+nplZWXFrazppK6uzu4iON6Y99+P6LzG3/1OLadOxbk0Ji6X8tes0Wd/9jMNaW8PHv4oO1t/WrZMrS6XtGtX4soTIV6T1qEurUE9xu706dMDvjalw5EkGedO2/3kbvbc4wFr165VRUVF8HlnZ6cKCwvldruVnZ0dv4KmAY/Ho7q6Os2YMSNsqx0i4/F49IfXX4/o3KtmzdKVCWw5kiTNni1VVurMSy8FN5rNnDJFV2dk6OrElqRfvCatQ11ag3q0TrvpBi1aKR2O8vLy1NbW1uPY8ePHNWjQoF6DjsvlksvlCjmemZnJC9Ui1GXs2i+/XL4xY2QcOxZ+fI9hSAUFGuR22zO2JzPTv0K2Q/CatA51aQ3qMXax1J+jZqtFq7S0NKRpsra2VpMmTeJFB2fLyJB340b/v89tBbVrmjzQn2RZkwvoh6PC0YcffqjGxkY1NjZK8k/Vb2xsVHNzsyR/l9jXvva14PnLly/XO++8o4qKCv3lL3/Rz3/+c23ZskWrV6+2o/iApXzz54dd1NDWafJAb5JpTS6gH47qVjt48KDcbnfweWBs0JIlS1RdXa3W1tZgUJKk4uJi7dq1S3fddZceeeQRjR49Wg899BDT+JE6ysv90/X37g2O79HUqbQYIbn0tudeYE0uwjySjKPC0fTp04MDqsOprq4OOTZt2jT94Q9/iGOpkFK8XucFjcCihkAy6mvPPbvW5AL64ahuNSCuaPYH+jaQMUP97bn3yZpc2rvXqlICMSMcAdI/mv2TZSsOINkM9OYh0j33otmbz04MKk8LhCOgv2Z/yd/szy9BpKtYbh4i3XMvmr357ELrctogHAE0+wO9i/XmIbDnXi8L78owpMJC/3nJLNqASAuToxGOgFRr9gesFOHNg/HSS+E/n5EhPfig/99OXZMr2oBIC5PjEY6AVGr2B6xmxc1Debmz1+SKpnWZ8YspgXAEpEqzPxAPVt08lJdLR45I9fXS1q3+x6am5A9GUuQBsaWF8YspgnAEpEKzPxAvEd48+KZM6f//CqzJdfPN/kenvKciDYjvvsv4xRRBOAIk5zf7A/HCzUPkrcsXXRTZ/8f4xaRHOAICnNzsD8RTut88RBoQz62f3jB+Mek5avsQIO7YigMIL9338QsExDvv7Nl1VlDgD0bl5f6xRAUF/rFH4cYdGYb/84xfTHqEIwBAZNL95qG/gBhoYbrpJn8QMgekdOmCTBGEIwAAItVfQIykhQlJj3AEAICV0r0LMgUQjpDevF5+gQGwXrp3QToc4Qjpq6YmfNP3gw/S9A0AaYyp/EhPLPEPAOgF4QjpJ9ZdxgEAKY1whP55vVJDg7Rtm//R6aEhmk0kAQBphzFH6FsqjsuxYpdxAEDKouUIvUvVcTlW7TIOAE6Taj0BcUI4QnipPC4n0k0kWeIfQCqpqZGKiiS3W1q0yP9YVOTcG904IhwhvFQel8Mu4wDSTar2BMQJ4QjhOWlczkCaidN9l3EA6SOVewLihAHZCM8p43JiGTDOEv8A0kE0PQGs6i2JcITeBMbltLSEv9swDP/n7RyXE2gmPrd8gWbiSFqAWOIfQKpzUk9AkqBbDeEl+7gcmokBIDJO6QlIIoQj9C6Zx+Wk8oBxALASM3SjRrca+pas43JoJgaAyAR6Am66yR+EzC3uydATkIQIR+hfMo7LoZkYACIX6AkIN4GlqooZuucgHMGZnDBgHACSSbL2BCQhx405evTRR1VcXKzBgwerpKREe/sYU9LQ0CDDMEI+3nzzzQSWGHGR7APGASAZBXoCbr7Z/8jvyLAcFY62b9+uVatW6Z577tGhQ4c0depUzZo1S83NzX1ed/jwYbW2tgY/Lr300gSVGHGVzAPGAQCO5ahutY0bN2rZsmW65ZZbJElVVVV68cUXtWnTJq1fv77X60aNGqURI0ZE9DW6urrU1dUVfN7Z2SlJ8ng88ng8Ay88gvVnaT3OmSPNni3jpZeCzcS+KVP8d0Mp+vOKSz2mKerSOtSlNahH68RSh4bPF27ARvLp7u5WVlaWfvnLX2r+/PnB43feeacaGxu1e/fukGsaGhrkdrtVVFSkjz/+WJdffrm++93vyu129/p1KisrtW7dupDjW7duVVZWljXfDAAAiKvTp09r0aJF6ujo0LBhw6K61jEtRydOnJDX61Vubm6P47m5uWprawt7TX5+vh5//HGVlJSoq6tLTz/9tL70pS+poaFB119/fdhr1q5dq4qKiuDzzs5OFRYWyu12Kzs727pvKA15PB7V1dVpxowZyszMtLs4jkU9Woe6tA51aQ3q0Trt7e0DvtYx4SjAOGfwrc/nCzkWcNlll+myyy4LPi8tLdXRo0e1YcOGXsORy+WSy+UKOZ6ZmckL1eu1ZJYDdWkN6tE61KV1qEtrUI+xi6X+HDMgOycnRxkZGSGtRMePHw9pTerLddddp7feesvq4qW+mhqpqEhyu6VFi/yPRUX+4wAA5/B6pYYGads2/yPbLIVwTDg6//zzVVJSorq6uh7H6+rqNHny5Ij/n0OHDimfhQGjE9jg9dztOgIbvBKQAMAZuNGNiKO61SoqKrR48WJNmjRJpaWlevzxx9Xc3Kzly5dL8o8Xamlp0VNPPSXJP5utqKhIV1xxhbq7u/XMM89ox44d2rFjh53fhrP0t8GrYfg3eJ07l/UyACCZBW50z/19HrjRZQmUIEeFo69+9atqb2/XD37wA7W2tmrixInatWuXxo0bJ0lqbW3tseZRd3e3Vq9erZaWFg0ZMkRXXHGFfvvb32r27Nl2fQvOE80Gr8m2xQgAwI8b3ag4KhxJ0ooVK7RixYqwn6uuru7xfM2aNVqzZk0CSpXC2OAVAJyPG92oOGbMEWzCBq8A4Hzc6EbFcS1HSDA2eAUA+8W6lAo3ulGh5Qh9Y4NXALCXFTPMAje6vawLKMOQCgu50f0E4Qj9Y4NXALCHVUupcKMbFbrVEJm5c6Xhw/0Lhkn+AXvTp/NGAoB4sXqGWeBG9847e4atggJ/MOJGN4hwhP7V1IS+maqr/XchvJkAID7iMcOsvNwfpsKNX7Joi6hUQDhC31g0DADsEa8ZZhkZoWEq3E1wQUHa3gQz5ijenLyHTX9NupK/SddJ3xMAOEWiZpixRVQIwlE82bmHjRWhLJomXQCAtRIxw4yb4LAIR/FiZxK3KpSxaBgA2CcRM8y4CQ6LcBQPdiZxK0MZi4YBgL3ivZQKN8FhEY7iwa4kbnUoY9EwALBfebl05IhUXy9t3ep/bGqyZqA0N8FhEY7iwa4kbnUoY9EwAEgOgRlmN99s7Rpz3ASHRTiKB7uSeDxCGatjA0Dq4iY4LMJRPNiVxOMVyuLZpAsAsBc3wSFYBDIeAkn8ppv8Qcg8BiieSTwQylpawo87Mgz/5wcSysItGgYASA19rZydhmg5ihc7kjjNowCAgYrXuCYHIhzFkx3dUTSPAgAQE7rV4s2O7iiaRwEAGDDCUapijBAAAANCOAIAAPbwepOyl4NwBAAAEq+mxr+rg3nx4oIC/8Qim8fHMiAbAIBU4PVKDQ3Stm3+x3js32kVOzdnjwDhCAAAp6upkYqKJLdbWrTI/1hUZHvICMvOzdkjRDgCAMDJkrwVJoRdm7NHgXCU7pzUDAsA6MkBrTAh7NqcPQqEo3TmpGZYAEAoB7TChLBrc/YoEI7SldOaYQEAoRzQChPCrs3Zo0A4SkdObIYFAISKtHXlrbfiW45oOGAfUMJROnJiMywAIFR/rTABlZXJ1SOQ5PuAOi4cPfrooyouLtbgwYNVUlKivf38Ad+9e7dKSko0ePBgXXzxxdq8eXOCSprEnNgMCwAIFWiFCdcTcK5k6xGwY3P2CDkqHG3fvl2rVq3SPffco0OHDmnq1KmaNWuWmpubw57f1NSk2bNna+rUqTp06JDuvvtu3XHHHdqxY0eCS55kHDAYDgAQofJyad26vs9J1h6BwD6gN9/sf0yCrUOkAYSjpUuXas+ePfEoS782btyoZcuW6ZZbbtGECRNUVVWlwsJCbdq0Kez5mzdv1tixY1VVVaUJEybolltu0Te+8Q1t2LAhwSVPMg4YDAcAiMKll0Z2Hj0CEYl6b7WTJ0+qrKxMhYWF+vrXv64lS5ZozLl9hnHQ3d2tV199Vf/3//7fHsfLysr08ssvh71m3759Kisr63Hshhtu0JYtW+TxeJSZmRlyTVdXl7q6uoLPOzs7JUkej0cejyfWbyNpGD/+sTIWLpQMQ4apOdb3SWDybtgg39mz0tmzln3NQP2lUj3agXq0DnVpHerSGgOtR+OiiyL6g37moovkS5OfUSyvxajD0Y4dO9Te3q5nnnlG1dXVuvfee/XlL39Zy5Yt09y5c8MGDiucOHFCXq9Xubm5PY7n5uaqra0t7DVtbW1hzz9z5oxOnDih/DDdRuvXr9e6MM2T9fX1ysrKiuE7SDIul/LXrNFnf/YzDWlvDx7+KDtbf1q2TK0ul7RrV1y+dF1dXVz+33RDPVqHurQOdWmNqOvR61VZdrYGt7crXJ+AT9JHOTmq6+yM2+/2c8uT/ec/a/D77+vjCy9U++WXJ7zL7PTp0wO+NupwJEnZ2dm68847deedd+rQoUP6+c9/rsWLF+uCCy7Qv/7rv2rFihW6NNImvigZ53QF+Xy+kGP9nR/ueMDatWtVUVERfN7Z2anCwkK53W5lZ2cPtNjJafZsqbJSZ156yd/Ump+vzClTdHVGhq6Ow5fzeDyqq6vTjBkz4hai0wH1aB3q0jrUpTViqUfj0UelhQvlk8L2CJz/yCOaPWeOlcUNX46dO5VRUSGjpeUfZRgzRt6NG+WbPz/uXz+g3XTjH60BhaOA1tZW1dbWqra2VhkZGZo9e7beeOMNXX755XrggQd01113xfLf95CTk6OMjIyQVqLjx4+HtA4F5OXlhT1/0KBBvQYdl8sll8sVcjwzMzM13/CZmdKXv5zgL5midZlg1KN1qEvrUJfWGFA9/su/SIMG+dexMy3XYhQUSFVVGpSIWWA1NdLChSGz54xjxzRo4cKETtOP5XUY9YBsj8ejHTt26Ctf+YrGjRunX/7yl7rrrrvU2tqqJ598UrW1tXr66af1gx/8YMCFCuf8889XSUlJSFNjXV2dJk+eHPaa0tLSkPNra2s1adIk3rwAgNRj5/T4FFpgOOqWo/z8fJ09e1Y333yzXnnlFV111VUh59xwww0aMWKEBcXrqaKiQosXL9akSZNUWlqqxx9/XM3NzVq+fLkkf5dYS0uLnnrqKUnS8uXL9fDDD6uiokK33nqr9u3bpy1btmjbtm2Wlw0AgKQQmB6faNEsMGxH+aIQdTj6yU9+ogULFmjw4MG9nnPhhReqqakppoKF89WvflXt7e36wQ9+oNbWVk2cOFG7du3SuHHjJPm7+cxrHhUXF2vXrl2666679Mgjj2j06NF66KGHdOONN1peNgAA0loKLTAcdThavHhxPMoRsRUrVmjFihVhP1ddXR1ybNq0afrDH/4Q51IBAJDmUmiB4ZgGZMNGXq+/afKTWWaaOjVpVhYFAKShwALDLS3hxx0Zhv/zDlhg2FHbh+ATNTVSUZHkdkuLFvkfi4qSa1NBAEB6CezzJoXuwBB4XlXliBt5wpHT1NRIN90UOuitpcV/nIAEALBLebl/uv65O2cUFCR0Gn+s6FZzkv6mSRqGf5rk3LmOSOYAgBRUXu7/O+TgoR+EIydJoWmSAIAUZtdyAhahW81JUmiaJAAAyYpw5CQpNE0SAIBkRThyksA0yd422jUMqbDQEdMkAQBIVoQjJ0mhaZIAACQrwpHTpMg0SQAAkhWz1ZwoBaZJAgCQrAhHTuXwaZIAACQrutUAAABMCEcAAAAmdKsBAIDE8nqTetws4QgAACROTY1/n1DzdlgFBf6lapJkxjXdagAAIDFqaqSbbgrdJ7SlxX+8psaecp2DcAQAAOLP6/W3GPl8oZ8LHFu1yn+ezQhHAAAg/vbuDW0xMvP5pKNH/efZjHAEAADir7XV2vPiiHAEAADiLz/f2vPiiHAEAADib+pU/6y0czdODzAMqbDQf57NCEcAACD+MjL80/Wl0IAUeF5VlRTrHRGOAABAYpSXS889J40Z0/N4QYH/eJKsc8QikAAAIHHKy6W5c1khGwAAICgjQ5o+3e5S9IpuNQAAABPCEQAAgAndarBGku+wDABApAhHiJ0DdlgGgIhxs5f2HNOt9v7772vx4sUaPny4hg8frsWLF+uDDz7o85qlS5fKMIweH9ddd11iChwPXq/U0CBt2+Z/jGZzvliu7YtDdlgGgIjU1EhFRZLbLS1a5H8sKuJ3WZpxTDhatGiRGhsb9cILL+iFF15QY2OjFi9e3O91M2fOVGtra/Bj165dCShtHMTyho3Xm91BOywDQL+42cMnHBGO/vKXv+iFF17Qz372M5WWlqq0tFT/8R//od/85jc6fPhwn9e6XC7l5eUFP0aOHJmgUlsoljdsPN/sDtphGQD6xM0eTBwx5mjfvn0aPny4rr322uCx6667TsOHD9fLL7+syy67rNdrGxoaNGrUKI0YMULTpk3Tfffdp1GjRvV6fldXl7q6uoLPOzs7JUkej0cej8eC7yZKXq8G3XGH5PMpZDcan08+w5DuvFNnZs8O7ROP5doIGEePRvQCOnP0qHym+rOlHlMI9Wgd6tI6Tq9LY/duDYrgZu9Mfb1806bFrRxOr8dkEksdOiIctbW1hQ00o0aNUltbW6/XzZo1SwsWLNC4cePU1NSk733ve/riF7+oV199VS6XK+w169ev17p160KO19fXKysra+DfxABlv/66prS09Pp5w+eT/v53/X7DBrV/9rOWXRtR2d55R1MiOG//O++o3dSdWVdXF/XXQijq0TrUpXWcWpdj9uzRpAjOa/zd79Ry6lTcy+PUekwmp0+fHvC1toajysrKsEHE7MCBA5IkI8wuvj6fL+zxgK9+9avBf0+cOFGTJk3SuHHj9Nvf/lblvcyiWrt2rSoqKoLPOzs7VVhYKLfbrezs7D7LGg/GJy1X/blu3Dj5Zs+27NqI3HCDfJs3S8eO+YPWOXyGIY0Zo2tXr5YyMuTxeFRXV6cZM2YoMzMz+q8HSaIeLURdWsfpdWl86lPSxo39nnfVrFm6Ms4tR06ux2TS3t4+4GttDUcrV67UwoUL+zynqKhIr732mv73f/835HPvvvuucnNzI/56+fn5GjdunN56661ez3G5XGFblTIzM+15oRYWRnTaoMJC6dzyxXJtJDIzpYce8o9dMoyeffWG4e/Ke/BBZQ4efM5lNtVliqEerUNdWsexdel2+5cgaWkJP+7IMKSCAg1yuxMyrd+x9ZhEYqk/W8NRTk6OcnJy+j2vtLRUHR0deuWVV/SFL3xBkvT73/9eHR0dmjx5csRfr729XUePHlV+fv6Ay5xwU6dG9IbV1KnWXhupwA7L4dY5qqpinSMAzpCR4V+brZebPUn+32msd5QWHDFbbcKECZo5c6ZuvfVW7d+/X/v379ett96qr3zlKz0GY48fP147d+6UJH344YdavXq19u3bpyNHjqihoUFz5sxRTk6O5s+fb9e3Er3AG1b6xxs0oL83bCzXRqO8XDpyRKqvl7Zu9T82NRGMADhL4GZvzJiexwsK/Mf5nZY2HBGOJOkXv/iFPvvZz6qsrExlZWX63Oc+p6effrrHOYcPH1ZHR4ckKSMjQ6+//rrmzp2rz3zmM1qyZIk+85nPaN++fRo6dKgd38LAxfKGTdSbPbDD8s03+x+5uwLgRNzsQQ6ZrSZJI0eO1DPPPNPnOT5TM+iQIUP04osvxrtYiVNeLs2dO7Al7WO5FgDSTeBmD2nLMeEIiu0Ny5sdAICIOKZbDQAAIBFoOQIAIBG8XoY3OAThCACAeKupCb/kyYMPMtg7CdGtBgBAPMVzA3DEBeEIAIB48Xr9LUbhFuINHFu1yn8ekgbhCACAeNm7N7TFyMznk44e9Z+HpEE4AgAgXlpbrT0PCUE4AgAgXiLdy9NJe36mAWarAQAQL4ENwHvrWuttA3CvV/qf/2Hav00IRwAAxMuvfy199FH4z/WyAXj+vn0adPvt/tlsAUz7Tyi61QCvV2pokLZt8z8yawSAFQJT+Nvbw39+5MiQDcCNnTt1zQ9/2DMYSUz7TzDCEdJbTY1UVCS53dKiRf7HoiJ+AQGITV9T+AOGDPFvCm66JqOiQpJknHsu0/4TinCE9MXCbADipb8p/JL/8+Yp/Hv3ymhpCQ1GAUz7TxjCEdITC7MBiKeBTOFn2n/SIBwhPbEwG4B4GsgUfqb9Jw3CEdITd2gA4ikwhd/opZPMMKTCwp5T+KdOlW/MGPU6SincNYgLwhHSE3doAOIpI8M/9V4KDUi9TOFXRoa8GzdKknyRXoO4IBwhPQ3krg4AolFe7p+qP2ZMz+MFBSFT+AN88+frwHe+I40eHfE1sB6LQCI9eL3+8UPm1WYffNA/K80weg7M5g4NiJ9w78VUfp+Vl/un60fxPbeWlupMZaUy9+9Pn3pKMoQjpL6aGv/MNPMA7MBqs889F/5zVVXcoQFW6+u9mMrvt4wMafr0+F8DyxCOkNoCaxmdO2U/sJbRc89JR46k150skCjmVqK33pIqK/t+L6ZyQIKjEI6Quvpby8gw/GsZzZ3LHRpgtXCtROGc+17kxgRJgAHZSF2sZYRU45R9AHtbfb43vBeRZGg5QupiLSOkEqeM14lkT7He8F5EkqDlCKmLtYyQKhy0D6Dx0kuRtxidi/cikgThCKmLtYyQCpy2D+BAWn94LyLJEI6QugayQi2QbJw2di7a1h/ei0hChCOktgGsUAskFYeNnfNNmdJ3i+25eC8iCTEgG6lvACvUOk66rTqcrOLxc3Da2LlAi21vq8/7fNK6ddKll/JaRdIiHCE9pPJqs06ZxZTq4vVzCIyda2kJP+7IMPyfT6bxOoEWW1afh0M5plvtvvvu0+TJk5WVlaURI0ZEdI3P51NlZaVGjx6tIUOGaPr06XrjjTfiW1AgkRw0iymlxfPn4NSxc+Xl/tXn6+ulrVv9j01NBCM4gmPCUXd3txYsWKDbbrst4mseeOABbdy4UQ8//LAOHDigvLw8zZgxQydPnoxjSYEEcdosplSViJ+DU8fOBVpsb77Z/5hsAQ7ohWO61datWydJqq6ujuh8n8+nqqoq3XPPPSr/5BfHk08+qdzcXG3dulXf/OY3w17X1dWlrq6u4PPOzk5JksfjkcfjieE7QKD+qMfYBOrP29CgzAhmMZ2pr5dv2rQElc5ZrHhNGrt3a1AsPwev17820CfjlHxTpoQPEXPmSLNnhz83Cd5TvL+tQT1aJ5Y6dEw4ilZTU5Pa2tpUVlYWPOZyuTRt2jS9/PLLvYaj9evXB4OYWX19vbKysuJW3nRSV1dndxFi5/Uq+89/1uD339fHF16o9ssvT/hd8Z/q6jQpgvMaf/c7tZw6FffyOFksr8kxe/YM+OeQv2+fPvuzn2lIe3vw2EfZ2Xr9llvUWlra+382bJh06pT04osDLHX8pMT7OwlQj7E7ffr0gK9N2XDU1tYmScrNze1xPDc3V++8806v161du1YVFRXB552dnSosLJTb7VZ2dnZ8CpsmPB6P6urqNGPGDGVmZtpdnAEzdu5URkWFjJaW4DHfmDHybtwo3/z5cf/6gXqcOGOGtHFjv+dfNWuWrqTlKCwrXpPGpz41oJ+DsXOnMh54IKQ7bvB77+maBx6Q99lnE/J6skqqvL/tRj1ap9100xEtW8NRZWVl2FYaswMHDmjSpEjuy8IzzhnA6PP5Qo6ZuVwuuVyukOOZmZm8UC3i6LqsqZEWLgz5g2YcO6ZBCxcmdPxHxvTpEc1iGuR2M9ajHzG9Jt3u6H8OXq/0f/5P2PONT3apH7R6tXTjjY772Tn6/Z1EqMfYxVJ/toajlStXauHChX2eU1RUNKD/Oy8vT5K/BSnftP7H8ePHQ1qTgIj0N/DWMPwDb+fOTcwftP7Wk5GScxZTqhnIzyGaVa9TdQkKIInZOlstJydH48eP7/Nj8ODBA/q/i4uLlZeX16Pftru7W7t379bkyZOt+haQTpJxG4dYZzF5vVJDg7Rtm/+RmW0DE+3PwWGrXgPpxjFjjpqbm/Xee++publZXq9XjY2NkqRLLrlEF1xwgSRp/PjxWr9+vebPny/DMLRq1Srdf//9uvTSS3XppZfq/vvvV1ZWlhYtWmTjdwLHStY/aANdAZzFI60Vzc/BaateA2nGMeHo+9//vp588sng86uvvlqSfxbZ9E+anQ8fPqyOjo7gOWvWrNFHH32kFStW6P3339e1116r2tpaDR06NKFlR4pItj9oXq/0P/8zsK0qAosWnttFGFi0MJnXzklmka7E7sRVr4E04phFIKurq+Xz+UI+ppt+Efl8Pi1dujT43DAMVVZWqrW1VR9//LF2796tiRMnJr7wSA2BP2i9Deg3DKmwMCF/0PL37dOgSy7xDwZetMj/WFQU2UrMLB5pjVi6JJ266jWQJhwTjgDbJckfNGPnTl3zwx/6Wx3MIt2qIhnHTjlNTY0/jA4knAY4ddVrIA0QjoBo2P0HzetVxifrcIW0X0Xa6pOsY6ecwsp91Nh/DEhKjhlzBCSNgQ6AtsLevT0WnwwRyRTwZBs75STxWM4h0nFKABKGcAQMhF1/0Kxo9WEw8MCxPhGQFuhWA5zEilafJBk75Uh0SQJpgXAEOMnUqfKNGaMw7T1+kc6Ys3vslFPRJQmkBcIR4CQZGfJ+ssmpL9ZWHwYDRy+JlnMAED+EI8BhfPPn68B3viONHt3zEwNp9QmMnbr5Zv8jXWl9o0sSSAuEI8CBWktLdebtt2n1sQNdkkDKY7Ya4FRMAY+O1xu6/MJA2bmcA4C4IxwBSH29bLJr/PjHkss1sP+TcAqkLLrVAKS2Pla0zli4UPn79tlTLgBJi3AEIHVFsMnuxC1b2GQXQA+EIwCpq58VrQ2fT1knTsh46aUEFgpAsiMcAUhdrGgNYAAIRwBSFytaAxgAwhGA1NXPitY+w9DpnBz5pkxJcMEAJDPCEYDUFcGK1n9atoz1iQD0QDgCkNr6WNHa++yzai0ttadcAJIWi0ACSH29rGjtO3tW2rXL7tIBSDKEIwDpIdyK1mfP2lIUAMmNbjUAAAATwhEAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPHhKP77rtPkydPVlZWlkaMGBHRNUuXLpVhGD0+rrvuuvgWFAAAOJpjwlF3d7cWLFig2267LarrZs6cqdbW1uDHLvZRAgAAfXDM3mrr1q2TJFVXV0d1ncvlUl5eXhxKBAAAUpFjwtFANTQ0aNSoURoxYoSmTZum++67T6NGjer1/K6uLnV1dQWfd3Z2SpI8Ho88Hk/cy5vKAvVHPcaGerQOdWkd6tIa1KN1YqlDw+fz+SwsS9xVV1dr1apV+uCDD/o9d/v27brgggs0btw4NTU16Xvf+57OnDmjV199VS6XK+w1lZWVwVYqs61btyorKyvW4gMAgAQ4ffq0Fi1apI6ODg0bNiyqa20NR70FEbMDBw5o0qRJwefRhKNztba2aty4cXr22WdVXl4e9pxwLUeFhYVqbW1VdnZ21F8T/+DxeFRXV6cZM2YoMzPT7uI4FvVoHerSOtSlNahH67S3tys/P39A4cjWbrWVK1dq4cKFfZ5TVFRk2dfLz8/XuHHj9NZbb/V6jsvlCtuqlJmZyQvVItSlNahH61CX1qEurUE9xi6W+rM1HOXk5CgnJydhX6+9vV1Hjx5Vfn5+wr4mAABwFsdM5W9ublZjY6Oam5vl9XrV2NioxsZGffjhh8Fzxo8fr507d0qSPvzwQ61evVr79u3TkSNH1NDQoDlz5ignJ0fz58+369sAAABJzjGz1b7//e/rySefDD6/+uqrJUn19fWaPn26JOnw4cPq6OiQJGVkZOj111/XU089pQ8++ED5+flyu93avn27hg4dmvDyAwAAZ3BMOKquru53jSPz2PIhQ4boxRdfjHOpAABAqnFMtxoAAEAiEI4AAABMCEcAAAAmhCMAAAATwhEAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgAnhCAAAwIRwBAAAYEI4AgAAMCEcAQAAmBCOAAAATAhHAAAAJoQjAAAAE8IRAACACeEIAADAhHAEAABgQjgCAAAwIRwBAACYEI4AAABMCEcAAAAmhCMAAAATwhEAAIAJ4QgAAMDEEeHoyJEjWrZsmYqLizVkyBB9+tOf1r333qvu7u4+r/P5fKqsrNTo0aM1ZMgQTZ8+XW+88UaCSg0AAJzIEeHozTff1NmzZ/XYY4/pjTfe0E9+8hNt3rxZd999d5/XPfDAA9q4caMefvhhHThwQHl5eZoxY4ZOnjyZoJIDAACnGWR3ASIxc+ZMzZw5M/j84osv1uHDh7Vp0yZt2LAh7DU+n09VVVW65557VF5eLkl68sknlZubq61bt+qb3/xmQsoOAACcxRHhKJyOjg6NHDmy1883NTWpra1NZWVlwWMul0vTpk3Tyy+/3Gs46urqUldXV/B5Z2enJMnj8cjj8VhU+vQUqD/qMTbUo3WoS+tQl9agHq0TSx06Mhz99a9/1U9/+lP9+Mc/7vWctrY2SVJubm6P47m5uXrnnXd6vW79+vVat25dyPH6+nplZWUNsMQwq6urs7sIKYF6tA51aR3q0hrUY+xOnz494GttDUeVlZVhg4jZgQMHNGnSpODzY8eOaebMmVqwYIFuueWWfr+GYRg9nvt8vpBjZmvXrlVFRUXweWdnpwoLC+V2u5Wdnd3v10PvPB6P6urqNGPGDGVmZtpdHMeiHq1DXVqHurQG9Wid9vb2AV9razhauXKlFi5c2Oc5RUVFwX8fO3ZMbrdbpaWlevzxx/u8Li8vT5K/BSk/Pz94/Pjx4yGtSWYul0sulyvkeGZmJi9Ui1CX1qAerUNdWoe6tAb1GLtY6s/WcJSTk6OcnJyIzm1paZHb7VZJSYmeeOIJnXde3xPtiouLlZeXp7q6Ol199dWSpO7ubu3evVs//OEPYy47AABITY6Yyn/s2DFNnz5dhYWF2rBhg9599121tbUFxxUFjB8/Xjt37pTk705btWqV7r//fu3cuVN/+tOftHTpUmVlZWnRokV2fBsAAMABHDEgu7a2Vm+//bbefvttFRQU9Picz+cL/vvw4cPq6OgIPl+zZo0++ugjrVixQu+//76uvfZa1dbWaujQoQkrOwAAcBZHhKOlS5dq6dKl/Z5nDkqSv/WosrJSlZWV8SkYAABIOY7oVgMAAEgUwhEAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgAnhCAAAwIRwBAAAYOKIjWftFNjM9uTJk8rMzLS5NM7m8Xh0+vRpdXZ2UpcxoB6tQ11ah7q0BvVonZMnT0oK3ZQ+EoSjfrS3t0uSiouLbS4JAACIVnt7u4YPHx7VNYSjfowcOVKS1NzcHHXloqfOzk4VFhbq6NGjGjZsmN3FcSzq0TrUpXWoS2tQj9bp6OjQ2LFjg3/Ho0E46sd55/mHZQ0fPpwXqkWGDRtGXVqAerQOdWkd6tIa1KN1An/Ho7omDuUAAABwLMIRAACACeGoHy6XS/fee69cLpfdRXE86tIa1KN1qEvrUJfWoB6tE0tdGr6BzHEDAABIUbQcAQAAmBCOAAAATAhHAAAAJoQjAAAAE8JRFP75n/9ZY8eO1eDBg5Wfn6/Fixfr2LFjdhfLcY4cOaJly5apuLhYQ4YM0ac//Wnde++96u7utrtojnTfffdp8uTJysrK0ogRI+wujmM8+uijKi4u1uDBg1VSUqK9e/faXSRH2rNnj+bMmaPRo0fLMAz96le/srtIjrR+/Xpdc801Gjp0qEaNGqV58+bp8OHDdhfLkTZt2qTPfe5zwYU0S0tL9bvf/S6q/4NwFAW3263//M//1OHDh7Vjxw799a9/1U033WR3sRznzTff1NmzZ/XYY4/pjTfe0E9+8hNt3rxZd999t91Fc6Tu7m4tWLBAt912m91FcYzt27dr1apVuueee3To0CFNnTpVs2bNUnNzs91Fc5xTp07pyiuv1MMPP2x3URxt9+7duv3227V//37V1dXpzJkzKisr06lTp+wumuMUFBTo3//933Xw4EEdPHhQX/ziFzV37ly98cYbEf8fTOWPwfPPP6958+apq6uL3ZNj9KMf/UibNm3S3/72N7uL4ljV1dVatWqVPvjgA7uLkvSuvfZaff7zn9emTZuCxyZMmKB58+Zp/fr1NpbM2QzD0M6dOzVv3jy7i+J47777rkaNGqXdu3fr+uuvt7s4jjdy5Ej96Ec/0rJlyyI6n5ajAXrvvff0i1/8QpMnTyYYWaCjo2NAmwMC0eru7tarr76qsrKyHsfLysr08ssv21QqoKeOjg5J4vdijLxer5599lmdOnVKpaWlEV9HOIrSd77zHX3qU59Sdna2mpub9etf/9ruIjneX//6V/30pz/V8uXL7S4K0sCJEyfk9XqVm5vb43hubq7a2tpsKhXwDz6fTxUVFZoyZYomTpxod3Ec6fXXX9cFF1wgl8ul5cuXa+fOnbr88ssjvj7tw1FlZaUMw+jz4+DBg8Hzv/3tb+vQoUOqra1VRkaGvva1r4meSb9o61KSjh07ppkzZ2rBggW65ZZbbCp58hlIXSI6hmH0eO7z+UKOAXZYuXKlXnvtNW3bts3uojjWZZddpsbGRu3fv1+33XablixZoj//+c8RXz8ojmVzhJUrV2rhwoV9nlNUVBT8d05OjnJycvSZz3xGEyZMUGFhofbv3x9Vc12qirYujx07JrfbrdLSUj3++ONxLp2zRFuXiFxOTo4yMjJCWomOHz8e0poEJNq3vvUtPf/889qzZ48KCgrsLo5jnX/++brkkkskSZMmTdKBAwf04IMP6rHHHovo+rQPR4GwMxCBFqOuri4ri+RY0dRlS0uL3G63SkpK9MQTT+i889K+EbOHWF6X6Nv555+vkpIS1dXVaf78+cHjdXV1mjt3ro0lQzrz+Xz61re+pZ07d6qhoUHFxcV2Fyml+Hy+qP5Wp304itQrr7yiV155RVOmTNGFF16ov/3tb/r+97+vT3/607QaRenYsWOaPn26xo4dqw0bNujdd98Nfi4vL8/GkjlTc3Oz3nvvPTU3N8vr9aqxsVGSdMkll+iCCy6wt3BJqqKiQosXL9akSZOCLZfNzc2MexuADz/8UG+//XbweVNTkxobGzVy5EiNHTvWxpI5y+23366tW7fq17/+tYYOHRps2Rw+fLiGDBlic+mc5e6779asWbNUWFiokydP6tlnn1VDQ4NeeOGFyP8THyLy2muv+dxut2/kyJE+l8vlKyoq8i1fvtz397//3e6iOc4TTzzhkxT2A9FbsmRJ2Lqsr6+3u2hJ7ZFHHvGNGzfOd/755/s+//nP+3bv3m13kRypvr4+7OtvyZIldhfNUXr7nfjEE0/YXTTH+cY3vhF8b1900UW+L33pS77a2tqo/g/WOQIAADBhoAcAAIAJ4QgAAMCEcAQAAGBCOAIAADAhHAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgAnhCAAAwIRwBAAAYEI4ApBW3n33XeXl5en+++8PHvv973+v888/X7W1tTaWDECyYONZAGln165dmjdvnl5++WWNHz9eV199tf7pn/5JVVVVdhcNQBIgHAFIS7fffrv+67/+S9dcc43++Mc/6sCBAxo8eLDdxQKQBAhHANLSRx99pIkTJ+ro0aM6ePCgPve5z9ldJABJgjFHANLS3/72Nx07dkxnz57VO++8Y3dxACQRWo4ApJ3u7m594Qtf0FVXXaXx48dr48aNev3115Wbm2t30QAkAcIRgLTz7W9/W88995z++Mc/6oILLpDb7dbQoUP1m9/8xu6iAUgCdKsBSCsNDQ2qqqrS008/rWHDhum8887T008/rZdeekmbNm2yu3gAkgAtRwAAACa0HAEAAJgQjgAAAEwIRwAAACaEIwAAABPCEQAAgAnhCAAAwIRwBAAAYEI4AgAAMCEcAQAAmBCOAAAATAhHAAAAJv8fxmf9o5CT9lgAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "a = 2.5\n", "b = 1.1\n", "m = 50\n", "theta = np.linspace(0, 2*np.pi, m)\n", "sigma = 0.1\n", "xi = a * (np.cos(theta) + sigma*np.random.randn(m))\n", "yi = b * (np.sin(theta) + sigma*np.random.randn(m))\n", "plt.plot(xi, yi, 'or'); plt.axis(\"equal\"); plt.grid()\n", "plt.xlim([-3, 3]); plt.xlabel('x'); plt.ylabel('y');" ] }, { "cell_type": "markdown", "id": "bb8405e5-978b-4c48-b76e-06857d5e8b32", "metadata": {}, "source": [ "Par ces points, on va essayer de faire passer \"au mieux' une ellipse de paramètres d'axe\n", "$\\hat a$ et $\\hat b$, d'équation \n", "\n", "$$\n", "\\frac{x^2}{\\hat a^2} + \\frac{y^2}{\\hat b^2} = 1\n", "$$\n", "\n", "ou, de manière équivalente, d'équations paramétrées\n", "\n", "$$\n", "\\begin{array}{l}\n", "x(\\theta) = \\hat a \\cos(\\theta), \\\\\n", "y(\\theta) = \\hat b \\sin(\\theta), \\quad \\theta\\in[0, 2\\pi].\n", "\\end{array}\n", "$$\n", "\n", "On pourrait chercher $(\\hat a,\\hat b)$ qui réalisent le minimum de\n", "\n", "$$\n", "\\min_{(a,b)} \\quad \\frac{1}{2} \\sum_{i=1}^m \\left[ \\frac{x_i^2}{\\hat a^2} + \\frac{y_i^2}{\\hat b^2}-1\\right]^2\n", "$$\n", "\n", "mais ce problème ne peut pas être mis sous la forme $\\min_u \\frac{1}{2} \\|Au-b\\|^2$. On considère plutôt le problème aux moindres carrés\n", "\n", "$$\n", "\\min_{u=(u_1,u_2)} \\quad \\frac{1}{2} \\sum_{i=1}^m \\left[ x_i^2\\, u_1 + y_i^2\\, u_2 -1\\right]^2\n", "$$\n", "puis on affectera $\\hat a = \\frac{1}{\\sqrt{u_1}}$ et $\\hat b = \\frac{1}{\\sqrt{u_2}}$. Ecrire ce problème de minimisation sous la forme\n", "\n", "$$\n", "\\min_u \\quad \\frac{1}{2}\\, \\|Au-b\\|^2\n", "$$\n", "\n", "où l'on déterminera la matrice $A$ et le second membre $b$. Résoudre les équations normales, calculer\n", "$\\hat a$ et $\\hat b$. Sur le même graphique, tracer le nuage de points $(x_i,y_i)$ et l'ellipse de pramètres\n", "d'axe $\\hat a$ et $\\hat b$." ] }, { "cell_type": "code", "execution_count": null, "id": "da6e3f75-6ced-432b-84af-f7dd03d90b9a", "metadata": { "tags": [] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "bdf77217-9360-404d-8649-25b1a7b1e007", "metadata": {}, "source": [ "#### Exercice 3. Moindres carrés avec régularisation de Tykhonov.\n", "\n", "Soit $L=3$. On considère le nuage de points $(x_i,y_i)$, $i=1,...,m$ avec $x_i\\in [0,L]\\ \\forall i$ tel que défini ci-dessous :" ] }, { "cell_type": "code", "execution_count": 5, "id": "528e7d74-76d4-4192-adc8-9805891a7647", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXw9JREFUeJzt3Xt4VNW9P/73JEwGgkwiCSRRbpFLIKHcVeIFbyVcKtLq+X2rqMdqyzmo0Et+/jgieISqRZ/allYqqEUp5iA+PcgRDjESvxpQDCgQrgEEGsBiEpoAGUxkMpnZvz/iDpnMZa89s/eevWfer+fBx0z2zN5ZWZn57LU+67NskiRJICIiIrKIpFhfABEREZEaDF6IiIjIUhi8EBERkaUweCEiIiJLYfBCRERElsLghYiIiCyFwQsRERFZCoMXIiIispRusb4Arfl8Pnz99dfo1asXbDZbrC+HiIiIBEiShIsXL+Kqq65CUlL4sZW4C16+/vpr9O/fP9aXQURERBH46quv0K9fv7DHxF3w0qtXLwDtP7zT6dT9fB6PB1u2bEFRURHsdrvu57MitpEYtpMytpEytpEYtpMyo9vI5XKhf//+HZ/j4cRd8CJPFTmdTsOCl9TUVDidTv4BhMA2EsN2UsY2UsY2EsN2UharNhJJ+WDCLhEREVkKgxciIiKyFAYvREREZCkMXoiIiMhSGLwQERGRpTB4ISIiIkth8EJERESWwuCFiIiILCXuitQRERHpxeuT8HnNOZy9eAl9e3XHdbm9kZzEffSMxuCFiIhIQNnBWizZVI3apksdj+WkdcczM/IxdWRODK8s8XDaiIiISEHZwVo8WrLHL3ABgLqmS3i0ZA/KDtbG6MoSE4MXIiKiMLw+CUs2VUMK8j35sSWbquH1BTuC9MDghYiIKIzPa84FjLh0JgGobbqEz2vOGXdRCY7BCxERURhnL4YOXCI5jqLH4IWIiCiMvr26a3ocRY/BCxERURjX5fZGTlp3hFoQbUP7qqPrcnsbeVkJjcELERFRGMlJNjwzIx8AAgIY+etnZuSz3ouBdA1eVqxYgVGjRsHpdMLpdKKwsBDvv/9+yOMrKipgs9kC/h05ckTPyyQiIgpr6sgcrHhgHLLT/KeGstO648+zxiKtRwre23sGlScauerIALoWqevXrx9eeOEFDBkyBADw17/+FTNnzkRVVRUKCgpCPu/o0aNwOp0dX/fp00fPyyQiIlI0dWQOJudn+1XYPd/cimc3s3Cd0XQNXmbMmOH39fPPP48VK1Zgx44dYYOXvn37Ij09Xc9LIyIiUi05yYbCwRkA2gvXPb52T0D9F7lw3YoHxjGA0Ylh2wN4vV787W9/Q3NzMwoLC8MeO3bsWFy6dAn5+flYtGgRbrvttpDHut1uuN3ujq9dLhcAwOPxwOPxaHPxYcjnMOJcVsU2EsN2UsY2UsY2EhNtO3l9EhZvPBSycJ0NwJJNh3Dr0AzL5sIY3ZfUnMcmSZKuk3MHDhxAYWEhLl26hCuuuAJr167F9OnTgx579OhRbNu2DePHj4fb7cZbb72FlStXoqKiApMmTQr6nMWLF2PJkiUBj69duxapqama/ixEREQAcKzJhuXVyYrHzc33Ymgac2BEtLS0YNasWWhqavJLHQlG9+CltbUVp0+fxoULF7B+/Xr85S9/wdatW5Gfny/0/BkzZsBms2Hjxo1Bvx9s5KV///5oaGhQ/OG14PF4UF5ejsmTJ8Nut+t+PitiG4lhOyljGyljG4mJtp027a9F8d8OKB73+//ne5gxyppTR0b3JZfLhczMTKHgRfdpo5SUlI6E3QkTJuCLL77AH//4R7z66qtCz584cSJKSkpCft/hcMDhcAQ8brfbDf3DNfp8VsQ2EsN2UsY2UsY2EhNpO+Wk9xQ+zuq/B6P6kppzGF7nRZIkv5ESJVVVVcjJsWbUSkRE8YmF62JL15GXp556CtOmTUP//v1x8eJFrFu3DhUVFSgrKwMALFiwAGfOnMGaNWsAAMuWLcOgQYNQUFCA1tZWlJSUYP369Vi/fr2el0lERKSKXLju0ZI9sAF+ibssXKc/XYOX+vp6PPjgg6itrUVaWhpGjRqFsrIyTJ48GQBQW1uL06dPdxzf2tqKJ554AmfOnEGPHj1QUFCAzZs3h0zwJSIiihW5cN2STf51XrJZ50V3ugYvq1atCvv91atX+309f/58zJ8/X8crIiIi0k6wwnXX5fbmiIvODKvzQkREFI86F64jYzB4ISKiuOT1SRwRiVMMXoiIKO6UHawNyEXhnkPxw/Cl0kRERHoqO1iLR0v2+AUuwOU9h8oO1sboykgrDF6IiChueH0SlmyqDrnnEAAs2VQNr48l+62M00ZERBQ3Pq85FzDi0pkEoLbpEnadOm/cRX2HOTjaYfBCRERx4+zF0IGL/3FuKG+rqB3m4GiL00ZERBQ3+vbqLnhc4J54wXh9EipPNOK9vWdQeaIxoukm5uBojyMvREQUN+Q9h+qaLgXNe7GhvQLuhIFX4oPD4V9Li9ESpRwcG9pzcCbnZ3MKSQWOvBARUdyQ9xwCELBpopo9h7QaLRHNwfm85pzQ61E7Bi9ERBRX5D2HstP8p5Cy07pjxQPjFEdNtFyxJJ6DI3YcteO0ERERxZ1o9hxSM1qitC2AeA6O2HHUjsELERHFpUj3HNJytEQ0B+e63N7qLjLBcdqIiIioEy1HS7TKwSF/DF6IiIg6kUdLQoUTNrSvOhIdLYk2B4cCcdqIiIioE3m05NGSPbABftM9kY6WRJODQ4EYvBAREXUhj5Z0rfOSHUVV3EhzcCgQgxciIqIgOFpiXgxeiIiIQuBoiTkxeCEiUom7AxPFFoMXIiIVuDswUexxqTQRkSDuDkxkDgxeiIgEaLnfDRFFh8ELEZEA7g5sDV6fhMoTjXhv7xlUnmhkMBmnmPNCRCSAuwObH/OREgdHXoiIBHB3YHNjPlJiYfBCRCRA6/1uSDtWzUfiFFfkOG1ERCRAj/1uSBtq8pHMUnCOU1zR4cgLEZEg7g5sTtHmIxk9AhLNFBdHa9px5IWISAXud2M+0eQjfXCoHs+/f9SwERClKS4b2qe4JudnB/QpjtZcxpEXIiKV5P1uZo65GoWDMxi4xFik+Uj7Gm2Yt26foUm+kS65Z0KyP12DlxUrVmDUqFFwOp1wOp0oLCzE+++/H/Y5W7duxfjx49G9e3dcc801WLlypZ6XSEQUFQ7jx56cjwQgIIAJlY/k9Ul492SSYUm+cj95XzDI6DzFZdWEZD3pOm3Ur18/vPDCCxgyZAgA4K9//StmzpyJqqoqFBQUBBxfU1OD6dOnY/bs2SgpKcH27dvx2GOPoU+fPrjnnnv0vFQiItU4jG8eU0fm4N8m5eL1T2ogdfkM/8Go9qm+znadOo8LraFHzLRM8g3WT5R0nuKyYkKy3nQNXmbMmOH39fPPP48VK1Zgx44dQYOXlStXYsCAAVi2bBkAYMSIEdi1axdeeuklBi9EZBiRXaPlYfyu97ryMD4TeI1VdrAWr22rCfh9SAD+d38tPj3egBfu/l7H7+TsRbfQ60ZbdDBUPwnFhvYE8M5TXCyQGMiwhF2v14u//e1vaG5uRmFhYdBjKisrUVRU5PfYlClTsGrVKng8Htjt9oDnuN1uuN2XO6HL5QIAeDweeDweDX+C4ORzGHEuq2IbiWE7KTOijT44VI/nSo+gznX5fSXb6cCi6cMxpSALQHtws3jjIYWky0O4dajx+TCJ2I/C/T5kF1o8mFOyB8vvHY0pBVno3SNZ6LUzUrtF3JYi19WZ3FMWTsuDz9sGn/fyNYiI5lqDMbovqTmP7sHLgQMHUFhYiEuXLuGKK67Ahg0bkJ+fH/TYuro6ZGVl+T2WlZWFtrY2NDQ0ICcn8C5m6dKlWLJkScDjW7ZsQWpqqjY/hIDy8nLDzmVVbCMxbCdlerXRvkYb3vhSTgW8HHTUuS5h7rq9eGSYD6MzJBxrsqHOFfrDr30Y343l75RhaFps8hASqR8p/T4uk7Do3b1w13hxwmVDarcktLQBgZky7cempwD/rN6B0sN6X1e7tBQJdw/ywXtqN0pPXX7cJwHpKcm40KrftYZjVF9qaWkRPlb34CUvLw979+7FhQsXsH79ejz00EPYunVryADGZvP/xUjfTV52fVy2YMECFBcXd3ztcrnQv39/FBUVwel0avRThObxeFBeXo7JkycHHRkitpEotpMyPdvI65Ow9HfbAASbTrDBBuD9+lTMv38SvAfrgOoDiq95odcgZIzIwoSBVxo2ApOI/WjT/lqh3wdgw4VW4PkDPXCuJfRdvu27/z539+iO0TY9r+uB6/tjakH4fmIfVI956/YBCFYgMfprDcboviTPnIjQPXhJSUnpSNidMGECvvjiC/zxj3/Eq6++GnBsdnY26urq/B47e/YsunXrhoyM4ElIDocDDocj4HG73W7oH67R57MitpEYtpMyPdpo14lGv6miruTRlP/6/B8439Iq9JolO79Cyc6vYpLEm0j9KCe9p6rjwwUuQHvOiRa/L9Hr+sGoq4Mm2vrnXvXAn2eNxbObD/sl70ZzrSK5XYBxfUnNOQwvUidJkl+OSmeFhYXYtGmT32NbtmzBhAkTEuaPkIhiQzTZ8dnN6sflmcSrL7nOi5rVPMH0dCRj5f3jccOQTE1GyuTrqmu6FDTvJVhyrizUSranf5CPK3umRF0g0eor5XSt8/LUU0/hk08+wcmTJ3HgwAEsXLgQFRUVuP/++wG0T/n867/+a8fxc+bMwalTp1BcXIzDhw/jjTfewKpVq/DEE0/oeZlERLruBp2otTiMItd5iTbcaHZ78Yt39qK8uk75YBXXBYjXnwHCF6R7fO0eNH3bGlWBxHgoeKdr8FJfX48HH3wQeXl5uOOOO7Bz506UlZVh8uTJAIDa2lqcPn264/jc3FyUlpaioqICY8aMwbPPPos//elPXCZNRLpTqtIarVCVU0kb8r5T6anRjdKfa27V9ANc7X5Yeheki5eCd7pOG61atSrs91evXh3w2C233II9e/bodEVERMGF2zVaRFF+FrZU1ysel0i1OIwm7zu1/KNjeHP7SVz49nJuS++edpxrFl+KG2p/oWiuSyS/RO+CdGpef8IA/Re9RIobMxIRfUe+S1ZbDRUAhmVdIRS86Dk9Re1B6C++Pwxzbx/qFyyMH3glbvntxyHzTzrTo2KtvB+WEr0L0ql7fQYvRESW0PUuueGiWyhJt/CaTKzfcyai5EzSXrBgQR5ZE2XkKJm88udY/UWh40WD4K4rijKvCFydG83rxwqDFyKiLjp/8Hl9Ev7yaY1iUDJxcEbIaadwyZlkHHlk7akNB4SmkCL5ABddftyZmr2P1ATBwV432+lAeqodTS0exSDb521TPEesMHghIgojXC5M16Ak1LSTVnVDElkkQUEwU0fm4PbhWZj4mw9xrqUVwSrWRjpKFsnyYzV7H6kJgkO9br3L3fGYUn+WtycwIwYvREQK1AQlapIzSYzWNUlSuiXh13flY+66vZqNkkWyUWe4lT/BZKd1x73XDoC7zYfKE40h+5XSiiIbgLRUO7p3S0ady5pBNoMXIiIBaoIS0eRMUqbX7t1TCrLwyDAfSutS/TfhjOADXCRYCLZ6SWnlT2dXOJLh8/nwhw+/7HgsVAAnsqLoQosH//XTcUhKslkyyGbwQkQkiEGJsSINCkSNzpAw//5JqPrHxag+wCNd3qwmIfgbtxffuP3ncUIFcKKv29DsxswxVwtfg5kweCEiIlPSu+YJoE1AGuny5mhX9IQK4ERf1+wrisLRtcIuERFRpPSueaKVSIMFuapzNIJVblaqFm1D+5STlZftM3ghIgrC65NQeaIR7+09g8oTjaYvlx6PrDKCEGmw0Hnvo2h1DuAi3VPJShi8EJkUPzxjp+xgLW568SPc9/oO/GLdXtz3+g7c9OJHltiwLp5YZQQhmmBh6sgcrNRgT6auAZzaPZWshjkvRCZk9e3qrUyv1S2knpoaO+EEqxGjtWhq/Fzek+k43txe47cnU05ad3zr8QoVlQv1uvG4bJ/BC5HJfHCoHvPW7Qt4o6rlh6fu9F7dInoN8fhhE6loC/+FuhFYOC1Pl2uNNFho35NpKObePiTg+eXVdREHcCIJyVbscwxeiEzEJwFLS4+ELFolAVjw7gFdPzwTmRGrW8LhiFtwkQYF4UbR5q3bh4eH2TBd42uNdvVSsOdrWbm5a6Byvrl9765gfe6OvMyIfw69MXghMpETLptfwaxgzrd4sPyjY/jF94cZdFWJI5arWzhdFVzXD9s7R10lFLiLjKK9ezIJ830Soss2MYYWU0CieyjJfe7le0dHe9m6YfBCZCIu5b3iAABvbj+JubcP5eiLxmK1usUM01VmpHYkqnOg03DRrVxlttWGXafO46ZhWULXo/X0itrXi2ZUR80eSnKfe/79I5g/IqLT6Y7BC5GJOAVvAS9869Ft6iKRyatblHaQ1jrhM9bTVWakdiRKzc7MnZ29GH6kM9zrRzOlZ+QUodo9lAC5z7lxwmXOYJlLpYlMZLBTQloPsQgm1oW54lGkS16jXdZulWJsRlEaiQLaR6LkdpYDHbWBCwD07eVQPCbU68uBlNol9Fq/Xldd++OOvzdG1DaA+Giw0TjyQmQiSTbgJ4UD8MePTigeG+vCXPFKbXKkFnfQVinGZhQ1I1HX5fZWPaoAfLezcoqECQOvDHuc1lN6ek8RBuuP6YI3RMGIjgYbjcELkck8ess1+OuO07jQEvyWR6+pC7qcg+Bu8+GlfxkN2ICGb9wh8xG0SrKN1XSVWakZiVKzM3NXdw/yKQYIWk/p6TlFGKo/dq4bI6q9zzkw2Nms+rlG4LQRkckkJ9nwwt3fC/q9eCntbUZdq+rev2onnvjbPji6JaFwcEbQqSI1UxvhJEI5dzXUjERFOpX20xsHYnSG8u9G6yk9vaYII8lrCUXuZQunDYdZuxyDFyITkkuGd920LV5Ke5tNJDkIau6gRcR7OXc11GwLEOlU2uaDdRBJTdJ6Sk+PKUKvT8Lq7TURj0B1leV0YMUD4zClQGwVVixw2ojIJLw+CceabNi0vxY56T0xOT87bkt7m0mkOQh63EHHczl3NdRsC6A05RaK6Eoaraf0tH49taus0nvY/aaRrky1Q5K6Ti2Zv78xeCEygbKDtVi88RDqXMlA9QEArKxqlEhzEPRKso22Qmu8CJU43btnCmaOuQppPVLg9UlhAx0lIitp1ARSInVbtNqvCVBXu0X251njkJRkw9mLl3CyoRl/+PBYwDH1LhapIyIFrKwaW5GOoDDJVn+dR6I+rK7Dhr1n0Njcije2n8Qb20/6BfjBAh0loitpRFagqVl1pkW5f7U5LnJ/nPhd/pbXJ+GmFz8KeiyL1BFRWKysGnuRjqBoeQdNoSUn2dD0bXvAohTgy4FOXdO3eHbzYZxvbg0TWKpbSRNuSi+SG5BopwjVrLIK1h/FRhxZpI6IgtA66ZPUU5Mc2hWTbPWnZlWXPOX2o3H98JsfjQQQevVWJCtp5NefOebqjhVo0aw66/x61+X2xuc15xQLHcoF6N5XUcguWH8UHXFkkToiCsDKqrEX7QgKk2z1FWlOktLUzB15mSg9Ffh6avcb0qJui+iUUyRbIDz9gxH4yY25AT+D6Igji9QRUQBWVjWHaHMQmGSrn2gC/HCBpccTOKQQSbXkaG9ARKec1CbnyjkuwQIXQDRny7xF6hi8EMUQkz7NQ/6g2/H3RlSeaAQgofCaTExkUBJTJxvEPjxDBfiigWWkifPR3ICI5rzdPjxLdXIuoLwKSmnEceG04fCe2i14VmPpmvOydOlSXHvttejVqxf69u2LH/7whzh69GjY51RUVMBmswX8O3LkiJ6XShQTrKwaW103sPvgYB2e+Ns+LP/4OJZ/fAL3r9qJm178KOqN8igyXp+Etz8/rXhcqJwkNeeJNG8lmpwp0SmntypPqpoq6pzj0rVy9H2v7+jo00o5WwlbpG7r1q14/PHHce2116KtrQ0LFy5EUVERqqur0bNnz7DPPXr0KJxOZ8fXffr00fNSiWJGfgNpr/Pi7nhczbJJEtP5DvRkQzPe/vy0X5sHwyXrsfN5zTnF3w8A3HvtgKgC/GjyVqLJmRKdcjp1rkXouH8tHIhpI3NUr4JSM7VmFroGL2VlZX5fv/nmm+jbty92796NSZMmhX1u3759kZ6eruPVEZnH1JE5uHVoBpa/U4ZrCsYgJ70nkz41FkmyI8Al67Ek+uE+KDPVkPOEOi7SnCnRKaeBvcV+vmkjczqCK7VlGKyWs2VozktTUxMAoHdv5eG9sWPH4tKlS8jPz8eiRYtw2223BT3O7XbD7b4cmbtcLgCAx+MxJGqUz2HmCDXW2EZifN42DE2TMHlEJux2O3zeNvi8sb4qc4m0L31wqB7z1u2LeNM6+c678vhZXG/y/KN4+nvLSBX7iMpI7ab65+3cTlqc5468TNw69GbsOnUeZy+60beXAxMGXgkA2HqkDju/K3cwMbd3x43J2H69kO10oN7lDps0e++Eq/H6J39XPG5sv14d17dTcDQpXJ82ui+pOY9NkiQtNqFUJEkSZs6cifPnz+OTTz4JedzRo0exbds2jB8/Hm63G2+99RZWrlyJioqKoKM1ixcvxpIlSwIeX7t2LVJTo4vGicj6fBKwZE8yLrQC0e7Z8q9DvRifachbJkHkdychPQV4Zpw3qt2P9TrPvkYb1v09CS1t/k9K7Sbh3mt8GJ0hYV+jDW98Kaefdj6uvZ89MkzdcbLdDTasOZaseI1m6tMtLS2YNWsWmpqa/NJGgjEseHn88cexefNmfPrpp+jXr5+q586YMQM2mw0bN24M+F6wkZf+/fujoaFB8YfXgsfjQXl5OSZPngy73aQL4mOMbSSG7aQskjbaWXMOD7yxS5PzlzwywRIjL/HUj+RRMyB4PsnL946OKLG0aztpfZ4PDtVj7nevF8ry717zg0P1eK70iF9+T06aAwunDfc7p+hxgHi/D9enje5LLpcLmZmZQsGLIdNG8+bNw8aNG7Ft2zbVgQsATJw4ESUlJUG/53A44HA4Ah632+2G/uEafT4rYhuJYTspU9NGjS1tUZ9PXrJeOKSvZXJe4qUf3TmmH7p1S45qH6Bw5HbS8jxen4TnSsOvrAWA50qPYNqoq3HnmH6YNupqxeJ4oscBQOGQvkJlGET6tFF9Sc05dA1eJEnCvHnzsGHDBlRUVCA3Nzei16mqqkJODrP8iUi9aAv8ccl67BlVxVir87SvklJOAq5zuTtWMIkmzao5Lp733tI1eHn88cexdu1avPfee+jVqxfq6uoAAGlpaejRowcAYMGCBThz5gzWrFkDAFi2bBkGDRqEgoICtLa2oqSkBOvXr8f69ev1vFQiilNKhQCVcMm6ORi1IkaL86jZzkPPrT+02L3arHQNXlasWAEAuPXWW/0ef/PNN/GTn/wEAFBbW4vTpy8XIWptbcUTTzyBM2fOoEePHigoKMDmzZsxffp0PS+ViOJUuDvQrnLSuuPpH4zAlT0d3KeIIqZmtE/vrT/ide8t3aeNlKxevdrv6/nz52P+/Pk6XRERJaJQd6A5ad1x77UDMCgzNW7e1Cn2rsvtjWxnd8Wpo2ynw5CtP6xYx0UJ9zYiooQQr3egZD7JSTYsvisfc0r2hD1u8V0FMel/anfONiMGL0SUMOLxDpTMaerIHKx8YByefPcALrT4F19LT7Xjhbu/F5Ock0h2zjYjBi9EREQ6CLdTeSxGOiLdOduMGLwQWVg8DP8SxbPkJBtuHJKJG4dkxvQ61O51ZHYMXogsKl6Gf4lIf9HsnG1GScqHEJHZyMO/Xd+M5OHfsoO1MboyIjKjaHfONhsGL0QWozT8C7QP/3p95thsjYhiT7SejN51Z7TCaSMii4m34d94FioniblKZDSlStPyXkdG1J3RAoMXohgI9uElKt6Gf60uVCASKifprtE52LivlrlKZKh42+uIwQuRwUJ9qC2clif0/JMNLULHWWX418rCBSivbasJuMOtbbqEV7fVBLyOFZeqkvXE015HDF6IDBSuzsK8dfvw8DAbwu3i5fVJePvz02GOaGdU2fFEFup3GSpACceKS1XJmuKl0jSDFyKDiNRZePdkEub7JNhDvMbnNecU90sBgB9f299yb0ZWEu53GSnmKpFR4qHSNIMXIoOIJNpeaLVhZ805pNjtQe+KRPNY/vrZKYzIcVpqGNhKlH6X0WCuEpEyBi9EBhH9UPr5O/vQ9G1bx9fZTgfuu24ABmX2RMNFt9BrXPjWwxwKHekZYDBXiUgZgxcig4h+KHUOXACgzuXGHz481vF1kg0QLeHCHAp96BFgWG2pKlEssUgdkUHkOgvhwwjlqEQ0cOmcQ0HaEvtdQvH7XY+z0lJVolhi8EJkELnOAhDuQ038g0v0SOZQaC/c79L23b9/n5SL7DT/EZqctO7490m5yOnyeHZad07xEanAaSMiA4Wqs5CeaseFFo+q1xJd6cIcCn2I1MyYP3VE0CWpoR4nIjEMXogMFqzOgs8n4f5VO1W/VnoPO5q+9cRFuW8rUqqZEWpJajwsVSWKJQYvRDHQ9cPL65OQ7XR8V8NF/A784RtzsezDL+Oi3LdVMRAhMh5zXohMIDnJhkXThwMQC11saM+fmHv7EKx4YFxAboWcQzE5PxuVJxrx3t4zqDzRyJ2miSgucOSFyCSmFGThkWE+lNalos4Vup5L11GVrlMXmT0dgA34v4fr8dSGgzjX3NrxXG4ASETxgMELkYmMzpAw//5JqPrHRZy9eAknG1rw9uen/bYECLaJmjx1UXawFk/8976Q1V+5ASARRSPULupGY/BCZDJdcyjm3j5E6M0i1EaBnXEDQCKKVKhd1GMxmsucFyKTk4OZmWOuRuHgjKABh5qNAlm8jojUkm+Ouo7qyqO5ZQdrDb0eBi9EcSCSjQJZvI4oMXh9UlSJ++FujuTHlmyqNnRBAKeNiOJAJIEIi9cRxT8tpnp2nTof9uao82iuUWUDOPJCFAfUBCLyMmsWryOKb1pN9ZwV3M3eyNFcBi9EJuH1STjWZMOm/bWqh3bVbhTI4nVkZtFOc5C2Uz19ezmEzmnkaC6njYgMIi8xrGv6FueaW9H7Cgeyne0jIOXVdVi88RDqXMlA9QEA6oZ25Y0CHy3ZE1Btt7Ngy6yJzMRMK1qsTCkPTs1Uz4SBVyInrTvqmi6ZZisSBi9EBgj2hiwLtSmj2posoTYK7N3Tjh+NuRrfz8/mBoAmZpb6GbEUark/6xOpJzqFI3JcuJujWI3m6hq8LF26FO+++y6OHDmCHj164IYbbsCLL76IvLy8sM/bunUriouLcejQIVx11VWYP38+5syZo+elEulGqf5KqN2kI6nJorRRIJkTRxuUpzlYn0gd0Skc0eNEdlE3kq7By9atW/H444/j2muvRVtbGxYuXIiioiJUV1ejZ8+eQZ9TU1OD6dOnY/bs2SgpKcH27dvx2GOPoU+fPrjnnnv0vFwizampvxJMJFn83CjQWjja0E7LaQ66nAen5VSPmW6OdA1eysrK/L5+88030bdvX+zevRuTJk0K+pyVK1diwIABWLZsGQBgxIgR2LVrF1566SUGL2RqwYb9I6m/EgxrssQnr0/Ck+8e4GgDtJ3mIP2mesxyc2RozktTUxMAoHfv0JFeZWUlioqK/B6bMmUKVq1aBY/HA7vd7vc9t9sNt/vyMi6XywUA8Hg88HiCD8drST6HEeeyqkRoow8O1eO50iN+GypmOx2YWpClyetnpHaL6/YTFW996eWPToScNgQujzZUHj+L6wXvkK3aRhmpYh9HWv0tWLWd1LgjLxMv3zs68L0pzYGF04bjjrzMsD+/0W2k5jw2SZIMWYMmSRJmzpyJ8+fP45NPPgl53LBhw/CTn/wETz31VMdjn332GW688UZ8/fXXyMnxHz5dvHgxlixZEvA6a9euRWpqqnY/AFEI+xpteONLuepA57sYKchjaklITwGeGedFnN94G84nASdcNrg8gNMODHZKhraxTwIWfpGMFq/ySf91qBfjM+N7ubBPApbsScaFViD43wz/FiIV674uqqWlBbNmzUJTUxOcTmfYYw0beZk7dy7279+PTz/9VPFYm82/VeX4quvjALBgwQIUFxd3fO1yudC/f38UFRUp/vBa8Hg8KC8vx+TJkwNGhahdPLeR1ydh6e+2AQhWxMnW8RYcyceO7bv/Pnf3aEzRaATH6rTqSx8cqsfSICNli6YPN6ytd9acQ8uOXULHFt18vaqRF6v+vdkH1WPeun0Agk1zaPu3YOV2MorRbSTPnIgwJHiZN28eNm7ciG3btqFfv35hj83OzkZdXZ3fY2fPnkW3bt2QkRE4z+ZwOOBwBBbQsdvthnZIo89nRfHYRrtONPp9AHalJmjpmZKM5lZvx9fRZvHH89LbaPpS2cFazFu3L+B3U+9yY966fYYlyDa2tAkdl55qR+GQvqp/d1b8e7tzTD9065Zs6IoWK7aT0YxqIzXn0DV4kSQJ8+bNw4YNG1BRUYHc3FzF5xQWFmLTpk1+j23ZsgUTJkxgByPT0TJ5cMldI3Cyeh+uKRiDnPSeEQUbcsBSXl2H/9n7Nc41t3Z8L9GW3gZjpuW4oktUH74hN26CThFmWtFC5qVr8PL4449j7dq1eO+999CrV6+OEZW0tDT06NEDQPu0z5kzZ7BmzRoAwJw5c7B8+XIUFxdj9uzZqKysxKpVq/D222/realEEdGyHHa2szvsaRKmj8qJKFAPVwgPSLylt8GYaTmu0lJWoH3UZe7tQ3S9DjMyy4oWMi9d9zZasWIFmpqacOuttyInJ6fj3zvvvNNxTG1tLU6fPt3xdW5uLkpLS1FRUYExY8bg2WefxZ/+9CcukyZTUtpTyIb2XIpsZ/hjctK6Y8LAKyO+jlAbsHUWq63rzcRMy3HlpaxA6JTuF+7+HkcciILQfdpIyerVqwMeu+WWW7Bnzx4drohIWyK1FBbfVQAAupXWVlMIL9ELfWlddTRaoaqWcoqPKDzubUQUJdGy2UrHyDUOvD4Ju040Cs/3R1IIL1ELfelRdTRazPEgUo/BC5EGRD6ARI7Z12jD0t9t81vBpHQXHkkgYuTW9WZixg3m5OvqOhIWz6vFiKLF4IVIIyJJhuGO+eBQ/XfF7vyXXisl2qoJRGIxsmA2ZttgLhhu1EgUHoMXIhPw+iQ8V3ok6PeUlvCKrFoBYjuyYCZen4S0HimYPyUP55pb0fuK9oRqs4xscKNGImUMXohM4POac99NFQX/8AyXaBtuKqQzM40sxEq4EQ0zBC5mqkNDZGa6LpUmIjHRLuGVp0Ky0/ynkHr3tOOnNw7C27Mn4tP/uD3hA5dgy8nlEY2yg7UxurLL1NShIUpkHHkhMgEtlvBy1UpoVhnRMFMdGqMwMZkiweCFSGORvBlfl9sb2U4H6lyXEGzqSDTRlpVJgzNTZd1wzFaHRm9MTKZIMXgh0lCkb8bJSTYsmj4cc9ftNdUS3nhhlRENM9ah0QsTkykazHkh0ki0ORVTCrLwyDAfspz+u6Rnp3XnG3mUrDKiEW7LgHgKYpWm8YDE3saClHHkhUgDWuVUjM6QMP/+Saj6x0XdcwASKdfASiMaVqhDEy2rTOOReTF4IdKAlm/GeuWtdA5WTja04O3PT3+XY9MunnMNzFpZN5R4T762yjQemReDFyINmP3NOFguTlfxnmtgtRGNeE6+tso0HpkXgxciDYi+yTZcdOO9vWcMvZMOlRjZlZmWDOsl3kc0rMJK03hkTgxeiDQgUqI/yQY8u/lwx9dGTNOEy8UJJhFyDeJ5RMMqrDaNR+bD1UZEGgi3SkTWdeGEEZVdlXJxQmGuAektVFVorq4jERx5IdJIqJyKJFtg4AIETtPoIdIghLkGZARO41GkGLwQaajrm3HDRbffVFFXnadpJgxwan49aoMQ5hqQ0TiNR5HgtBGRxuQ345ljrkZmL4fyE6DfNI2ciyNyH8tcAyKyCgYvRDqK9ZJQkVwcGXMNiMgqOG1EpCM1S0J93jZdriFkfROnA/ddNwCDMntaPtcgkaoFExGDFyJdqVkS6vPqdx3xnBgZTzsTMwgjEsPghUhnZqns2jUx0uuTUHmi0dIflPG0M3E8BWFEemPwQmQAs418xMMHpVabYZpBPAVhREZgwi6RQTqvQiocnBHTwOXRkj0BxeuMKJqnJTWbYZqZUhAGtAdh3mDFgogSFIMXogQSTx+UZt8MU1S8BGFERuK0EVECUfNBafbCYZEuQzdbUqzVgzCztSclBgYvRAnE6h+UnUWyM7EZc31iXQsoGmZsT0oMnDYiSiBW/qDsKlwBvmDVgs2a66NUBdmG9oDAbFs2aNme8sq39/aeQeWJRktMW1JsMXghSiDX5fZGtjN0YGLWD8pQRHcmNnOuj9ogzAy0bM+yg7W46cWPcN/rO/CLdXtx3+s7cNOLH1kmcZxig9NGRAmkvLoOl9qCV8Mz6welEpFl6GbP9TFLLSBRWrUnl4hTpHQNXrZt24bf/va32L17N2pra7Fhwwb88Ic/DHl8RUUFbrvttoDHDx8+jOHDh+t4pUTxL9QHhSw91Y6ld3/Pkh8WSjsTWyHXx2y1gMKJtD07J/dmXuHA4o2H4qJODxlP1+ClubkZo0ePxsMPP4x77rlH+HlHjx6F0+ns+LpPnz56XB5Rwgg3zC9zdEvC5Pxsw67JSFbJ9VEKwsKRA4PaC834e5MNXp8Eu8bXJ4ukPYMl94YT69EwMjddg5dp06Zh2rRpqp/Xt29fpKena39BRAlKaZgfAOpc7rj9oIhkZZKVBAYGyXjrhY/xyI25mHv7UM1HLtS2p9KoXzhWWPlGxjNlzsvYsWNx6dIl5OfnY9GiRUGnkmRutxtut7vja5fLBQDweDzweDy6X6t8DiPOZVVsIzF6tlPthWbh4zwep/KBMRJNGy2clod56/aF3CBz4bQ8+Lxtum6QqYcPDtVj3rp9AYFB07dt+MOHx/Dm9pN4bmY+phRkaXpe0fb0eKSQ00MiMlK76fI3wfclZUa3kZrz2CRJMiS93mazKea8HD16FNu2bcP48ePhdrvx1ltvYeXKlaioqMCkSZOCPmfx4sVYsmRJwONr165FamqqVpdPcc4nASdcNrg8gNMODHZKiKdp9mNNNiyvTlY8bm6+F0PT4neZ6r5GG949mYQLrZd/uekpEu4e5MPoDOv93D4JWLInGRdagcC1SrL2n+uRYdr/jCLtKdr3AklITwGeGeeNq79FCq2lpQWzZs1CU1OTX+pIMKYKXoKZMWMGbDYbNm7cGPT7wUZe+vfvj4aGBsUfXgsejwfl5eWYPHky7Ha9Zpitzext9MGhejxXegR1rsv9KNvpwKLpwzW/Ww1Hz3by+iTc+rttqHe5wwzzO/Bx8SRTJ0dq0UZen4Rdp87j7EU3+vZyYMLAK039M4ezs+YcHnhjl9CxOTr9fpXac9P+WhT/7YCq15Sf/fK9o3X7GzT7+5IZGN1GLpcLmZmZQsGLKaeNOps4cSJKSkpCft/hcMDhcAQ8brfbDe2QRp/PiszYRmUHa4MOude73Ji3bl9Mlmrq0U52AIvvKsCjJXtCDvM/M6MA3R0pmp5XL9G0kR3ATcOMC0r11NjSJnxsbZMbVf+4qHlOk1J75qT3VP2aRi4RN+P7ktkY1UZqzmH64KWqqgo5OdZbuknmp1RoK96WakZaS4R715iX2tVRsUh+FU3ufelfRqOh2c0+RkJ0DV6++eYbHD9+vOPrmpoa7N27F71798aAAQOwYMECnDlzBmvWrAEALFu2DIMGDUJBQQFaW1tRUlKC9evXY/369XpeJiUosxcu04PaWiJq965hoGMsOTAQXX4ci6XgcgXh8KN++bhxaKbh10bWpWvwsmvXLr+VQsXFxQCAhx56CKtXr0ZtbS1Onz7d8f3W1lY88cQTOHPmDHr06IGCggJs3rwZ06dP1/MyKUFZoXCZHkRriaitfspN+ownBwZzSvaEPS7WS8GtVkGYzE/X4OXWW29FuHzg1atX+309f/58zJ8/X89LIupglcJlsaB2So1l3mNn6sgcrHxgHJ589wAutAQuNTXLtg9WqiBM5seNGSlhWXU3XyOomVIz86aHiWLqyBzsXjQZv/r+MKT18E967LpJpV5EdoaWR/1mjrkahYMzGLhQxEyfsEukF9G5+K5vsImQ16FmSi0Rc4fMKDnJhl98fyj+/eaBWP5OGa4pGIOc9J6G9E9OGZLRGLxQQlM7F58ob9JqptQSNXfIrJKTbBiaJmH6qBxDlrdyypBigcELJTzRufhEepNWs3fN5zXnhF4zEXOH4l2ilRsg82DwQgTlFTiJ9iatZkot3jc9pNBEpwx3nGhEUpItrqdayVgMXiihieavJGJeh+iUWqS5Q6Q9r0/Czppz2N1gQ0bNORQO6atru4tOBT6+dg8ufHt5JVQ8TrWSsRi8UMJSk7+SqHkdolNqrOMRe/79ORlrju3SPUgQnQrsHLgA8TnVSsZi8EIJSU3+itcnoeGiO/BFgjBrXkc0K6REi9qxjkfsRJuPFWn/UJoyDCUep1rJWAxeKOGoyV8pr64LGE0Ixsx5HUaukBINdEg70eZjRdM/Ok8ZqhWPU61kHBapo4Qjmr+y/KPjeLRkj1DgApgzr0O+I+/6M8h35GUHa2N0ZaQVNflYXWnRP6aOzMG/TcpVfd2yeJtqJWMweKGEI/pm+eb2GqGhcKMqmKrFyreJIdJ8LK36h9cnYeO+yINgs061krkxeKGEE2mSYTALpw/Hp/9xu+kCFyC6O3Kyjkj36NKqfyi9TiiJvP0GRY/BCyUckT2N0nuIVSb988cnUF5dp9m1aSlRV0glmkj36NKqf3wYQf8381QrWQODF4sQ2fSMxMhJhgAC3vDlrx++cZDQa1341mPa3BHump0YRPpzsCBBi/5RdrAWq7afFLzSy8w61UrWwdVGFpAo++kYSakuyeT8bKz74ivhJaBmXPLJyreJI5I6O9H2DzlnRokNQJbTgd/9nzFo+MbNJfSkCQYvJpdI++kYTakuiegSULMu+WTl2/iiVItF7s+Vx89iyyc7UXTz9WEr7EbbP0RzXSQAi+8qwI1DMsV+UCIBDF5MTHQ1gNnu+K0kXF0S+W72yfUHhJJ3zZg7wsq38UF09DU5yYbrc3uj8bCE6wVGN6LpH6L9/ZEbB7GfkeYYvJiYyJ2NGe/448nUkTno1d2O+/+yU/HYYLkB0VS21Qor31qb3qOvkfYP0ZyZyfnZEV8bUSgMXkxM9M6mvLqOwYuOJl6TEVFugJlylVj51pqM2s08kv5xXW5vZDu7o84V/H2KOVWkJ642MjHRO5v39n7N1Uc6imQ1ByvbkhbMXKunvLoOl9q8Qb/HnCrSG4MXE7sutzd691SuN9LY3MpCYzqTcwOy0/wDymBLPlnZlrQiOvq6/XiDof1JDs4vtATPBUtPtXMxAemK00Ymlpxkw4/GXC1URyHYm5wZ8i3iiWhugJq7ZU7liPP6JOysOYfdDTZk1JwLu5ImXoiOvi7/+DjW7/kHnpmRjzvy9F3VEy44lzm6JTHXhXTF4MXkvp+fLRS8dH2TM1O+RTwRyQ1gZVvt+ffnZKw5tish+rNSLZbO5CnJl+8dres1iSwkqHO5GZyTrjhtZHKRlP5mvkVssbKtthK5P4fLt+pKDm6ef/8I9JxBYnBOZsDgxeTUJosy3yL2It1rhgKxP4fOtwqmfUrSjRMu/abTGJyTGTB4sQA1yaJmXp2QKCLda4YCsT+3mzoyB5/+x+2Ye9tgoeNdyjUVI8bgnMyAOS8WIZosyiFdc0ikyrZ6JoazP1+WnGTDjUP6YPnHJxSPdYptih7xdXDbCYo1Bi8WIpIsyiFd80iEyrZ6J4azP/sT20zRgcHOZl2vI5GCczInBi9xhjsJm0s8V7Y1YtNQLXY+jqfgUWTUY+G04fCe2q37tSRCcE7mxeDFpCJ90+WQLhnByLL1Sv356R+MCPq3Eq/lApRGPe7Iy0TpKWOuJZ6DczI3XYOXbdu24be//S12796N2tpabNiwAT/84Q/DPmfr1q0oLi7GoUOHcNVVV2H+/PmYM2eOnpepKS3u9KJ90+WQLunNyEJ84frzXaNz8Ov/rUady335cacDM8dchde21eg6KhRL4UY9PB4ds3WJTELX4KW5uRmjR4/Gww8/jHvuuUfx+JqaGkyfPh2zZ89GSUkJtm/fjsceewx9+vQRen6saXGnp9VQPId0SU9GJ9LK/bny+Fls+WQnim6+Hq5LPjy2dk/AsXUuN17dVhP0dbQcFYo1jnpQItM1eJk2bRqmTZsmfPzKlSsxYMAALFu2DAAwYsQI7Nq1Cy+99JLpg5cPDtVj3rp9UQUdWg/F882N9BKLRNrkJBuuz+2NxsMSJgy8EhNfrIjodbg9A5H1mSrnpbKyEkVFRX6PTZkyBatWrYLH44HdHrj+z+12w+2+PGTscrkAAB6Px5DhU4/HA58E/Kb0iELQcQi3Ds0IG3TsFByKrzx+FtdbKOFW/j2YbTjb65Ow69R5nL3oRt9eDkwYeGVM78TN2k7BjO3XC9lOB+pd7rCrXsb266XpzyO/1mfHzobcFFBU7YVmeDxOLS7LVPTuR2b7u4mUlf7eYsXoNlJzHlMFL3V1dcjKyvJ7LCsrC21tbWhoaEBOTuDIxdKlS7FkyZKAx7ds2YLU1FTdrrWzEy4b6jvNuXclV71c/k4ZhqaFrgS6u8EGIFnxfFs+2YnGw9arKFpeXh7rS+iwr9GGd08m4ULr5Tfd9BQJdw/yYXRGbNvWTO0UzvRsG95wyXUuO394SZAATMtqwQdl7+ty7nUf74HI30o4fz+0F6X/qNLmgkxIj35k5r+bSFnl7y2WjGqjlpYW4WNNFbwAgM3mH8FLkhT0cdmCBQtQXFzc8bXL5UL//v1RVFQEp1P/uyqPx4Pdb38odOw1BWMwfVToqaOMmnNYc2yX4usU3Xy95UZeysvLMXny5KCjZ0b74FA93qwMnOJrarXhzS+T8fK9ozGlICvoc/VktnZSMh3AuEP1eK70iF/CbE5adyycNlyXNpTbaFBuLnDmdESvIY8Kzf3xJEuOGCjRqx+Z9e8mUlb7e4sFo9tInjkRYargJTs7G3V1dX6PnT17Ft26dUNGRvC5aYfDAYfDEfC43W43rEOKVrPMSe8Z9poKh/QVqmlROKSvJd90jfydhOL1SXj+/aNhp/ief/8opo26OmZtbIZ2EnXnmH6YNupqQxPD9zXasOEfkW3GeLlcQAG6O1K0uygT0rIfWeHvJlJW+nuLFaPaSM05TBW8FBYWYtOmTX6PbdmyBRMmTDB15xrslATm/5ULw7FGi/5El/iu3l6DzF6OhF+lJbL038jE8A8O1eONL5MAKM+N90xJRq/udtS5WC4gWkYujScSoWvw8s033+D48eMdX9fU1GDv3r3o3bs3BgwYgAULFuDMmTNYs2YNAGDOnDlYvnw5iouLMXv2bFRWVmLVqlV4++239bzMqCXZgEXTh2Peun1RBx2halqk9bDj4RsHYXJ+tqbXnmhEl+4+u/lwx//HQ2GzSJityJvXJ+G50iPCx//u/4xmuQCNcI8pMhtdd5XetWsXxo4di7FjxwIAiouLMXbsWPznf/4nAKC2thanT1+et87NzUVpaSkqKiowZswYPPvss/jTn/5k+mXSADClIEt452cl8g6yv/r+MKT3aB9xuvCtB3/48BhuevEjlB2MbMicgJMN6vd8kZe7J1K7y/WGut5tx7ItPq85911uTfjgo3dPO1Z+9zcnjwrNHHM1CgeHX+1HoXGPKTIbXUdebr311o6E22BWr14d8Ngtt9yCPXsCC09ZgZaF4cqr67Dswy/jtkJoLHh9Et7+XH2SZzwVNhNhVOl/tUTv6p++s4B/GxrjnmlkNrqOvCQiLe70lD48gPYPD6/PmksTY+Xynbt6nef0452a/AYjid7VZzt59681OR8PCBz3Yj4exQKDFxMy64eH1WkxH58Ic/pmyW/w+iRUnmjEe3vPoPJEI8YPvBLZTgcQNKxv/xDNSbC7/65tpOcNjZyPp8XUOFG0TLXaiNqZ5cMj3mgxH58Ic/pmyG8IlSx85/ey8ZftJ7kaD6HbaOG0PE1eP9hKs65T45lXOAAJaGh2o/JEIxOiyTAMXkzIDB8e8Uhp3j6cRJrTj3V+Q7jNSVdtP4Xbc3yobk7130k6wVaEhWujeev24eFhNkyP8vXDrTQrHJyBsoO1eOJv+0yzGo0SC6eNTEj+8Ah1/5KIw+NaEJm3D/e9RLmrj2V+g0i+185/JsHn8/l/L8zCgHgj0kbvnkwKOoUkMs0kstLMjKvRKLEweDEhLT88jJwTt4Jw8/YrHxiHlZzTBxC7/AaRfK/mNhvOfuNfpK7e5U6YD02RNrrQasOuU+f9Hi87WIubXvwI972+A79Ytxf3vb4joPSCSGC0eOMhLN7IBQUUW5w2MqlQxerUDI+brciYWSgtaWdhs3ZaLv0XFWkeVyItZxfPibs8rVZ2sBZzSgJLUHQtvSASGCmt2GO1XTICgxcTi+bDI9ycOOvEhC9pb2S5e7Mzui2iyeNKlA9N8Zy49j3fvD4JT757IOgxXYM+LRcBcEEB6YnBi8lF8uFh1iJjRhPZl4fMJZqkalm8f2iKJFSnpUiYMPBKAMDyj47jQkvovaA6B31aLgLI7Bm4YS6RVpjzEodYJ0Zsfp/MJ1y+l6h4X4UnkhN39yAfkpNs8PokvLrthNDrnr14SWixQLbTgWxn6GNk/+/f9vHvjXTD4CUOJXqdmHArIeaU7MGvNx1i8rKJhUoWzknrjvQe3cAideETql++dzRGZ7S30Y6/N6Kl1Sv0mn17dRcKjBbfVYDFdykHmPUurjwi/XDaKA4lcp0YkdUSb2w/iTe2n2TyskmIFEOTH39//xnMXbeXReoQOifO521D6an2YypPNAq9Vk9HckfQJ7pYYMUD47B446GQCbyJNEVNxmPwEodiXWQslpSmzDpj8nLsiRRD62xKQRYeGeZDaV1iF6mTBcuJ8/kNtIiNLt4yNNMvuBBZLDB1ZA56dbfj/r/sDPm6iZJETcZj8BKH5KHfR0v2JNwdqpqpMN4ZxlakK+JGZ0iYf/8kVP3jIpOxFRRek4nlHyvnvNx//aCAx0QWCzR8I7bRabxOUVPsMOclTiXqJmpqp8ISIXnZjKLdOV2L3dsTwcTBGUhPtYc9Jj3VjokRjook8hQ1xRZHXmLAqCW8SkO/8biUONKltrwzNJaaFXGcbohccpINL9z9vaAF6mQv3P29iP/uE3mKmmKLwYvBjK56G2roN16r74abMguHd4bGSvQVcUaaOjIHK4Mk12Y7HVh8V0FUf++JPEVNscXgxUBmqXprluvQS6jVEsHwzjA2ON1gLD23etBiKxMitRi8GMQsVW/Nch166/xm/WF1HVZtPxlwDO8MY4fTDcbTc6uHWOyDRYmNCbsGMUvVW7NchxHkN+unZxRg5QPjkJNgyctmpuXO6WQOTKImI3HkRSNen4RdJxpD3nWIzt3XNX2LyjCvE61EzTXgnaH5cLqBiCLF4EVQuJU5+xptWPq7bX7JcF2TX0Xn7p9+7xC+cbeFfB2la1FysqFF6Lh4zDXgbtHmw6DSvOJxNSLFDwYvAsKtzGlr8+KNL5MA+Bdr6pr8el1ub2Q7HSFLacs6By7BXieaVUJen4S3Pz+t+PNmOx3MNSDDMKg0n3hdjUjxgzkvCpQ2+Vv0XnXQ53UttJWcZMN91w1Qff7Or1O6P/S1iGyA9nnNOdS5lKeD7rtuAO+wiBJUuPc8brRIZsHgJQyRKqAXvvUg1N6qXZNfB2X2jOg65NdZ9N7BiCuSAuJ5LJFeJxFZW7SVj4mMwuAlDDWb/IUjBw3R5pGca24N+T2RVUKsrUFE4STSakSyNgYvYWi14uZY/TeoPNGI8QOvRE5a9xDjNNoId81ybY1Q57ehfV6b+S5EiSlRVyOS9TB4CUN8BCL8EOryj4/jvtd34Jbffoy7Rrcnu6kJYGwAevcMv7maLNw1s7YGEYXD0VmyCgYvYYiMVKSnduv4fyV1TZfw2rYa/Nuk3IDdnuWdX0MFFc/NHKnJqIlcWyPLyYJtROSPo7NkFQxewhAZqXjurgI8MsyHLKdD8fXk8ZmN+2qx9f+7DW/Pnog/3jsGb8+eiN2LJmPlA+MCgho5qJg+6iqNR038R4skiQl4RImOo7NkFazzokCpCugdeZnwnpIw//5JqPrHRWw//k8s//hEyNeTE952nzofUNtCqWCXFhVJQ23KWO9yx8WmjEQUHVY+JiswJHh55ZVX8Nvf/ha1tbUoKCjAsmXLcPPNNwc9tqKiArfddlvA44cPH8bw4cP1vtSgwgUVHo8HwOVCW9EmvCkV7IqmImmibMpI8aVzpdeM1G7gKl39sfIxmZ3uwcs777yDX/7yl3jllVdw44034tVXX8W0adNQXV2NAQNCF207evQonE5nx9d9+vTR+1LDEq0CakTCW6QVSdUsg2TFUzKDYJVe01OSYR9UjzvH9IvhlRFRLOkevPz+97/HT3/6U/zsZz8DACxbtgwffPABVqxYgaVLl4Z8Xt++fZGenq735WlOTnira7oUcg1Skg043xx+mwA9cBkkmYXIvjmhpjgvtALz1u1Dt27JCT2FoefeQ9wegMxO1+CltbUVu3fvxpNPPun3eFFRET777LOwzx07diwuXbqE/Px8LFq0KOhUEgC43W643ZcDAZfLBQDweDwdUzp6ks/R+VwLp+Vh3rp9IZ/jk4DH11bhZZ8PUwqydL9GWUaq2K+7vqkFl9ytmr0RBmsjCpQo7fTBoXo8V3rEb5+vbKcDi6YP7/h78PokLN54KMQNQHu/XLLpEG4dmpGQUxnh2vD2Ye0rgSLtRx8cqse8dfsC2l7eHuDle0cb+r6ll0T5e4uG0W2k5jw2ScdlJl9//TWuvvpqbN++HTfccEPH47/5zW/w17/+FUePHg14ztGjR7Ft2zaMHz8ebrcbb731FlauXImKigpMmjQp4PjFixdjyZIlAY+vXbsWqamp2v5AKuxttGH1l0mQQi46lJCeAjwzzguj3nt9ErBkTzIutAJKi7vTUyTcPciH0RlMMCDt7Gu0fbeRKeDfB9v72SPD2vvcsSYbllcnK77e3HwvhqYlVh8VbcNIKL9HGP++RYmjpaUFs2bNQlNTk1/aSDCGJOzabP69XJKkgMdkeXl5yMvL6/i6sLAQX331FV566aWgwcuCBQtQXFzc8bXL5UL//v1RVFSk+MNrwePxoLy8HJMnT4bdfrmQXEbNObz55a4wz7ThQivQJ38irjewZoJ9UH3HqFC4t7emVhve/DJZk7usUG1E/uK9nbw+CUt/tw1dd2BvZ4MNwPv1qZh//yR4D9YB1QcUX/OagjGYPipxpjFE2rC0LhXf692MKUXq+9HOmnO4sMN871t6iPe/Ny0Y3UbyzIkIXYOXzMxMJCcno66uzu/xs2fPIitL/ANx4sSJKCkpCfo9h8MBhyOwxordbje0Q3Y9X2NLm9DzGlvaDL3OO8f0Q7duyQHz2V3Jq4+ef/8opo26WpOheaN/J1YVr+2060Sj3zRHV+0J425U/eMictLFNgfNSe8Zl20Vikgb1rncOOGyRdSPzPq+pad4/XvTklFtpOYcuhapS0lJwfjx41FeXu73eHl5ud80kpKqqirk5Fjr7srMZbanjszBp/9xO57+wYiwx3ETNtKSmoRxVnoNTrQNXZ1SB7w+CZUnGvHe3jOoPNEYdkdoM79vEXWm+7RRcXExHnzwQUyYMAGFhYV47bXXcPr0acyZMwdA+7TPmTNnsGbNGgDtq5EGDRqEgoICtLa2oqSkBOvXr8f69ev1vlRNKa06sqG96JP85qvnyoFgkpNsyOylXBUY4Ooj0oaaD0a50uujJXtgQ0A9aAC2hKz0KtqGzu9uYNWuGlL7vkUUK7oHLz/+8Y/R2NiIX//616itrcXIkSNRWlqKgQMHAgBqa2tx+vTpjuNbW1vxxBNP4MyZM+jRowcKCgqwefNmTJ8+Xe9L1VS4N9+uZbZjtSyRd1lkJPmDMdx0JQCcb24FELrSa3oK8NzdoxNyya5YcOHAYGez4qqhYNW01bxvEcWSIXsbPfbYYzh58iTcbjd2797tl3i7evVqVFRUdHw9f/58HD9+HN9++y3OnTuHTz75xHKBi0x+8w21X9HUkTkdtSy6vqHLbzBlB2t1uz69h+bl4epN+2txrMkWdria4l9ykk1xqhIAnt1c3dFX5ClOeR+wkkcm4Jlx3rhYqhsJkb2HFk5rr0T+XOmRkNW0gfZq2sH+JkXet4hijXsb6Sxcme1Yl+vX8y4rcDQpGf/9u21YfFcB3/wSSNfp0LTUFMXndK3y3LmitMfjQelhXS/Z9ET2W/vjfptAcnToatrcHoDMjsGLAUKV8zdDuX49NmHj5o8EhCjt30NsNQHzrMJT2m/NJVjrK1w7R7oNCZERGLzEUKzL9ct3xe42H176l9GADWj4xh3VXVasR5PCXRfvIo0TsrT/t2KfqsyzUhYuuHAKrjhlO5NVMXiJoVgmzIZLEo7kbksODrYf/2fMR5M6X8/Zi5dwsqEFb39+GnUu7tNihHABrBKuZtHGYKeEbKcD9S43Vw1RXGLwEkOxWpYY6q443CoEpddTKnrXlZ7TAiLXE+nPSsqUpkND4WoW7STZgEXTh2Peun1cNURxyZDVRhScyMoBrd9glKZ1gNCrEIIJtVpKiV7D1aLXE8nPSmJEA9Ou+S9czaKtKQVZXDVEcYsjLzGmR8JsOFomCUcyPaDncLXa6zFqCivRiAamf541DklJNuYh6YirhiheMXgxASPfYLRMElY7PaD3cHWk0xVc2aItpelQoD3naOLgDH6IGoCrhigeMXgxCaPeYLRMElb7oZ+d5sAzM/Sr8xJpEMIVF9rqXD8olLtG5zBwIaKIMeclwZxvDl24SiZaVVf0Q/+xW3IxN9+Lj4snmWK7A1mibu5nhKkjc/Bvk3JDfv+1bTW6Vo8movjG4CWBeH0Snt2sXJ706R+ITeuIbi/w89uHYGiapPudttL1dL02gCsu9OL1Sdi4L3xwwmRpIooUg5cEIpoTcmVP5RLuQGxWS0V6PV1xxYW+1CSGExGpxZyXBKJHRV+R1VIej2Ctcg2EvB6nA/ddNwCDMntyxYUBYl09mojiG4OXBKJXRV+zLcc02/UkolhWjyai+MfgJYHoWdHXbMsxzXY9iSZW1aOJKDEw5yWBmC1HheIX+xoR6YnBS4x4fRIqTzTivb1nUHmi0bBVF3JOCEuGk97Y14hIL5w2ioFwOzob8YbOnBAyCvsaEemBwYvBtN7ROVLMCSGjsK8RkdY4bWQgrXd0JiIiSkQMXgwUi8JdscqtISIi0gunjQxkdOGuWObWeH1SR55DRmo3MGZKTJ37AfNdiEgrDF4MZGThrljm1pTu/xqL3juIc82XK+umpyTDPqged47pp8s5yXxinZhORPGL00YGEt3IMNrCXbHMrVlaWo3H1lb5BS4AcKEVmLduX8idhDm9FV/k4LnrNKkcPHNHaSKKBkdeDCQX7nq0ZA9sgF9woWXhLjW5NVquAindX4tXt9WE+K4NEtqDpsn52X4/I+/Q44tS8GxD8H5ARCSKIy8GM6JwlxG5NV1HSlrbfFj03kHF53VNSOYdenzx+iSs3l7DHaWJSFcceYkBvQt36Z1bE2ykpHdPe8BUUShy0MQ79PgSrF+Ewx2liShSDF5iRM/CXVpsihdqlUioRGDRwAW4HDTFanqLtBeqX4TDHaWJKFIMXgxi9JLRe6/tjz98eCzgcZHcmlA5KE//IB/Pbg4+UiKqd097R9Bk9NJx0ke4EbRguKM0EUWLwYsBjExIVRq6z1Y4b7gl1o+t3RPFlbVPBD03cyQAoPJEI47VXxR6Ju/QzU1pBK0z7ihNRFpg8KIzI+utKA3d/+r7wzD39iEhPzREllhH42c3DkRSkg03vfiR0Icd79CtQc3ImFLwTEQkwpDVRq+88gpyc3PRvXt3jB8/Hp988knY47du3Yrx48eje/fuuOaaa7By5UojLlNzRtZbURq6twFY98XpsK+h5g5ajd497Xh4mA9j+qcHXVkUDO/QrUN0ZOzpH4zAp/9xOwMXIoqa7sHLO++8g1/+8pdYuHAhqqqqcPPNN2PatGk4fTr4B2lNTQ2mT5+Om2++GVVVVXjqqafw85//HOvXr9f7UjVn5F5GWpxLj9ySp38wAp/NvxWjekt4rvSI8AiOlkvHSV+ixRd/cmMuA1Ei0oTuwcvvf/97/PSnP8XPfvYzjBgxAsuWLUP//v2xYsWKoMevXLkSAwYMwLJlyzBixAj87Gc/wyOPPIKXXnpJ70vVnJEJqVqcS4/cksxeDiQn2XDCZUOdy614/NzbhuDt2RN5h24hcvFFAAEBDEfQiEgPuua8tLa2Yvfu3XjyySf9Hi8qKsJnn30W9DmVlZUoKirye2zKlClYtWoVPB4P7Ha73/fcbjfc7ssfii6XCwDg8Xjg8Ygv342UfI5g58pIFWvejNRuUV+rFuca268Xsp0O1LvcYZZYO7BgSh6e+d/DON+ifM3y+VyCP941mT0wtl8vVB4/i7MX3ejby4EJA69MiA++cH3J7O7Iy8TL947Gc6VH/ILU7DQHFk4bjjvyMjX5uazcRkZhG4lhOykzuo3UnEfX4KWhoQFerxdZWVl+j2dlZaGuri7oc+rq6oIe39bWhoaGBuTk+N+NL126FEuWLAl4nS1btiA1NTXKn0BceXl5wGM+qX1DwgutQOA9KQBISE8B/lm9A6WHozu/Vueanm3DGy55QK7z60iQAIzp1YJdVVW4dwCw5ngSmtuUz1d+GHDaxYKPv/7ffXh6w340t10+Pj1Fwt2DfBidkRj7HQXrS1bxH/nACZcNLg/gtAODnc3wntqN0lPansfKbWQUtpEYtpMyo9qopaVF+FhDVhvZbP4fXJIkBTymdHywxwFgwYIFKC4u7vja5XKhf//+KCoqgtPpjOayhXg8HpSXl2Py5MkBo0IAYB9Uj3nr9gEItpeRDc/dPRpTCrICnhcJLc41HcC4Q/UBd9DpqXZAsuH9f3TaKbpHN6CtLcQ+TZfP5/F44NtSjiynA2dDjOrI9p0LnMlsarXhzS+T8fK97a/n9UnYdep83I3MKPUlYhuJYBuJYTspM7qN5JkTEboGL5mZmUhOTg4YZTl79mzA6IosOzs76PHdunVDRkZglVWHwwGHwxHwuN1uN7RDhjrfnWP6oVu35IDaK3osGdXqXHeO6Ydpo67uKKp3sqEFyz78MiDoaPq2DQCQlmrHhU5TSMHOl2QDnp4+HPPW7QsIdpTIWwU8//5RJCUl49nN8b2Jo9F914rYRsrYRmLYTsqMaiM159A1eElJScH48eNRXl6OH/3oRx2Pl5eXY+bMmUGfU1hYiE2bNvk9tmXLFkyYMMGyHUzvvYz0OJe8fYHXJ+GmFz8Ku/9Q925J+K+fXY+Gb9xhzzelIAsrHhinav+bzueqDVEoT4+aOUREZF66TxsVFxfjwQcfxIQJE1BYWIjXXnsNp0+fxpw5cwC0T/ucOXMGa9asAQDMmTMHy5cvR3FxMWbPno3KykqsWrUKb7/9tt6Xqis99zLS81wiS7DrXG4k2WyYOeZqxdfrHFxtP/5PLP/4RNTXyE0ciYgSi+7By49//GM0Njbi17/+NWprazFy5EiUlpZi4MCBAIDa2lq/mi+5ubkoLS3Fr371K/z5z3/GVVddhT/96U+455579L5UCkKP5d5ycKVlXRlu4khElDgMSdh97LHH8NhjjwX93urVqwMeu+WWW7BnTzT76JBWRGu/hDvO65NwrMmGTftrkZPes2NaSY+6MmcvXjJ8E0wiIjIW9zaisOTqqXVNl8LUfgm9/1DZwVos3ngIda5koPoAgMsJtpPzs8O+diTKq+vxwvtH4jqhl4go0RmytxFZVzTVU+WNIrtW1pUTbMur60K+dmc5ad3xyqyxYUvQy/53f21Ajo58vrKDtQrPJiIiK2DwQoqmjszBigfGITvNf5on3P5DoptSTs7PDvraGT1T8MiNgzq2Cpg+6io8MyM/ohEarTfBJCKi2OK0EQlRuwRbzUaRoq89dWQOfvX9ofjDh8dUXz8TeomI4geDFxKmZgm22lVKoq89KLOn0OtGe11ERGRenDYiXWixSkmL47V+PhERxR6DF9KFvEopVIKtDe2JuKFWKUX6uqFEej4iIjIfBi+kKa9PQuWJRvzv/q9x77UDAKhfpRROuNVPoURzPiIiMh/mvJBmyg7WBuxb1L4bNXDh2/AbN6ohr34S3SNJj00wiYgodhi8kCbkmi4BO0+3eCABmNbPi8mFY/wq7EZDXqG0/KNjYVcf/er7wzD39iEccSEiiiOcNqKoKdV0sQGoPJuE6SOzUTg4Q9NAYt0XX4X8ng3Aui9Oh/w+ERFZE4MXippITZcLrTbsOnXe8PPKtV2IiCh+MHihqInXdHErH6TLeVnbhYgonjB4oaiJ13RxxOi8rO1CRBRPGLxQ1ERquqSnSJgw8ErDz8vaLkRE8YfBC0VNZOfpuwf5NF/xE82O10REZF0MXkgT4Xaefvne0Ridoc9uzpHseE1ERNbGOi+kmVC7Q/u8bSg9Zfx5OeJCRBSfGLyQpoLtDu3zxua8REQUnzhtRERERJbCkRcyBa9P4rQPEREJYfBCMRdsQ8ccbqZIREQhcNqIYkre0LFrmf+6pkt4tGQPyg7WxujKiIjIrBi8UMwobegIAEs2VcPr02eZNRERWRODF4oZbqxIRESRYPBCMfNhdZ3QcdxYkYiIOmPwQjFRdrAWq7afFDqWGysSEVFnXG1EhpNzXZTY0F7mnxsrEhFRZxx5IcMp5brIJHBjRSIiCsTghQwnmsPyyI2DWOeFiIgC6Bq8nD9/Hg8++CDS0tKQlpaGBx98EBcuXAj7nJ/85Cew2Wx+/yZOnKjnZZLBRHNYJudn63wlRERkRboGL7NmzcLevXtRVlaGsrIy7N27Fw8++KDi86ZOnYra2tqOf6WlpXpeJhnsutzeyEnrjlCTQTa0V9hlrgsREQWjW8Lu4cOHUVZWhh07duD6668HALz++usoLCzE0aNHkZeXF/K5DocD2dm8645XyUk2PDMjH4+W7IEN8CtSJwc0zHUhIqJQdAteKisrkZaW1hG4AMDEiRORlpaGzz77LGzwUlFRgb59+yI9PR233HILnn/+efTt2zfosW63G263u+Nrl8sFAPB4PPB4PBr9NKHJ5zDiXFYVrI3uyMvEy/eOxnOlR1Dnuvz7y05zYOG04bgjLzPh2pR9SRnbSBnbSAzbSZnRbaTmPDZJknSpvf6b3/wGq1evxpdffun3+LBhw/Dwww9jwYIFQZ/3zjvv4IorrsDAgQNRU1ODp59+Gm1tbdi9ezccDkfA8YsXL8aSJUsCHl+7di1SU1O1+WFINz4JOOGyweUBnHZgsFMCB1yIiBJPS0sLZs2ahaamJjidzrDHqh55CRUsdPbFF18AAGy2wE8hSZKCPi778Y9/3PH/I0eOxIQJEzBw4EBs3rwZd999d8DxCxYsQHFxccfXLpcL/fv3R1FRkeIPrwWPx4Py8nJMnjwZdrtd9/NZEdtIDNtJGdtIGdtIDNtJmdFtJM+ciFAdvMydOxf33ntv2GMGDRqE/fv3o76+PuB7//znP5GVlSV8vpycHAwcOBDHjh0L+n2HwxF0RMZutxvaIY0+nxWxjcSwnZSxjZSxjcSwnZQZ1UZqzqE6eMnMzERmZqbicYWFhWhqasLnn3+O6667DgCwc+dONDU14YYbbhA+X2NjI7766ivk5LDeBxEREem4VHrEiBGYOnUqZs+ejR07dmDHjh2YPXs27rzzTr9k3eHDh2PDhg0AgG+++QZPPPEEKisrcfLkSVRUVGDGjBnIzMzEj370I70ulYiIiCxE1zov//Vf/4Xvfe97KCoqQlFREUaNGoW33nrL75ijR4+iqakJAJCcnIwDBw5g5syZGDZsGB566CEMGzYMlZWV6NWrl56XSkRERBah68aMvXv3RklJSdhjOi926tGjBz744AM9L4mIiIgsjnsbERERkaUweCEiIiJLYfBCRERElsLghYiIiCxF14TdWJATgNVU6ouGx+NBS0sLXC4XCx2FwDYSw3ZSxjZSxjYSw3ZSZnQbyZ/bIrsWxV3wcvHiRQBA//79Y3wlREREpNbFixeRlpYW9hjdNmaMFZ/Ph6+//hq9evUKu4eSVuS9lL766itD9lKyIraRGLaTMraRMraRGLaTMqPbSJIkXLx4EVdddRWSksJntcTdyEtSUhL69etn+HmdTif/ABSwjcSwnZSxjZSxjcSwnZQZ2UZKIy4yJuwSERGRpTB4ISIiIkth8BIlh8OBZ555Bg6HI9aXYlpsIzFsJ2VsI2VsIzFsJ2VmbqO4S9glIiKi+MaRFyIiIrIUBi9ERERkKQxeiIiIyFIYvBAREZGlMHgR8MorryA3Nxfdu3fH+PHj8cknn4Q9fuvWrRg/fjy6d++Oa665BitXrjToSmNHTRtVVFTAZrMF/Dty5IiBV2ysbdu2YcaMGbjqqqtgs9nwP//zP4rPScR+pLadEq0vLV26FNdeey169eqFvn374oc//CGOHj2q+LxE60uRtFOi9aUVK1Zg1KhRHQXoCgsL8f7774d9jpn6EYMXBe+88w5++ctfYuHChaiqqsLNN9+MadOm4fTp00GPr6mpwfTp03HzzTejqqoKTz31FH7+859j/fr1Bl+5cdS2kezo0aOora3t+Dd06FCDrth4zc3NGD16NJYvXy50fCL2I0B9O8kSpS9t3boVjz/+OHbs2IHy8nK0tbWhqKgIzc3NIZ+TiH0pknaSJUpf6tevH1544QXs2rULu3btwu23346ZM2fi0KFDQY83XT+SKKzrrrtOmjNnjt9jw4cPl5588smgx8+fP18aPny432P//u//Lk2cOFG3a4w1tW308ccfSwCk8+fPG3B15gNA2rBhQ9hjErEfdSXSTonel86ePSsBkLZu3RryGPYlsXZK9L4kSZJ05ZVXSn/5y1+Cfs9s/YgjL2G0trZi9+7dKCoq8nu8qKgIn332WdDnVFZWBhw/ZcoU7Nq1Cx6PR7drjZVI2kg2duxY5OTk4I477sDHH3+s52VaTqL1o2glal9qamoCAPTu3TvkMexLYu0kS8S+5PV6sW7dOjQ3N6OwsDDoMWbrRwxewmhoaIDX60VWVpbf41lZWairqwv6nLq6uqDHt7W1oaGhQbdrjZVI2ignJwevvfYa1q9fj3fffRd5eXm44447sG3bNiMu2RISrR9FKpH7kiRJKC4uxk033YSRI0eGPC7R+5JoOyViXzpw4ACuuOIKOBwOzJkzBxs2bEB+fn7QY83Wj+JuV2k92Gw2v68lSQp4TOn4YI/HEzVtlJeXh7y8vI6vCwsL8dVXX+Gll17CpEmTdL1OK0nEfqRWIveluXPnYv/+/fj0008Vj03kviTaTonYl/Ly8rB3715cuHAB69evx0MPPYStW7eGDGDM1I848hJGZmYmkpOTA0YQzp49GxCByrKzs4Me361bN2RkZOh2rbESSRsFM3HiRBw7dkzry7OsROtHWkqEvjRv3jxs3LgRH3/8Mfr16xf22ETuS2raKZh470spKSkYMmQIJkyYgKVLl2L06NH44x//GPRYs/UjBi9hpKSkYPz48SgvL/d7vLy8HDfccEPQ5xQWFgYcv2XLFkyYMAF2u123a42VSNoomKqqKuTk5Gh9eZaVaP1IS/HclyRJwty5c/Huu+/io48+Qm5uruJzErEvRdJOwcRzXwpGkiS43e6g3zNdP4pJmrCFrFu3TrLb7dKqVauk6upq6Ze//KXUs2dP6eTJk5IkSdKTTz4pPfjggx3H//3vf5dSU1OlX/3qV1J1dbW0atUqyW63S//93/8dqx9Bd2rb6A9/+IO0YcMG6csvv5QOHjwoPfnkkxIAaf369bH6EXR38eJFqaqqSqqqqpIASL///e+lqqoq6dSpU5IksR/J1LZTovWlRx99VEpLS5MqKiqk2trajn8tLS0dx7AvRdZOidaXFixYIG3btk2qqamR9u/fLz311FNSUlKStGXLFkmSzN+PGLwI+POf/ywNHDhQSklJkcaNG+e33O6hhx6SbrnlFr/jKyoqpLFjx0opKSnSoEGDpBUrVhh8xcZT00YvvviiNHjwYKl79+7SlVdeKd10003S5s2bY3DVxpGXYXb999BDD0mSxH4kU9tOidaXgrUNAOnNN9/sOIZ9KbJ2SrS+9Mgjj3S8Z/fp00e64447OgIXSTJ/P7JJ0ncZN0REREQWwJwXIiIishQGL0RERGQpDF6IiIjIUhi8EBERkaUweCEiIiJLYfBCRERElsLghYiIiCyFwQsRERFZCoMXIiIishQGL0RERGQpDF6IiIjIUhi8EBERkaX8/zM4a/D1twTCAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "m = 200\n", "# Data : \n", "L = 3.0\n", "xi = np.linspace(0.0, L, m) + 0.05*np.random.randn(m)\n", "yi = 0.8 + 0.5*xi-np.sin(3*np.pi*np.cos(xi)) + 0.3*np.random.randn(m)\n", "plt.plot(xi, yi, 'o'); plt.grid()" ] }, { "cell_type": "markdown", "id": "d2de79e1-2d70-42c5-a0df-51879d797d5d", "metadata": {}, "source": [ "On réutilise la fonction $\\phi$ de l'Exercice 1. Soit $n\\leq m-1$ et $\\mu$ un scalaire positif ou nul. On considère la fonction de régression affine par morceau\n", "\n", "$$\n", "f(x) = \\sum_{j=0}^n u_j\\, \\phi(x - \\frac{j}{n}L)\n", "$$\n", "\n", "et le problème de minimisation régularisé\n", "\n", "$$\n", "\\min_{\\mathbf{u}}\\quad \\frac{1}{2} \\sum_{i=1}^m \\Big[ f(x_i)-y_i \\Big]^2\n", "+ \\frac{1}{2} \\mu \\sum_{i=1}^m u_i^2.\n", "$$\n", "\n", "Réécrire le problème de miminimisation sous la forme\n", "\n", "$$\n", "\\min_{\\mathbf{u}}\\quad \\frac{1}{2} \\| A\\mathbf{u} - \\mathbf{y}\\|_2^2\n", "+\\frac{1}{2}\\mu\\, \\|\\mathbf{u}\\|_2^2.\n", "$$\n", "\n", "Pour la mise en oeuvre, on considèrera $\\mu = 10^{-15}$ et $n=80$.\n", "Codez la matrice $A$ et résoudre les équations normales associées au problème de minimisation.\n", "Afficher la valeur de $\\|\\mathbf{u}\\|_2$." ] }, { "cell_type": "code", "execution_count": null, "id": "e129c9fc-4425-4d0b-b790-4a4efc320cf3", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "d8783e39-f30d-40a9-9654-c12c3668f954", "metadata": {}, "source": [ "Codez une fonction\n", "\n", "```\n", "def f(x, u):\n", " # ...\n", "```\n", "\n", "qui implémente la fonction de régression $f$ avec le vecteur de coefficients $\\mathbf{u}$. " ] }, { "cell_type": "code", "execution_count": null, "id": "4dbd14b5-9d6b-4c6a-be16-d724a67d018f", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "6c11b67a-f0bd-4876-8a04-f40cbd3b23d3", "metadata": {}, "source": [ "Sur le même graphique, tracez le nuage de points ainsi que la fonction de régression.\n", "Affichez aussi la valeur du conditionnement en norme 2 de la matrice $A^T A + \\mu I$." ] }, { "cell_type": "code", "execution_count": null, "id": "b9aec732-27f6-42b8-a016-1527fb24416c", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 418, "id": "db9a7765-8143-4524-888a-e93cbccd24e9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "84072032.97441417\n" ] } ], "source": [ "print(np.linalg.cond(A.T@A+mu*Id))" ] }, { "cell_type": "markdown", "id": "5fbcf543-73f1-41ae-b301-15f7c2197e8f", "metadata": {}, "source": [ "Recommencez en testant différentes valeurs de coefficient de régularisation $\\mu$ :\n", "$\\mu = 10^{-4}$, $10^{-3}$, $10^{-2}$, $5.10^{-2}$, $10^{-1}$, $1$ et 10. Affichez à nouveau le nuage de points et la fonction de régression résultante. Observez l'influence du coefficient de régularisation sur la régression." ] }, { "cell_type": "code", "execution_count": null, "id": "47624213-36c1-427a-9bab-32feff2fcd34", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "dd3e5800-d44d-41fa-94a3-09a11fa6946b", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "d9b440a3-22b5-422d-b773-fa920821e940", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a7070eff-ee1f-4504-8274-9718413617f7", "metadata": {}, "source": [ "#### Exercice 4 (optionnel) - Identification de paramètres d'un système masse-ressort\n", "\n", "On considère ici un système masse-ressort linéaire avec amortissement, $m_0$ est la masse, $k$ est le coefficient de raideur du ressort et $a$ est le coefficient d'amortissement. Le mouvement du point matériel est gouverné par l'équation différentielle du second ordre:\n", "\n", "$$\n", "m_0\\, \\ddot{x}(t) + a\\, \\dot x(t) + k\\, x(t) = 0\n", "$$\n", "\n", "où $x(t)$ représente la position de la masse au temps $t$.\n", "L'équation caractéristique correspondante est:\n", "\n", "$$\n", "r^2 + \\frac{a}{m} r + \\frac{k}{m} = 0\n", "$$\n", "\n", "avec le discriminant $\\Delta$ fonction de $m_0$, $k$ et $a$ :\n", "\n", "$$\n", "\\Delta = \\left(\\frac{a}{m_0}\\right)^2 - 4 \\frac{k}{m_0}.\n", "$$\n", "\n", "On suppose $\\Delta<0$, si bien que la masse à un mouvement oscillatoire amorti. Sans conditions initiales, la solution générale est de la forme\n", "\n", "$$\n", "x(t) = A\\, \\exp(-\\frac{a}{2m}t)\\, \\cos(\\frac{\\sqrt{-\\Delta}}{2} t) \n", "+ B\\, \\exp(-\\frac{a}{2m}t)\\, \\sin(\\frac{\\sqrt{-\\Delta}}{2} t) \n", "$$\n", "\n", "avec deux constantes réelles arbitraires $A$ et $B$. \n", "\n", "Dans la suite, on considère $m=1$, $a=0.4$ et $k=11$ et on suppose $A=1$ et $B=0$ correspondant à $x(0)=1$ et $\\dot x(0) = -0.2$. La solution est donc\n", "\n", "$$\n", "x(t) = \\exp(-\\frac{a}{2m}t)\\, \\cos(\\frac{\\sqrt{-\\Delta}}{2} t).\n", "$$\n", "\n", "Sur brouillon, calculez $\\dot x(t)$ et $\\ddot x(t)$. Ecrire une fonction \n", "\n", "```\n", "x, xdot, xdotdot = xressort(t)\n", "``` \n", "\n", "qui retourne les valeurs $x(t)$, $\\dot x(t)$ et $\\ddot x(t)$ au temps $t$. Afficher le résultat de\n", "```xressort(0)```." ] }, { "cell_type": "code", "execution_count": null, "id": "eb6d8244-4b12-4736-ad05-a26597942ac4", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "23051076-08f4-4714-a60c-6af738f1633a", "metadata": { "tags": [] }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "54492e61-049f-442c-b8b8-4a8fdec24246", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "83dcf94b-8685-4779-a328-2223f0afe4fa", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ff41dcdc-59ba-4892-9a14-4e0fd7962852", "metadata": {}, "source": [ "L'exercice consiste à étudier si, à partir de la connaissance de la position $x(t^n)$, la vitesse $\\dot x(t^n)$ et l'accélération $\\ddot x(t^n)$ à différents instants $t^n$, **connaissant la masse $m$**, on est en mesure d'identifier les caractéristiques mécaniques du ressort, c'est-à-dire la **raideur $k$** et l'**amortissement $a$**. Pour chaque instant $t^n$, on a donc\n", "\n", "$$\n", "\\ddot x(t^n) = - \\frac{a}{m}\\, \\dot x(t^n) - \\frac{k}{m}\\, x(t^n)\n", "$$\n", "\n", "Les mesures de $x(t)$, $\\dot x(t)$ et $\\ddot x(t)$ ont sujettes à des incertitudes de mesures représentées par un bruit gaussien. Considérez le code python ci-dessous pour la génération de données. Dans ce cas,\n", "\n", "$$\n", "\\ddot x(t^n) \\approx - \\frac{a}{m}\\, \\dot x(t^n) - \\frac{k}{m}\\, x(t^n).\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "313da49f-942b-4cc1-be18-6935656d34b0", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "77e697c7-ec3c-4bd1-a9f7-7e1af8cbfd12", "metadata": {}, "source": [ "On considère le problème aux moindres carrés\n", "\n", "$$\n", "\\min_{u=(k,a)}\\quad \\frac{1}{2} \\sum_{n=1}^N \\left[ \\ddot x(t^n) + \\frac{a}{m} \\dot x(t^n)\n", "+ \\frac{k}{m} x(t^n)\\right]^2.\n", "$$\n", "\n", "Ecrire le problème sous la forme matricielle\n", "\n", "$$\n", "\\min_{u=(k,a)} \\quad \\frac{1}{2} \\|A u - b\\|^2\n", "$$\n", "\n", "en déterminant la matrice $A$ et le second membre $b$. Résoudre les équations normales et affichez les valeurs\n", "estimées $(\\hat k, \\hat a)$ et $k$ et $a$." ] }, { "cell_type": "code", "execution_count": null, "id": "46704a0e-50c5-40dc-a4db-91ffb1951cbe", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "26e72adf-65a2-4d61-b3fa-e239c3af5b99", "metadata": {}, "source": [ "Sur le même graphique, tracer la solution exacte $x(t)$ et $\\hat x(t)$ correspondant à la solution de paramètres mécaniques estimés $\\hat k$ et $\\hat a$." ] }, { "cell_type": "code", "execution_count": null, "id": "7d0ddf4b-9ca2-4e96-9ecb-7c92bd24db97", "metadata": { "tags": [] }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "bbbebb41-a357-4531-9f74-96aaeb955c4a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" } }, "nbformat": 4, "nbformat_minor": 5 }