{ "cells": [ { "cell_type": "markdown", "id": "363da173-3471-4c86-8279-422b91f8c6f9", "metadata": {}, "source": [ "# RANSAC Tutorial\n", "\n", "This notebook contains a simple demonstration of the RANdom SAmpling Consensus (RANSAC) algorithm. While old, it is provided as a demonstration of a method that attempts to work around noisy data.\n", "\n", "* First, read the classic paper, [here](https://dl.acm.org/doi/pdf/10.1145/358669.358692).\n", "* Second, read the online course notes, [here](https://mphy0026.readthedocs.io/en/latest/calibration/ransac.html).\n", "\n", "The aim of this notebook, is to illustrate the main points of the Fischler, Bolles paper, just fitting data to a line. We also have an additional notebook, [here](https://mphy0026.readthedocs.io/en/latest/notebooks/RANSAC.html), that shows how this is applied to pivot calibration.\n", "\n", "**Note:** The code in this notebook is meant to be simple and easy to read, not efficient." ] }, { "cell_type": "code", "execution_count": 98, "id": "bfeb9810-e21b-4e94-8024-80d91049abba", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 99, "id": "6e313a41-a142-4e49-af65-ff6256a1ce96", "metadata": {}, "outputs": [], "source": [ "# Jupyter notebook sets the cwd to the folder containing the notebook.\n", "# So, you want to add the root of the project to the sys path, so modules load correctly.\n", "import sys\n", "sys.path.append(\"../../\")" ] }, { "cell_type": "code", "execution_count": 100, "id": "628c89d9-7e6b-4c48-8ab7-57de6cab8798", "metadata": {}, "outputs": [], "source": [ "# All imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "8aea9e80-631f-4290-83f6-b9c9a7522d67", "metadata": {}, "source": [ "# Zero-mean Gaussian noise.\n", "So, imagine we had some data, that should fit to a line, but its corrupted by a small bit of noise." ] }, { "cell_type": "code", "execution_count": 101, "id": "23420348-f4d9-4889-a926-8e57a55b0984", "metadata": {}, "outputs": [], "source": [ "# Reference data for initial line-fitting test, fitting data to y=mx + c.\n", "number_of_samples = 100\n", "ref_m = 1\n", "ref_c = 5\n", "sigma = 2" ] }, { "cell_type": "code", "execution_count": 102, "id": "9b8c5257-1a63-4bfe-87c6-41e1ff2ccd56", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "m=1.0037316008860027, c=4.9589812721009\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUnFJREFUeJzt3Qd4VFX6BvB3ZtJ7AxKkRUCKSFMpgrq0BenK6tKUtiwqCKEq+gcLCAIKCAgISlGautIRlCYIht6kd0FICOm9zdz/cy7OMDeZJJNkkmnv73myYc69mQwjS17OOd93VJIkSSAiIiKyIWprvwAiIiKivBhQiIiIyOYwoBAREZHNYUAhIiIim8OAQkRERDaHAYWIiIhsDgMKERER2RwGFCIiIrI5LrBDOp0Od+/eha+vL1QqlbVfDhEREZlB9IZNSUlB5cqVoVarHS+giHBStWpVa78MIiIiKoHbt2+jSpUqjhdQxMyJ/jfo5+dn7ZdDREREZkhOTpYnGPQ/xx0uoOiXdUQ4YUAhIiKyL+Zsz+AmWSIiIrI5DChERERkcxhQiIiIyOYwoBAREZHNYUAhIiIim8OAQkRERDaHAYWIiIhsDgMKERER2Ry7bNRGRERExaPVSThyIx4xKZmo6OuBZuFB0Kht9zw7BhQiIiIHt+NsFD7cch5RSZmGsTB/D7zfrT46NQiDLeISDxERkYOHkzdWnVCEEyE6KVMeF9dtEQMKERGRAy/rfLjlPCQT1/Rj4rq4z9YwoBARETkgrU7CioM38s2cGBOxRFwXe1NsDfegEBEROcGek8KIjbO2hgGFiIjIAfecSMX4GlHVY2u4xENEROQEe05MUf1dzSNKjvVycnIwadIkHD16FNbEgEJEROQgjtyIN3tZR98BRZQa6/uh3LhxA8899xymTp2KPn36ICUlBdbCgEJEROQgYoqxlyTU3wNf9G0Cf083bDp1Bx/OXYLGjRvj0KFD8vVr167h3wOGWq3Ch3tQiIiIHERFM/eSTOpSD2H+npiy7Tzu3E9A/M4vkXZ2l+Ietbs3jurC0XrGHqs0dOMMChERkYNoFh4k7ylRFbHnRIST4WtO4Oals4haMSpfOHGv8jjCBs+Hd93WVmvoxoBCRETkIDRqlTzbIeQNKSqj2ZOPtpxF4uH1iP52HHIT7hrdpIZ/636o1GcaXPwqWrWhGwMKERGRA+nUIAyL+jeV95gYE4/FuC49Eae/ehuJvy4DdLmG6xq/CqjUdzoCWvWBSq2xekM37kEhIiJywJDSoX5ovtOLf/l5B/r3exWZCXGK+73qtEZQpxHQePjYTEM3BhQiIiI7p9VJ+cKIWO5pWTNYvp6VlYVxY8dg7ty5iq9TubojsN0w+DTsAJWqoJ0r1mnoxoBCRETkYG3tw/w9DJU3Fy9elHuanDp1SvF1rhXDUaH7BLgGVy3ye6j+XiIybuhW1hhQiIiIHKytfXRSJl7/9jhe9L6EJTMmIz09XXG9R///4GSlzlC7uBXZddZUQ7fywIBCRETkYG3ttZmpiPv5C8y9+JtiPCQkBCtWrECXLl1MzrwEeLnKnxPTcwxjYubEGn1QGFCIiIgcqK195l8XELtlFrTJMYrxDh06YOXKlQgLCyt0I63+ufPuZylvxS4z3r9/P7p164bKlSvLG2o2btyouC5JEiZPniy/AZ6enmjfvj2uXLmiuCc+Ph79+vWDn58fAgICMGTIEKSmppb+d0NEROSAMyWR1+LkdvTis74XSd6KGkmnReLBtbi35m1FONG4uGDWrFnYsWOHIZwYrv29kbZH40fkz+KxqTFrKPYMSlpaGho1aoTBgwfjpZdeynd95syZmDdvnpzSwsPD5RMRO3bsiPPnz8PD48HuXxFOoqKisHPnTvnUxEGDBuG///0v1qxZY5nfFRERkYNvgK1oVFGTm3wfsVs/Q9bts4qvdwkMw5fLvsHgnu1hb1SSmPIo6RerVNiwYQN69uwpPxZPJWZWxo4di3HjxsljSUlJqFSpkrzm1bt3b1y4cAH169eXj3F+6qmn5HtEquvcuTP++usv+euLkpycDH9/f/m5xSwMERGRs2yAVf39WRz0N2XbBVw7sgdxO+ZBl6lcifBu0A71e41C5OSuVpsFKc3Pb4t2khXHNEdHR8vLOnrihTRv3hyRkZHyY/FZLOvow4kg7ler1Th8+LDJ5xX12+I3ZfxBRETkjBtgpb8/f7DhJAKOr8D9jdMU4UTl5omQbuNQoctofPSvp2wmnFh1k6wIJ4KYMTEmHuuvic8VK1ZUvggXFwQFBRnuyWv69On48MMPLflSiYiI7G4DrF5WzA2c2jwTOXG3Ycwt7DGEdJ+AatVrWKXyxumqeCZOnIgxY8YYHosZlKpVi24sQ0REZI9iCmgpL7ZSpJzYioS9ywBtjmLLxavDRqH74FGoHORrtcobmw0ooaGh8ud79+4pdgqLx40bNzbcExOjLH3Kzc2VK3v0X5+Xu7u7/EFEROQMKppoKa9NT0LcT3ORce2oYlzs3fz222/Rtm1bOBKL7kERVTsiZOzevVsx2yH2lrRs2VJ+LD4nJibi+PHjhnv27NkDnU4n71UhIiJy9pLi6KQMBHm7GTbEZtw8hajlb+ULJ926dcfp06cdLpyUaAZF9Cu5evWqYmOs6O8v9pBUq1YNERERmDp1KmrXrm0oMxbpTl/pU69ePXTq1AlDhw7F4sWL5TLjESNGyBU+5lTwEBEROUtJsaTNReJvq5B8+Eej7bGigYkr3nz7AyyYOtGsQ/6cIqAcO3YMbdq0MTzW7w0ZMGCAXEo8YcIEuVeK6GsiZkpat24tlxHre6AIq1evlkNJu3bt5OqdXr16yb1TiIiInJGpkuKchCjEbpmJ7Chls1PPijXw2aKv8cZLjjdrYrE+KNbCPihERORIyzqtZ+xRzJykntuL+F8WQsrOUNz7Yt+BWPnlAvj6eMMeFefnt11U8RARETlDSbEuKx3xOxch7dxexT1qDx98PPsLvPPGa3AWDChEREQ2UFKcFXUZsZtnITcxSnHdvWoDhHQdh3ot28GZMKAQERFZUYi3G5IO/w+J+78FdNqHF1Rq+LfuC/8WL0Ol1ihKj8WykC2cOFyWGFCIiIisRByc+/4bryHx112KcY1fRVToPh7uj9STS41D/R+EkKIOELTnzrFl2geFiIiIzLNt2zY0bNgQu3crw4lX3WdRefB8QzgRRPgQMyT6ap+8bfCjkzLlcXHdUTCgEBERlaPMzEyMGjUKXbt2RWxsrGHc3dMTj/YaJ5+lo3Z/UKUjZk4W9W8qz4yYc4CguC7ucwRc4iEiIionFy5ckBuTnjlzRjHepEkTrF27FrVqP1bg3pIjRRwgKGKJuC7ua1kzGPaOAYWIiKiMiZZjS5YsxaiICGRlKnub9B78BnoNG4d4F3/5cUHhIqaAAwRLep+tY0AhIiKnVR7VMAkJCej271dxcOc2xbiLdwCqvjgOkRUaI3L9hSI3u1Y0cYCgKebeZ+sYUIiIyCmVdTWMCD9Lvt+GiW8NRVJstOKaR/iTCOkSAZ13oMnNrov+3ndiTIQn8frEPaZ2meSt9rF33CRLREROp6yrYbaduo3qHQbgzb49lOFE7YLANkNQ8eX3ockTTora7KpRq+TwJOSd48lb7eMIGFCIiMiplHU1zMpfjqJX1064s+dbQNIZxl2CHkHoq5/Cr9mLUKnUZm12zUvMqojZFTFTYsy42sdRcImHiIicSllWw3z//Q8YMmAQtJlpinHvJ9ojqP0wqN08zX6umAI2u4oQ0qF+KDvJEhEROZKyqIZJS0tDREQEvvrqK8W4ys0LwR2Hw7v+88V+nRUL2ewqwogjlBIXhgGFiIiciqWrYU6dOoU+ffrg4sWLinG3ynUQ0m08XANCi/X6HG2za0kxoBARkVMpbTWMvjT5XnIGfl3/DRbN/BDZ2dmKZ/Br+QoCWvWBSlO8H7OOuNm1pBhQiIjIqeirYUS1jogAUjECgr40+a+70Yj7aS4yrh9TPrdvCEK6joVHtSeKfB0BXq7y58T0HMNYqAMe+ldSDChEROR09NUwefugFBYQ9KXJ6TdOIm7bbGjTEhTXPWu3QPALo6Dx9C30ew9pVQPt64caZmgcfbNrSTGgEBGRUypONYxY1nl/w2nE712G5CPrFddULm4IbDcUPo06QaUqOFwU1ATO0Te7lhQDChEROS1zq2F+3HsUpxaORHb0FcW4a0h1+fRhtwrVC/zaEW1qoVWtEM6OFBMDChERUSGH/H3zzTd4/Y03kZ2Rrrjm27QLAv4xGGpX90Kfo3YlH86SlAADChERkQnJycl44403sGbNGsW42sMXwZ1Hwat2C6c6vK+8MaAQERHlcejQIfTt2xc3btxQjHtUa4jgrmPg4htS5HOwn0npMKAQERH9TavVYubMmZg0aZL8az2NRoNXR0zAXrcWUKk1JvunGGM/k9JjQCEiIrunb55WmnLdO3fu4NVXX8XevXsV4+Hh4Vi7di2aN29u6INiXJrMfiZlgwGFiIjsmqnQUFBJb0E2b96MwYMHIy4uTjEulnkWLlwIf3//QkuTBfYzsSyVJLYo2+HGJfGHJSkpCX5+ftZ+OUREZCX65ml5f5Dpo4FoxlZYSMnIyMD48ePxxRdfKMY9PL0w/qNZmDz6dbho1GXwyp1TcjF+fvNdJyIim1++ibwWh02n7sifxWP9uJg5MfWvbOnvj3d+/AMHr8YavsbYuXPn5GWbvOHELbQWgvrPwTex1fHszL1yCKLyxxkUIiKyy+Ubf0839Fl6yKznMV7yET/2vvzyS4wePRqZmQ+fV/Br9hICnnsVKs2DfSX6s3pGt6+NGiHeXL4px5/fDChERGSXyzeDW9XA1wdvmvVc+q+Z2S0caz/7P2zYsEFx3dUnEIGdR8MzvGmRz1Xc/S30EAMKERHZNbEk03rGHsXMSd7AEejtivi0h5UzRcm69Qfif5qN7KT7ivGWz7fH7ScGQuMdYNbzmLu/hfLjHhQiIrJroiKmoHAiiH9Zi3Di6+FiCAwF3qvTInH/t4he+64inLi5uWHOnDkYP3uZ2eFE/70FsfRkam8LWQbLjImIyOaIcl1zpGTmFno9N+keYjfPQtbdi4rxOnXqyL1NmjRpIm+8LS4RS0SAEkGK5+yUDQYUIiKyOZY4vybtwn7E7VgAKVt5yF+3l/th7fIv4e3tLT8Wm17FvpLopMwiO8SWNEhR8XGJh4iIbI4+NBSnVibQywX+nq7QZWcg9qe5iN08UxFOVO7eqPHKexg0cQbORGcalmdERY7Y9CrfU8zXyYMAyw4DChER2ZyShIaE9Fx0qJiOqJURSPtjl+Ka+yP1UXnQfEjhLTH6u1NyebLYhKvvcSI2u4pNr6JFvTnEaxIBigcBlh1W8RARkV31QTFFknRIOboJKb99g9xco8oelRr+Lf8N/1a95UP+iqrGMT7T52ZsOubuuvzg+Yv4OjIPy4yJiMhhiNCw4uANTNl2wfT1tATEbpuLzBvHFeOVwh7BiI/mYn2UP+LTsk1+rQgbYtbkwNttTTZfs8Q5P1Syn9/cJEtERDZNBIeBrcLx1YEb+TayZlw/jthtc6BLT1R8Ta9evbB06VJcjNfhq0K6zRZVjVPQ4YDsJFv2GFCIiMhu9qSIzrIiGui0OUjc9w2Sjyo7wnp6euLzzz/Hf/7zH6hUKsT8eafU1Tjie7OUuPxxkywREdkF/UZW/5xYRH87Ll84adiwIY4dO4ahQ4fK4aQ4VTasxrE9DChERGTzJxcLYstk1NEduPrlcGTfu6b4upEjR+Lw4cOoX/9B5Y+55cqsxrFdXOIhIiKbUNiG1JZVvfD6669j3bp1iq8JDg7G8uXL0a1bN7OWhkxV44jr3FNie1jFQ0RENn1yceadC5D2zsO9O7cV19q1a4dvvvkGlStXNuv5WY1jfaziISIiuyGWcUR4kEwd8nfoByQdWCMeGMZdXFwwdepUjB8/Hmq1eTsVWI1jfxhQiIjI5k4uzk2ORezWT5F1+6xi/NFHH5UP+WvWrFmxvw+rcewLAwoREVlV3hLf9Mu/I277POgyUxXjz3d5CZvXLOfSvpNgQCEiolIzbhFf3OUTfYmvLicLCXu+Quqp7YrrKjdPBP3zTUyfO5HhxIkwoBARUamUdgOqCDP+GVG4sGYKcmJvKa65hdVGhW4TULVGOEuBnQz7oBARUamrb/LuIREt6cW4/rTggohC0i8XL8KFxSPyhRO/5v9CWL+ZcA0MYymwE+IMChERWbT6RhBjIk6I66J6xlS4uBdzH//q+xoO7N6hGNd4ByK461h41mjMUmAnxoBCREQlYqr6xtyD+GYs+x8mjX4dOclxivGnn22HDz79AjluPiwFdnIMKEREVCKFHbBn7ODV+4agkZOTg1eHj8N3S+cr+7pqXBHUZjBimnaF2ssPPThj4vQYUIiIqETMPWBvwd5r+PHEHQxr4oNFH4zCkcOHFdddgqqgQo8JcKv4aJHLQuQ8uEmWiIhKpKiD+IxdjdyBwT3a5AsnPo06IWzgXDmc5F0WIudm8YCi1WoxadIkhIeHw9PTEzVr1sSUKVPkndp64teTJ09GWFiYfE/79u1x5coVS78UIiIqQ/qD+ISCQoouKx2x2+YgdsunkLIzDONqd2+E9JyI4E4joHb1KPHyETkuiweUGTNmYNGiRViwYAEuXLggP545cybmzxfrjQ+Ix/PmzcPixYvl47G9vb3RsWNHZGbyDyQRkT1U70Rei8OmU3fg7+mGL/o2Rah//pCRFXUFUStHIe3sbsW4e5XHETZ4PrzrtCr18hE5LovvQfn999/Ro0cPdOnSRX5co0YN+dyEI0eOGGZP5s6di//7v/+T7xPEaZSVKlXCxo0b0bt3b0u/JCIiKuOmbJO61MP5qBQs2HsVkqRD8pENSNz/LaDLffjFKjX8W/VB2PN9kGE0bEzMxIiww6ZsZPEZlGeeeQa7d+/G5cuX5cenT5/GgQMH8MILL8iPb9y4gejoaHlZR08cvdy8eXNERkaafM6srCz5iGbjDyIisp2mbMPXnISrRoXc1HjEfP8+En9drggnGr8KqNT3EwS0KjycCGzKRmUyg/LOO+/IAaJu3brQaDTynpSPP/4Y/fr1k6+LcCKIGRNj4rH+Wl7Tp0/Hhx9+yP9iREQ23JRt0bc/4N6PnyI3LVFx3atO6wd7TTx8Cv0eYuaETdmozALK999/j9WrV2PNmjV4/PHHcerUKURERKBy5coYMGBAiZ5z4sSJGDNmjOGxCEBVq1a14KsmIqLCDgEUvUwKasom5eYgft8KpBzbpBhXubojsN0w+DTsAJWq4BmRAE9XfNGvKVo8GsyZEyq7gDJ+/Hh5FkW/l+SJJ57An3/+Kc+CiIASGhoqj9+7d0+u4tETjxs3bmzyOd3d3eUPIiKy7n6TvHLibuP+5lnIibmuGPcKq4mALuPgGlz0PyYTM3KgVqkYTqhs96Ckp6dDrVY+rVjq0el08q9F+bEIKWKfivGMiKjmadmypaVfDhERWXC/iZ4oeEg5/QuiVkbkCydi1vz+tbP43zsv47WW1c36fiwrpjKfQenWrZu856RatWryEs/Jkycxe/ZsDB48WL4upvnEH96pU6eidu3acmARfVPEElDPnj0t/XKIiMiC+00EXWYq4nYsQPqlA4pxF+8ArF+3Ct26PqjibFnzQanwN5F/Fvk9WVZMZR5QRL8TETjefPNNxMTEyMFj2LBhcmM2vQkTJiAtLQ3//e9/kZiYiNatW2PHjh3w8OAfUCIiWz4EMPOv83LTNW1yjGLco0YTfL18Obr9o5HJbrOi0sdU4GFZMRVEJRm3eLUTYklIlCYnJSXBz8/P2i+HiMihiAZso9adUoxJOi2SIr9H0sG14sHDC2oXVOs4GIs+mYzODR8pdLlIfh6jcf2Ok0X9m7Jyx0kkF+PnN8/iISIixfJObEqWYiw3+T7urX0XSQdWK8JJSOXq+PrH7bi+dXGB4UQQ4UOEkLzdZsVjhhMqCE8zJiKiAqt20i4dRPz2edBlpSnuDWn6T1zZ8wMC/M2bxRYhRJxQLJaPxIZYsedELOuwcocKwoBCRORkPU1MBQT9Mox+CUaXk4mE3V8h9fQOxXOo3DwR3HE4vp02xuxwoie+V8uawRb7/ZBjY0AhInLiM3RE51Yxs2FctZMdcx2xordJ3G3Fc7iF1UH9fu9h+oD2XJahMseAQkTk4PLOjuiJyhoxHtG+thxc5N4mJ7YiYe8yQJtjdKcKfi1fxsyPP8J/nn+MyzJULhhQiIic/Ayd5QdvQpuehLif5iLj2lHFPRqfIAR3HQvP6o0QGujDcELlhgGFiMhJe5roQ0rUhaOI2zYb2tR4xTXPWs0R/MJIaLz85cdspkbliQGFiMiBFdZCXtLmIvG3VUg+/KOyQ4nGFUFt/wOfJp3l7t9spkbWwIBCROTAFTtX7qWYvJ6TEIXYLTORHXVFMe4aUg0h3SfArUIN+bF+QUdspuXyDpUnBhQiIic7hTj13F7E/7IQUnaGYrzLvwfgfv1XcC/94WyKmDkR4YRVO1TeGFCIiJygYkfQZaUjfucipJ3bqxj39Q/EyuVf48UXXyy0VwpReWJAISJygoqdrKjLcm+T3MQoxfgTT7XETxu+R5UqVeTHbKZGtoIBhYjIgSt2JEmH5MPrkfjbt4BOaxjXaDSY/P4HaPvv/+J4bA5uZ8VxtoRsCgMKEZGDVuzkpsYjbutsZP6pPJm4YuWqeGfGF/j+theWLTuar7Ms95uQLeBpxkREDsK4T0n6taOIWjYiXzjxqvssRi/ciHln1flmW/SdZcU+FiJr4wwKEZGd029sjU7KQIAbcP2nL5FyfIviHpWrB4LaD0PNVl2w8VxCoZ1lxT4WcT4Pl3vImhhQiIgcpKQ4J/Y27m+egZz7NxX3uFWqiQrdJ8A16BH0bV4dc3Ype5/k6yyblCkHHm6WJWtiQCEistMZk53no7Hs4E35kL/U0z8jYfdSSLlZint9n+6JwOcGoHKwr7y/JCtXV+oOtETlgQGFiMiOm7BpM1MRv30e0i//rrhP7RWA8F7jMGv0QIT6PexnEnktzqzvw3N3yNoYUIiI7LQJW+bts4jd8hm0KfcV93mEP4mQLhHI9Q6Uw4nxUo0IKqJaR2yINbUPhefukK1gFQ8RkZ01YZN0WiT+thr31r6rDCdqFwS2GYKKL78PjXegyaUaMYsilnqEvFtgee4O2RIGFCIiO2rClpsUg3trJiLp97UiqRiuuwQ9gtBXP4VfsxehUqkLXaoRfU4W9W8qz5QYE4/FOPugkC3gEg8RkR0QMyFpFw8gbsd8SFlpimveT7SXS4jVbp5mL9WIECJKiXnuDtkqBhQiIhtf2tl37jamvTMGsT//qLimcvNCcKcR8K73nHLczKUanrtDtowBhYjIhjfFTvhyMy6smYrc+L8U19wr10VI9/Fw8a+U7+vEzAlb1pO9Y0AhIrJB2/+4i34RHyBh33JAm2t0RQX/Z/4N/1Z9oFJrFF8zpFUNtK8fyqUacggMKEREVm64lncPSFT0PfT914tIvHxEcb/GNwQhXcfCo9oTinEe8keOiAGFiMgGGq7pg0bX4BjMfi8CibExivs9H2uJ4E4jofH0NYyNaFMLrWqFcMaEHBIDChGRlRuuCZI2B+c3LsOhI+sV96pc3BDYbih8GnWCSqUMIbUr+XCTKzksBhQiIis1XNPLib+D2C2zkB19VXGva4UaCOk2Hm4Vqpt8LrajJ0fGgEJEZIWGa4I45C/t7B7E71wEKUfZ8TWkeQ94tXoNKlf3fM/BdvTkDBhQiIjKkb71vC4rDXE/L0T6hX2K62pPPwS/MAojBv4byw/elMeMZ1vYjp6cBQMKEVE5EssyWXcuyks6uUn3FNfcqzVESNcxcPENkbu8ihmSvBtp2eOEnAUDChFROdFqtdjz3ZeIXjMJ0D08RwcqNQKe7Q+/5r2gVmsMyzdihoTt6MlZMaAQEZWDO3fu4LXXXsOePXsU46ITbEj3CXCvXEdevhHLOb2froqtZ+4aAgkrdcgZMaAQEZWxzZs3Y/DgwYiLi1OMBzdqC682r0Pt7iU/9vdylT/P2XXFcA+bsJGzengmNxERWVRGRgZGjBiBHj16KMKJj48PVq5ciejjO/HdiLb4vHdjjG7/GJLSc5CYnqN4juikTLlniuidQuRMOINCRFQGLesvXjiP3r174+zZs4p7a9ZriGnzl6JXm6cNpwmLr289Y4+iWkdPjImlH7FZVuxH4f4TchYMKEREFmxZL3qbaC7twp0dXyI7S9nbxK/ZS8h57lVM2Hkfc47sMSzdGPdGQQEhRVwX93E/CjkLBhQiIgu1rNdmJCNu+zxkXDmkuE/jHYjgLmPgGd4k39LNov5NkZVrVNFjRg8VImfAPShERBZoWZ956w9ELXsrXzgJeKwZwgbNV4QTQf914jlCfPJ3izWFre3JmXAGhYioBPTLMpJOi6QDa5AU+b2y56vGBYH/GATfJ7vnO+Qv79KN+IWo1hGzKqb2obC1PTkjzqAQEZWAWG4RnWDvrX4bSZHfKcKJS1AVhL06G35P9SgwnBiLTcuS96MIee9ma3tyVgwoREQlcHrfT7i77C1k3b2oGPdp+E+EDZgLt0qPmv1cYulGbJYV+1HETIkx8ViMsw8KORsu8RARFaOEuH4FN4yOGIXly5cr7lO7eyOo01vwrtva7OfOu3QjQghb2xM9wIBCRGRmCXFW9FUkbp2FzLg7ivvcq9RHSLdxcPGraPZzF7R0o++NQuTsGFCIiIooIZYkHVKObkLCvpWALtdwj1qtRp9hEbhauSOiU5UdYIvCU4mJCseAQkRUSAmxNi0BsdvmIPPGCcU9bv4V8POm/+Efzz9nWAY6ePU+Fuy9VuTzT+pSDwNbhXPphqgQDChERAWUEGdcPy6HE116ouK6V51W8n6TY5mV4H4tznDisPj844k7RZYLM5wQFY1VPEREedyJS0L8nq8Q88P7inCicnGXg0lIj3eg8fDBgr1X0WfpIfkcHbEkJEIHy4WJLIMBhYjob2KpZt3OQ4jo2xUpRzcqrrlWqCGXD/s26pivt4nxicMsFyayDJUkTrayM8nJyfD390dSUhL8/Pys/XKIyAFKiG/cT8P8xUtxY8sCSDnKM298n+wmd4VVubgV+Fz65ZsDb7eVZ0hMnXDMmRNydsnF+PnNPShEBGcvIdZlpiLu5y+QfvE3xT1qTz8Ed46AV61mRT5f3hOHWS5MVDoMKETkFIxnNG7GpmPurstyqMj86wJit8yCNjlGcb9H9UYI7joWLj7FO/+GJw4TWQYDChE5XcM1QT7k79AP8kF/kHQPb1ZrEPDsq/Br/hJUqofb9Ho2royNp+4W+b144jCRDW+SvXPnDvr374/g4GB4enriiSeewLFjxwzXxbaXyZMnIywsTL7evn17XLlypSxeChE5yexI5LU4bDp1R/4sHuvHPtpyDq+vOqEIJ7nJsbi37j0k/bZKEU5cAsIQ2m8m/Fv8SxFOhOcfqyCfOFzQLhIxLq7zxGEiG51BSUhIQKtWrdCmTRts374dFSpUkMNHYGCg4Z6ZM2di3rx5WLlyJcLDwzFp0iR07NgR58+fh4cH//VBRKWbHQnwcpU/J6bn7+6afjkScdvnQZeZohj3frwNgjq8AbW7l8nvE+rvKZcIi2odEUaMqwtYQkxkB1U877zzDg4ePIjfflNuNtMT365y5coYO3Ysxo0bJ4+J3byVKlXCihUr0Lt37yK/B6t4iChvO/qi6HKykLDnK6Se2q4YV7l5Iuifb8Ln8TZmVeeYCkRi5oRt64lsvIpn8+bN8mzIyy+/jH379uGRRx7Bm2++iaFDh8rXb9y4gejoaHlZR0+82ObNmyMyMtKsgEJEZNyOvijZ928idvNM5MTeUoy7hT2GkG7j4RpoOliYmhnhicNE5cPiAeX69etYtGgRxowZg3fffRdHjx7FyJEj4ebmhgEDBsjhRBAzJsbEY/21vLKysuQP4wRGRM5N346+MGLGNvXkNsTv+RrQGi/3qODXohcCWveHSuNS7AP9WEJMZIcBRafT4amnnsK0adPkx02aNMHZs2exePFiOaCUxPTp0/Hhhx9a+JUSkT0rqpxXm54k7zXJuHpYMa7xCUJwlzHwrNE439eE+rmjT7NqqBHizZkRIkcLKKIyp379B2dR6NWrVw8//vij/OvQ0FD587179+R79cTjxo3z/4UhTJw4UZ6RMZ5BqVq1qqVfOhHZkcLKeTP/PIPYrZ9CmxqvGPes1QzBL4yCxsvfsNF1dPvaDCREzhBQRAXPpUuXFGOXL19G9erV5V+Lqh0RUnbv3m0IJCJwHD58GG+88YbJ53R3d5c/iIj0RJgQm1ONTw6WtLlIPLgGyZE/KOtsNK4IbDMYvk27Gs7RKWj5hogcNKCMHj0azzzzjLzE88orr+DIkSNYsmSJ/CGIvxwiIiIwdepU1K5d21BmLCp7evbsaemXQ0QOSn9ysL7sNzsxGrGbZyE7SvkPJNfgqgjpPgFuFcPlx0Na1UD7+qGcLSFyxsMCt27dKi/LiP4nIoCI5Rl9FY8gvuX7778vh5bExES0bt0aCxcuxGOPPWbW87PMmIj0RNnvW1MX4NqGOZCyMxTXfBp3QmDb/0Dt6sFSYCIbUJyf3zzNmIjsVkpKCkaMGIFvvvlGMS4aQy5ZshSPNPkHS4GJbAhPMyYihz/07/dDhzH3vRH4688biuvPPfccVq1axY30RHaOAYWI7Go554NNZ3Fp1zok7v8G0OUarqnVannp+L333oNGo7Hq6ySi0mNAISK7CSdDF+9C7LY5yLx5UnFN41cBMxYsxdhXe1jt9RGRHZxmTERk6WWdiE+X4+7yt/KFE6+6z6LyoPn4318+8n1E5Bg4g0JENre/xHhjqzjmosdrb+DS/1Yo7lW5uiOw3TD4NOwgty8Qbe/F17IFPZFjYEAhIptg6pRgj7Qo/Pm/6ciMvq6417Xio6jQfbzc46Q47e+JyH4woBCRTYQT0XBNMj7k78wvuLVrCaTchweFCr5P9UDg8wOhcnEtVvt7IrIvDChEZPVlHTFzog8n2sxUxO+Yj/RLBxX3qb38EdJ5NDxrPpXvOVR/t64XS0JE5BgYUIjIqsS+Ef2yTuZf5xC75VNok+8r7vGo0QQhXcZA4xOY7+v1rddEl1g2YiNyHAwoRGTVDbHbz0ZB0mmR9Pt3SPp9HSDpHt6kdkHg86/B9+meUKlMFx3y0D8ix8SAQkTlUo1jPLthvCE2NzlGnjXJ+uu84jlcAivLh/y5h9Yq8PtM6lIPA1uFc+aEyAExoBBRmVfjGB/UZ7whNu3iAXm/iS4rTfEc3g3aI6jDMKjdPE1+D/2eE4YTIsfFgEJEZVaNoxedlCmPf9G3CaZsuwBtdiYS9ixF6umfFfep3LwQ3HE4vOs/X+D34J4TIufAgEJEFlm+yVuNY0yMiSjxf5vOIvrGJdzfNBO58X8p7nGrXAch3cbDNSDUMBbg9aCUODE9xzDGPSdEzoEBhYhKtHfEePlG7AURY8bjeekkCTf3/Q8Jvy4HtA8P+RPRxa/lywho1RcqzYO/kl5rWR0vNAgzlA0X9pqIyDExoBBRgeEj1M8dfZpVQ40Qb9yMTcfcXZfzzZCI+99cozwfJy9tehLifpqLjGtHFeManyCEdB0Hj+oNFeMinBi3rGf7eiLnw4BCRAXvHUnOwpxdV0r13Bk3TiJu22xo0xIU4561miP4hZHQePkbxthwjYj0GFCInFxhe0dKQ9LmIPG3VUg+/KPygsYVQW3/A58mneVD/vS4+ZWIjDGgEDk5406ulpKTcBexW2YhO0o5++IaUk3ubeJWoUa+r+HmVyIyxoBC5OQsfQJw6tk9iN+5CFJ2hmJczJgEthkCtau7YjzA0xVf9GuKFo8Gc+aEiAwYUIicnKVOANZlpcvBJO3cXsW42sNX3mvi9VhLk1+XmJEDtUrFcEJECgwoRE5ObEgVpcKimVpJ96FkRV1G7OZZyE2MUox7V2+IgM5j4OIXUq6zOERk/0yfvkVETkPMXIi9H0Jx5zAkSYekQ/9D9KrxynCiUiPg2Vcx5ct1RYYTS87iEJHjYEAhInlj6qL+TeWNqoXRB5jR7Wvj/bahqHjgUyTuWwHotIZ7NP6V8PiwuVi7cAZGdqgrz84UFHzEuLjOsmIiyotLPERkCCkd6ocauraKxmxrj9xCdHJmvkqb3JvHMWjQIMTGxiqeo/U/u2PyjLlo27CGYU+JuF/0WBGPjJeQWFZMRIVRSZJk6fYHZS45ORn+/v5ISkqCn5+ftV8OkdO0vm8Y5oWJ77yN+fPnK+7z9vbGggULMGDAAEVvE3NPOCYi55BcjJ/fnEEhogKJmQ19m/nz58/jmZZ9cObMGcU9TZs2xdq1a/HYY4+ZPTvDM3WIqCgMKERUKDHJumTJEowePRoZGcreJmPHjsW0adPg5uZWrLBDRFQUBhQiKlB8fDyGDh2K9evXK8YrVaqElStXomPHjlZ7bUTk2BhQiMik/fv3o1+/fvjrr78U4506dcKKFSvkkEJEVFZYZkxECrm5uZg8eTLatGmjCCeurq6YPXs2tm3bxnBCRGWOMyhEZPDnn3+ib9+++P333xXjderUkTfCNmnSxGqvjYicC2dQiEj2/fffo1GjRvnCyZAhQ3D8+HGGEyIqV5xBIXJyaWlpGDVqFL7++mvFuOhVIKp3XnnlFau9NiJyXgwoRE7s5MmT6NOnDy5duqQYf+aZZ7B69WrUqFHDaq+NiJwbl3iInJBOp8OcOXPQokULRThRq9XyBtl9+/YxnBCRVXEGhcjJxMTEYODAgdi+fbtivEqVKvKsyXPPPWe110ZEpMeAQmRHZ+GUtj38L7/8gtdeew337t1TjL/44ov46quvEBTEU4WJyDYwoBDZKHMP2DMnxGRnZ+Pdd9/FZ599phh39/DE4LHvo//AIfAPCCyH3xURkXl4mjGRjYaTN1adQN7/c+pjx6L+TeWQYk6IuXLlirwRVpQKG/MKDYd/l/FwC6lm8uuIiKz585sBhcjGiBmR1jP2KEJH3pAS6u+BSV3qY/iagkPMwn5NcO/4LxgxYoRcSmzMt2lXBLYZDJWLW4Hhh4jImj+/ucRDZGPEck1B4UQQgURc/79NZ/OFE/11KSsNr/bvj7gzexXXgoODEdI5ApmVm5j8OhFSxIxMh/qhpdrrQkRUWiwzJrIxYi+JOeLTsk2OZ925iDvLR+YLJ23btsWyTb+aDCd5w48ISURE1sSAQmRjxEbXkpB0WiRFfo/o1ROgTXpYpePi4oLp06fLFTySd6BFQxIRUVnhEg+RjRFVOGLDanRSpsklHLHwEujtivi0HMNYbkosYrd+hqxbfyjuDQqtghkLvsagFzvISzbmhp+ShiQiIkvhDAqRjRFBQlTTCHl3gegfT+3RQA4x4nH6lUOIWvZWvnDiXf8f8O49G1OP5sibbkXFjz78FLS7RIyL6+I+IiJrYkAhskGiikZU04hqHWPisRjv3LAy3ukQjrhfFuH++qnQZaYY7lG5eSK4y2iEdBsHtbuXPCZmY15fdQIL9lzBCw1CDRtijekfi3DEDbJEZG0sMyayYQU1YTt37hx69+6Ns2fPKu53C62NkO7j4RpYucjnFhlEZ/T/fvZBIaKyxjJjIgchwkjLmsGGx+LfE4sWLcKYMWOQmancyPqPl4fgevWuUGlczXpufTgZ0qoG2tcPLXUbfSIiS+ISD5GdiIuLw0svvYQ333xTEU5CQ0PlCp2Idz80O5zoiTjy09lohhMisjkMKERWXL6JvBaHTafuyJ/F44L8+uuvaNSoETZu3KgY79y5M86cOYMOHTqUqPKGfU+IyFZxiYfIhg8CzMnJwYcffohp06bJyzt6bm5umDlzJkaOHAmVSmVWeXJh2PeEiGwNZ1CIrHQQYN529iJYiHFxXbhx4waee+45fPzxx4pwUrduXRw+fBijRo0yhJOiypOLwr4nRGRrGFCIypFYxhEzJwWdoSOI66tXr0Hjxo1x6NAhxT1Dhw7F4SNHkeFb1eTSUEHlyQVh3xMislVc4iGyoYMAtdkZ+GPNXPQ/u0sxHhAQgKVLl8Knbit0XHC40KUh8Vkc9qcvT74Zm465uy7L14yDEfueEJEtY0AhKkeF7fXIir6K2M0zkZtwVzH+7LPPYtWqVTif7CovAeWdfdEvDYmZE31IyVueXCfUJ9+eFzHLwr4nROS0SzyffPKJvE4eERFhGBMlksOHD5ePfvfx8UGvXr1w797Dw82IHJWpvR6SpEPykfWI/nacIpyo1Wp88MEH2LNnDx6pUtWspaGCKoFECDnwdlusHdoCn/duLH8WjxlOiMgpZ1COHj2KL7/8Eg0bNlSMjx49Gtu2bcMPP/wgd5QbMWKE3N/h4MGDZflyiKwub6WNNjUBsdtmI/PmScV91apVw+rVq9G6dWv5sdhrUtjSkHG5sPHMibG8sypERE45g5Kamop+/frJ6+aBgQ+PeBftbb/++mvMnj0bbdu2xZNPPonly5fj999/z7chkMjRGFfaZFw/jrvL38oXTlp36IpTp04ZwklxyoBZLkxEjqLMAopYwunSpQvat2+vGD9+/Ljc28F4XJRNin8xRkZGmnyurKwsuX+/8QeRvWpTOwhN7m5GzA/vQ5eeaBhXu7pj1AezsP/nzYpQX5wyYJYLE5GjKJMlnnXr1uHEiRPyEk9e0dHRcpMpUZVgrFKlSvI1U6ZPny43qyKyd5cuXUKfPn1w8qRy1qR2vQZY/8N3aPD4g9mVvIpqwqb6e9Mry4WJyFFYfAbl9u3bcgMpsX7u4WGZf81NnDhRXhrSf4jvQWRPcrU6vPfJPDRq3CRfOBH/fzlz4miB4aSoJmwsFyYiR2TxgCKWcGJiYtC0aVO4uLjIH/v27cO8efPkX4uZkuzsbCQmPpzaFkQVjzj0zBR3d3f5WGbjDyJ78b/fL6JSk7aYNnEUsjIzDON+gUHYsmUL5s6da1aYL6gJm3hsXGJMROQILL7E065dO/zxxx+KsUGDBsn7TN5++21UrVoVrq6u2L17t1xerJ/2vnXrFlq2bGnpl0NkVbNXbcaE4UOhTY5RjHtUbwy/rmPgUuPJYj1f3iZsYs8JTyImIkdk8YDi6+uLBg0aKMa8vb3lnif68SFDhmDMmDEICgqSZ0PeeustOZy0aNHC0i+HqMyJ3iN5AwMkHaZ+/DE+EHundLqHN6s1CHjuNfg1exFqlVruXSICR3ECBsuFicgZWKWT7Jw5c+QmVGIGRVTodOzYEQsXLrTGSyGy+KnEQVIycnbPw9njyrJ5l8AwhHSbAPew2mb3LiEiclYqyfiYVDshyoxFgzexYZb7UcjapxIb/x8o/fLviNs+D7rMVMW93g3aIqj961C7e+V7HtHZtUfjR8rhFRMR2c/Pb57FQ2SBU4l1OZlI2PMVUk/tUNyncvNEcMfh8K7/D7N6l5haLuL+EiJyRgwoRMUkQsSKgzcMyzrZMTcQu3kWcuJuKe5zC3sM1f/1DrK9Kpp8nry9S0wtF+U9qZiIyFkwoBAVg3GIEKujKSe2ImHvMkCbY3SXCn4teiGgdX9ka0z/Xyxv7xJTy0UFnVRMROQMGFCIzGQcIrTpSYjb/jkyrh5R3KPxCUJwlzHwrNG40OcKNZoZybtcZEz6O8yUpNqHiMieMaAQmbH/wzhEZPx5GnFbP4M2NV7x9Z61miH4hVHQePkX+D0CPF3xRb+maPFosOG5xfcs7UnFRESOhgGFyIz9HyIc3I1PReKBVUg+9OPfseFvGlcEthkM36ZdoVIVPsORmJEDtUqlmAnhScVERPkxoJDTM2f/x80b1xG9egKyoy4r7nENroaQ7uPhVjHc7O+XN2jwpGIiovwYUMipmbP/460p83F323xkpyl7m/g0fgGBbYdA7Vq84JA3aPCkYiKi/BhQyKkVtv9Dl5WO+F2LkXZ2j2Jc7eGD4E4j4VXnGUWIqOTnLv/qXnLxgob+pGIxWyPuMf5anlRMRM7K4qcZE9nT7MnBq7Emr2VFXUbUylH5wol71QaoPGh+vnAifND9cXzQvb5izNygwZOKiYiU2OqenJKpTbGCJOmQfGQDEvd/A+i0hnGNRoP3338fTbsPwtSfLhXaTK00DdfYSZaIHFlyMX5+M6CQ0yloU2xuajzits1B5s2TinG3gErYteVHPNu6ldkhgkGDiCg/nsVDVMxNsRnXjiJ22xzoMpIV4151n8XKr5fg2WfqGsZE0CiqH4k59xARUcEYUMipN8VKudlI+HU5Uo5vUdyncnVHeLcRWPDBGLzwRGUrvFIiIufGgEJOxbgHSU7sbdzfMhM5MTcU97hVqolRU+dj+uBOXJYhIrISBhRyKmI/iNh2lXrmFyTsWgIpN0tx3fepHgh8fiBebNOM4YSIyIoYUMip1A5QIfWnWYg/u18xrvYKQEiX0fB69Ek2RSMisgEMKOQ0Dhw4gH79+iH+1i3FuEeNJgjpMgYuPoHyYzZFIyKyPgYUcni5ubn4+OOP8dFHH0Gn0xnGVRoX+D/7Gvya9YRKpZZnTszpVUJERGWPAYUckr4PybnL17Bg8iicPnZIcb127dpYtXoNtEHh7FVCRGSDGFDIrpjTAE3fyfXq4V2I3zEfuqw0xfVBgwZh3rx58PHxKedXT0RE5mJAIbthTgt5cc+wZb8jfs9SpJ7+WfH1KjcvTJj6KT4Z/4bZ35MdYYmIrIMBhey6PX10UqY8Lg7U61A/FBOWbMHd1VOQG/+X4j63ynVQodt47MutIYcOc0JGac7UISKi0uFpxmS37ekF/dgHm89hzORp+GPh8DzhRAW/lv9GaN8ZcAkIlcOGmBExNxDlPUxQH4jEdSIiKjucQSG7a0+fV25aIs78MBeHrx9TjGt8ghHSdSw8qjcssJtsSQKRmHsR18WMDZd7iIjKBgMK2bzCAkXGjZOI2zYb2rQExbhn7RYIfmEkNJ75T8sUe0lKE4hESNHPxPBAQCKissGAQjbPVKCQtDlI3P8tko+sV4yrXNwQ2PY/8Gn8AlQq5eyGeGROl9iiZliKex8RERUf96CQzROBQmxO1ceNnIS7iF41IV84adCgARZ+vwN+TTpDbSKcmNsltqgZluLeR0RExceAQjZPBAoRLMQhf2lndyNqxShkR19R3NO6ez9MWroBjRo2xBd9m8ozJcbEY1HpY071Td5AlJcYF9d5Xg8RUdnhEg/ZjMJ6jjxTzRvhZ5fh1582KL5G4+mLqj3H4Ha1pzFh4yV5TISHSV3qIdDbvUT9S/SBSFTriK+QSjgTQ0REJaeSxD9L7UxycjL8/f2RlJQEP7/8myDJ/hTWcyQw7Rb69u2L69evK76mVqPmyGz1BjS+IYpxfWwwd8akJK+JfVCIiMr25zcDCll9dqSgJmyQdEg6/CNSDq6GNjfXMKzRaPDBhx9iK55GdEqOye+n3xB74O22pZrpYCdZIiLr/PzmEg9ZXN4f6glp2ZiyzfRMhOglYqrnSG5KHOK2fYbMP88oxsPDw7FmzRpIFWrjq6XKAwDLohRYhBGWEhMRlT8GFLIoU8sipug7ska0r53v3vSrhxH30+fQZSQrxvv06YNFixbJ6XvTqTtmvR6WAhMR2ScGFLKYApdqCunIuvzgzYdjudlI+HU5Uo5vUdyrcvVAtzf+DyPeGgYf3wdTgiwFJiJybAwoZBGFtYcviLg3MePBHpLs2FuI3TwTOfcfBhbBLbQWQrqNx2nPR9D3q8OKpSHxazETY+p7mtuUjYiIbBP7oJBFFNUeviBij3bWHz8jeuXofOHEr9lLCO0/C65Bj+RbGtp5PloOKkLeLassBSYisn8MKGQRJdnroc1IQezG6Yj+aT6k3CzDuNo7ABVf+QiBbQZDpXFVfI1+tkR/WJ8oJTanKZuY4Ym8FifvXRGfxWMiIrJdXOKhUhM/7GNTHgYMc2Te+gOxWz+DNiVWMe7x6JMI6TwaGu8Asyp0RAgRQaWwUmD2MyEisj8MKFQuVTt6kk6LpINrkRT5vdznxEDjgsDnB8H3qW7wdnNFeo7W7FmbwkqBC9q4q18qKm0zNyIiKhsMKFQuVTtCbtI9xG75FFl3LijGXYIeQYXuE+BWqab82JxwYk6FTmEbd/VVRPqlIu5VISKyLdyDQmVetSOWU3oH3ULS6tH5wolPw38ibMDnhnCiF+DpWurD+orauGu8VERERLaFMyhUoi6xB6/eN2tZZ3ybajj+3RzMWL5cMa5y90ZwxxHwrvesya8b1Cocc3ddLtVhfeZu3GUzNyIi28OAQmW23yQr+iqmDRuJu38qD/nzqfY4AjqPhca/YoH9S0a0rYU6oT75vl9oMTa3spkbEZH9YkAhy3eJlXRIObYZCb+uAHQPD/lTq9X4v//7PzR78T8Yse7BGTuFzY6YU6FTGHEvm7kREdknBhSy6H4TbVoCYrfNReaN44rxKlWqYPXq1XjuuefkxxoXF7NmR0pzWJ/4WvF8IliVZqmIiIjKHwMKWaxLbMaNE4jdNhu6tETF+EsvvYSlS5ciKOjhTEVpZ0fMJb6PKCUuzVIRERGVPwYUKvUmUkmbg8R93yD56AbFuKenJ+bOnYuhQ4dCpcofPEozO1Ic5RWGiIjIchhQqFSbSHPi7yB2yyxkR19VjD/xxBNYt24d6td/cF6OtZVXGCIiIstgHxQye7OpKs8hf6l/7EbUilH5wslbb72FI0eO2Ew4ISIi+8OAQmZvNhVESNFlpcmzJnE/zYGU83D5Jzg4GJs3b8a8efPg4cHSXSIiKjkGFCrWZlPvpOu4u3wk0i/sV1xv27Ytzpw5g27dulntNRIRkePgHhQyq3NsVGIafl7zJS4unQGt9uFZOS4uLpgyZQrGjx8PjUZj1ddKRESOgwGFiuwce/v2X4jd9hmybv2huP7oo49i7dq1aNasmdVeIxEROSYGFCq0c2za5UjEbZ8HXWaK4nqbLi9h45rl8PPzs9prJCIix2XxPSjTp0/H008/DV9fX1SsWBE9e/bEpUuXFPdkZmZi+PDh8qZKHx8f9OrVC/fu3bP0S6FSLOu8v/4kYn9ZiPsbPlaEE5WbJ0K6jEFGqzfh7eNr1ddJRESOy+IBZd++fXL4OHToEHbu3ImcnBz885//RFpamuGe0aNHY8uWLfjhhx/k++/evSt3G6WyDR2R1+Kw6dQd+bN4XJC1Px/EyflvIPXkT4pxt9DaCBv4ObwbtJW7soq9KURERGVBJYmGFmXo/v378kyKCCLiHJakpCRUqFABa9aswb/+9S/5nosXL6JevXqIjIxEixYtinzO5ORk+Pv7y8/FJYaSnUIcZqLVu/ijsGjRIoweMwbZWVmK5/Br3gsBz/aHSuNqGPu8d2P0aPxIOf0uiIjI3hXn53eZlxmLFyHoz2E5fvy4PKvSvn17wz1169ZFtWrV5IBiSlZWlvybMv6g4u0lyXuWjjjhV4yL60JsbKy8HCdmv4zDicY7EBVfmYLAfwxShJOiOswSERHZ7CZZnU6HiIgItGrVCg0aNJDHoqOj4ebmhoCAAMW9lSpVkq8VtK/lww8/LMuX6nSnEEt/N10T113uXcCA116Vl9qMedZ8GsGdI6Dx8leMq/4+bE90mCUiIioLZTqDIv41fvbsWflMltKYOHGiPBOj/7h9+7bFXqMzn0Ks0+bi/JYl+GeH9opwIgLk6+98hIq9JsPFRDgRxPIQD9sjIiK7m0EZMWIEtm7div3796NKlSqG8dDQUGRnZyMxMVExiyKqeMQ1U9zd3eUPstwpxDmJ0YjdPAvZUcoKK7HcJgJlo0aNTO5dCTWxd4WIiMjmA4rYaCkOi9uwYQN+/fVXhIeHK64/+eSTcHV1xe7du+XyYkGUId+6dQstW7a09Mtxqm6vIpCIfSFi6UXMbhS0RyTt/D7E/bwAUnaGYtynUSdMX7oAjRrVlB+LENKhfqjJ5yYiIrKrgCKWdUSFzqZNm+ReKPp9JWLXrqenp/x5yJAhGDNmjLxxVuziFYFGhBNzKnhIyeQsh587+jSrhmpBXgjydkNCWra850SXlY74XV8i7exuxXOo3b0R9MJI+NRphRm7bqLbk48aQoj43LJmcLn/voiIyLlZvMxYpTL9r+vly5dj4MCBhkZtY8eOldukiwqdjh07YuHChQUu8eTFMmNlhY45/wGzoq8idvMM5CY8qNrRc6/yOEK6jYWLX0XD2NqhLRhKiIjI4orz87vM+6CUBQaUB8s6rWfsKXQTrCBJOiQf2YjE/d8AutyHF1Rq+LfqA/+Wr0ClVh7y91rL6nihQRiXc4iIyGo/v3kWj4NW6Aja1ATEbpuNzJsnFeMavwoI6TYOHlUeN/l130T+KX+YauZGRERUHsq8URuVf4WOkHHtGO4uH5EvnATUb43Kg+YXGE4Ka+ZGRERUXhhQ7FRBFTpSbg7idy9FzP8+gC79QRdfQeXijqBObyFi+iJoPHwM/UwKo1/7E5twCzu7h4iIyNIYUOyU2B8ilmCMg0ZO3G1EfTsWKcc2Ke51rRiOsAFz4duoI/75eBgW9W8q9zMxh4glPBiQiIjKG/eg2CmxeVXsDxFLMJAkpJzZiYTdX0LKUR7y5/tkdwT+YyDULm6G9vTia/X9TbafjZL3m5R2SYmIiMiSGFDKsXGapYnNq7O618Sw119H/B/7FNfUXv4I6Rwhn6djqj29cX8TcwIKDwYkIqLyxIBSDo3Tyqoa5uDBgxjXrx/i/1QGDL+aTeHbKQIuPkFFtqfXLxWJDbGmdpnwYEAiIrIGBpRyaJymr4YRez8sEVK0Wi2mTZuGDz74QD4xWs/FxUUejxg9Bsf+TDRrBsd4qUjcYfzaeTAgERFZCxu1lVPjNP1MxIG325r9w97UUtHdO3+hf//+8iGMxmrVqiUfMfD000/b/MwPERE5p2Q2arO9xmnG1TBi70dR+1RMBQb320dxd+vnSE1OVDz3gAEDMH/+fPnso5LiwYBERGRLGFAsxNwqF3FfUbMVeZeKdDmZSNjzFVJP7VA8lwgkixcvRt++fS3ye+DBgEREZCsYUCzE3CqXm7HpmLvrcoH7VL7o2wRTtl0wXM+OuYHYzbOQE3dLcX/tBo0xbd5XePH5phb6HRAREdkO7kGx8B6UwqphKvm5y7+KTi54tsXXQ4OUTC3Ef5aUE1uRsHcZoM1RPJNfi38hoHU/qDQu3CdCREQO+fObnWQtRF8NI+TdtaF/3KdZtULDiSDCiTY9Cfd//AgJu75UhBONTxAq9p6KwOcHyOFE4Hk5RETkiBhQLEjMYphqIy8ei/EaId5FPkfGzVOIWv4WMq4dVYx71mqGsEHz4Vm9kWKc5+UQEZEj4h4UCyusGibyWlyBXydpc5F4YBWSD/2o7EaicUVQ2yHwadIFKpXKrAohIiIie8eAUgZMVcOI2Q2dTkKApysSM4z3lAA5CVGI3TIL2VGXFeOuIdUQ0n0C3CrUMOv78rwcIiJyFAwo5cBUWbFe6rm9iP9lIaTsDMW4T5POCGwzBGpXsbHWPDwvh4iIHAUDipXa3+uy0hG/cxHSzu1VjKs9fBD8wkh4PfaMYtzXwwWpmbk8L4eIiJwCN8mWIbGsI2ZO8oaKrKjLiFoxKl84ca/aAGGDFijCiQgfopR4xksNDY+N8bwcIiJyRJxBKcf295KkQ/Lh9Uj87VtApzWMq9Ua+LXqA/8WLwNqjcnwIVcIqZvmWyoq7KRiIiIie8WAUoaMN63mpsYjbutsZP55SnGPxq8iPv58CRo91azI8MHzcoiIyFkwoJQh/abV9KtHEPfTXOgykhXXveo+i+COw/Hcs63lqh9zwgfPyyEiImfAgFKGGoZ5IWv/V7gfuVExrnJ1R1D71+HzRHuEBXgaNrcyfBARET3AgFJGLly4gN69eyP6zBnFuFulmgjpNh5uwVXkx9zcSkRElB8DioWJQ/6++uorjBo1ChkZyt4mvk/3ROBzA6ByceXmViIiokIwoFiglFi/b8RDl4klH7+N9etFu/qHKlasiGXLVyCoTjNubiUiIjIDA0opAsnN2HSsPXJLPqE48/ZZxG75DNqU+4r7O3bsiJUrV6JSpUpWe81ERET2hgGllC3rJZ0WSb+vQ9Lv34kHhnEXF1fMmPEJIiIioFazHx4REVFxMKCUomV9blIMYrd8iqw75xX3ugRWRr2+kzAqYhjUXMYhIiIqNgaUApZvjPeJmGpZn3bxAOJ3zIcuK03xHN5PtEdQ+2FIdvOUn4tlw0RERMXHgFLI8k3Y35U2/p5uhnFddiYSdi9B6plfFF+vcvOSm65513/eZCdZIiIiMh8DSiEnDkcnZcrjg1vVkB9n37uO+5tnIjf+L8V97pXrIrjbOLgGhJrsJEtERETF4/QBpaAThwUxJnaQrD/5F5KPbkLCvuWANtfoDhX8W74C/9Z9ocpzyJ/oc6LvEEtERETF4/QBJe+Jw3nlpiXi0g9zkXH9mGJc4xOMkG5j4VGtoWLc+ARi9jkhIiIqGacPKIXtE8m4cRKx2z6DLi1RMe5ZuwWCXxgJjadfvq9hh1giIqLSc/qAYmqfiKTNQeL+b5F8ZL1i3M3dA490GgZtnfZQqR7MjoT6uaNPs2qoEeLNDrFEREQW4vQBRQQKUa0jNsSKPSc58Xfk3ibZ0VcU9zVo0ABr165FvfqPmyxFJiIiIstx+oAiwoVYknn92+NIO7sHcTsXQcpRLvt06z0Q3y1bCE9PT/kxe5sQERGVLfZgF+flPB6K1il7EfvTHEU4cfH0xeTPl2Hz2uWGcEJERERlz+kDiiQ9KDCeNOq/8PLyMow3bdEaly+cQ6cu3bHp1B1EXouTS5KJiIio7Dn1Eo8+nIgNr3Xq1MG8efMwbNgwTJkyBQ07v4q+ay6Z7CzLCh0iIqKypZL0P6XtSHJyMvz9/ZGUlAQ/v/ylviUl3orLly/jRo6fyc6y+q2wi/o3ZUghIiIqw5/fTr/EY0zMpNSq/VihnWUFcZ3LPURERGWHAaWYnWVFLBHXxX1ERERUNhhQ8jD3BGKeVExERFR2GFDyMPcEYp5UTEREVHYYUAroLFtQb1gxLq7zpGIiIqKyw4BSQGdZIW9I4UnFRERE5YMBxQRRQixKicXJxMbEY5YYExERlT2nbtRWGBFCOtQP5cGAREREVsCAUggRRngwIBERUfnjEg8RERHZHAYUIiIisjkMKERERGRzGFCIiIjI5lg1oHzxxReoUaMGPDw80Lx5cxw5csSaL4eIiIicPaB89913GDNmDN5//32cOHECjRo1QseOHRETE2Otl0RERETOHlBmz56NoUOHYtCgQahfvz4WL14MLy8vLFu2zFoviYiIiJw5oGRnZ+P48eNo3779wxeiVsuPIyMj892flZWF5ORkxQcRERE5LqsElNjYWGi1WlSqVEkxLh5HR0fnu3/69Onw9/c3fFStWrUcXy0RERGVN7voJDtx4kR5v4peUlISqlWrxpkUIiIiO6L/uS1Jkm0GlJCQEGg0Gty7d08xLh6Hhobmu9/d3V3+yPsb5EwKERGR/UlJSZFXRGwuoLi5ueHJJ5/E7t270bNnT3lMp9PJj0eMGFHk11euXBm3b9+Gr68vVKqSH94ngo4IOeK5/Pz8Svw8ZB6+3+WL73f54vtdvvh+2+f7LWZORDgRP8dtdolHLNkMGDAATz31FJo1a4a5c+ciLS1NruopithQW6VKFYu9FvFm8w94+eH7Xb74fpcvvt/li++3/b3fRc2cWD2g/Pvf/8b9+/cxefJkeWNs48aNsWPHjnwbZ4mIiMj5WHWTrFjOMWdJh4iIiJyLU5/FIzbeik62xhtwqezw/S5ffL/LF9/v8sX32/Hfb5VkTq0PERERUTly6hkUIiIisk0MKERERGRzGFCIiIjI5jCgEBERkc1x6oDyxRdfoEaNGvDw8EDz5s1x5MgRa78kuycOdnz66aflLr8VK1aUOwVfunRJcU9mZiaGDx+O4OBg+Pj4oFevXvmOPaCS+eSTT+TuyhEREYYxvt+WdefOHfTv319+Pz09PfHEE0/g2LFjhuui7kD0dwoLC5Ovi1Par1y5YtXXbK/EobKTJk1CeHi4/F7WrFkTU6ZMUZzjwve7dPbv349u3brJnV3F3x0bN25UXDfn/Y2Pj0e/fv3kBm4BAQEYMmQIUlNTS/nKHnxzp7Ru3TrJzc1NWrZsmXTu3Dlp6NChUkBAgHTv3j1rvzS71rFjR2n58uXS2bNnpVOnTkmdO3eWqlWrJqWmphruef3116WqVatKu3fvlo4dOya1aNFCeuaZZ6z6uh3BkSNHpBo1akgNGzaURo0aZRjn+2058fHxUvXq1aWBAwdKhw8flq5fvy79/PPP0tWrVw33fPLJJ5K/v7+0ceNG6fTp01L37t2l8PBwKSMjw6qv3R59/PHHUnBwsLR161bpxo0b0g8//CD5+PhIn3/+ueEevt+l89NPP0nvvfeetH79epH6pA0bNiium/P+durUSWrUqJF06NAh6bfffpNq1aol9enTp5SvTJKcNqA0a9ZMGj58uOGxVquVKleuLE2fPt2qr8vRxMTEyH/o9+3bJz9OTEyUXF1d5b9o9C5cuCDfExkZacVXat9SUlKk2rVrSzt37pSef/55Q0Dh+21Zb7/9ttS6desCr+t0Oik0NFSaNWuWYUz8N3B3d5fWrl1bTq/ScXTp0kUaPHiwYuyll16S+vXrJ/+a77dl5Q0o5ry/58+fl7/u6NGjhnu2b98uqVQq6c6dO6V6PU65xJOdnY3jx4/LU1XG5/uIx5GRkVZ9bY4mKSlJ/hwUFCR/Fu97Tk6O4r2vW7cuqlWrxve+FMQSTpcuXRTvq8D327I2b94snx/28ssvy0uYTZo0wdKlSw3Xb9y4IR/dYfx+i3NHxBIy3+/ie+aZZ+RDZC9fviw/Pn36NA4cOIAXXnhBfsz3u2yZ8/6Kz2JZR/z/Qk/cL36mHj582H5b3VtLbGysvLaZ99wf8fjixYtWe12ORpxQLfZCtGrVCg0aNJDHxB92cZq1+AOd970X16j41q1bhxMnTuDo0aP5rvH9tqzr169j0aJF8mGn7777rvyejxw5Un6PxeGn+vfU1N8tfL+L75133pFP0RWhWqPRyH9vf/zxx/J+B4Hvd9ky5/0Vn0VYN+bi4iL/o7S0/w2cMqBQ+f2r/uzZs/K/eKhsiKPPR40ahZ07d8qbvansQ7f4l+K0adPkx2IGRfwZX7x4sRxQyLK+//57rF69GmvWrMHjjz+OU6dOyf/oERs6+X47Pqdc4gkJCZHTeN5KBvE4NDTUaq/LkYhDILdu3Yq9e/eiSpUqhnHx/ooltsTERMX9fO9LRizhxMTEoGnTpvK/WsTHvn37MG/ePPnX4l86fL8tR1Qy1K9fXzFWr1493Lp1S/61/j3l3y2WMX78eHkWpXfv3nK11KuvvorRo0fL1YIC3++yZc77Kz6Lv4OM5ebmypU9pf1v4JQBRUzHPvnkk/LapvG/jMTjli1bWvW12Tuxz0qEkw0bNmDPnj1yeaAx8b67uroq3ntRhiz+gud7X3zt2rXDH3/8If/LUv8h/oUvpsD1v+b7bTliuTJv2bzYH1G9enX51+LPu/hL2fj9FksUYi2e73fxpaeny3sZjIl/XIq/rwW+32XLnPdXfBb/ABL/WNITf/eL/0Zir0qpSE5cZix2Iq9YsULehfzf//5XLjOOjo629kuza2+88YZckvbrr79KUVFRho/09HRF2asoPd6zZ49c9tqyZUv5gyzDuIpH4Ptt2VJuFxcXufz1ypUr0urVqyUvLy9p1apVirJM8XfJpk2bpDNnzkg9evRg2WsJDRgwQHrkkUcMZcaiFDYkJESaMGGC4R6+36WvADx58qT8ISLB7Nmz5V//+eefZr+/osy4SZMmcun9gQMH5IpClhmX0vz58+W/uEU/FFF2LGq4qXTEH3BTH6I3ip74g/3mm29KgYGB8l/uL774ohxiqGwCCt9vy9qyZYvUoEED+R84devWlZYsWaK4LkozJ02aJFWqVEm+p127dtKlS5es9nrtWXJysvxnWfw97eHhIT366KNyz46srCzDPXy/S2fv3r0m/84W4dDc9zcuLk4OJKJHjZ+fnzRo0CA5+JSWSvxP6eZgiIiIiCzLKfegEBERkW1jQCEiIiKbw4BCRERENocBhYiIiGwOAwoRERHZHAYUIiIisjkMKERERGRzGFCIiIjI5jCgEBERkc1hQCEiIiKbw4BCRERENocBhYiIiGBr/h8tzHXi7jFY6AAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Now:\n", "# 1. create samples based on a line.\n", "# 2. Add random gaussian noise in y to simulate measurement error.\n", "# 3. Fit a line to the data.\n", "x_array = np.arange(1, number_of_samples)\n", "y_array = []\n", "for x in x_array:\n", " value = ref_m * x + ref_c\n", " value += np.random.normal(0, sigma)\n", " y_array.append(value)\n", " \n", "m, c = np.polyfit(x_array, y_array, deg=1)\n", "print(\"m=\" + str(m) + \", c=\" + str(c))\n", "plt.scatter(x_array, y_array)\n", "plt.plot(x_array, c + m * x_array, color=\"k\", lw=2.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "31bd8076-f2fa-4ba0-9a2d-dfbe64482c4a", "metadata": {}, "source": [ "# Outlier Data\n", "Try re-running, the same cell over and over. You can see that the fit of the line is reasonable, as long as the noise is random, zero-mean, and fairly small sigma (<1.0).\n", "\n", "But what if, occasionally, we get a rediculous spike in noise, i.e. representing some kind of a device or measurement failure? Try running the cell below, and you get quite vastly different results!" ] }, { "cell_type": "code", "execution_count": 107, "id": "29d75f37-5645-4856-b9cf-8c2a861c1013", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "m=1.3792878631632493, c=8.194608938811692\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQMxJREFUeJzt3Ql0VFW2x/+dhCTMYyRhBgWBCIrMiBMQQEFk7EbbAW2X/qXB14Lt+LRt225xeKtt7Vb5P/+v1fcUcUBAQFFmBJkHZRaUSUgYJcxTUv+1j97bVaHG1HSr6vtZqzqpujeVS5lO/XLOPvukuVwulwAAADhIerwvAAAAoCwCCgAAcBwCCgAAcBwCCgAAcBwCCgAAcBwCCgAAcBwCCgAAcBwCCgAAcJwKkoBKS0tl7969Uq1aNUlLS4v35QAAgCBob9hjx45J/fr1JT09PfkCioaTRo0axfsyAABAOezevVsaNmyYfAFFR06sf2D16tXjfTkAACAIR48eNQMM1vt40gUUa1pHwwkBBQCAxBJMeQZFsgAAwHEIKAAAwHEIKAAAwHEIKAAAwHEIKAAAwHEIKAAAwHEIKAAAwHEIKAAAwHESslEbAAAIX0mpS5ZvPyz7j52WutUqSudmtSUj3Rl73BFQAABIQTPXF8oz0zZKYfFp+7F6NSrK0wPy5YY29STemOIBACAFw8nId1d7hBNVVHzaPK7H442AAgBAik3rPDNto7i8HLMe0+N6XjwRUAAASCHLtx++YOTEncYSPa7nxRMBBQCAFLL/2OmInhctBBQAAFJI3WoVI3petBBQAABIIZ2b1TardXwtJtbH9bieF08EFAAAUkhGeppZSqzKhhTrvh6Pdz+UkALKG2+8IZdffrlUr17d3Lp16yaff/65ffz06dMyatQoqVOnjlStWlWGDh0q+/bt83iOXbt2Sf/+/aVy5cpSt25defjhh+X8+fOR+xcBAAC/tM/JG7e3l7wantM4el8fd0IflJAatTVs2FCef/55adGihbhcLnnnnXdk4MCBsmbNGrnssstkzJgxMmPGDPnoo4+kRo0aMnr0aBkyZIgsXrzYfH1JSYkJJ3l5efL1119LYWGh3HnnnZKZmSnPPfdctP6NAACgDA0hvfPzHNtJNs2lSSMMtWvXlpdeekmGDRsmF110kUyYMMF8rjZv3iytW7eWJUuWSNeuXc1oy0033SR79+6V3Nxcc8748ePl0UcflQMHDkhWVlZQ3/Po0aMmABUXF5uRHAAA4HyhvH+XuwZFR0MmTpwoJ06cMFM9q1atknPnzklBQYF9TqtWraRx48YmoCj92LZtWzucqL59+5oL3rBhg8/vdebMGXOO+w0AACSvkAPKunXrTH1Jdna23H///TJ58mTJz8+XoqIiMwJSs2ZNj/M1jOgxpR/dw4l13Drmy7hx40zism6NGjUK9bIBAEAyB5SWLVvK2rVrZdmyZTJy5EgZMWKEbNy4UaLp8ccfN8NB1m337t1R/X4AACDBdjPWUZLmzZubzzt06CArVqyQV155RYYPHy5nz56VI0eOeIyi6CoeLYpV+nH58uUez2et8rHO8UZHa/QGAABSQ9h9UEpLS02NiIYVXY0zZ84c+9iWLVvMsmKtUVH6UaeI9u/fb58za9YsUyij00QAAAAhj6DoVMuNN95oCl+PHTtmVuzMnz9fvvjiC1Mbcs8998jYsWPNyh4NHQ888IAJJbqCR/Xp08cEkTvuuENefPFFU3fy5JNPmt4pjJAAAIByBRQd+dC+Jdq/RAOJNm3TcNK7d29z/OWXX5b09HTToE1HVXSFzuuvv25/fUZGhkyfPt3UrmhwqVKliqlh+fOf/xzKZQAA4CglpS7H9hNJVGH3QYkH+qAAAJwSImauL5Rnpm2UwuJ/7/6re9lou3gndGRN1PfvkItkAQBINNEKEfq8I99dLWX/0i8qPm0ed0rb+ETEZoEAgKRmhQj3cOIeIvR4eUdkNPR4m4awHtPjeh5CR0ABACStaIYInS4qG3rKPr8e1/MQOgIKACBpRTNEaC1LJM+DJ2pQAABJKxohwiq23brvWFDna0EuQkdAAQAkrWDDQbDneSu29UXXB+XV+Hm1EEJHQAEAJC0NB7paRwtiXWGGCF8rdryxFi/rKiH6oZQPNSgAgKSl4UBDgkoLI0T4K7b1RkMPS4zDwwgKACCpaUjQsFB2aiYvhD4ogYptLaN7NJfuzXPoJBsBBBQAQNLTENI7P6/cnWSDLaJtkVtVul1SJ8yrhSKgAABSgoaR8oaHSBfbIjBqUAAACLLY1td4iz6ux1mxEzkEFAAAYlRsi+ARUAAACKHYVotr3bFiJzqoQQEAIEbFtggeAQUAgBgV2yJ4TPEAAADHIaAAAADHIaAAAADHIaAAAADHoUgWAOJEN6BjNQjgHQEFAOJg5vrCCzavqxfC5nVAsmOKBwDiEE5Gvrv6gt1xi4pPm8f1OJDqCCgAEONpHR05cXk5Zj2mx/U8IJURUAAghrTmpOzIiTuNJXpczwNSGQEFAGJIC2IjeR6QrAgoABBDulonkucByYqAAgAxpEuJdbWOr8XE+rge1/OAVEZAAZAStOh0yfeHZOraPeZjvIpQtc+JLiVWZUOKdV+P0w8FqY4+KACSXjg9R6LRTE2/5xu3t7/gmvLogwLY0lwuV8KtZTt69KjUqFFDiouLpXr16vG+HAAJ0HOk7C86K2JoUPAVCKLdTI1Oskg1R0N4/2YEBUBK9xx5YvI6OXWuVPKqewYEX8HGaqbmL9gES79Xt0vqhPUcQLIioABI2Z4j6vCJczLmg7UeoyO98/P8BhuNMHpcz2PEA4gOimQBJK1Qe4lYoyP/nLuVZmpAnBFQACStUHuJWCMmby3eEdT5NFMDooeAAiBle474CilHTp0L6lyaqQHRQ0ABkLT9S/z1HAmkZqVMmqkBcUSRLADHisQyX189RwK5u3sz+fvs70wYcY9ENFMDYoM+KACSrn+Jv54jRcWn5NkZm+SnE2e9rtJJ+6Vh2qJHe8qsjUVR7YMCpJqj9EEBkMz9S8qzzNe950ilrAwTfgKNjmgI0e9BMzUg9qhBAZBw/UvCXeZrTfvoSIk7vV92ZMYKNgPbNTAfCSdAbDCCAsBxLduDXb4bzjJfRkcAZyOgAAhLNParCXb5brjLfGk1DzgXUzwAwi5kLTsdY3Vk1ePR6F/CMl8g+RFQAERtIz49HkrfkmD6l7DMF0gNIQWUcePGSadOnaRatWpSt25dGTRokGzZssXjnOuvv17S0tI8bvfff7/HObt27ZL+/ftL5cqVzfM8/PDDcv78+cj8iwCkXCErgBSvQVmwYIGMGjXKhBQNFE888YT06dNHNm7cKFWqVLHPu/fee+XPf/6zfV+DiKWkpMSEk7y8PPn666+lsLBQ7rzzTsnMzJTnnnsuUv8uAFFGISsAxwSUmTNnetx/++23zQjIqlWr5Nprr/UIJBpAvPnyyy9NoJk9e7bk5uZKu3bt5Nlnn5VHH31U/vSnP0lWVlZ5/y0AYohCVgCOrUHRTnCqdm3PQrX33ntPcnJypE2bNvL444/LyZMn7WNLliyRtm3bmnBi6du3r+kut2HDBq/f58yZM+a4+w1AfFHICsCRAaW0tFQefPBB6d69uwkilt/85jfy7rvvyrx580w4+b//+z+5/fbb7eNFRUUe4URZ9/WYr9oXbY1r3Ro1alTeywYQIRSyAnBkHxStRVm/fr0sWrTI4/H77rvP/lxHSurVqye9evWS77//Xi655JJyfS8NOmPHjrXv6wgKIQWIP18b8WkhK/vVAIh5QBk9erRMnz5dFi5cKA0bNvR7bpcuXczHbdu2mYCitSnLly/3OGffvn3mo6+6lezsbHMD4DwUsgKI+xSPbnys4WTy5Mkyd+5cadasWcCvWbt2rfmoIymqW7dusm7dOtm/f799zqxZs8yuhvn5Pw8XA0gs7FcDIK4jKDqtM2HCBJk6darphWLVjGhdSKVKlcw0jh7v16+f1KlTR7799lsZM2aMWeFz+eWXm3N1WbIGkTvuuENefPFF8xxPPvmkeW5GSQAAgEpz6bBIkLTpmjdvvfWW3HXXXbJ7925TEKu1KSdOnDB1IoMHDzYBREdILDt37pSRI0fK/PnzTf+UESNGyPPPPy8VKgSXl7QGRUORriJyf14AAOBcobx/hxRQnIKAAgBA4gnl/Zu9eAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgONUiPcFIHWVlLpk+fbDsv/YaalbraJ0blZbMtLT4n1ZAAAHIKAgLmauL5Rnpm2UwuLT9mP1alSUpwfkyw1t6sX12gAA8ccUD+ISTka+u9ojnKii4tPmcT0OAEhtBBTEfFpHR05cXo5Zj+lxPQ8AkLoIKIgprTkpO3LiTmOJHtfzAACpi4CCmNKC2EieBwBITgQUxJSu1onkeQCA5ERAQUzpUmJdreNrMbE+rsf1PABA6iKgIKa0z4kuJVZlQ4p1X4/TDwUAUhsBBTGnfU7euL295NXwnMbR+/o4fVAAADRqQ1xoCOmdn0cnWQCAVwQUxI2GkW6X1In3ZQAAHIgpHgAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAAkNgBZdy4cdKpUyepVq2a1K1bVwYNGiRbtmzxOOf06dMyatQoqVOnjlStWlWGDh0q+/bt8zhn165d0r9/f6lcubJ5nocffljOnz8fmX8RAABIrYCyYMECEz6WLl0qs2bNknPnzkmfPn3kxIkT9jljxoyRadOmyUcffWTO37t3rwwZMsQ+XlJSYsLJ2bNn5euvv5Z33nlH3n77bfnjH/8Y2X8ZAABIWGkul8tV3i8+cOCAGQHRIHLttddKcXGxXHTRRTJhwgQZNmyYOWfz5s3SunVrWbJkiXTt2lU+//xzuemmm0xwyc3NNeeMHz9eHn30UfN8WVlZAb/v0aNHpUaNGub7Va9evbyXDwAAYiiU9++walD0G6jatWubj6tWrTKjKgUFBfY5rVq1ksaNG5uAovRj27Zt7XCi+vbtay56w4YN4VwOAABIEhXK+4WlpaXy4IMPSvfu3aVNmzbmsaKiIjMCUrNmTY9zNYzoMesc93BiHbeOeXPmzBlzs2iYAQAAyavcIyhai7J+/XqZOHGiRJsW5+qQkHVr1KhR1L8nAABIsIAyevRomT59usybN08aNmxoP56Xl2eKX48cOeJxvq7i0WPWOWVX9Vj3rXPKevzxx810knXbvXt3eS4bAAAkY0DReloNJ5MnT5a5c+dKs2bNPI536NBBMjMzZc6cOfZjugxZlxV369bN3NeP69atk/3799vn6IogLZbJz8/3+n2zs7PNcfcbAABIXhVCndbRFTpTp041vVCsmhGddqlUqZL5eM8998jYsWNN4awGiQceeMCEEl3Bo3RZsgaRO+64Q1588UXzHE8++aR5bg0iAAAAIS0zTktL8/r4W2+9JXfddZfdqO2hhx6S999/3xS26gqd119/3WP6ZufOnTJy5EiZP3++VKlSRUaMGCHPP/+8VKgQXF5imTEAAIknlPfvsPqgxAsBBQCAxBOzPigAAADRQEABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACJH1AWLlwoAwYMkPr160taWppMmTLF4/hdd91lHne/3XDDDR7nHD58WG677TapXr261KxZU+655x45fvx4+P8aAACQmgHlxIkTcsUVV8hrr73m8xwNJIWFhfbt/fff9ziu4WTDhg0ya9YsmT59ugk99913X/n+BQAAIOlUCPULbrzxRnPzJzs7W/Ly8rwe27Rpk8ycOVNWrFghHTt2NI/94x//kH79+sl//dd/mZEZAACQ2qJSgzJ//nypW7eutGzZUkaOHCmHDh2yjy1ZssRM61jhRBUUFEh6erosW7bM6/OdOXNGjh496nEDAADJK+IBRad3/vd//1fmzJkjL7zwgixYsMCMuJSUlJjjRUVFJry4q1ChgtSuXdsc82bcuHFSo0YN+9aoUaNIXzYAAEjkKZ5AbrnlFvvztm3byuWXXy6XXHKJGVXp1atXuZ7z8ccfl7Fjx9r3dQSFkAIAQPKK+jLjiy++WHJycmTbtm3mvtam7N+/3+Oc8+fPm5U9vupWtKZFV/y43wAAQPKKekD58ccfTQ1KvXr1zP1u3brJkSNHZNWqVfY5c+fOldLSUunSpUu0LwcAACTjFI/2K7FGQ9T27dtl7dq1poZEb88884wMHTrUjIZ8//338sgjj0jz5s2lb9++5vzWrVubOpV7771Xxo8fL+fOnZPRo0ebqSFW8AAAED8ul8ussv3oo4/khx9+kEmTJsXtWtJcejUh0FqSHj16XPD4iBEj5I033pBBgwbJmjVrzCiJBo4+ffrIs88+K7m5ufa5Op2joWTatGlm9Y4GmldffVWqVq0a1DVoDYoWyxYXFzPdAwBAGDQGrFy5Uj788EMTTHbu3Gkf05DSrFkziZRQ3r9DDihOQEABAKD89K1fSy00kGgw2bFjh9fzdDWuzoTE4/074qt4AACAM0PJ6tWr7ZESLdHwR+tCGzduLPFCQAEAIIlDyZo1a+xQolM2/nTu3Fl+/etfy7Bhw6RJkyYSTwQUAACSLJSsXbvWDiW6YMWfTp06ya9+9Stza9q0qTgFAQUAgCQIJd98840dStxX23qj281YoSSSRbCRREABACBBQ8m3335rh5KtW7f6Pb99+/Zm+kZDiTZRdToCCgAACRRK1q9fb4eSLVu2+D3/yiuvtEOJbjuTSAgoAAA4PJRs2LDBDiWbN2/2e367du3sUKKNUhMVAQUAAAfa4BZKNm3a5PfcK664wq4pufTSSyUZEFAAAHCIjRs3mlDy4YcfBgwlbdu2NYFER0tatmwpyYaAAgBAHGkQsUKJBhR/2rRpI8OHDzfBJBlDiTsCCgAAMaZ1JFabeS16DRRKfv1LTUmrVq0kVRBQAACIAV1xY4WSdevW+T03Pz/fDiX6eSoioAAAECXfffedHUq0Z4k/rVu3tkPJZZddJqmOgAIAQARpwzQNJXrTlvP+6JSNhhK9EUo8EVAAAAiTtpa3QoluzuePFrdaIyVaX5KWlhaz60wkBBQAAMpBdwa2pm9Wr17t91ztTWKFEl0eTCgJjIACAECQtm/fboeSVatW+T1Xu7haoUQbqRFKQkNAAQDAjx07dtihZOXKlX7P1f1urFCiLecJJeVHQAEAwEso+fjjj00oWbFihd9zdWdgq9CVUBI5BBQAAERk586ddihZvny533ObNWtmj5S0b9+eUBIFBBQAQMratWuXHUqWLVvm99ymTZvaoaRDhw6EkigjoAAAUsru3bvtULJ06VK/5zZp0sQOJR07diSUxBABBQCQ9H788Uc7lCxZssTvuY0bN7ZrSggl8UNAAQAkpT179phQoitwFi9eHDCU6CiJhpJOnToRShyAgAIASBp79+61Q8miRYv8ntuwYUM7lHTp0oVQ4jAEFABAQissLJRJkyaZ6RsNJS6Xy+e5DRo08Agl6enpMb1WBI+AAgBIOEVFRSaU6EjJwoUL/YaS+vXr26Gka9euhJIEQUABACRMKPnkk0/MSEmgUFKvXj0ZNmyYCSVXXXUVoSQBEVAAAI61b98+O5QsWLAgqFCioyXdu3cnlCQ4AgoAwFH2799vQolO38yfP19KS0t9npuXlydDhw6V4cOHE0qSDAEFABB3Bw4ckMmTJ5uRknnz5vkNJbm5ufb0jYaSjIyMmF4rYoOAAgCIi4MHD3qEkpKSEp/n1q1b156+ueaaawglKYCAAgCImUOHDtmhZO7cuX5DyUUXXWSmb3Sk5NprryWUpBgCCgAg6qFkypQpJpTMmTPHbyjJycnxCCUVKvA2lar4Lw8AiLjDhw97hJLz58/7PLdOnTomlOj0zfXXX08ogcFPAQAgIn766ScTSnT1zaxZswKGkiFDhphQ0qNHD0IJLsBPBACg3I4cOSJTp041IyUaSs6dO+fz3Nq1a8vgwYPN9I2GkszMzJheKxILAQUAUK5QoiMlX375pd9QUqtWLTuU9OzZk1CCoBFQAAABFRcX26Hkiy++8BtKatas6RFKsrKyYnqtSA4EFACAV0ePHpVPP/3UTN9oKDl79qzfUDJo0CBTU1JQUEAoQdgIKAAAj1Aybdo0E0pmzpzpN5TUqFHDjJQQShANBBQASHHHjh3zCCVnzpzxG0qskZLevXsTShA1BBQASNFQMn36dBNKPv/8c7+hpHr16jJw4EBTU6KhJDs7O6bXitgqKXXJ8u2HZf+x01K3WkXp3Ky2ZKSnxfgqCCgAkDKOHz8uM2bMMKHks88+k9OnT/s8t1q1anYo6dOnD6EkRcxcXyjPTNsohcX//tmoV6OiPD0gX25oUy+m10JAAYAkduLECY9QcurUKZ/nVq1a1SOUVKxYMabXiviOjMxcXygj310trjJfV1R82jz+xu3tYxpSQg4oCxculJdeeklWrVolhYWFZtMnnY+0uFwuefrpp+XNN980a+V1K+w33nhDWrRo4dEC+YEHHjBznunp6abF8SuvvGL+zwEACD+UaBjRUKLhJFAoGTBggAwfPlz69u1LKEmiaZgSt3N2HDwp7y/fJUVHvY+M6Lk6clI2nCh9TJ9Zj/fOz4vZdE+F8vzgX3HFFfLb3/7WtCku68UXX5RXX31V3nnnHWnWrJk89dRT5od+48aN9g/+bbfdZsKN1XXw7rvvlvvuu08mTJgQmX8VAKSYkydPeoQSve9LlSpV5KabbjKh5IYbbpBKlSrF9FoRXr3HzCCmYbydI35GRmpUyvJ7roYUPa7X2u2SOhILaS4d8ijvF6eleYyg6FPVr19fHnroIfnDH/5gN/fJzc2Vt99+W2655RbZtGmT5Ofny4oVK6Rjx47mHK0a79evn/z444/m64NZBqeV5PrcWrwFAKlIQ4gWuGrzNB2RDiaU6PTNjTfeSChJoHqPkjIjIX+f/d0FIx1WpNGwobxN1XijX5dXo6I8ckMrGfPB2oDnv3JLOxnYroGUVyjv3xGtQdm+fbsUFRWZ9fAWvZAuXbrIkiVLTEDRj9rQxwonSs/XqZ5ly5aZNfVlaXW5e4W5/gMBIBXpdI17KNFRbV8qV65sQokuCdY/AvU+nCHYeo+ZQYyEKOt5Hv34W8nISA8qnFhfp899+LjvVVzudJQnViIaUDScKB0xcaf3rWP6sW7dup4XUaGC2UTKOqescePGyTPPPBPJSwWAhAolOtJshRJdjeOLjoxYIyWEEmcKtt6jtFRk1ITgRkIsxad97yDtT+0qWWb0RgOSy89Ii05BxUpCrOJ5/PHHZezYsR4jKI0aNYrrNQFANOkSYCuUaLv5QKGkf//+ZqREP+p0DpxLp2uCqfd4cur6kMJJOPJqVDJTSzp6o2HE5WX6SI/Hsh9KRANKXl6e+bhv3z6pV+/f82d6v127dvY5+/fv9/i68+fPm5U91teXpevvWYMPIBVCie55Y4USbabmiy460BESLXTVj6yCTBxaSxKMwyd8bzMQKe4jIxo+dGqp7JRSXjL0QdFVOxoy5syZYwcSHe3Q2pKRI0ea+926dTPLj3WZcocOHcxjc+fOldLSUlOrAgCpROvrvvzyS7P6RncLDiaU6PSNjpQQShKLVey6dZ/v/8axlOZlZERDiC4lTshOsjrMuG3bNo/C2LVr15oaksaNG8uDDz4of/nLX0zfE2uZsa7MsVb6tG7d2ixru/fee2X8+PFmmfHo0aNNAW0wK3gAIBlCibZZsEKJv8J/HT12DyXa4RWJIVAfEn+qZGXIibMlUb0+XyMjGkZitZQ4ogFl5cqV0qNHD/u+VRsyYsQIs5T4kUceMVXl2tdER0quvvpqM4/q3vznvffeM6GkV69edqM27Z0CAMkeSnT6RkOJLrP0F0r0DzmtKbn55psJJTHsQ9KhSS1ZtfMnj9EDFYleJaE4EYVwklc9W27t3Fia5lSJ68hITPqgxAt9UAAkgrNnz9qhZMqUKX5Die4KrKFER0q0syu/26LPW4jQ9+tSt3fFmpUzzccjJ88FvTeNryXEkZL2SxHr73s1l3e+3ilHTv372sqqXSVTnrrpMsmr7oxAErc+KACQ6jSUzJ492w4lOpLsL5Rop21rpER/cSM2fIUI93BSNpgEszeNvyXE/lTNzpDjZ0pCnpppXa+6uRbxsfLmucFtY17cGikEFAAIk9bS6eIArSnR7tr+QklmZqYdSnRjPkJJ7JU3RASzN02gJcS+BBNORvdoLt2b53iMhGj4cNLKm0gioABAOUOJrkDUUKIjJdoqwV8o0d2BdfpGR0q0mzbip7whIpi9aYJdQlweLXKrei1eddLKm0gioABACKFk3rx59kiJv1CiHbJ79+5tQomOlNSqVSum15pqQtl0L5IhouxzhdoKXq+wVpVMOXzCdx1JMM/tlJU3kURAAQA/tJHk/PnzTSj55JNP5NChQwFDiTV9o+0XEJ9iV38rViK5n0zZ59Lv469lvDsrPv1lYBt5dsYmR7WZdwICCgB4CSULFiwwha6TJk2SgwcP+g0l2jJBR0q03xOhJPqC2d236OgZeXn2Vq8rb0IJEb74Cg0agny1jC/LvU4kPT3NUW3mnYBlxgDwSyhZuHChPVJy4MABn+dmZGSYXdgJJc5veGax3trddwn2tvqlPM8ViVEdX19TLwmKXcv7/k1AAZCySkpKTCixRkrK7hNWNpT07NnTTN8MHjxYcnJyYnqtqSrchmfeRj0WPdrTBINo9UEpT11MOF+TSAgoAOAnlHz11Vf2SIluZuqLdrrWUKIjJYSS2ItWw7P37+1qF5RGq5MsvKNRGwCUCSWLFi0yIyUff/xxwFBy/fXXm1AyZMgQueiii2J6ranOCgxFxadM4Wg0/oJ2X3njbfWLt9UwybZCJhEQUAAkbShZvHixHUqKior8hpLrrrvODiV169aN6bUmo2CmKsqe89OJs/LsjMhM5/gTyVU8iB4CCoCkUVpaKl9//bWZvtFQUlhYGFQo0emb3NzcmF5rqhWylq3biGRtSbBSdbluoiKgAEiKUGKNlOzdu9fnuWlpaSaUaKGrjpTk5eXF9FqTVTBhQ5f03v/uahlT0EKKT52Tfy3eEdNrTOXluomKgAIgIUPJkiVL7FCyZ88ev6HkmmuuMaFk2LBhhJIIrxgJtpDVOu7emyRSrN4hGn50Ca+3EZxk2Jsm1RBQACRMKFm2bJmZvtFgEiiUXH311Wb6ZujQoVKvHm9KZfnqufFU/9ZSq0p2UKEl3E33IsVb+BjdszkrbxIcAQWAY2kXBCuU6EjJ7t27/Z7vHkrq168fs+tMNL5GPTSs/G7CGo/H/DUXi8Sme+URqOFZsu5Nk2oIKAAcF0qWL19uh5Jdu3b5Pb979+729E2DBg1idp2JOH2jln5/SB6btC7oUQ9/LeOjuXNvWbWrZMpTN10medUZDUkVBBQAjgglK1assKdvAoWSq666yg4lDRs2lFQXzHJdbx1Ry0OLXXX0Rdu8x2K5rhVDnhvclvqRFENAARC3ULJy5Uo7lOzcudPv+V27djXTNxpKGjVqFLPrdLpgl+uGG0ws1sjLE5PXyX/2y5faVbJMIIpWHQrFramLgAIgpqFk1apVJpBoMNmxY0fAUGKNlDRu3Dhm15koIyazNhbFfLmu5fCJc/LQR99E9DlDLdJFciOgAIh6KFmzZo09UvLDDz/4Pb9Lly52KGnSpImkmkBLf+PR4CzcQlZd9vv32d+Zx72NtNzTvakU5OcRRuCBgAIgKqFk7dq1dij5/vvv/Z7fqVMne/qmadOmkqq8hQ9vb/TxXtZbnkLWlnlVvS5rZvoGvrCbMYCI0F8l33zzjT19s23bNr/nd+zY0Q4lzZo1k1QXrZ17Y819p+BoNIZDYmM3YwAxCyXffvutPVKydav/LqEdOnQw0zd6u/jii2N2nU4X74ZnNStlymu3tZfik+fC3qzP39JjepMgFAQUACGHknXr1plQordAoeTKK680IyUaSi655BJJBpEeCYhXwzPrip8f2la6N88xn/dtk+d3079A2CkYkUJAARBUKFm/fr09UrJly5agQolO3zRv3lxSoUV82VqKUEJMLBqeeeuD4m0Jb9lRDqtlfFHxKXl2xiafS4rZKRiRRkAB4NOGDRvskZLNmzf7Pbddu3b29E2LFi0kWVfQeKsTcW9epm/2wYaYaI06+Fquq0Id+XEPLJWyMsy/09qcz8JOwYgGimQBXBBKrELXTZs2+T23bdu2Mnz4cBNKLr30UklkgUKFhperX5jrcyrGGkF4qn++jJpwYYix3ratEOPOem4NOuH8Qo7Fct1QwxdQ3vdvAgoAE0SskZKNGzcGDCVWTUnLli0lGfgaGXEPFTUqZcmtby4N+FzaWfXwibMBC1K7XlzH6+iMCvWXcqwDAqtxUF4EFAAB6ZSNNVKi9SX+tGnTxg4lrVq1kmQS7MjIIze0kjEfrI3Y9/UWKoJpwmZNr4wpaOF3N1/AiVhmDMArLW61QomuxAkUSqyaktatW0uyCrSCRsOAHj98/ExEv2/ZuhWlH3vn+19Fw940SBUEFCDJfffdd/bqG+1Z4k9+fr49UqKfp4JgV9Do1I2OeviqE9Hxi1pVMs0eNcFw/fI1OmKiocQaAfG1iobpFKQaAgqQhLQ3iRVKtLurPzplYxW6XnbZZZJqgl1Bk1ejkhm58LeK5S8D25iluMEWu1qjMxpAfDUwo7kZUhUBBUgS2lreCiW6D44/WtyqIyV601CSlpZ8f5EHW8ipjwcaGbH6e+jX65TMBfvluE27pKeneQ0x/sSiDwqQaAgoQALTTfismhLdMdgfXQZsTd/oSpxkDCXlWQqroSPQyIh7f4+ydSIafjo0qSWrdv4kU9fuMfdf+037kFrG030VuBCreIAE88MPP9ihZPXqn5el+qIN0zSQaDC5/PLLkzqUhLJk2FuBaXn7e/j6Om2UpkuTtSfKkVPe61Ks0ZlFj/akrgQp4SjLjIHksn37dhNK9LZy5Uq/52preWuk5Iorrki6UOJv6ibYJcO+AkGo/T2CCUPKW3+TQIEJSEYsMwaSwI4dO+xQsmLFCr/n6iZ8VijRlvPxDiXRauQVaJQj2CXDvopSQylI9bcDsfsKHQ1DgepWAFyIgAI4yM6dO+Xjjz820zfLly/3e+7FF19sT9/o5nzxDiXRboUezD44Z86XxqwoNZQw5K1uheXCgH8EFCDOdu3aZYeSZcuW+T23WbNm9khJ+/btHRNKQt1ML9SRl0CjFeqJyevk1s6NY1aUGmzIsc5juTAQGgIKEAe7d++2Q8nSpf73d2nSpIkdSjp27Oi4UBLqlId7U7JgR1602DTQihhtkPbavO/9nuO+ZDhcwYYcVugA5UNAAWLkxx9/tEPJkiVL/J7buHFju0+Jk0NJeaY8Xp71nXRvnuMxxRFo5OW33ZuGfX3elgyHI5T+KQBCR0ABomjPnj12KPn6668DhhKrpqRTp04JEUrKM+Xxz3nbzC2veraZkmlcu7Lpvupv5GXy2j1hX1+ki1JD7Z8CIDQEFCAKoWTSpEkmlCxevNjvuY0aNZJhw4aZVvOdO3dOuFASzlRG0dEz8vLsrQHPc/0yfaN74fx04mzQ3Vkto3s0v2DEJlI07LBCB4gOAgoQAXv37vUIJf7aCzVs2NAeKdFQkp6eLsmwFDjQlEe4BrWrL28t3hFSC3nVIrdqVItTWaEDJEhA+dOf/iTPPPPMBft+bN682Xx++vRpeeihh2TixIly5swZ6du3r7z++uuSm5sb6UsBoqqwsNAOJYsWLfIbSho0aGBCid66du0at1AS6lLgUEKMvymPSNAQoN+/7HU7oUiVFTpAgoyg6OZjs2fP/vc3qfDvbzNmzBiZMWOGaT6l3eRGjx4tQ4YMCTgUDjhBUVGRCSX687tw4UK/oaR+/fpm+kZHSrp16xb3UBLqUmB/LdxrVcn2Glp8TXmEo+xmfdZoRVHxKVO74mvahyJVILFFJaBoIMnLy7vgcW1t+z//8z8yYcIE6dmzp3nsrbfektatW5ullvqXJeA0+/btk08++cSMlCxYsMBvKKlXr549UnLVVVc5JpRYgu0nsmz7YTOdUpaGjt9N8NyUsOzIi/uUx+JtB+SfAZb++uOt2NR9tKJSVgZFqkCSikpA2bp1q/nrsWLFiuYvx3HjxpkVCqtWrZJz585JQUGBfW6rVq3MMV126Sug6FSQ3tx7+QPRtH//fo9QUlrqu0OphnEdKdFQcvXVVzsulLiHk7cXbw+qn4i3cOKLtyZsVojQ0YtJq/eUuy4lULEpRapA8op4QOnSpYu8/fbbpu5E5+i1HuWaa66R9evXm+HxrKwsqVmzpsfXaP2JHvNFA07ZuhYg0g4cOGCHkvnz5wcMJUOHDjXTN927d5eMjAyJl2DqRLxN10SKvyZs5alLqV0lU5666TLJqx5csSlFqkByivpuxkeOHDGdMP/2t79JpUqV5O677/YYDVG6kqFHjx7ywgsvBD2Cossz2c0YkQglkydPNqFk3rx5fkOJBmkrlOhISTxDSSjFrr5qTqLh/Xu7ei0WDSYgsbsvkPyOOmk3Yx0tufTSS2Xbtm3Su3dvOXv2rAkt7qMoOsfvrWbFkp2dbW5AJBw8eNAjlJSUlPg8t27dunYo0ZFAJ4SSYIpd7393tYwpaOG3CVosm7WVHeXYcfCkvL98lxQdZVoGQJwCyvHjx+X777+XO+64Qzp06CCZmZkyZ84c80tfbdmyxWyWprUqQLQcOnTIDiVz5871G0ouuugis7JMQ8m1117rsQotkYpdg2mCFmn+lvSWXYo7umdzpmUA+BTx37x/+MMfZMCAAWZaR5tXPf300+avzltvvdUM69xzzz0yduxYqV27thneeeCBB0w4YQUPohFKpkyZYkKJhmJ/oSQnJ8cOJdddd13cQkmwfUcC7XsTa+VZ0kvvEAD+VIjGhmgaRvTNQf8S1bl6XUKsn6uXX37ZrHLQERT3Rm1AJBw+fNiEEu1Tor14zp8/7/PcOnXqmFCiq2+0BireIyXBNk8LZd+bYFWrWEGOnz5frqkglvQCSMgi2XgX2SD5/fTTTzJ16lQzUjJr1iy/oURH7jQcxzqUBBoZCVTIek/3plLwSydVpcuFtbYkUiMfT/XPl1ETVpvHXH6+vzZFe3ZGcCEKAMJ5/yagICFpobV7KNH+Ov5CiftIidZBxVKgjqxWR9TDJ84GfK6alX++9iMnff97g1V21Uw02t8DgDsCCpI2lHz66acmlHz55Zd+Q0mtWrVk8ODBJpT06tUr5qHEEsslvqEifABI6WXGQDj0h9g9lOgydV906bqGEi10jWUo8famrpZ+f0gem7TOUeEkUBM0ClcBOAUBBY5M2FYo+eKLL/yGEk3i7qFEOxXHkrdpkUhOw4RqWPsGprW8+Nib5rnBbakVAZAQCChwTCiZNm2aWX0zc+bMC7oNlw0lgwYNMtM32vwv1qEk0PRNPIKJ5ZpLL5KC/Fz2pgGQ8AgoiJtjx47J9OnTzUjJ559/7jeU6FzlzTffLMOHDzehJBqdhUOpv/DXKC2e9Lp1ioa9aQAkOgIK4hJKdKTks88+8xtKqlWrJgMHDjTTN3369Inqdgeh9CBJhEZp1JIASHQEFEgstjuYMWOGGSnRUHL69Gm/oURHSqxQUrGi79bpsdzTpmlOFY+RiEg3SitbyDqwXQN5a/GOoHYAplEagGREQEFUnDhxwiOUnDp1yue5VatWtUOJdhaORSgp7542Vv+Sg8d8j/yUV9lC1i7NagdVgEt9CYBkRB8URDSUaBjRUKLhJFAo0T2brFBSqVKluCwFjlRH1mDVrJQpr93WXopPnguqI6uv66a+BEAiog8KYubkyZMeoUTv+1KlShU7lNxwww1RDSXub+w7Dp6U95fvkqKj8VsKbMWH54e2le7Nc8znfdsELmT1VUtCfQmAZEdAQcg0hOiqGy101aXB/kJJ5cqV7VBy4403+g0lkepi6q3gtaxoBpNgp2EoZAUA3wgoCIpO17iHEp3O8RdKbrrpJjuU6P1w96sJNrTEu7W8Xu9d3ZuZz5mGAYDyI6DAbyjRpmlWKNHVOL7oyEj//v1NKOnXr5+ZzgmWr1ChYeV3E9Z4POavTsPadC+eRVU51bLtIMLoCACUHwEFHnQJsLaX15oSbTcfKJRoGNFQouEklFBS3oZnuvRXw4y/HXjjSUdLAADhI6DAhBLdiM8KJdpMzRddAqyhRNvM6zSOrsYJR6gNz6wgo5vwbS48Jq/M2RqVERNfe9oE2ygNABAeAkqK0g6uViiZOnWq31CiHVy1lkTbzAcTSkJZGlvehmdHTp2Tv8/5d3+SSLGCxgvDrvC6p42vr1E0SgOAyCGgpFgomTVrlh1KdD16oFCi0zcaSrTDayR397VqSZw0JVI2aOgUUtk9bX46cfaC/iU0SgOAyKNRW5I7e/asRyjR18xfKOl8TU/p2ONG6ddvgPS4vElIq1FCWUFjPcNrv7nSFLZqbUksfxD9haZAQSNSy6EBINUcDeH9m4CSpKFk9uzZZvXNlClT5MiRIz7PzcrKMk3TWnbrI/NONZIDZzLK9Saub9pXvzA3pHoSazrlqf75MmrCavNYtH4Y86pny62dG3vsqaMIGgAQO3SSTUHnzp2TOXPmmJGSyZMnBwwl2l5ep2+0idqS3Se9jnx4a2ZWdhVNOLv76vfTr6lVJcs8X6RX4+ime0/ddJnkVfcdPlgKDADOREBJ8FAyd+5cO5T89NNPPs/NzMw0oURX3+jGfDVr1nRb5rsy6JEL91U01SpmSteL64S9u69+re7eG0y9R3k23QMAJB4CSgKGknnz5tmh5PDhw35DSUFBwc+rbwbcLN/9VGre/DcdKpEOVUtl1c6fZPG2A+UatdBVNLf9f8siUuxqfa231u+6X83S7w+ZKSD9nsGgaBUAEh8BJQGcP3/eI5QcOnTI57kVKlSQgoLe0qFHP2ndtYdc3CDPjETc9P+u9ggiOttRGoGCD2vKR4tdNayEUuwaTO8QDS3dW+SYTfb0+yhXmefQ+2MKWnjUl1BLAgCJjSJZB4eS+fPn26Hk4MGDAUJJgZm+qdayq/xtYWHMO6tqvYdO07y1eIcdGvyx4kPZWpby7NfDaAkAJAZW8SRwKFmwYIFZfTNp0iS/oSQjI8OEElPoevNA2VYsMmtjkfxr8Q6Jt7KjM+Es6S2LJb4AkLhYxZNgoWThwoVmpOSTTz6RAwcO+A0lvXr1MqFk0KBBUqdOHTOqcPObax2zF42ywsk93ZtKQX5eRJf0eqtTAQAkHwJKHJSUlJhQYo2U7N+/328o6dmzpx1KatWuY97oF+0+LTvWbJW/z/4uar1DvI18BEujx2fri+SJ/v9u/06wAAAEi4ASw1Dy1Vdf2aFk3759Ps9NT083oWTYsF9J4/bXy9nMn4s/l+89K8++GVoztPIY3aO5dG+eY498hLqKxr3HiYYpggkAIFQElCiHksWLF5vpm48//thvKJG0dLmiUzfpO2Cw5F9VIEdKK8tby3dJ0aTIb4gXaFXNmN6Xeky/+FtFE0g4/VEAAKmLgBJhpaWldijRkZLCwkLfJ6elS3ajNlKl1dVS+dKr5EiVmvLBcRH58sdYXvLPlxJgR14tZi1Pt1cnbQYIAEgcBJQw6aqSpd8flIWLFsmaBTNl8azpsnfvXt9fkJb2Syi5Ripf2k0yqtQSJ6y0Caa5mfvuvkXFp8wmf9pjxVXOHicAAPhCQAljpOTl96bJf43/Xznw7QIpOe67eVpaWppce+210rFHf/ngQD3JqBrbUKJLep/q31pqVcm2V9F0aFLLdJINdVWN+yqaSlkZZtqnbN+TQKMxAAAEQkAJMZQsXbrUTN+8+/6Hcmh/od9Qcs0115jVN0OGDJF69erJ1LV75OOJa2N2ve7LfKOxUZ6vaR9azQMAwkVACSKULFu2zC50/fFHf/UhaZLdMN/UlDTt2EPm/vXXHsEgVvUYseyu6j7tQ/M0AECkEFC80Oa6GkomfvCBTPzgI9lXuMfv+RpKKre8Wiq3vEoqVMsxjx1y/dyYzH2UQt+4Q92vJhAn7EVD8zQAQKQRUNxCyZKly+Sf//N/MnvGVDlQFCCUNMiXyq08Q0lZn6//eQrIvZPqjW3yTDv6YParCQbTKQCAZMRePL+Ek/94+kX557OP+T0vu0FrM1Lyq2FD5Yud58PqyFqe3YS9FbsynQIASBTsxROiLzYUyZSDeV6PZdVv+fOS4JbdpUL1i34+P4Rw4qtVfNn9anS57rMzyhSbVs+WWzs3jtvUDQAA8ZLyAUX7mOgqlIzqOWaE5MyeTZJVT0PJ1VK5lYaSulH73mX3q+nbhmJTAABUygcUDQTWqEWtgv9HMipVkwo1cmPyvcvuV0OxKQAAP0v5gOK+V0x2XvO4XwMAABBJlxTnhL1inHANAAA4ScoHFKs3SSQqPWpXyQrpefRc/d7sVwMAgKeUDyha96F9RFR5Q4oVNP4ysE3Qz8N+NQAAODSgvPbaa9K0aVOpWLGidOnSRZYvXx6X67D2lNGmZ+50ma92aH3llnYypuBSEyrS/ASNfpd7fx7tg2L1QrGfu0ZFcy4N1gAAcFCjtg8++EDuvPNOGT9+vAknf//73+Wjjz6SLVu2SN26dWPaqM19ybG/Zb4z1xdesDGet31vvD2PYgkxACCVHQ3h/TtuAUVDSadOneSf//ynvSlfo0aN5IEHHpDHHnssLgElEiEGAAAkaCfZs2fPyqpVq+Txxx+3H0tPT5eCggJZsmTJBeefOXPG3Nz/gfFCrxIAAJK0BuXgwYNSUlIiubmeDdH0flFR0QXnjxs3ziQu66YjLQAAIHklxCoeHWnR4SDrtnv37nhfEgAAiKK4TPHk5ORIRkaG7Nu3z+NxvZ+Xd+GmfdnZ2eYGAABSQ1xGULKysqRDhw4yZ84c+zEtktX73bp1i8clAQAAB4nbXjxjx46VESNGSMeOHaVz585mmfGJEyfk7rvvjtclAQCAVA8ow4cPlwMHDsgf//hHUxjbrl07mTlz5gWFswAAIPXErQ9KOOLZBwUAAET//TshVvEAAIDUQkABAACOE7calHBYs1Lx7CgLAABCY71vB1NdkpAB5dixY+YjHWUBAEjM93GtRUm6IlntmbJ3716pVq2apKWlhZXkNORoZ1qKbaOP1zu2eL1ji9c7tni9E/P11sih4aR+/fpmD76kG0HRf1TDhg0j9nz6YvMDHju83rHF6x1bvN6xxeudeK93oJETC0WyAADAcQgoAADAcVI6oOgGhE8//TQbEcYIr3ds8XrHFq93bPF6J//rnZBFsgAAILml9AgKAABwJgIKAABwHAIKAABwHAIKAABwnJQOKK+99po0bdpUKlasKF26dJHly5fH+5IS3rhx46RTp06my2/dunVl0KBBsmXLFo9zTp8+LaNGjZI6depI1apVZejQobJv3764XXMyef7550135QcffNB+jNc7svbs2SO33367eT0rVaokbdu2lZUrV9rHdd3BH//4R6lXr545XlBQIFu3bo3rNSeqkpISeeqpp6RZs2bmtbzkkkvk2Wef9djHhdc7PAsXLpQBAwaYzq76u2PKlCkex4N5fQ8fPiy33XabaeBWs2ZNueeee+T48eNhXtnP3zwlTZw40ZWVleX617/+5dqwYYPr3nvvddWsWdO1b9++eF9aQuvbt6/rrbfecq1fv961du1aV79+/VyNGzd2HT9+3D7n/vvvdzVq1Mg1Z84c18qVK11du3Z1XXXVVXG97mSwfPlyV9OmTV2XX3656/e//739OK935Bw+fNjVpEkT11133eVatmyZ64cffnB98cUXrm3bttnnPP/8864aNWq4pkyZ4vrmm29cN998s6tZs2auU6dOxfXaE9Ff//pXV506dVzTp093bd++3fXRRx+5qlat6nrllVfsc3i9w/PZZ5+5/vM//9P1ySefaOpzTZ482eN4MK/vDTfc4LriiitcS5cudX311Veu5s2bu2699dYwr8zlStmA0rlzZ9eoUaPs+yUlJa769eu7xo0bF9frSjb79+83P/QLFiww948cOeLKzMw0v2gsmzZtMucsWbIkjlea2I4dO+Zq0aKFa9asWa7rrrvODii83pH16KOPuq6++mqfx0tLS115eXmul156yX5M/xtkZ2e73n///RhdZfLo37+/67e//a3HY0OGDHHddttt5nNe78gqG1CCeX03btxovm7FihX2OZ9//rkrLS3NtWfPnrCuJyWneM6ePSurVq0yQ1Xu+/vo/SVLlsT12pJNcXGx+Vi7dm3zUV/3c+fOebz2rVq1ksaNG/Pah0GncPr37+/xuipe78j69NNPpWPHjvKrX/3KTGFeeeWV8uabb9rHt2/fLkVFRR6vt+47olPIvN6hu+qqq2TOnDny3XffmfvffPONLFq0SG688UZzn9c7uoJ5ffWjTuvo/y8ser6+py5btiys75+QmwWG6+DBg2ZuMzc31+Nxvb958+a4XVey0V2ntRaie/fu0qZNG/OY/rBnZWWZH+iyr70eQ+gmTpwoq1evlhUrVlxwjNc7sn744Qd54403ZOzYsfLEE0+Y1/w//uM/zGs8YsQI+zX19ruF1zt0jz32mNlFV0N1RkaG+b3917/+1dQ7KF7v6Arm9dWPGtbdVahQwfxRGu5/g5QMKIjdX/Xr1683f/EgOnTr89///vcya9YsU+yN6Idu/UvxueeeM/d1BEV/xsePH28CCiLrww8/lPfee08mTJggl112maxdu9b80aMFnbzeyS8lp3hycnJMGi+7kkHv5+Xlxe26ksno0aNl+vTpMm/ePGnYsKH9uL6+OsV25MgRj/N57ctHp3D2798v7du3N3+16G3BggXy6quvms/1Lx1e78jRlQz5+fkej7Vu3Vp27dplPrdeU363RMbDDz9sRlFuueUWs1rqjjvukDFjxpjVgorXO7qCeX31o/4Ocnf+/Hmzsifc/wYpGVB0OLZDhw5mbtP9LyO9361bt7heW6LTOisNJ5MnT5a5c+ea5YHu9HXPzMz0eO11GbL+gue1D12vXr1k3bp15i9L66Z/4esQuPU5r3fk6HRl2WXzWh/RpEkT87n+vOsvZffXW6codC6e1zt0J0+eNLUM7vSPS/19rXi9oyuY11c/6h9A+seSRX/3638jrVUJiyuFlxlrJfLbb79tqpDvu+8+s8y4qKgo3peW0EaOHGmWpM2fP99VWFho306ePOmx7FWXHs+dO9cse+3WrZu5ITLcV/EoXu/ILuWuUKGCWf66detW13vvveeqXLmy69133/VYlqm/S6ZOner69ttvXQMHDmTZazmNGDHC1aBBA3uZsS6FzcnJcT3yyCP2Obze4a8AXLNmjblpJPjb3/5mPt+5c2fQr68uM77yyivN0vtFixaZFYUsMw7TP/7xD/OLW/uh6LJjXcON8OgPuLeb9kax6A/27373O1etWrXML/fBgwebEIPoBBRe78iaNm2aq02bNuYPnFatWrn++7//2+O4Ls186qmnXLm5ueacXr16ubZs2RK3601kR48eNT/L+nu6YsWKrosvvtj07Dhz5ox9Dq93eObNm+f1d7aGw2Bf30OHDplAoj1qqlev7rr77rtN8AlXmv5PeGMwAAAAkZWSNSgAAMDZCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAMBxCCgAAECc5v8HSU1sgwvWeFYAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Now:\n", "# 1. create samples based on a line.\n", "# 2. Add random gaussian noise in y to simulate measurement error.\n", "# 3. With a certain probability (fail_probability), add errors that are (fail_magnitude) times worse.\n", "# 4. Fit a line to the data.\n", "fail_probability = 0.10\n", "fail_magnitude = 100\n", "\n", "x_array = np.arange(1, number_of_samples)\n", "y_array = []\n", "for x in x_array:\n", " value = ref_m * x + ref_c\n", " \n", " noise = np.random.normal(0, sigma)\n", " likelihood = np.random.uniform(0, 1)\n", " \n", " value += noise \n", " if likelihood < fail_probability:\n", " # Deliberately adding large positive numbers.\n", " spike = sigma * fail_magnitude\n", " value += spike\n", " y_array.append(value)\n", " \n", "m, c = np.polyfit(x_array, y_array, deg=1)\n", "print(\"m=\" + str(m) + \", c=\" + str(c))\n", "plt.scatter(x_array, y_array)\n", "plt.plot(x_array, c + m * x_array, color=\"k\", lw=2.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8e29d08c-b824-489d-9f2a-20e958fd3260", "metadata": {}, "source": [ "# The RANSAC algorithm in English\n", "\n", "1. Pick the minimum number of datapoints to solve your problem (e.g. fit a line, compute 3D/2D pose)\n", "2. Fit the model (e.g. for a line, that means, compute slope and intercept)\n", "3. For a given distance threshold, count how many data samples 'agree with' that model, i.e. it's like they are voting for that model, i.e. a consensus. In terms of a line, it just means, 'they are close to the line'.\n", "4. If the number of votes is greater than another threshold, then compute the model using that consensus set.\n", "5. If the number of votes is less than this threshold, go to 1.\n", "6. If after a certain number of iterations, no solution has been found, compute the model with the largest set so far, or fail." ] }, { "cell_type": "code", "execution_count": 105, "id": "7ad67a03-8c5e-4f0b-8051-19425f8428a4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "m=1.0000658856550042, c=5.3643822445749, votes=88\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPlhJREFUeJzt3Qt0VOW5//FncptAIAkBcoGEi6BACGpFhRS1VUFUtN7aU61atCxdUvFUsNbSI3osp8XantNqj8o5Xa3av6Kn/k/RilX/FAS8RG6KmnAREIGQG7dcuOQ2s//reZMZMmGS7CSTZM/M97PWdJjZO5mdkTK/vO/7PK/LsixLAAAAHCSmry8AAACgNQIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwnDgJQ16vV0pKSmTgwIHicrn6+nIAAIAN2hu2pqZGhg0bJjExMZEXUDSc5OTk9PVlAACALti/f79kZ2dHXkDRkRPfD5icnNzXlwMAAGyorq42Awy+z/GICyi+aR0NJwQUAADCi53lGSySBQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjhOWjdoAAEDP8Hgt2bDniFTU1Er6wES5cHSaxMb0/r53BBQAAGC8XVgqj72xVUqrapueEJGslER59NpcuTIvS3oTUzwAAEA0nMx98eOAcKLKqmrN83q8NxFQAACIch6vZUZOrCDHfM/pcT2vtxBQAACIchv2HDlt5KQljSV6XM/rLQQUAACiXEVNbUjPCwUCCgAAUS59YGJIzwsFAgoAAFHuwtFpplqnrWJifV6P63m9hYACAEAU8HgtKdh9WF7fcsDct1zwqn1OtJRYtQ4pvsd6vDf7odAHBQCACPe2jf4mev/sbeeddl5mH/VBcVmW1Xs1QyFSXV0tKSkpUlVVJcnJyX19OQAAOL6/idXqed9YiIaSluGjJzvJdubzmxEUAACitL+Jq7m/yYzcTH8I0fv8MYOlr3VqDcqzzz4rZ599tkk9esvPz5e33nrLf7y2tlbuvfdeGTx4sAwYMEBuuukmKS8vD/ge+/btk1mzZkn//v0lPT1dHnzwQWlsbAzdTwQAABzb36RHAkp2drY8/vjjsnnzZtm0aZNcdtllct1110lRUZE5Pn/+fHnjjTfk1VdflbVr10pJSYnceOON/q/3eDwmnNTX18uHH34oL7zwgjz//PPyyCOPhP4nAwAgylU4sL9Jr61BSUtLk1//+tfy7W9/W4YOHSrLli0zf1bbt2+XCRMmSEFBgUydOtWMtlxzzTUmuGRkZJhzli5dKg899JAcPHhQEhISbL0ma1AAAOiYVuvc8oePOjzv5bum9sq0Tmc+v7tcZqyjIa+88oocP37cTPXoqEpDQ4NMnz7df8748eNlxIgRJqAovZ80aZI/nKiZM2eaC/aNwgRTV1dnzml5AwAA4dffxK5OB5TPP//crC9xu91yzz33yPLlyyU3N1fKysrMCEhqamrA+RpG9JjS+5bhxHfcd6wtS5YsMYnLd8vJyensZQMAEHViHdjfpMcCyrhx42TLli2yfv16mTt3rsyePVu2bt0qPWnhwoVmOMh3279/f4++HgAAvdkorSdd2dzfRPuZtKSPW5cYO0mny4x1lGTs2LHmz5MnT5aNGzfKk08+Kd/97nfN4tfKysqAURSt4snMzDR/1vsNGzYEfD9flY/vnGB0tEZvAABEaqO0nnRlXpYpJe6p/iaObHXv9XrNGhENK/Hx8bJq1Sr/sR07dpiyYl2jovRep4gqKir856xcudIslNFpIgAAIrVRWuty37KqWvO8Hu8Nsc39Ta47d7i5d3I46fQIik61XHXVVWbha01NjanYWbNmjbzzzjtmbcicOXNkwYIFprJHQ8d9991nQolW8KgrrrjCBJHbb79dnnjiCbPu5OGHHza9UxghAQBEmq40SkMXAoqOfHz/+9+X0tJSE0i0aZuGkxkzZpjjv/3tbyUmJsY0aNNRFa3QeeaZZ/xfHxsbKytWrDBrVzS4JCUlmTUsP//5zztzGQAARFyjNCd0b3US9uIBAKCH6ILYH72ypcPznrz5XDP1Eumqe6MPCgAAaJ8uRg3ledGEgAIAQA8J50ZpfY2AAgBADwnnRml9jYACAEAPCtdGaWHXqA0AAER+o7S+RkABAEREvxGnf/j7GqXBHgIKACCs9XUbefQM1qAAAMJ2szyntJFH6DGCAgAIy1EO2shHNkZQAAA9qqdGOTrTRh7hh4ACAOgxHY1yKD3elekeXRAbyvPgLAQUAECP6clRDtrIRzYCCgCgx/TkKAdt5CMbAQUA0GN6cpSDNvKRjYACAOgxPT3KQRv5yEWZMQCgx/hGObRaR8OI1QOjHLSRj0wuy7K63ymnl1VXV0tKSopUVVVJcnJyX18OAKADdHtFZz+/GUEBAPQ4RjnQWQQUAECvYLM8dAaLZAEAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOMQUAAAgOPE9fUFAIBTeLyWbNhzRCpqaiV9YKJcODpNYmNcfX1ZQFQioACAiLxdWCqPvbFVSqtq/c9lpSTKo9fmypV5WX16bUA0YooHQNTTcDL3xY8Dwokqq6o1z+txAL2LgAJAon1aR0dOrCDHfM/pcT0PQO8hoACIarrmpPXISUsaS/S4ngeg9xBQAEQ1XRAbyvMAhAYBBUBU02qdUJ4HIDQIKACimpYSa7VOW8XE+rwe1/MA9B4CCoCopn1OtJRYtQ4pvsd6nH4oQO8ioACIetrn5NnbzpPMlMBpHH2sz9MHBeh9NGoDgOaQMiM3k06yQDiOoCxZskQuuOACGThwoKSnp8v1118vO3bsCDjnm9/8prhcroDbPffcE3DOvn37ZNasWdK/f3/zfR588EFpbGwMzU8EoNdpj5CC3Yfl9S0HzH249gzRMJI/ZrBcd+5wc084AcJkBGXt2rVy7733mpCigeJnP/uZXHHFFbJ161ZJSkryn3fXXXfJz3/+c/9jDSI+Ho/HhJPMzEz58MMPpbS0VL7//e9LfHy8/PKXvwzVzwWgl9AiHkBPcFmW1eVfdQ4ePGhGQDS4XHLJJf4RlHPPPVd+97vfBf2at956S6655hopKSmRjIwM89zSpUvloYceMt8vISGhw9etrq6WlJQUqaqqkuTk5K5ePoAQtYhv/Y+Ib9yB9RsAuvr53a1FsvoCKi0tsPzupZdekiFDhkheXp4sXLhQTpw44T9WUFAgkyZN8ocTNXPmTHPRRUVFQV+nrq7OHG95A9C3aBEPwJGLZL1er9x///0ybdo0E0R8vve978nIkSNl2LBh8tlnn5mREV2n8te//tUcLysrCwgnyvdYj7W19uWxxx7r6qUC6OMW8bqeAwB6JaDoWpTCwkJ5//33A56/++67/X/WkZKsrCy5/PLLZffu3TJmzJguvZaOwixYsMD/WEdQcnJyunrpAEKAFvEAelKXpnjmzZsnK1askHfffVeys7PbPXfKlCnmfteuXeZeF8eWl5cHnON7rMeCcbvdZq6q5Q1A36JFPADHBBRdT6vhZPny5bJ69WoZPXp0h1+zZcsWc68jKSo/P18+//xzqaio8J+zcuVKEzpyc5u6OQJwPlrEA3BMQNFpnRdffFGWLVtmeqHomhG9nTx50hzXaZzFixfL5s2b5auvvpK//e1vpoRYK3zOPvtsc46WJWsQuf322+XTTz+Vd955Rx5++GHzvXWkBEB4oEU8AMeUGWvTtWCee+45ueOOO2T//v1y2223mbUpx48fN+tEbrjhBhNAWk7L7N27V+bOnStr1qwx/VNmz54tjz/+uMTF2VsSQ5kx4Bz0QQFgV2c+v7vVB6WvEFAAZ9FSYlrEAwjl5zd78QAIWYt4AAgVdjMGAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOQ0ABAACOw148CHtsVAcAkYeAgrD2dmGpPPbGVimtqvU/l5WSKI9emytX5mX16bUBALqOKR6EdTiZ++LHAeFElVXVmuf1OAAgPBFQELbTOjpyYgU55ntOj+t5AIDwQ0BBt2kIKNh9WF7fcsDc90Yo0DUnrUdOWtIr0ON6HgAg/LAGBWG5BkQXxIbyPACAszCCgrBcA6LVOqE8DwDgLAQUhOUaEC0l1pGatoqJ9Xk9rucBAMIPAQVhuQZE+5zoNJJqHVJ8j/U4/VAAIDwRUBC2a0B0jcuzt50nmSmB0zj6WJ+nDwoAhC8WySKs14BoCJmRm0knWQCIMAQUdGsNiC6IDbbKxNU8ktEba0A0jOSPGdzjrwMA6D1M8aBLWAMCAOhJBBR0GWtAAAA9hSkedAtrQAAAPYGAgm5jDQgAINSY4gEAAI5DQAEAAI7DFE8f0RbwrNsAACA4AoqDdgBeNGuCDEpyE1oAAFGPgNJHOwC3bm6mYeWHyz4JeE5Di/YSoVwXABBtWIPikB2Ag9EurRpmNNQAABBNCCgO2gG4NV+Q0VCj4QYAgGhBQOlFXdnZV2OJhhoNNwAARAsCSi/qzs6+XQk3AACEKwJKH+wA7OrlcAMAQLghoDhkB+C26HkaajTcAAAQLQgoDtkBOBhfiNFQQz8UAEA0oQ+KQ3YAPnq8Xha/Gdi8TUMMfVAAANGIgOKgHYBn5gWGFjrJAgCiFQHF4aEFAIBoxBoUAADgOAQUAAAQ3gFlyZIlcsEFF8jAgQMlPT1drr/+etmxY0fAObW1tXLvvffK4MGDZcCAAXLTTTdJeXl5wDn79u2TWbNmSf/+/c33efDBB6WxsTE0PxEAAIiugLJ27VoTPj766CNZuXKlNDQ0yBVXXCHHjx/3nzN//nx544035NVXXzXnl5SUyI033ug/7vF4TDipr6+XDz/8UF544QV5/vnn5ZFHHgntTwYAAMKWy7KsLu9Cd/DgQTMCokHkkksukaqqKhk6dKgsW7ZMvv3tb5tztm/fLhMmTJCCggKZOnWqvPXWW3LNNdeY4JKRkWHOWbp0qTz00EPm+yUkJHT4utXV1ZKSkmJeLzk5uauXDwAAelFnPr+7tQZFX0ClpTV1Od28ebMZVZk+fbr/nPHjx8uIESNMQFF6P2nSJH84UTNnzjQXXVRUFPR16urqzPGWNwAAELm6HFC8Xq/cf//9Mm3aNMnLyzPPlZWVmRGQ1NTUgHM1jOgx3zktw4nvuO9YW2tfNHH5bjk5OV29bAAAEMkBRdeiFBYWyiuvvCI9beHChWa0xnfbv39/j78mAAAIs0Zt8+bNkxUrVsi6deskOzvb/3xmZqZZ/FpZWRkwiqJVPHrMd86GDRsCvp+vysd3Tmtut9vcAABAdOjUCIqup9Vwsnz5clm9erWMHj064PjkyZMlPj5eVq1a5X9Oy5C1rDg/P9881vvPP/9cKioq/OdoRZAulsnNbdrpFwAARLe4zk7raIXO66+/bnqh+NaM6LqQfv36mfs5c+bIggULzMJZDR333XefCSVawaO0LFmDyO233y5PPPGE+R4PP/yw+d6MkgAAgE6XGbtcwTeue+655+SOO+7wN2p74IEH5OWXXzbVN1qh88wzzwRM3+zdu1fmzp0ra9askaSkJJk9e7Y8/vjjEhdnLy9RZgwAQPjpzOd3t/qg9BUCCgAA4afX+qAAAAD0BAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAAAI/4Cybt06ufbaa2XYsGHicrnktddeCzh+xx13mOdb3q688sqAc44cOSK33nqrJCcnS2pqqsyZM0eOHTvW/Z8GAABEZ0A5fvy4nHPOOfL000+3eY4GktLSUv/t5ZdfDjiu4aSoqEhWrlwpK1asMKHn7rvv7tpPAAAAIk5cZ7/gqquuMrf2uN1uyczMDHps27Zt8vbbb8vGjRvl/PPPN8/9/ve/l6uvvlp+85vfmJEZAAAQ3XpkDcqaNWskPT1dxo0bJ3PnzpXDhw/7jxUUFJhpHV84UdOnT5eYmBhZv3590O9XV1cn1dXVATcAABC5Qh5QdHrnz3/+s6xatUp+9atfydq1a82Ii8fjMcfLyspMeGkpLi5O0tLSzLFglixZIikpKf5bTk5OqC8bAACE8xRPR26++Wb/nydNmiRnn322jBkzxoyqXH755V36ngsXLpQFCxb4H+sICiEFAIDI1eNlxmeccYYMGTJEdu3aZR7r2pSKioqAcxobG01lT1vrVnRNi1b8tLwBAIDI1eMBpbi42KxBycrKMo/z8/OlsrJSNm/e7D9n9erV4vV6ZcqUKT19OQAAIBKneLRfiW80RO3Zs0e2bNli1pDo7bHHHpObbrrJjIbs3r1bfvKTn8jYsWNl5syZ5vwJEyaYdSp33XWXLF26VBoaGmTevHlmaogKHgAAoFyWZVmdeSt0Lcmll1562vOzZ8+WZ599Vq6//nr55JNPzCiJBo4rrrhCFi9eLBkZGf5zdTpHQ8kbb7xhqnc00Dz11FMyYMAAW9ega1B0sWxVVRXTPQAAhInOfH53OqA4AQEFAIDw05nPb/biAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjkNAAQAAjhPX1xcAAAC6x+O1ZMOeI1JRUyvpAxPlwtFpEhvj6rGv6w0EFAAAwtjbhaXy2BtbpbSq1v9cVkqiLJo1QQYludsMH2193aPX5sqVeVnS11yWZVkSZqqrqyUlJUWqqqokOTm5ry8HAIA+8XZhqcx98WOx80HeMny09XW++PLkd3JluByRyZMn99nnNwEFAIAw5PFactGvVgeMgLTHFz6e/t7XZPGb2wK+zlt3QuoObJPa4iKp218o9aU7JUa8UllZKQMGDOiTz2+meAAACEMb9hyxHU6U1RxSHn69UA4ePCR1xUVSu7/I3NeXfylieQPO94hIQUGBzJgxQ8KiimfdunVy7bXXyrBhw8Tlcslrr70WcFwHZB555BHJysqSfv36yfTp02Xnzp0B5xw5ckRuvfVWk55SU1Nlzpw5cuzYse7/NAAARImKGvvhpLH6kBzfukYOvfOfUvjkXVL8+1vl4PJfSs2m16W+bNdp4aTlZ35f6fQIyvHjx+Wcc86RH/zgB3LjjTeedvyJJ56Qp556Sl544QUZPXq0LFq0SGbOnClbt26VxMREc46Gk9LSUlm5cqU0NDTInXfeKXfffbcsW7YsND8VAAARLn1g02dqazpQ0Hi0xD86olM2jVXlNr+rS+LTR0liTp48dMcNctc/zZK+0q01KDqCsnz5crn++uvNY/1WOrLywAMPyI9//GPznM4zZWRkyPPPPy8333yzbNu2TXJzc2Xjxo1y/vnnm3Pefvttufrqq6W4uNh8fUdYgwIAiETByn5VsFJg/xqUyhNSf3Cv1O4vlLr9RVJbXCje45ViS0ysJGSOlcTsieIeMUkSh0+Q2MQBkpmSKO8/dFnIS477bA3Knj17pKyszEzr+OiFTJkyxcxjaUDRe53W8YUTpefHxMTI+vXr5YYbbgjlJQEAEBaClf2m9o8395UnGvzPZQyIk++OqpcvCzfLkdVrZP9nm8Rbd9zWa7ji3DJwxAS55orLZFXlYHFnjRNXwqmRGF8c0Wqfvu6HEtKAouFE6YhJS/rYd0zv09PTAy8iLk7S0tL857RWV1dnbi0TGAAATtaZJmhtlf1WnmgQb0Ot1Jd80TRCUlwk+0q2y4aGU5+J7XG5k8yoiDsnTxJzJoo7c6wsnT3FX2rcOhBlOqgPSlhU8SxZskQee+yxvr4MAABs6UwTNA0yeq4vnOhoSF3x1lNrSEp3ingbbb1uUkqaXDB1mow75wL58ES6VPXLEldMbNDX1/sZuZnR0Uk2MzPT3JeXl5sqHh99fO655/rPqaioCPi6xsZGU9nj+/rWFi5cKAsWLAgYQcnJyQnlpQMAEBJtjYaUVdWa55+97Tx/SNBw8tSKjbJrwyqzmFVDSUPFnuai4I7FJqc3jYxkT5R+OXmSPXqM/OOnl/vXqHQUPvRx/pjB4kQhDShataMhY9WqVf5AomFC15bMnTvXPM7PzzeNXzZv3uzvULd69Wrxer1mrUowbrfb3AAAcCJfGCirOmmaoFnt9CH52f95VyomeuTlv70ja9auk9pD+22/TlxatqmwcedMNMEkLjlwyURZdZ25Dg0dTg4fPRJQtF/Jrl27AhbGbtmyxawhGTFihNx///3yb//2b3LmmWf6y4y1MsdX6TNhwgS58sor5a677pKlS5eaMuN58+aZBbR2KngAAHD6dE5Aye+RYv90ja4j+ar6oMy2841dMZKQPtqMjphQkp0rsUmpIe2PElEBZdOmTXLppZf6H/umXmbPnm1KiX/yk5+YXina10RHSi666CJTRuzrgaJeeuklE0ouv/xyU71z0003md4pAACE83SO5fVIw8GvmgKJTtkUbxXvCbslv3HizjqraXQkW6dtJkiMOylk/VHCDXvxAADQxWmdr//iHdm7o6m6xhdIrPoTtr7eFZ8o7mHj/dM1CVnjJCa+68sZXM1VOD3RvyRU2IsHAAAb7CwkbXnOwFiP1JfukFff+H+y9r335MuiLWI12iv51dEQDSOpo88WK3OCJGSMEVds+x/DqUH6oATjpP4loUJAAQBEZRg5erxeFr/Zfinwqx9sk0X//Vcp3f6x2enX7Fvj1W30OhabNOhU/5GcPEkYMkJcrhi5f/qZ8tt/BO5RF8yiWRPkjmmjzZ87um4n9S8JFaZ4AAARPerR1od60O997KgJImMa98n+bZul5MsvbJf8xqVmNi9obQokcalZZkuY1uFHe49oi3otO7a6OFXj6UQTOCdhigcAELUN0NqrqmlJfz/3VFe02MOmSBqPHDDHDtm4lvjBI8Q9Iq95QauW/A457Zy0pHhZdM1EyUwODBF6vbq4Vh9ZXZiqiQ3zEmI7GEEBAIR1AzTfx7g2QFPBzlH6cddweL+/3FdDiafGThRpLvnNOKNpysYEklyJ7Z/S9uktrqmtaZfOdJuNFJ35/CagAAAcz79zbxujIhoIMpK1AsYlZdW1/pLf+oo9/h1+9d570uZebrFa8juuOZDkinu4lvz2t329doNGuE7VdBVTPACAPhXqD179Xu1N2ehv2qVHjkld2RdNgUT7kBzQkt+T9kt+h0/wd2nVfiSuuIROX+e8S8fKtLFDbP+80TBV01UEFABASIV66kLDzge7Tp+K8daflLoD25v7jxRJXckOEU/75bg+MYkDTzVE0x4kWvLbvKleV/gWts6fcVZEj4D0JgIKAKBPNsoLpr3qG8/JGrPLb1MgKZT6st06j2PrumIHpDWX/Da1jI9vLvkNhUjsQeIEBBQAQEhouNCRk/Y2ytPjWmYb7IO89chLY83h5gWtTfvYaAt5u+IGZZ3aw0ZLflMyAkp+u0Mv3dvih4zEHiROQEABAISErXUiVbX+3XZbjpbsOXhcfvO/75nqGl+FTWNlqe3Xjh86SmZe9g3JO3+qvLinv8QNHNxh9xKddtJmaIOS3LJya5n86YOvgpb96uP508+UUUOSzKjO5JGDZPPeo1GzsLWvEFAAAN3W1jqRYPSD/e+fHZCFf3xLird/0hxICsVz7Ij9kt/Msc3rR/JkZO7XZPF3p/pHML7RxhoYXxgJFiw0MOnj1l/X1ugIC1t7HmXGAIBuVePYaYxmSn7Ld5uRkeF1X8nWjzeIt7bG3kXExot72Dj/dI3+OSahnzk079IxMn/GuA47ydod5Yi2st/eRpkxAKDXurYGWxRrNdZLXekX/umaupLt/pLfox28riuhX4uS3zxxZ54prrimTfNamzZ2aNAA0dXyXcp+nYOAAgA4bQTBtyZDglTj3PPix2ZNxoi0/rL4zW0mnHjrTkjdgW1N5b4aSkq/EPE02nrNmH7JzSW/TT1IEtJHd1jy6yvr1REORCYCCgDA/lRN8/1v/rYpoMKmvvxL+yW/A4c0BRLTpTVP4gZnd6rChrLe6EBAAYAo1XK9xVeHTsjv/vFFm5UvjdWHpK5YK2x0hKRIGg7vs/06Q4ePkm9deblccsklMmDkJJn/9wPdKvmlrDc6EFAAIAq1N1qitRONR0v8oyM6ZdNYVW7zO7skPn1Uc0O0pk6tf5l/lX9dh4aiX394xEwV2anQ6Kj6BpGLgAIAUb6+xLK8pgmab3REQ4nneEdLWZvFxDaV/PoDSa7EJA5oc52IBgsd/dCFtXZ6jhBGohcBBQCibMTE8jRKfdmuUwtai7eKt+64re+jG+glDBsvibqhnvYhGTZeYhISTz+vnXUiOjWjLe/t9hxBdCKgAICDBevLoTrq1dHy674oPiz//tKbpwKJlvw21Nl6fZc7SRKHT2ja4VerbLLGiis2eMlvSx2FDX1eW97TcwRtIaAAQBitE0nt3xQOKk80tNmb5H8Ldsi//Nf/Sun2T5rWkJTuFPHaLPntn+ofHdFpm/ihI23v8puWFC+Lrpkomcn2wgY9R9AeAgoAOFBbDdBaBhOf4pIyuX3RKhnr2S/F2z6W4l3bWq3uaFtscnpAIIlLGx60wmbOtFGS3C/BVPpIkLUj6pc3TGJ6BiFDQAGAMNoVWDVWVfina3Rha+ORYvO8nZ1w4tKymzu0aiCZKHHJ6e2e33p0ZlzmANaOoFcQUACgj7S170vLXYFNye+RA00t401jtELxVB/sRMnv6KaGaKbKJldik1LbObvjKhrWjqC3EFAAoA/CyNHj9bL4zdP3ufnZlWfJ6oLNUr3praYREq2wOVFp70Vi4sSddZa4c3Kb2sZnT5AYd5Lta7Q7EsLaEfQGAgoAdJKdHW9bd2l9ecM+KasO0hTN09BU8ru/SCr2F8p1/7pVrPoTtq7DFZ9oynx9beMTss6SmHh3p38eXV8yPTeTkRA4CgEFALpZWZOZ7JZbLhzhnxYJNjri422olbqSHf71I/UlO8RqtFfyq6MhvnJfXT+SkDFGXLFd/2e89foSwEkIKADQzcqasuo6+e0/dgb9Gk/tMdMIrWm6psiMlojXY+v1YpMGidusH2mqsjElv66YLl07XVoRbggoABCCyhr/eceOnqqwKS6ShoqvbJf8xqVmNpf7TjTBJC4167SS32B9UOyg0gbhhoACAO3wrSX5YNfB06ZstMLGU13RVGGzv8gEEq24sSt+yIhTgSQ7T+KSh7R7/rxLx8j8GePMn1vuqcOeNohEBBQAsLnexJT8Hi6W2uJCfyjx1NjpPqKpIcasGdFSX3/Jb/+UTl3PtLFD/SFDq2j0psGDviSIRAQUAFFbadNROLnnzxulrmJP8+hIUyDxnqy29w1iteR3XFOFjW6qN1xLfvt36WcJtiuwD31JEKkIKACiosdI60qbYJvunTMsST7evEnWrlsnv37hNan+qlCs+pP2S36HT/B3adV+JLrzb3e1tyuwD31JEIkIKAAiruw3mNaVNrrY1FN3Uip2fe7v0FpfqiW/9hafxiQOPDU6oj1IMs6wvamer8R30awJMijJ3X6wYroGUYqAAiDiyn7b4jlZ4y/5LfWV/FpeW18bOyCtueRX149MlPghOZ0q+Q02ghNsRGRmHtM1gCKgAAjLNSJ2yn4baw43j440lf02HNpr+/vHDcry7/BrSn5TMk4r+Z136ViJj405rUtssNERuz8f0zVAEwIKAEdMy3S2q2nLDfX8FTaVZf7+I7qgtbGy1PY1xQ8ZGdAULT0zy0y5WO0sWp0/4ywTKOZdNpZRDyDECCgAenV0pM1urFW1cs+LH5/Wu0O1/F6TRw6SzXuPyt8/OyD1B78yUzZNJb+F4jl2xH7Jb+ZY//oR3VQvtl9ywCnXnztMnmujx0jrRauMegChR0AB0GujI1oO29a0jO+51gtZfV1TLa9H6su/lPriQjlppmyKxFtbY++iYuPFPWycf/2Ie/h4iUno1+6X6LXSYwToOwQUACHT0ejIt88b3mG1jY/VWC9lO5o7tOoIScl2+yW/Cf3EPVwbojWPkGSeKa64+E73HNGREXqMAH2DgAIgJNpbtOp77v9+3HYbeG/dCak7sM2/fqSudIeIp9HWa8f0Sw4s+U0f3amSXx+mbwDnIKAACInWi1Y74jlRdWr9iJb8ln9pv+R34JCmQKKLWnUPm8HZARU2ratovjp04rRKm2Cb7jF9AzgHAQVASGgQaE9j9SF/QzQdIWk4vM/2944bNKyp5HdEUx+S2OT000p+1ffzR8pVeVlBp2GCVdoopm8AZyKgAAgJ/YAPKPk9WtLUf8RM2RRKY1W5ze/kkvihI/39R3TaJnbAIFtfqeGkremYtqZqmL4BoiSg/Ou//qs89thjAc+NGzdOtm/fbv5cW1srDzzwgLzyyitSV1cnM2fOlGeeeUYyMjJCfSkAeonX65XEmmJxbX1HKnZuaSr5PX7U3hfHxDaV/PoqbHSX38QBIdtMD0B46pERlIkTJ8o//vGPUy8Sd+pl5s+fL2+++aa8+uqrkpKSIvPmzZMbb7xRPvjgg564FCAq2enS2p1Org0NDfLxxx/LunXrzE3//3v0qL1AohvoJQwb72+I5h42Xv5p6hj53+YFtHbb1vu/n43N9ACEnx4JKBpIMjMzT3u+qqpK/vjHP8qyZcvksssuM88999xzMmHCBPnoo49k6tSpPXE5QFRpqw9JR4tG2+vkeuLECVm/fr289957JpAUFBSY5+xwJfSXxOxcf5dWHS1xxcYHjHz86tvnyPTcjNOuO9hCVs0g3hYphoWtQGTqkYCyc+dOGTZsmCQmJkp+fr4sWbJERowYIZs3bza/eU2fPt1/7vjx480x/QevrYCiU0F686muru6JywYitg+Jfuj/cNkn7X6t9irRr332tvMkP6e/GRVZu26dvLXyXdn22SfSaHOX36FDh8rFF18iIyZOlhG554kMGilPvbvbHGuvI6sGjGA9R9rqJMvCViCyhTygTJkyRZ5//nmz7qS0tNSsR7n44oulsLBQysrKJCEhQVJTUwO+Rtef6LG2aMBpva4FQOc3z2vza49XNvcfKZQb/7xN6sq/NOtK7EhIzZCLLrpIbv7WTPP/df3/fusKmwnDU2x1ZLW7kJWFrUDkC3lAueqqq/x/Pvvss01gGTlypPzlL3+Rfv3aby3dloULF8qCBQsCRlBycnJCcr2Ak7ReF9KZ0YLO9CFprK7w7/Cr941Him1fY1xadnOFjfYhmSjxyenypYjk5J8n48cHn2Zpa3SEkQ8AfVZmrKMlZ511luzatUtmzJgh9fX1UllZGTCKUl5eHnTNio/b7TY3INrWjrReb9HeOpG2+pCYkt8jB/wN0TSQeKorbF6VSxIyzjCVNf6S36TAEVC9PI0Zeu0aQtoKHXRkBeCogHLs2DHZvXu33H777TJ58mSJj4+XVatWyU033WSO79ixQ/bt22fWqgDRqq21Iy3DSet1Iq1Diq8PiW6q13Bwr3+H39rireI9UWnvQmLixJ11lrhzck2HVt3lN8ad1OGX6WVqsNIREkIIAEcGlB//+Mdy7bXXmmmdkpISefTRRyU2NlZuueUWU1Y8Z84cM12TlpYmycnJct9995lwQgUPolVn1o74zvnZ8s/lZINXMpMTZVJWkix7c4289946Ofr2KqnZVyRW3XFbr+2Kd5tdft05k5oqbLLGSUy8u8e6yQJAnwWU4uJiE0YOHz5sVvPr4jktIdY/q9/+9rcSExNjRlBaNmoDonW9yQe7DnZqDxtvfa0c+GqL3PnO82a6pr5kh1iNp6rc2qOjIbp2xLSN1031MsaIKzauR7rJAkB3uCydoA4zukhWR2O0r4qOwgCRsN6kLZ7aY2ZTvabpmiKpL9sl4vXYep3YpEH+/iMaSrSFvMsVI6Hm62fy/kOXsfAVQEg+v9mLB3DIehMfz7Gj/pJfXUei60ns9leNS81sHh2ZKBlnnSuLvz9DslL6ydHj9bL4zVZlvsluueXCETIirb8sfnObOacrv63QyRVATyCgAN3UmZbxrdeb6ACmVtT4dvjVYKIVN3bFDxnhn67R+7jkIf5jugpFw4lv0erMvLbLfPslxJrQpI/aCynBOrvSyRVATyCgAN0IH3bayrfsZ/L+zgrZt3un1BY3jY7U7d8qnpqD9i7GFWPWjGjJb1MgyZXY/im2F622V+ar4UIrg05rptY8yjJqSFKbnV3pZwKgJ7AGBWhDW+HDN1rQ0VSNr+S3vmKP1BcXycnmURLvSZtbNcRqye+45v4jueIeriW//Tv1M7x819ROlf12ZwNBAOgIa1CAbmorfPj6kDz9va+ZdRutj1uNDVJXttPfobXuwFax6k/aek1XfKIJITo6csu3rpB/uvpS+fpZWWbkpazqZKfWifgWrfpGPOyimRoApyCgAJ3oS+J77qG/fiY1tR7x1p+UugPbmzu0Fkp96RdiNdbbep2YxIFN7eJ1dCS7aZffmJhYEyx+v+BUNYwvMNhdJ8KiVQCRgICCqBBs6kIF2/emvb4knpM1puT3aMuSX8vepnqxA9KaFrSOmGTWj+gC15Ylvx0Fi7bWibRuh8+iVQCRgICCiA8kXx06IS9v2Cdl1bXtVqO0/qBXjTWHzeiIbw+bhoNf2b6GuEFZpypscvIkLiXjtF1+W7ITLIJtuteZDQUBIFwQUBCVTdBaBpOWoaaxqtzff0QXtDZWltp+3fghIwOaosUNtLeWY96lY2Xa2CG2g0WwdSKsGwEQaQgoiCh2Kmt8LMsrDYf2+9ePaDDxHDtiv+Q3c0zThnqm5HeCxPbrXEWZbyHr/BlnMeIBAK0QUBA1m+6Zkt/y3f4dfk3Jb22NvW8eG2821fM1RHMPHy8xCf26fK0sZAWA9hFQ4Gid6cuh57Wc1tFqmrrSL/zTNXUHtonVYG9TPldCP3EP14ZoE5s218s8S1xxTetW2qMN2rR7a+u28ixkBYDOIaAgbBultba37JCc/HKzfx8bDSfiabT1WjH9kptLfnXKZqIkpI8WV0ys//gPvzlGXtm4v80+JL7pmjumjTYBqnVbeRayAkDnEFAQlo3StNz2/Mx4ef/992XdunXm9sknn4jXa7Pkd+CQpkDSPGUTPzgnaIWNL3g8cMU4OTs7JWgfkmDTNSxkBYDuIaCgU3qjFXpba0kaqw81lfzuL5Qbn9sqJyt0l1974tKGS6KuHWkOJbHJ6e2W/AYLHm3uV8N0DQCEHAEFPTbl0lUagEoqT0rj0RJ//xENJVoCbI9L4tNH+UdHNJjEDhgUcIadPijBgkewPiRM1wBA6LFZILo15eL7WNaRhfZCSkcjLzo18+lnn8uLy9+Wt1e9Kzu2bBDP8aO2ri02Lk4uOP98ufjii+Ub3/iGnBw0Rv597YFO78rLOhEAcM7nNwEFHdJwcdGvVrfZ+My3TuP9h07tH9PRyEvGgDi55YxGaTiwVd577z15d806OVZTZet6XHEJktBc8mtGSbLGyX/94OsBAYldeQHAedjNGCHVuny3NU24elzP04WgrdvM/+4fX4inoU7qS3c0T9cUyb6SbbKhoc7W67vcSZI4fIK/S6tuqueKPVXyq7FDA5BOvbS3SBUAED4IKFGiOyMK+jV2z/ONlhyo0D1stkltcXMPktKdIl6bJb/9U5t2+G0eIYkfOjKg5LejgAQACH8ElCjQ3cWtGmg64jleKa8vXy7/982VZmFrfcUe+7v8Jg89taBVN9VLG95hhU13ghQAwPkIKFG6uFXDyj0vfixzpo2S6bmZbY6o6MiL12tJar94qTx5quKlsbrCX12jIyQNR4rlOZvXFJeW3dyhNc9U2MSlpEso2AlSAIDwQECJ4r1p1B8/+MrcdERF27QPSnL7p4G0a6q2bDclv0cONLWMby779VRX2N9UL3100/41plNrrsQmBZb82pGWlNBhF1dfZQ4AIPwRUCJ4vckHuw62u7i1JT3vh8s+8W+q13Bwr3+HX91Yz3ui0t6Lx8SJO+vM5jAyUdzZuRLjTuryz+ILH4tm5cq9y+x1cQUAhD8CShSsN+mI5WmQ+rJd/imbWt1Ur+64ra/t37+/5Ofny+i88+WN8hRJGHaWxMSHZqqlZfgwXVxj6OIKANGCgBIB1ThKn1u5tUz+9MFXHX4Pb32t1JVs90/X1JfsEKvRXsmvjoaY3X2zJ8qv590ss791qcTHx5vrKvzVarNXTqga67QOH3RxBYDoQUAJ89GRYC3bW/PUHpO64q3NgaTQjJaI12PrNXW9iH/9iK/k1xVjjg0dM9GEE3NejMuEibY209PH86efaTq5+nqjqGBhpr2Fu/Q3AYDoQEAJ82qcYMHEc+yo1Oqmes2BpKFCR1XsjWvEpWQ0j5A0NUWLGzSszZLf1lUzndlMb1zmgF7Z1wcAEJ4IKGFejaM7FWhFTdOC1iITTLTixq74wSP8oyM6UhKXPKTDr2mvasbuNAzTNQCA9hBQwqzVvAaSxsPFTSMkuqBVS35rDtr7Rq4YGTLyLKkb0rSPjVbYxPZP6dS12KmasTsNw3QNAKAtBBSH83g8UrBho1RvfL1pyqa4SLwn7G2qJ7Fa8nuWaYjWL3ui5Ew4Vz569FqzmLazlT4+VM0AAHoDAaWbQr1rbl1dnWzatEnWrVtndvn94IMPzO6PdrjiE8VtNtVrnrLJOsvs/Ou7msXfOc9cW7DpFV9TttZrQlo3b2MaBgDQGwgoPbDHTWc+1I8dOyYFBQUmjGgoWb9+vdTW2hvZiEkcaKZpfG3jE9LPEFdsnK1Rj2DTKzPzWBMCAHAGl6WLGsKMjiikpKRIVVWVJCcnO6qqJpiWoWV3camUbN8iZV98Ih+8/75s3rzZTOPYETsgzb+hno6SxA8Z4S/5DaajfXYAAHDq5zcjKD20x41PY81h2bWtSL776m/NotaGQ3ttv05calZzIMkVd84kiUvNPK3kN1gfFMp1AQDhjoDShXDy/Ad7gi4wNRU2lWX+/iNa9ttYWWr7e+fl5cnFF18sK4+kycnBZ0rcwOAlv2lJ8bLomomSmRzYSZapGQBApCCgtLPgdfLIQbJ579E2F5JallcaDu1vEUgKxXPsiP1dfjPHyOAx58jTC26VSy65WAYPHiwFuw/L3//wUbv/YY4cbzDhpOUaEsp1AQCRhIDSzoJXHYTwtpjH0V1+68u/bN7hV/uQFIm3tsbeC8TGi3vYuKYdfnUNyfDxEpPQzxzKPHuqCSdKw5Adds8DACAcEVDaWfDqaaiXutIvmjq06ghJyXax6k/a+p6uhH7iHu6rsJko7kwt+W1aL9Je2GjdPr4tds8DACAcRX1Aabng1Vt3QuoObPOPjtSV7hDxNNr6PjH9kpv6j2Q3VdgkpI8WV0xsp8OGrh/RRa5t7QrcXpt5AAAiRdQHFF8b+RNffCgHX3tc53FsfV3swCEBe9jED85pc1O9tgQLGx3tCtxRm3kAACJB1AcU3/RK/NBR7YaTuLThzetHmkJJbHL6aYFEe51kpfQ7rSNrMO2Fjc7sCgwAQCSK+oDim17RniPaCK2pCscl8UNHNjdEyzPBJHbAoA5HQu6YNtqEjdYdWYO1ke8obLDbLwAgmkV9QGm55iP1m3dKjDupaZffxAG2vj7YSEio2siz2y8AIFq13Se9Fzz99NMyatQoSUxMlClTpsiGDRt6/Rp8az7UwImXSv+xF9oOJ76REJ2O6WjaxRc2rjt3uLlnJAQAAAeOoPzP//yPLFiwQJYuXWrCye9+9zuZOXOm7NixQ9LT03v1Wtpa89G6Dwq7+wIAEOGbBWooueCCC+Q///M/zWOv1ys5OTly3333yU9/+tM+2Sywo06yhBEAACJ4s8D6+nqzi+/ChQv9z8XExMj06dOloKBA+kqwNR+sAQEAoPf1SUA5dOiQeDweycjICHheH2/fvv208+vq6sytZQIDAACRq08Xydq1ZMkSMyTku+lUEAAAiFx9ElCGDBkisbGxUl5eHvC8Ps7MzDztfJ0K0vkq323//v29eLUAACAqAkpCQoJMnjxZVq1a5X9OF8nq4/z8/NPOd7vdZjFNyxsAAIhcfVZmrCXGs2fPlvPPP18uvPBCU2Z8/PhxufPOO/vqkgAAQLQHlO9+97ty8OBBeeSRR6SsrEzOPfdcefvtt09bOAsAAKJPn/VB6Y6e6oMCAACc8fkdFlU8AAAguhBQAACA44Tlbsa+WSkatgEAED58n9t2VpeEZUCpqakx9zRsAwAgPD/HdS1KxC2S1Z4pJSUlMnDgQHG5XN1KchpytPEbi217Hu937+L97l28372L9zs832+NHBpOhg0bZvbgi7gRFP2hsrOzQ/b9aP7Wu3i/exfvd+/i/e5dvN/h9353NHLiwyJZAADgOAQUAADgOFEdUHSPn0cffdTco+fxfvcu3u/exfvdu3i/I//9DstFsgAAILJF9QgKAABwJgIKAABwHAIKAABwHAIKAABwnKgOKE8//bSMGjVKEhMTZcqUKbJhw4a+vqSwt2TJErngggtMl9/09HS5/vrrZceOHQHn1NbWyr333iuDBw+WAQMGyE033STl5eV9ds2R5PHHHzfdle+//37/c7zfoXXgwAG57bbbzPvZr18/mTRpkmzatMl/XOsOHnnkEcnKyjLHp0+fLjt37uzTaw5XHo9HFi1aJKNHjzbv5ZgxY2Tx4sUB+7jwfnfPunXr5NprrzWdXfXfjtdeey3guJ3398iRI3LrrbeaBm6pqakyZ84cOXbsWDevrOnFo9Irr7xiJSQkWH/605+soqIi66677rJSU1Ot8vLyvr60sDZz5kzrueeeswoLC60tW7ZYV199tTVixAjr2LFj/nPuueceKycnx1q1apW1adMma+rUqdbXv/71Pr3uSLBhwwZr1KhR1tlnn2396Ec/8j/P+x06R44csUaOHGndcccd1vr1660vv/zSeuedd6xdu3b5z3n88cetlJQU67XXXrM+/fRT61vf+pY1evRo6+TJk3167eHoF7/4hTV48GBrxYoV1p49e6xXX33VGjBggPXkk0/6z+H97p6///3v1r/8y79Yf/3rXzX1WcuXLw84buf9vfLKK61zzjnH+uijj6z33nvPGjt2rHXLLbd088osK2oDyoUXXmjde++9/scej8caNmyYtWTJkj69rkhTUVFh/tKvXbvWPK6srLTi4+PNPzQ+27ZtM+cUFBT04ZWGt5qaGuvMM8+0Vq5caX3jG9/wBxTe79B66KGHrIsuuqjN416v18rMzLR+/etf+5/T/wZut9t6+eWXe+kqI8esWbOsH/zgBwHP3Xjjjdatt95q/sz7HVqtA4qd93fr1q3m6zZu3Og/56233rJcLpd14MCBbl1PVE7x1NfXy+bNm81QVcv9ffRxQUFBn15bpKmqqjL3aWlp5l7f94aGhoD3fvz48TJixAje+27QKZxZs2YFvK+K9zu0/va3v8n5558v3/nOd8wU5te+9jX5wx/+4D++Z88eKSsrC3i/dd8RnULm/e68r3/967Jq1Sr54osvzONPP/1U3n//fbnqqqvMY97vnmXn/dV7ndbR/1/46Pn6mbp+/fpuvX5YbhbYXYcOHTJzmxkZGQHP6+Pt27f32XVFGt11WtdCTJs2TfLy8sxz+pc9ISHB/IVu/d7rMXTeK6+8Ih9//LFs3LjxtGO836H15ZdfyrPPPisLFiyQn/3sZ+Y9/+d//mfzHs+ePdv/ngb7t4X3u/N++tOfml10NVTHxsaaf7d/8YtfmPUOive7Z9l5f/Vew3pLcXFx5pfS7v43iMqAgt77rb6wsND8xoOeoVuf/+hHP5KVK1eaxd7o+dCtvyn+8pe/NI91BEX/ji9dutQEFITWX/7yF3nppZdk2bJlMnHiRNmyZYv5pUcXdPJ+R76onOIZMmSISeOtKxn0cWZmZp9dVySZN2+erFixQt59913Jzs72P6/vr06xVVZWBpzPe981OoVTUVEh5513nvmtRW9r166Vp556yvxZf9Ph/Q4drWTIzc0NeG7ChAmyb98+82ffe8q/LaHx4IMPmlGUm2++2VRL3X777TJ//nxTLah4v3uWnfdX7/XfoJYaGxtNZU93/xtEZUDR4djJkyebuc2Wvxnp4/z8/D69tnCn66w0nCxfvlxWr15tygNb0vc9Pj4+4L3XMmT9B573vvMuv/xy+fzzz81vlr6b/oavQ+C+P/N+h45OV7Yum9f1ESNHjjR/1r/v+o9yy/dbpyh0Lp73u/NOnDhh1jK0pL9c6r/Xive7Z9l5f/VefwHSX5Z89N9+/W+ka1W6xYriMmNdifz888+bVch33323KTMuKyvr60sLa3PnzjUlaWvWrLFKS0v9txMnTgSUvWrp8erVq03Za35+vrkhNFpW8Sje79CWcsfFxZny1507d1ovvfSS1b9/f+vFF18MKMvUf0tef/1167PPPrOuu+46yl67aPbs2dbw4cP9ZcZaCjtkyBDrJz/5if8c3u/uVwB+8skn5qaR4D/+4z/Mn/fu3Wv7/dUy46997Wum9P799983FYWUGXfT73//e/MPt/ZD0bJjreFG9+hf8GA37Y3io3+xf/jDH1qDBg0y/7jfcMMNJsSgZwIK73dovfHGG1ZeXp75BWf8+PHWf//3fwcc19LMRYsWWRkZGeacyy+/3NqxY0efXW84q66uNn+X9d/pxMRE64wzzjA9O+rq6vzn8H53z7vvvhv032wNh3bf38OHD5tAoj1qkpOTrTvvvNMEn+5y6f90bwwGAAAgtKJyDQoAAHA2AgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAAHAcAgoAABCn+f/j4ROAEL0O1gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "consensus_threshold = 5\n", "size_threshold = number_of_samples * 0.50\n", "number_of_iterations = 10\n", "most_votes = 0\n", "best_m = None\n", "best_c = None\n", "# Now solve the line-fitting using RANSAC.\n", "for i in range(number_of_iterations):\n", " \n", " # Need a minimum of 2 samples to fit a line, but use 3, to average out the baseline noise.\n", " indexes = np.random.randint(0, len(x_array), 3)\n", " new_x_array = [x_array[x] for x in indexes]\n", " new_y_array = [y_array[y] for y in indexes]\n", " m, c = np.polyfit(new_x_array, new_y_array, deg=1)\n", "\n", " # Now for each model, measure consensus set.\n", " votes = 0\n", " tmp_x = []\n", " tmp_y = []\n", " for j in range(len(x_array)):\n", " diff = y_array[j] - (m * x_array[j] + c)\n", " distance = np.sqrt(diff * diff)\n", " if distance < consensus_threshold:\n", " votes += 1\n", " tmp_x.append(x_array[j])\n", " tmp_y.append(y_array[j])\n", " if votes > size_threshold:\n", " m, c = np.polyfit(tmp_x, tmp_y, deg=1)\n", " if votes > most_votes:\n", " best_m = m\n", " best_c = c\n", " most_votes = votes\n", "\n", "print(f\"m={best_m}, c={best_c}, votes={most_votes}\")\n", "plt.scatter(x_array, y_array)\n", "plt.plot(x_array, best_c + best_m * x_array, color=\"k\", lw=2.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "458b0f32-0057-4442-a5ff-d3c3a6f8dca5", "metadata": {}, "source": [ "# How many samples?\n", "\n", "Given:\n", "\n", "* *w*: probability that any given point is within error tolerance of model\n", "* *z*: certainty that at least one of our subsets contains error free data\n", "* *k*: number of iterations to try\n", "* *n*: number of points to fit our model.\n", " \n", "$$\n", "(1 - w^n)^k = (1 - z)\n", "$$\n", "\n", "then\n", "\n", "$$\n", "k = \\frac{log(1-z)}{log(1 - w^n)}\n", "$$\n", "\n", "So, if we think that 20 percent of the data is bad, the probability of picking a good data point is w = 0.8. If we want 99 percent confidence that we will end up with a good model then, z = 0.99. So, in a line fitting problem, n = 3, so the number of iterations k we should attempt should be:\n", "\n", "$$\n", "k = \\frac{log(1-0.99)}{log(1 - 0.8^3)}\n", "$$\n", "so,\n", "$$k=6.4$$\n", "i.e. we would expect around 7 iterations to work.\n" ] }, { "cell_type": "markdown", "id": "2194f1fb-9d73-4b28-8079-227eb635031e", "metadata": {}, "source": [ "# Final Thoughts\n", "\n", "So, the RANSAC method can be used for these types of model fitting problems, where you know you could have quite large, randomly occuring errors. Note, that this is not Trimmed Least Squares (TLS). In TLS, we compute the model, throw away for example 20% of the worst fitting data, and re-fit. This can more more reliable than normal least squares, but can still fail. \n", "\n", "For RANSAC, given the inlier distance threshold, and the number of votes threshold, we can tune the method to be fairly reliable. The above method shows how to estimate the required number of iterations, however, I suspect that in practice, most people will just empirically test. It would depend on how long each model fitting takes, and how quickly you want an algorithm to either keep searching, or fail early. \n", "\n", "See also: Applying RANSAC to pivot calibration, [here](https://mphy0026.readthedocs.io/en/latest/calibration/ransac.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "21a7db3c-4d76-4869-8ceb-0feabc136074", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "mphy0026", "language": "python", "name": "mphy0026" }, "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.10.16" } }, "nbformat": 4, "nbformat_minor": 5 }