From ae0a8d5aa42b93ca452f93906e69a202a444841a Mon Sep 17 00:00:00 2001 From: Jessica Lundquist Date: Sun, 14 Nov 2021 23:02:00 +0000 Subject: [PATCH] updates to lab 7 and HW 7 --- modules/lab7/lab7-1.ipynb | 12 +- modules/lab7/lab7-3.ipynb | 406 ++++++++++++++++++++++++++++++++++---- modules/module7.md | 29 ++- 3 files changed, 397 insertions(+), 50 deletions(-) diff --git a/modules/lab7/lab7-1.ipynb b/modules/lab7/lab7-1.ipynb index c860920d..9b5338ad 100644 --- a/modules/lab7/lab7-1.ipynb +++ b/modules/lab7/lab7-1.ipynb @@ -155,6 +155,8 @@ "\n", "Take one step through the Markov chain. What fraction of people live in the city after one step?\n", "\n", + "This can be defined as the fraction of people who live in the city and stay in the city + the fraction of people who lived in the suburbs but move into the city.\n", + "\n", "$p_{city,1}$ = 0.95 * 0.6 + 0.03 * 0.4" ] }, @@ -182,6 +184,8 @@ "source": [ "What fraction of people live in the suburbs?\n", "\n", + "Similar to above, this is the fraction of people who lived in the city and moved to the suburbs plus the fraction of people who live in the suburbs and stay in the suburbs.\n", + "\n", "$p_{suburb,1}$ = 0.05 * 0.6 + 0.97 * 0.4" ] }, @@ -313,7 +317,7 @@ } ], "source": [ - "# three steps all at once by taking our Pmarkov matrix to the 3rd power, then multipying by our original populaiton percentage vector\n", + "# three steps all at once by taking our Pmarkov matrix to the 3rd power, then multipying by our original population percentage vector\n", "percent3 = np.dot(percent0, np.linalg.matrix_power(Pmarkov,3))\n", "print(percent3)\n", "\n", @@ -543,7 +547,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 22, "metadata": {}, "outputs": [ { @@ -560,7 +564,7 @@ ], "source": [ "# This counts the transitions from each state to the next and marks that count\n", - "S = sparse.csr_matrix((np.ones_like(data[:-1]), (data[:-1], data[1:])), dtype=np.float)\n", + "S = sparse.csr_matrix((np.ones_like(data[:-1]), (data[:-1], data[1:])), dtype=float)\n", "\n", "# This converts those counts to matrix form\n", "tm = S.todense()\n", @@ -704,7 +708,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.4" + "version": "3.8.12" } }, "nbformat": 4, diff --git a/modules/lab7/lab7-3.ipynb b/modules/lab7/lab7-3.ipynb index 78c7fd60..1d65d864 100644 --- a/modules/lab7/lab7-3.ipynb +++ b/modules/lab7/lab7-3.ipynb @@ -7,7 +7,7 @@ "# Lab 7-3: MCMC Rating Curves\n", "\n", "\n", - "(Steven Pestana, 2019 | Derived from CEE599_bayesian_rating_curves_Lab6_v2.m)" + "(Steven Pestana, 2019 | Derived from CEE599_bayesian_rating_curves_Lab6_v2.m by Jessica Lundquist)" ] }, { @@ -52,7 +52,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -77,7 +77,7 @@ "dict" ] }, - "execution_count": 4, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -95,7 +95,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -343,7 +343,7 @@ " [array(['6/21/06 3:47'], dtype='0\n", " [9/26/08 16:30]\n", " 0.1805\n", - " 0.0778599\n", + " 0.07786\n", " \n", " \n", " 1\n", " [9/19/08 0:02]\n", " 0.2197\n", - " 0.0719143\n", + " 0.071914\n", " \n", " \n", " 2\n", @@ -449,15 +449,15 @@ "" ], "text/plain": [ - " date_of_obs h1 Qobs1\n", - "0 [9/26/08 16:30] 0.1805 0.0778599\n", - "1 [9/19/08 0:02] 0.2197 0.0719143\n", - "2 [9/10/08 21:35] 0.2406 0.143829\n", - "3 [9/10/10 17:52] 0.2407 0.168177\n", - "4 [8/20/07 18:22] 0.2565 0.243489" + " date_of_obs h1 Qobs1\n", + "0 [9/26/08 16:30] 0.1805 0.07786\n", + "1 [9/19/08 0:02] 0.2197 0.071914\n", + "2 [9/10/08 21:35] 0.2406 0.143829\n", + "3 [9/10/10 17:52] 0.2407 0.168177\n", + "4 [8/20/07 18:22] 0.2565 0.243489" ] }, - "execution_count": 7, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -477,7 +477,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -486,7 +486,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -520,13 +520,13 @@ " 0\n", " 2008-09-26 16:30:00\n", " 0.1805\n", - " 0.0778599\n", + " 0.07786\n", " \n", " \n", " 1\n", " 2008-09-19 00:02:00\n", " 0.2197\n", - " 0.0719143\n", + " 0.071914\n", " \n", " \n", " 2\n", @@ -551,15 +551,15 @@ "" ], "text/plain": [ - " date_of_obs h1 Qobs1\n", - "0 2008-09-26 16:30:00 0.1805 0.0778599\n", - "1 2008-09-19 00:02:00 0.2197 0.0719143\n", - "2 2008-09-10 21:35:00 0.2406 0.143829\n", - "3 2010-09-10 17:52:00 0.2407 0.168177\n", - "4 2007-08-20 18:22:00 0.2565 0.243489" + " date_of_obs h1 Qobs1\n", + "0 2008-09-26 16:30:00 0.1805 0.07786\n", + "1 2008-09-19 00:02:00 0.2197 0.071914\n", + "2 2008-09-10 21:35:00 0.2406 0.143829\n", + "3 2010-09-10 17:52:00 0.2407 0.168177\n", + "4 2007-08-20 18:22:00 0.2565 0.243489" ] }, - "execution_count": 9, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -572,33 +572,135 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "---\n", - "### PART 1: Brute force parameter estimation\n", + "### Calculating rating curves\n", "\n", - "First, to simplify overall process, I'm pretending we know h11 (the transition between the two slopes in the rating curve) h11=0.54 EXACTLY, (based on work by former CEE M.S. student Gwyn Perry.)" + "At stream gauges all over the world, we measure a timeseries of water height at a fixed point, and we generally develop an empirical rating curve to determine discharge, which is volume of flow per unit time, in units of of $m^3/s$ (cubic meters per second or cms) or $ft^3/s$ (cubic feet per second of cfs).\n", + "\n", + "In developing this rating curve, we are trying to solve for a, b, and c in the equation\n", + "$Q = a(h-b)^c$" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAFzCAYAAAAKZcKfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk8UlEQVR4nO3dfXicdZ3v8c+3CTITJDwspQS7NKBBbYBWjU9H14pgF7ogtp4ttS3Lcrwsoe76sNZL9uxBPe3ltawH3F33pLOJyoGePUh8CI8HRI+ulh6NboACLQ/CSUVrh1AsC9hkWtN8zx8zSSfpTHJnOvc998y8X9eVK5l77rnny+Qq/fR3f3+/n7m7AAAA4mBOpQsAAAAYRzABAACxQTABAACxQTABAACxQTABAACxQTABAACx0VjpAoI45ZRTvLW1tdJlAACAMnjwwQdfcPe5hZ6rimDS2tqqgYGBSpcBAADKwMyeLfYct3IAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAEBsEEwAAKgj6XRaS5Ys0SOPPKIlS5boueeeq3RJkxBMAACoI5s2bdK2bdu0Zs0abdu2TRs3bqx0SZOYu1e6hhl1dHQ4e+UAAFC6ZDKpTCZT9PlEIqGRkZFIajGzB929o9BzjJgAAFAHBgcHtXr1ajU0NEw6nkwmtWbNGu3atWvi2Pjtnkrc5iGYAABQB8466yzdeuutOnTo0KTjIyMjam5u1mmnnTZxbPx2TyVu8xBMAACoA4ODg5o/f/6kEZNTTz1VV1555cTISDKZlJkplUppbGxMqVRKZqZkMhlZnQQTAADqQEtLiy655BK5uxKJhObMmaMPfehDuvnmm9XX1yfp8O2epqYmSVJTU9MRt3nCRjABAKBODA0NqbOzU/39/ers7Dyih6SlpUXNzc3KZDJKJBLKZDJH3OYJW2Nk7wQAACpqfGREkrq6ugqeMx5e1q1bp56eHqXT6ajKk8R0YQAAEDGmCwMAgKpAMAEAALFBMAEAALFBMAEAoAqVsjprJVd0DYpgAgBAFSplddZKrugaFLNyAACoIsU245tuE75SXhMmZuUAAFAjSlmdNQ4rugZFMAEAoIqUsjprHFZ0DYpgAgBAlZlpaflyvaYS6DEBAACRoscEAABUBYIJAACIDYIJAACIDYIJAACIjdCCiZn9oZn9q5k9YWY7zewTueMnm9n3zezp3PeTwqoBAABUlzBHTEYlfdrd3yjpHZI+ZmYLJV0r6Qfu3ibpB7nHAAAA4QUTd0+7+0O5n1+R9ISk10i6TNItudNukfTBsGoAAADVJZIeEzNrlfQmST+TNM/d01I2vEg6tchr1pnZgJkN7N27N4oyAQBAhYUeTMzs1ZK+I+mT7v5y0Ne5e4+7d7h7x9y5c8MrEAAAxEaowcTMjlE2lPwvd+/LHR4ys5bc8y2Sng+zBgAAUD3CnJVjkr4u6Ql3/3LeU3dJujL385WS7gyrBgAAUF0aQ7z2uyRdIekxM9ueO/afJV0v6Ztm9hFJv5L0pyHWAAAAqkhowcTdt0myIk9fENb7AgCA6sXKrwAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAASEqn01qyZImee+65SpdS1wgmAABI2rRpk7Zt26aNGzdWupS6Zu5e6Rpm1NHR4QMDA5UuAwBQg5LJpDKZzBHHE4mERkZGKlBR7TOzB929o9BzjJgAAOra4OCgVq9eraamJklSU1OT1qxZo127dlW4svpEMAEA1JWpvSQtLS1qbm5WJpNRIpFQJpNRQ0ODLr/8cvpNKoBgAgCoK4V6SYaGhtTZ2an+/n51dnbqgQceoN+kQugxAQDUhSC9JPSbRIMeEwBA3QvSS0K/SeURTAAAdaFQL0lzc7NOO+20WZ2DcBFMAAB1Y2ovSaHm1iDnIDz0mAAAgEjRYwIAAKoCwQQAAMQGwQQAAMQGwQQAAMQGwQQAAMQGwQQAAMQGwQQAAMQGwQQAUNWm7hY81fbt23XiiSfq0UcfjbgylIJgAgCoaoV2C863du1avfTSS1q9enXElaEUrPwKAKhKM+0EbGZFX5v/d186ndaqVavU29vLnjgRYeVXAEDNGRwcVFtb28TjqTsBP/zww1qwYMGk17S2tuqRRx6ZdGymERdEixETAEDVKTZa0tDQoNHR0YnH7e3tevzxxyc93rFjx7TXGB9xQXgYMQEA1JTBwUGtXr1ac+Zk/xpLJBJqa2vT0qVLJ5334osvqr29Xb29vWpvb9e+ffuOuEZTU5OkI0dcUBmNlS4AAIDZamlpUXNzs6RsKDl48KAuvPBCbd68edJ5e/bsmfh55cqVBa+RyWSUSCSUyWTU3NxMn0mFMWICAKhKQ0ND6uzsVH9/vzo7O4tOFw77GigvekwAAECk6DEBAABVgWACAABiI7RgYmY3mdnzZrYj79gXzOw3ZrY997UsrPcHAADVJ8wRk5slXVTg+N+7++Lc170hvj8AoAZMtxfOTPvkoPqEFkzcfaukfTOeCADANKZbmZVVW2tPqLNyzKxV0j3ufk7u8Rck/bmklyUNSPq0u78403WYlQMA9We6lVklsWprFYvTrJyUpNdKWiwpLenGYiea2TozGzCzgb1790ZUHgAgLqZbmZVVW2tXpMHE3Yfc/ZC7j0n6qqS3TXNuj7t3uHvH3LlzoysSABAL063MyqqttSvSYGJmLXkPl0vaUexcAACmW5mVVVtrU2g9Jmb2DUnvlXSKpCFJn889XizJJf1S0tXunp7pWvSYAABQO6brMQltEz93/3CBw18P6/0AAED1Y+VXAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGwQTAAAQGzNu4mdmp0p6l6TTJY1I2iFpwN3HQq4NAADUmaLBxMzOl3StpJMlPSzpeUkJSR+U9Foz+7akG9395QjqBAAAdWC6EZNlkj7q7r+a+oSZNUq6RNL7JX0npNoAAECdKRpM3P0z0zw3KumOMAoCAAD1a8bmVzP7hJk1W9bXzewhM1saRXEAAKC+BJmV859yfSRLJc2VdJWk60OtCgAA1KUgwcRy35dJ+h/u/kjeMQAAgLIJEkweNLPvKRtM7jez4yUxVRgAAJTdjOuYSPqIpMWSBt192Mz+QNnbOQAAAGU1YzBx9zEzG5X0ntw04XGPhlcWAACoR0FWfr1J0nmSdurwLRyX1BdiXQAAoA4FuZXzDndfGHolAACg7gVpfv2pmRFMAABA6IKMmNyibDh5TtIBZacKu7ufF2plAACg7gQJJjdJukLSY2KaMAAACFGQYPIrd78r9EoAAEDdCxJMnjSzWyXdreytHEmSuzMrBwAAlFWQYJJUNpDkb9zHdGEAAFB2QRZYY5VXAAAQiRmnC5vZLWZ2Yt7jk3KLrgEAAJRVkHVMznP3fx9/4O4vSnpTaBUBAIC6FSSYzDGzk8YfmNnJCtabAgAAMCtBAsaNkn5iZt9Wtul1paQvhloVAACoS0GaX7eY2YCk9ym76usKd3889MoAAEDdKRpMzOzV7v47ScoFkSPCSP45AAAAR2u6HpM7zexGM3uPmR03ftDMzjKzj5jZ/ZIuCr9EAABQL4qOmLj7BWa2TNLVkt6Va3r9vaSnJP1vSVe6+3PRlAkAAOrBtD0m7n6vpHsjqgUAANS5INOFAQAAIkEwAQAAsUEwAQAAsREomJjZu83sqtzPc83szHDLAgAA9SjIJn6fl/RZSX+dO3SMpH8JsygAAFCfgoyYLJf0AUn7Jcnd90g6PsyiAABAfQoSTA66uyu7T47yF1sDAFS/dDqtJUuW6LnnWJoKlRckmHzTzLolnWhmH5X0fyR9NdyyAABR2bRpk7Zt26aNGzdWuhRAlh0MmeEks/dLWqrsJn73u/v3wy4sX0dHhw8MDET5lgBQ85LJpDKZzBHHE4mERkZGKlAR6oWZPejuHYWeCzQrx92/7+6fcfcNUYcSAEA4BgcH1dbWNvG4qalJa9as0a5duypYFerdjCMmZvaKcv0leV6SNCDp0+4+GFJtExgxAYDyKjZa0tDQoNHR0QpUhHpytCMmX5b0GUmvkTRf0gZle0xuk3RTuYoEAERncHBQq1ev1pw52b8GEomE2tratHTp0gpXhnoXJJhc5O7d7v6Ku7/s7j2Slrl7r6STQq4PABCClpYWNTc3S8qGkoMHD+rCCy/UvfeybysqK0gwGTOzlWY2J/e1Mu+5mTtnAQCxNDQ0pM7OTvX396uzs5PpwoiFID0mZ0n6R0nvVDaI9Ev6lKTfSHqLu28Lu0h6TAAAqB0l95iYWYOka9z9Unc/xd3n5n5+xt1HogglAIBgWCgNtWDaYOLuhyS9JaJaAABHgYXSUAuC3Mq5UVKbpG8pt1+OJLl7X7ilHcatHAAojoXSUG2OdrrwyZJ+K+l9ki7NfV1SvvIAAEdjfOpvU1OTJBZKQ3VrnOkEd78qikIAAKUZn/qbyWSUSCSUyWTU3Nys0047beKcdDqtVatWqbe3d9JxIG5mHDExs4SZfczMNpvZTeNfAV53k5k9b2Y78o6dbGbfN7Onc99ZBwUAymCmqb/0n6BaBOkx+ZakJyWtlrRR0hpJT7j7J2Z43Xsk/U7SFnc/J3fsS5L2ufv1ZnatpJPc/bMzFUmPCQCUhv4TxNHR9pi8zt2vk7Tf3W+R9CeSzp3pRe6+VdK+KYcvk3RL7udbJH0wwPsDAEpE/wmqTZBg8vvc9383s3MknSCptcT3m+fuaUnKfT+12Ilmts7MBsxsYO/evSW+HQDUtyD9J0CcBAkmPblekOsk3SXpcUlfCrUqSe7e4+4d7t4xd+7csN8OAKpO0AXVWHoe1WTGHpOjurhZq6R78npMnpL0XndPm1mLpB+5++tnug49JgBwpPXr16u7u1tXX321Nm/eXOlygMCm6zGZcbqwmR0r6UPK3r6ZON/dS2ntvkvSlZKuz32/s4RrAEBdm9rQmkqllEqlaGhFTQhyK+dOZZtWR5Vd+XX8a1pm9g1JP5X0ejPbbWYfUTaQvN/Mnpb0/txjAMAs0NCKWjbjiImk+e5+0Wwv7O4fLvLUBbO9FgDgMBpaUcuCjJj8xMxmnB4MAIgODa2oVUWbX83sMUmu7KhKm6RBSQckmSR39/OiKpLmVwAAakepza9s1AcAACJV9FaOuz/r7s9KalF2Gfnxx/skcSMTAMok6HokQD0I0mOSUnbPm3H7c8cAAGXABnvAYUGCiXleI4q7jynYbB4AwDSSyaTMTKlUSmNjY0qlUjIzJZPJSpcGVEyQYDJoZh83s2NyX59QthEWAHAUpq5HMmfOHK1YsYL1SFDXggSTTkn/QdJvJO2W9HZJ68IsCgDqQf56JA0NDRobG9NTTz3FeiSoa6HulVMuTBcGUKvGA8lULC+PWjbddOEZR0zM7Etm1py7jfMDM3vBzNaWv0wAqD+7d+9meXkgT5BbOUvd/WVl1zXZLelsSZ8JtSoAqBMsLw9MFiSYHJP7vkzSN9x9X4j1AEDdYXl54LAg037vNrMnJY1IWm9mcyVlZngNAEDZxdNWrVql3t7eoqMgfX19Ez93dXVFVRoQSzOOmLj7tZLeKanD3X+v7AJrl4VdGADUAhZPA2Znuk383ufuPzSzFYWed/e+QsfDwKwcANUmmUwqkzlycJnZNkDps3KW5L5fWuCLDf4AYBpTF09jtg0QTNEeE3f/fO77VdGVAwC1gdk2QGmKBhMz+6vpXujuXy5/OQBQO8Zn26xbt049PT1Kp9OVLgmIvelm5Ryf+/56SW+VdFfu8aWStoZZFADUAmbbALM33a2c/ypJZvY9SW9291dyj78g6VuRVAcAAOpKkAXWzpB0MO/xQUmtoVQDAADqWpAF1v6npJ+b2e2SXNJySbeEWhUAAKhLMwYTd/+imd0n6Y9yh65y94fDLQsAANSjICMmcveHJD0Uci0AAKDOBekxAQBMkU6ntWTJEjbcA8qMYAIAJWAPHCAcBBMANSmsEY1kMikzUyqV0tjYmFKplMxMyWSyrO8D1CuCCYCaVO4RjfGg09/fzx44QIgIJgBqSlgjGuNBp7u7mz1wgBAFmpUDANVicHBQGzZs0B133KHh4WE1NTVp+fLluuGGG0q6XjKZVCaTmXicSqUkSQ0NDerv72cPHKDMGDEBUFPGd/UdGRnRnDlzNDIyclQjGoODgwVv3ezevVuLFi1SV1fXpD1xABwdggmAmjM0NKSFCxfK3bVw4cKjaoAdDzrcugGiYe5e6Rpm1NHR4QMDA5UuA0AVmHrrZVwikdDIyEhJ11yxYoVaWlq0bt069fT0aNeuXdq/f796e3sJKEAJzOxBd+8o9BwjJgBqSrFbL0cza6avr09dXV0Tt25aW1tZwwQICcEEQE0J89YLa5gA4SOYAKg5Q0ND6uzsVH9/vzo7O8u2yFoYozEAJmO6MICakz9Lpqurq+A56XRaq1atOqJPpNhxiUZYIAqMmACoS8VWhp1pxdiwRmMAZDErB0BdKTZrp5ijmc0DoDBm5QBATrE+ke3bt9M/AsQAwQRAXSnWJ7Jo0SL6R4AYIJgAqDvF+kToHwEqjx4TAHVjuhk3AKJDjwkAaOYZNwAqjxETADUvjP1zAJSOERMANS+dTmvJkiUF+0JYsRWoHgQTADVhuts0rNgKVA+CCYCqFnRjPWbcANWBYAKgqky9ZRP0Nk1fX5+6urq0aNEidXV1TdpPB0B8EEwAVJWpt2y4TQPUFoIJgKow3S0bbtMAtYNgAiBy082gKWZwcFBtbW0Tj8dv2fT39+u3v/2trrvuOm7TADWAYAIgcrNd6CyZTOr000/X008/PXFseHhYt912m7q7u1k0DaghLLAGIDKlLnSWTqe1YcMG3XbbbRobG1MikdCBAwdU6P9fLJoGxB8LrAGIhVIXOhtvcJWywePgwYO64oorWDQNqEEEEwCROZoZNOMNrnfffbfmzZunvXv3MhsHqEEEEwCRKnUGzfg6JH19fRoaGlJrayuzcYAaRI8JgKrARnxA7aDHBEDVYyM+oD40VuJNzeyXkl6RdEjSaLHUBADjWOEVqA+VHDE5390XE0qA6lPKAmmlvGYqekqA2leRHpPciEmHu78Q5Hx6TIB4Wb9+vbq7u3X11Vdr8+bNob0GQG2arsekUsFkl6QXJbmkbnfvme58ggkQD6U0oNK0CmCqODa/vsvd3yzpYkkfM7P3TD3BzNaZ2YCZDezduzf6CgEcoZQGVJpWAcxGRYKJu+/JfX9e0u2S3lbgnB5373D3jrlz50ZdIoACpjagjoyM6Ic//OGsXpPJZNTQ0KDLL7+cHhEAR4g8mJjZcWZ2/PjPkpZK2hF1HQBKk9+AunDhQqXT6Rk30JvatPrAAw+w8R6AgiLvMTGzs5QdJZGy05VvdfcvTvcaekyAeCm1b4R+EwBSzHpM3H3Q3RflvtpnCiUA4qfUvhH6TQDMhJVfAczabBc7G1/DxMxYJA3AtAgmAEoym8XONm3aNNFTwiJpAKbDJn4AQkNPCYBCYtVjAqB+0FMCYLYIJgBCw8Z7AGaLYAIgVPSUAJgNekwAAECk6DEBAABVgWACAABig2ACAABig2ACAABig2ACAABig2ACAABig2ACAABig2AC1InxHX5Z4AxAnBFMgDqRv8MvAMQVK78CNY4dfgHEDSu/AnWMHX4BVBOCCVDj8nf4PfbYYzU8PKzGxkZ2+AUQSwQToA6M7/D7gQ98QJK0devWClcEAIU1VroAAOG77777JvWZ7Nq1S2ZGnwmA2GHEBKgD9JkAqBYEE6AO5PeZJBIJZTIZNTc302cCIHYIJkCdGO8z6e/vV2dnJwutAYgl1jEBqkg6ndaqVavU29vLaAeAqsU6JkCVmrqMPKu3Aqh1BBMgRooFkfnz58vMlEqlNDY2plQqJTNTMpmscMUAUF5MFwZiJD+IHDp0qOh5TU1NWr58uW644YYIqwOA8DFiAoQo6I6+yWRy0ojI1FDS1NSktra2ibVHmFUDoFYRTIAQBe0JKbTOyNQgMjo6qmuuuYZZNQBqGrNygBCUsqPvNddco56eHr3qVa/SwYMHtWDBAl188cVat26denp6lE6n1dfXF3bpABA6ZuUAEStlpdWp64wsXrxYXV1dWrRokbq6ugglAOoCza9ACEpZaTU/eHR1dUVRJgDEDiMmQEhYaRUAZo8eEwAAECl6TIAKCDpVGABwGMEECAnLxwPA7HErByizUqYKA0A94VYOEKFSpgoDALIIJkCZlTJVGACQRTABQhDmVOGgTbU03wKoRvSYAFVm/fr16u7u1tVXX63Nmzcf9XkAELXpekwIJghdOp3WqlWr1NvbG/ntjEq+d7kFbaql+RZA3NH8ioqq5LTZWpqyG7SpluZbANWMYILQJJNJmZlSqZTGxsaUSqVkZkomkzX93mEJ2lRL8y2AakYwQWhK+Zf79u3bdeKJJ+rRRx+dOFZKE2etjhoEbaplnx4A1YrdhetMlD0XpfzLfe3atXrppZe0evVq7dixQ9Lk2zHXXXddoPprddQg6A7E7FQMoFoxYlJnou65KPQv9/wRkPGfzUxmpp07d0qSdu7cOXEs/3bM6aefrq1btwaqn1EDAKg+zMqpMcVGRCo1U6NQPfnTWPfv368tW7bokksu0WOPPaZnn3124rXz58/Xueeeqx//+McaHh4ueH1mmgBA9WFWTh0pNiJSqZ6L/HoKNaRu2bJFknTPPfdMCiWSdMIJJ2jBggXKZDI69thjJUmNjY2R1g8AiBbBpEbMNAsl6p6LQvVkMhnNmTNnIhwV09vbq/b2du3bt2/idszPfvYztbe3a3R0tKZ6RgAAkxFMakSQEZEwei6KzZgpVs/atWsnAspUbW1tSqfTWrlypXbs2KE9e/aor69PXV1dWrRokc4++2ytX7+enhEAqGHMyqkR4yMiIyMjmjNnjkZGRo4YUQhjpkb+rZr8Zc+LjdB0d3drbGys4LVGR0enHQFhpgkA1D5GTGrI0NCQFi5cKHfXwoULQxlRGB8hSSQSMy5gVmiEZvfu3ZNGUhoaGrRs2TJdeeWVWrx4cdnrBQBUF0ZMasTUWTc7d+7Uzp07lUwmZ5y1Mj5z5itf+Yo+/vGPT7tGyPgIydq1azU6Oqo77rhDw8PDampq0vLly3XDDTdMnFtshCN/JOXgwYNasGABm8wBACQxYlL1xkcw+vv7S551Mx421qxZU3SNk6nNrFu2bNGtt96q4eHhWTejsr4IAKAYRkyq3Hio6O7unugxkaTh4WE1NjZOGxQKjbJIUiqVUiqVmrRGyODgoDZs2DBphOTkk0/WBRdcoE996lPq6elROp0OVDO9IgCAYggmMTKb5eKnhopUKnXEOVu3bp32GuNh4/bbb590uyeZTGrFihWTbssUama99NJLJ27BEDAAAOXArZyIFdqkbtxsloufOh23kF27dk27o+542Dhw4IAaGhokZZtRDxw4UPC2DLdgAABhq+sl6cu5oV3Qa51zzjnauXOn2tvbJzapK3W5+MbGRh06dKjo8w0NDbrsssvU1dVVtKYVK1aopaVFTz75pIaGhjRv3jy94Q1vUDqdnnTLBQCAcpluSfqK3Moxs4sk/aOkBklfc/frK1FHsTU4wriWmU16PL5JnSTt2bPniP6NqTNcClm6dKmeeeYZ/frXvy4YbA4dOqR58+YFXhsEAIBKi/xWjpk1SOqSdLGkhZI+bGYLo6xhpuXbw7jWww8/rAULFkw61traqkceeaTk5eLvvfdeXXDBBTp48KASiYQk6fjjj9fKlSu1cuVKnXnmmdxuAQBUlUr0mLxN0jPuPujuByXdJumyKAso54Z2Qa+1ePFiHXfccZOOHXfccTrvvPMkld6/kf+69evX68ILL1Rvb696e3s1ODjIiAgAoKpU4lbOayT9Ou/xbklvn3qSma2TtE6SzjjjjLIWUM4N7WZzrRdffFHt7e363Oc+p40bN2rfvn0Tz5U6hZaptwCAWlKJERMrcOyIDlx373H3DnfvmDt3btmLKOcMk6DX2rNnj3bs2DFpkzoAAHBY5LNyzOydkr7g7n+ce/zXkuTuf1vsNWHNygEAANGbblZOJUZM/k1Sm5mdaWavkrRK0l0VqAMAAMRM5D0m7j5qZn8h6X5lpwvf5O47o64DAADET0XWMXH3eyXdW4n3BgAA8cWS9AAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYIJgAAIDYi3yunFGa2V9KzFSzhFEkvVPD9we+g0vj8K4vPv/L4HZTXAncvuENvVQSTSjOzgWKbDSEa/A4qi8+/svj8K4/fQXS4lQMAAGKDYAIAAGKDYBJMT6ULAL+DCuPzryw+/8rjdxARekwAAEBsMGICAABig2CSx8wuMrOnzOwZM7u2wPNrzOzR3NdPzGxRJeqsVTN9/nnnvdXMDpnZf4yyvnoQ5HdgZu81s+1mttPMfhx1jbUswP+DTjCzu83skdznf1Ul6qxVZnaTmT1vZjuKPG9m9pXc7+dRM3tz1DXWA4JJjpk1SOqSdLGkhZI+bGYLp5y2S9ISdz9P0iZxz7FsAn7+4+f9naT7o62w9gX5HZjZiZI2S/qAu7dL+tOo66xVAf8MfEzS4+6+SNJ7Jd1oZq+KtNDadrOki6Z5/mJJbbmvdZJSEdRUdwgmh71N0jPuPujuByXdJumy/BPc/Sfu/mLuYb+k+RHXWMtm/Pxz/lLSdyQ9H2VxdSLI72C1pD53/5UkuTu/h/IJ8vm7pOPNzCS9WtI+SaPRllm73H2rsp9pMZdJ2uJZ/ZJONLOWaKqrHwSTw14j6dd5j3fnjhXzEUn3hVpRfZnx8zez10haLumfI6yrngT5M3C2pJPM7Edm9qCZ/Vlk1dW+IJ//f5f0Rkl7JD0m6RPuPhZNedDs/55ACRorXUCMWIFjBacsmdn5ygaTd4daUX0J8vn/g6TPuvuh7D8YUWZBfgeNkt4i6QJJSUk/NbN+d/9F2MXVgSCf/x9L2i7pfZJeK+n7ZvaAu78ccm3ICvz3BEpHMDlst6Q/zHs8X9l/lUxiZudJ+pqki939txHVVg+CfP4dkm7LhZJTJC0zs1F3vyOSCmtfkN/BbkkvuPt+SfvNbKukRZIIJkcvyOd/laTrPbvOwzNmtkvSGyT9PJoS616gvydwdLiVc9i/SWozszNzzWSrJN2Vf4KZnSGpT9IV/Aux7Gb8/N39THdvdfdWSd+WtJ5QUlYz/g4k3Snpj8ys0cyaJL1d0hMR11mrgnz+v1J2tEpmNk/S6yUNRlplfbtL0p/lZue8Q9JL7p6udFG1hhGTHHcfNbO/UHa2R4Okm9x9p5l15p7/Z0mfk/QHkjbn/tU+yqZO5RHw80eIgvwO3P0JM/uupEcljUn6mrsXnFqJ2Qn4Z2CTpJvN7DFlbyt81t3Z8bZMzOwbys52OsXMdkv6vKRjpInP/15JyyQ9I2lY2REslBkrvwIAgNjgVg4AAIgNggkAAIgNggkAAIgNggkAAIgNggkAAIgNggmAo2Zmn8ytaxL2+7SY2T2zfM0NZva+sGoCUF4EEwDl8ElJoQcTSX8l6auzfM0/Sbo2hFoAhIB1TAAEZmbHSfqmsktxNyi74Nc8STdIekrZ5erPN7OUpLcqu5/Ot93987nXL5P0ZUkvSHpI0lnufknuuv8k6VxlF378grvfWeD9ByW90d0PmNmfS/pgro5zJN0o6VWSrpB0QNIyd9+Xe92Dkv7E3Z8r+4cCoKwYMQEwGxdJ2uPui9z9HEnfdfevKLtfyPnufn7uvL/JrYp8nqQlZnaemSUkdSu7z9S7Jc3Nu+7fSPqhu79V0vmS/lsurEwwszMlvejuB/IOnyNptaS3SfqipGF3f5Okn0rK3/n4IUnvKscHACBcBBMAs/GYpAvN7O/M7I/c/aUi5600s4ckPSypXdJCZTebG3T3XblzvpF3/lJJ15rZdkk/kpSQdMaUa7ZI2jvl2L+6+yvuvlfSS5LuzquzNe+85yWdHui/EEBFsVcOgMDc/Rdm9hZl9wv5WzP7nrtvzD8nN7KxQdJb3f1FM7tZ2aBRaMv4iZdJ+pC7PzXNOSO56+TLHz0Zy3s8psn/f0vkXg8g5hgxARCYmZ2u7O2Sf1G2r+TNuadekXR87udmSfslvZTbAffi3PEnJZ1lZq25x5fnXfp+SX9pud0xzexNBd7+F5o8CjIbZ0tis0GgCjBiAmA2zlW2/2NM0u8lXZM73iPpPjNL55pfH5a0U9KgpP8rSe4+YmbrJX3XzF6Q9PO8626S9A+SHs2Fk19KuiT/jd19v5n9PzN7nbs/E7RgMztG0uskDcz6vxZA5JiVAyAyZvZqd/9dLnx0SXra3f9+Fq9fLukt7v5fZvmaN7v7dbOvGEDUuJUDIEofzTW47pR0grKzdAJz99uVHU2ZjUZlpxIDqAKMmAAAgNhgxAQAAMQGwQQAAMQGwQQAAMQGwQQAAMQGwQQAAMQGwQQAAMTG/wefMx5KW1RhDgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "h11 = 0.54" + "# First, let's plot all of the data we just read in.\n", + "plt.figure(figsize=(9,6))\n", + "\n", + "plt.xlabel('stage (m)')\n", + "plt.ylabel('discharge (cms)')\n", + "plt.plot(df.h1,df.Qobs1,'k*')" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAFzCAYAAAAUrPIsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAApHUlEQVR4nO3df5hcdXn38c+dCTAzyIJKaFaWJPwIXmWtCXGlcGGJCuURSh9NHouYYJHWLhustbZpG8tl1Wx/PC3Rp5e6WdnWH7SgjcoSEWkFvayQp2xgAwESEMEN0TxMQiwIwewI2b2fP+bMMtmd2T2TnZlz5sz7dV1zZeb8mHMfh5Pcfr/39/s1dxcAAECczIk6AAAAgMlIUAAAQOyQoAAAgNghQQEAALFDggIAAGKHBAUAAMTO3KgDqMaJJ57oixYtijoMAABQA9u2bfuZu88rt6+pEpRFixZpeHg46jAAAEANmNnuSvvo4gEAALFDggIAAGKHBAUAAMQOCQoAAIgdEhQAABA7JCgAACB2SFAAAEDskKAAAIDYIUEBAACxQ4ICAAAkSblcTueee67OO+887d27N9JYSFAAAIAkad26ddq6dauGhoa0fv165XI5LV++PJJkxdy94Rc9Ul1dXc5aPAAA1FYmk1E+n6+4f82aNdq4cWPNr2tm29y9q9w+WlAAAGgxk1tGZmqs6O/vl5kpk8k0IjxJESYoZpY2s/vM7CEz22lmn4wqFgAAWklvb6+2bNmi9evXS5J27dqlM844Y8pxqVRKkpTNZrV69Wrt2rWrYTHObdiVpvqlpLe7+4tmdpSkLWb27+4+FGFMAAAk1uSunP7+fvX391c8fmxsTOl0Wvl8Xm1tbZo/f34jwpQUYQuKF7wYfDwqeDVPQQwAAE1mZGREq1atUjablfRKy8j27dvV0dEx0WKSSqWUTqd11VVXaWhoSD09PQ0vlI20BsXMUma2XdIzku5y961RxgMAQJK1t7erra1N+Xz+sJaRJUuW6LLLLpO7K51Oy9119dVX68tf/rKWLFmivr4+DQ4ONjTWSBMUdx9z96WSOiSdY2ZvmHyMmXWb2bCZDe/fv7/hMQIAkCT79u1TT0/PlJaRStujEpthxmb2cUm/cPcNlY5hmDEAAMkRy2HGZjbPzE4I3mckXSTph1HFAwAA4iPKUTztkm40s5QKidLX3P32COMBAAAxEVmC4u4PSzo7qusDAID4YiZZAABiKsq1cKJGggIAQExMTkh6e3t1zz33aNmyZS2XpERZgwIAAEoUp6Dv6OjQ2NjYxPZcLqf29nal02mNjo5GGGHjkKAAABCxmVYTLsrn88pkMi2RpNDFAwBAA+VyOZ177rk677zzJrptbrvtNh111FFKp9OSClPQn3LKKYedN3fu3IYv2BclEhQAABqot7dXW7du1dDQ0MRqwh/5yEf08ssvHzYF/c9//vOJc1KplA4dOtTwBfuiRBcPAAANUK4bp9xqwsVjDhw4MLFtbGxMqVSqpQplaUEBAKDOcrmclixZoksuuWRixWBJMjNlMpnDju3o6NAll1wyZcXhPXv2NHzBviiRoAAAUGe9vb26//77tXv37sNG5xRXDy51/PHHa+HChVNWHG6Vrp0iEhQAAOokk8nIzNTf36/x8XE9+uijU4557rnn1NnZqU2bNqmzs1PPPvts7FYWjkJsVjMOg9WMAQDNJJfLae3atdq8ebMOHjw40b0zNjambDarFStWaMOGDS3XOlIUy9WMAQBIuvb2drW1tU1014yNjWlsbKylu27CIkEBAKCOSrtrTj31VJ166qkt3XUTFl08AAAgEnTxAACApkKCAgAAYocEBQAAxA4JCgAAiB0SFAAAEDskKAAAIHZIUAAAQOyQoAAAgNghQQEAALFDggIAaHnbt2/XCSecoIcffjjqUBAgQQEAtLwrr7xSzz//vFatWhV1KAjMjToAAACiYmaHfd65c+fEtmZaqy6JaEEBACRaLpfT8uXLy64c/OCDD2rhwoWHbVu0aJEeeuihRoWHCkhQAACJ1tvbqy1btmj9+vVT9i1dulTHHnvsYduOPfZYvfGNb2xUeKiABAUAkEiZTEZmpv7+fo2Pj6u/v19mpkwmc1irynPPPafOzk5t2rRJnZ2devbZZ6MOHSJBAQAk1L333qt58+Ypk8lIkrLZrFavXq1du3Yd1qry9NNPa8eOHbr88su1Y8cOPf300xFHDkmyZioC6urq8uHh4ajDAAA0gWuvvVb9/f2SpHQ6rZdeeklmprGxsSnHptNpjY6ONjrElmdm29y9q9w+RvEAABIlk8kon88fti2fzyuVSuniiy/Wq1/9am3evFkHDx5UNpvVihUrtGHDhoiiRSV08QAAEmVkZESrVq1SNpuV9ErXzp49e3THHXeora1N+Xxe6XRa+XxebW1tmj9/fsRRYzISFABA05lu6HB7e/u0Sci+ffvU09OjoaEh9fT0lP0ORI8uHgBA0+nt7dU999yjZcuW6YEHHpjSAlJMQrq7uzUwMKBcLjexb3BwcOJ9X19fw2JGdSiSBQA0jXL1JRJFrs1quiJZungAALE1uStnZGREc+ZM/acrn89PDCdGMpCgAABia/IssO3t7Vq9evVhx8ydO3difhMkBwkKACB2ppsF9sUXX1RnZ6fMTKlUSocOHWIkTgKRoAAAYqfSUOFdu3ZpcHBQZ555ptasWaNt27bp2muvZSROAjGKBwAQOzMNFWYkTvLRggIAiCXmK2ltDDMGAEQql8vpiiuu0KZNm6gjaTGxHGZsZqeY2ffN7DEz22lmH44qFgDAkZluRtewJo/UAaRou3gOSfpTd/9VSedK+qCZnRVhPACAKh1pcpHL5ZRKpSqO1AEiS1DcPefuDwTvD0h6TNLJUcUDAAhvumHAYfT29srdtXjx4rIjdYBYFMma2SJJZ0vaWmZft5kNm9nw/v37Gx4bAGCq4jDguXMLg0HDTpZWmti4u5544gkdPHhQkjQ6Osp8JpgQeYJiZq+SdIukP3b3Fybvd/cBd+9y96558+Y1PkAAwBSnnXaavvKVr+jQoUOSpEOHDunmm2/WqaeeOu15k+c3mTNnjo4//nhJ0llnncVIHUyINEExs6NUSE5udvfBmY4HAMTDyMiIOjo6lEqlJEmpVEodHR0VW1CKxbRmNjG/iSSNj4/r+eeflyTt3LlTt956KzUokBTtKB6T9AVJj7n7p6OKAwBQvfb2dl122WVyd6XTabm7fvu3f3uie2by6J7SYtri/CZ33XWXFi9ePJHkUIOCUlHOJHu+pPdJesTMtgfb/tLd74guJABAWMVEo7u7WwMDA8rlchP7iglJR0eHxsbGJrb39/dLktLptPr6+nThhRfqxz/+cdnZYtHaIktQ3H2LJIvq+gCA2Sk33Xwmk5noviknm81qxYoV2rBhg6Tpkxy0NtbiAQDUzL333quLL75YL774okZHR5XNZnXyySfrySef1DHHHMOaOggt8lE8AIDkGBgY0P79+zU6OjrRbXPo0CGtWbOGNXVQFVpQAABVKbd2TrmunXw+r1QqpaVLl060jtBKgrBoQQEAVKXc9PaT5zcpjsjZs2fPYd04QFgkKACAikqHC083vX17e/vE/CaMyEEtkKAAACpat26d7r77bq1bt65iK0lx3pLiiBxqTVAL5u5RxxBaV1eXDw8PRx0GADS1cjUkk003XNjMdMwxx+ill17SNddco40bN9YzXCSYmW1z965y+2hBAYAWU66GZLLp/s/rWWedRSsJ6o4WFABoEZVaRdLptEZHRw/blsvldMEFF+jJJ5+s+H3lzgOqQQsKAGDGGpJS7e3tEysVF5ZOewVr5qARSFAAoEVUO9LmJz/5iaSp3T0HDx5khA7qjgQFAFpINSNt9uzZc1iLSyqV0qWXXqqrrrqK2hPUHTPJAkCTCjMaZ7pjZ5rVdXKLy0svvaSFCxcyagcNQQsKADSpMKNxjuTYUsxtgqgwigcAmkw1o3GqORZoNEbxAECCTB6Nk8lkdNJJJ2nr1q0zHpvNZrVy5UotXbqU1hDEGgkKADSZybUho6OjeuaZZ/T5z39+xmPz+bwef/xx3XfffVV39wCNRIICAE1o3759MrPDum9KF++bfGxPT4/Gx8c1Pj6unTt3TlnsD4gbEhQAaEKDg4P66U9/GmritcHBQfX19empp54KPVEbEDUSFABoUmEnXsvlclq+fLnMrKqJ2oAokaAAQBMLMwy4dIgxw4bRLBhmDAAJxRBjxB3DjAGgBVWzOCAQNyQoAJBQ1S4OCMQJCQoAJBg1J2hW1KAAAIBIUIMCAACaCgkKAACIHRIUAAAQOyQoAAAgdkhQAABA7JCgAACA2CFBAYCIFBfxY24SYCoSFACISOkifgAOx0RtANBgLOIHFDBRGwDECIv4ATMjQQGABmMRP2Bmc2c6wMxOknS+pNdJGpW0Q9Kwu4/XOTYASKziIn7d3d0aGBhQLpeLOiQgViomKGb2NknrJL1G0oOSnpGUlvQuSaeb2TckfcrdX2hAnACQKIODg5IKI3l27NihTZs2RRwREC/TtaBcKukP3P0nk3eY2VxJl0n6TUm31Ck2AEi80pE8GzdujDocIDYYxQMAEWAkDzDLUTxm1mVmHzGz681svZldbmavqX2YANA6GMkDTK9igmJm7zezByR9VFJG0uMq1KG8RdJdZnajmS1oTJgAkCyM5AGmN10NyrGSznf3sm2NZrZU0mJJU2pUwjKzL6pQy/KMu7/hSL8HAJoRI3mAyiKtQTGzCyS9KOlfwiQo1KAAAJAcs61B+QczazOzo8zse2b2MzO7shaBufvdkp6txXcBQBRY8A+ojzAzyV4czHVymaQ9ks6U9Gd1jaqEmXWb2bCZDe/fv79RlwWAUFjwD6iPMAnKUcGfl0r6qrs3tMXD3Qfcvcvdu+bNm9fISwNAWblcTqlUSmam/v5+jY+Pq7+/X2amTCYTdXhAIoRJUL5lZj+U1CXpe2Y2T9LUwfsAkFCTu3F6e3vl7lq8eDHDhIE6mXEtHndfZ2Z/L+kFdx8zs19Iemf9QwOAeCh243R0dGhsbGxi+xNPPDHxnmHCQG2FWSwwJek3JC0Kprgv+vRsL25mX5X0VkknmtkeSR939y/M9nsBoBYqzfZaNGfOHJ1++unq7+/X4OAgw4SBGgrVxSPp/ZJeK+m4ktesuft73b3d3Y9y9w6SEwCNEHbkTbnZXhcvXiwzUzqdliRddNFFuvDCC9XX1zexACCA2ZuxBUVSh7u/se6RAECDhF2gr9xsr4cOHdKaNWuYXA2osxknagvqT77n7nc2JqTKmKgNwGwcyQJ9K1euVHt7+2EJCS0lQG1MN1FbmARlhaSbVOgOelmSSXJ3b6t1oDMhQQEwG7lcTmvXrtXmzZt18OBBZbNZrVixQhs2bKC4FYjAdAlKmC6eT0k6T9IjHuW8+AAwSyzQBzSPMEWyT0jaQXICIAmKC/QNDQ2pp6eHKeqBmArTxfNlSadJ+ndJvyxud/dZDzOuFl08AAAkx2y7eHYFr6ODFwAAQF2FmUn2k40IBAAAoGjGGhQzu8vMTij5/Goz+05dowIAAC0tTJHsPHf/efGDuz8n6aS6RQQAAFpemARlzMwWFD+Y2UJJjOgBAAB1E6ZI9jpJW8zsB8HnCyR11y8kAADQ6mZsQXH3/5C0TNImSV+T9CZ3pwYFQNMJu0gggOhVTFDMbFHxvbv/zN1vd/dvufvPgv1mZh0NiBEAaqJ0kUAA8VZxojYz+7oKCcw3JW2TtF9SWtIZkt4m6UJJH3f3uxoTKhO1ATgyR7JIIID6m26itootKO7+O5I+Jun1kvok3aNCsvIBSY9LensjkxMAOFIjIyNatWqVstmsJCmbzWr16tXatWtXxJEBqGTaIll3f1SFIlkAaFosEgg0nzDDjAGg6bFIINBcZlwsME6oQQEAIDmOqAYFAAAgKmHW4jEzu9LM/ir4vMDMzql/aAAAoFWFaUHZKOk8Se8NPh9QYVQPAABAXYSZ6v7X3X2ZmT0oFRYLNLOj6xwXAABoYWFaUF42s5SCBQLNbJ6k8bpGBQAAWlqYBOUzkm6VdJKZ/Y2kLZL+tq5RAUCTCbvOD+sBAeGEWSzwZkl/LunvJOUkvcvdv17vwACgmYRd54f1gIBwZpwHxcxeU2bzAXd/uT4hVcY8KADiJuw6P6wHBEw123lQHlBhocAfSXoieL/LzB4wszfVLkwAaD5h1/lhPSCgOmESlP+QdKm7n+jur5V0iaSvSbpWhSHIACJETUO0wq7zw3pAQHXCJChd7v6d4gd3v1PSBe4+JOmYukUGIBRqGqIXdp0f1gMCwgtTg3KnpO9J+rdg03sk/aakd0i6392X1TXCEtSgAK+gpgFAs5ttDcoqSR2SNgevU4JtKUmX1yZEANVqppoGuqEAVGvaBCWYoO0f3f1D7n528PqQu+9395fc/ckGxQlgkmaqaaAbCkC1pk1Q3H1M0jymtgfiKe41DZlMRmam/v5+jY+Pq7+/X2amTCYTdWgAYi5MDcoNkpZJuk3SL4rb3f3T9Q1tKmpQgOaSy+W0du1abd68WQcPHlQ2m9WKFSu0YcOGWLb0AGis6WpQwiwW+HTwmiPpuFoGBiDZmqkbCkC8zJiguPsnGxEIgGQqdkN1d3drYGBAuVwu6pAANIEwXTzzVFiLp1NSurjd3d9e39CmoosHAIDkmO0w45sl/VDSqZI+KekpSffXLDoAkdu+fbtOOOEEPfzww1GHAgCSwiUor3X3L0h62d1/4O6/J+ncOscFoIGuvPJKPf/881q1alXUoQCApHAJSnHV4pyZ/ZaZna3CxG0AmpyZycy0c+dOSdLOnTsntjG5GoAohUlQ/trMjpf0p5LWSvpnSR+pa1QAGuLBBx/UwoULD9u2aNEiPfTQQ2UnVyNpAdAoMxbJxglFskDtdXZ26tFHH534bGYq9/dCOp3W1VdfrRtuuEHXXHONNm5kMXMAszOrIlkzm2dmf2lmA2b2xeKrRoG9w8weN7MnzWxdLb4TQHWee+45dXZ2atOmTers7NS8efOmrPEzZ84c5fN5ZoQF0DBhuni+Kel4Sd+V9O2S16wE6/z0SbpE0lmS3mtmZ832e4HZarVujKefflo7duzQ5Zdfrh07dmjfvn1TJle78sorm2ZhQgDJECZBybr7X7j719z9luKrBtc+R9KT7j7i7i9J+jdJ76zB9wKzwsJ2U9f4OXDgADPCAmioMBO1/bWk/3L3O2p6YbN3S3qHu38g+Pw+Sb/u7n846bhuSd2StGDBgjft3r27lmEAEzKZjPL5/JTt6XRao6OjEUQULytXrlR7e/thM8IODg5GHRaAJnZENShmdsDMXpD0YUm3m9momb1Qsn3WcZXZNiVbcvcBd+9y96558+bV4LJAeSMjI4nqxqh1V9Xg4KD6+vq0ZMkS9fX1kZwAqKuKCYq7H+fubcGfc9w9U/K5rQbX3iPplJLPHSosSghEImkL203uqmq12hoAzS3MKJ4VwTwoxc8nmNm7anDt+yUtNrNTzexoSVdIuq0G3wtMqPYf5cm1F834j3kmk5GZTRlxc8opp7R8bQ2A5hGmBmW7uy+dtO1Bdz971hc3u1TSP0pKSfqiu//NdMczDwqqde2117bcvB25XE5r167V5s2bdfDgwYrHUVsDIGqzXSyw3DFzZxdSgbvf4e5nuvvpMyUnQDUqtSK0wrwdk7uqzEyLFy+eqK2ZM2eOVq5c2bS1NQBaQ5gEZdjMPm1mp5vZaWb2fyRtq3dgwGwkreC1kkpdWKVdVWvWrNHLL7+sfD6vVCql8fFxPf74401bWwOgNYRpCfmQpI9J2qTCyJs7JX2wnkEBs5W0gtdKSgthS7uwSkfY9PX1TSQmRcVFAenmARBXM7aguPsv3H1d0Ed0jqS/c/df1D80YHaSUPBaSbVdWHv27GmJFiUAyRFmFM9XzKzNzI6VtFPS42b2Z/UPDZidJM/bUW0XVqu0KAFIjjA1KGe5+wuS3iXpDkkLJL2vnkEBmN6RJBxJblECkDxhalCOMrOjVEhQPufuL5vZ9GOTAdRdMeEonXp+OpPrUgAgzsIkKDdIekrSQ5LuNrOFkmox1T2AWSDhAJBkMyYo7v4ZSZ8p2bTbzN5Wv5AAAECrq5igmNmV7n6Tmf1JhUM+XaeYAABAi5uuSPbY4M/jKrwANAEWCQTQjCq2oLj7DcGfn2xcOABqrdJkbgAQZxUXCzSzz5TdEXD3P6pLRNNgsUAgvEwmo3w+P2U7s8cCiIsjXSxwW/BKS1om6YngtVTSWI1jBFBjrbIeEYBkqpiguPuN7n6jpMWS3ubun3X3z0q6UIUkBUAd1KpmhNljATSzMDPJvk6HF8W+KtgGoA5Ka0Zmi9ljATSrijUoEweYXS3pE5K+H2xaLukTQetKQ1GDgiSjZgRAqznSGhRJkrt/SdKvS7o1eJ0XRXICJB01IwDwiooJipktKr53973u/s3gtTfYb2bW0YAYgSPSbPN/UDMCAK+YrgXlejO7xcx+18w6zewkM1tgZm83s15J/1fSrzYoTqBqtazlaBRqRgCgYNoaFDM7S9JqSedLapc0KukxSd+W9A13n9phXkfUoCAMajkAoDkccQ2Kuz/q7te5+1vd/fXuvtTd3+vuNzU6OQHCopYDAJrfjKsZm9nKMpufl/SIuz9T+5CA2aGWAwCa34wJiqTfl3SeXhlm/FZJQ5LONLP17v6vdYoNOGLFWo7u7m4NDAwol8tFHRIAoAph5kH5lqQPuPu+4POvSOqX9AFJd7v7G+oeZYAaFAAAkmNW86BIWlRMTgLPSDrT3Z+V9HItAgQAACgVpovnHjO7XdLXg8/vlnS3mR0r6ef1CgwAALSuMAnKByWtlPQWSSbpRkm3eKFv6G11jA0AALSoGRMUd3cz2yLpJUku6T6fqXAFAABgFmasQTGzyyXdp0LXzuWStprZu+sdGAAAaF1himSvk/Rmd7/K3X9X0jmSPlbfsIDka7a1ggCgkcIkKHMmTcj23yHPAzCNZlwrCAAaJcw8KNdLeqOkrwab3iPpYXf/izrHNgXzoCAJWCsIAApmNQ+Ku/+ZpAEVkpQlkgaiSE6ApGCtIACYWZhhxnL3WyTdUudYgJbAWkEAMLOKLShmdsDMXijzOmBmLzQySCDuqi14La4VNDQ0pJ6eHgplAWCSGWtQ4oQaFMRFLpfTFVdcoU2bNmn+/Pm69tprdcMNN+iaa67Rxo0bow4PAJrCdDUoJCjAESgmJGamsbGxKfspeAWAmU2XoISqQQFQUGkETlE2m9WKFSu0YcOGBkYFAMnDfCZAFcqNwFm8eLHMjIJXAKghEhSgCuVG4Bw6dEhr1qyh4BUAaogaFDS1ycWqjbBy5Uq1t7eru7tbAwMDyuVyGhwcbMi1ASBJKJJFYjF6BgCaFwkKEofp4gGg+c1qqvt6MLPfMbOdZjZuZmUDA6bDdPEAkGxRFcnukLRS0t0RXR9NjuniASDZIklQ3P0xd388imsjOZguHgCSK9IaFDP7T0lr3b1iYYmZdUvqlqQFCxa8affu3Q2KDgAA1FMkM8ma2XcllWtvv87dvxn2e9x9QNKAVCiSrVF4AAAgxuqWoLj7RfX6bgAAkGzMJAsAAGInqmHGK8xsj6TzJH3bzL4TRRwAACCeIlnN2N1vlXRrFNcGAADxRxcPAACIHRIUJF4ul9Py5cuZJwUAmggJChKvt7dXW7Zs0fr166MOBQAQEosFIrFYUBAA4i12iwUCjcCCggDQvEhQkFgsKAgAzYsEBYnGgoIA0JyoQQEAAJGgBgUAADQVEhQAABA7JChINCZpA4DmRIKCRGOSNgBoTiQoSIzS1pJMJiMzU39/v8bHx9Xf3y8zUyaTiTpMAEAIJCioq0Z2sZS2ljBJGwA0t7lRB4BkK00aNm7cWJdrTJ7Svr+/X/39/UqlUnJ3JmkDgCZECwrqopFdLJVaSy6++GImaQOAJkULCupiZGREa9eu1ebNm3Xw4EFls1mtWLFCGzZsqPm1Kk1pf9NNN00c09fXV/PrAgDqhxYU1EWj18FhSnsASBZaUFA3xaShu7tbAwMDyuVydbvW4ODgxHtaSwCg+bEWDwAAiARr8QAAgKZCggIAAGKHBAUAAMQOCQoAAIgdEhQAABA7JCgAACB2SFAAAEDskKAAAIDYIUEBAACxQ4ICAABihwQFAADEDgkKAACIHRIUTCuXy2n58uXau3dv1KEAAFoICQqm1dvbqy1btmj9+vVRhwIAaCHm7lHHEFpXV5cPDw9HHUZLyGQyyufzU7an02mNjo5GEBEAIGnMbJu7d5XbRwsKyhoZGdGqVauUzWYlSdlsVqtXr9auXbsijgwA0ApIUFBWe3u72tralM/nlU6nlc/n1dbWpvnz50cdGgCgBZCgoKJ9+/app6dHQ0ND6unpoVAWANAw1KAAAIBIUIMCAACaCgkKpsU8KACAKJCgYFrMgwIAiAI1KCiLeVAAAPUWuxoUM7vezH5oZg+b2a1mdkIUcaAy5kEBAEQpqi6euyS9wd3fKOlHkj4aURyogHlQAABRiiRBcfc73f1Q8HFIUkcUcWB6zIMCAIhK5DUoZvYtSZvc/aYK+7sldUvSggUL3rR79+5GhgcAAOpkuhqUuXW86HcllesPuM7dvxkcc52kQ5JurvQ97j4gaUAqFMnWIVQAABAzdUtQ3P2i6fab2VWSLpN0oUfdjAMAAGKlbgnKdMzsHZL+QtJydz8YRQwAACC+ohrF8zlJx0m6y8y2m9nnI4oDAADEUCQtKO5+RhTXBQAAzYGp7gEAQOyQoAAAgNghQQEAALFDggIAAGKHBAUAAMQOCQoAAIgdEhQAABA7JCgAACB2SFAAAEDskKAAAIDYIUEBAACxQ4ICAABihwQFAADEDgkKAACIHRIUAAAQOyQoAAAgdkhQAABA7JCgAACA2CFBAQAAsUOCAgAAYocERVIul9Py5cu1d+/eqEMBAAAiQZEk9fb2asuWLVq/fn3UoQAAAEnm7lHHEFpXV5cPDw/X7PsymYzy+fyU7el0WqOjozW7DgAAmMrMtrl7V7l9Ld2CMjIyolWrVimbzUqSstmsVq9erV27dkUcGQAAra2lE5T29na1tbUpn88rnU4rn8+rra1N8+fPjzq0qlFHAwBIkpZOUCRp37596unp0dDQkHp6epr2H3jqaAAASdLSNShJQB0NAKBZUYOSYNTRAACSiASlySWpjgYAgCISlARISh0NAABF1KAAAIBIUIMCAACaCgkKAACIHRIUAAAQOyQoAAAgdkhQAABA7JCgAACA2CFBAQAAsUOCAgAAYocEBQAAxA4JCgAAiB0SFAAAEDtNtRaPme2XtLvK006U9LM6hNMMWvXeue/W0qr3LbXuvXPfybHQ3eeV29FUCcqRMLPhSgsRJV2r3jv33Vpa9b6l1r137rs10MUDAABihwQFAADETiskKANRBxChVr137ru1tOp9S61779x3C0h8DQoAAGg+rdCCAgAAmkwiEhQz+x0z22lm42ZWtsLZzF5vZttLXi+Y2R8H+z5hZv+vZN+lDb2BIxTmvoPjnjKzR4J7Gy7Z/hozu8vMngj+fHVjIp+9kL/5KWb2fTN7LDj2wyX7kv6bv8PMHjezJ81sXcn2pvzNw8Sd0Gc81O+V0Gc8zG+emGe80jNbst/M7DPB/ofNbFnYc5uWuzf9S9KvSnq9pP+U1BXi+JSkvSqMv5akT0haG/V91Ou+JT0l6cQy2/9B0rrg/TpJfx/1PdXy3iW1S1oWvD9O0o8knZX03zz47/vHkk6TdLSkh0ruuyl/82rjTtAzHuq+E/qMzxh7Up7x6Z7ZkmMulfTvkkzSuZK2hj23WV+JaEFx98fc/fEqTrlQ0o/dvdpJ32LlCO57sndKujF4f6Okd806qAYJc+/unnP3B4L3ByQ9JunkRsRXLyF/83MkPenuI+7+kqR/U+G3lpr3N6827kQ845r979Wsv7cUIvYEPePTPbNF75T0L14wJOkEM2sPeW5TSkSCcgSukPTVSdv+MGg2+2IzNYOG5JLuNLNtZtZdsv1X3D0nFR50SSdFEl0DmNkiSWdL2lqyOam/+cmSflryeY9e+Uu7WX/zauNOyjMe9r6T+IxXFXuTP+PTPbMzHRPm3KbUNAmKmX3XzHaUeVWVKZrZ0ZL+p6Svl2zul3S6pKWScpI+Vau4Z6tG932+uy+TdImkD5rZBXUKt6Zq+Ju/StItkv7Y3V8INif5N7cy22I/XI9nnGe8VZ7xMsI8s5WOacrnPYy5UQcQlrtfVKOvukTSA+6+r+S7J96b2T9Jur1G15q1Wty3uz8d/PmMmd2qQpPg3ZL2mVm7u+eCpsJnZnutWqrFvZvZUSr8xXWzuw+WfHeSf/M9kk4p+dwh6engfWx/8+nu28yqiTsxz3jY+07iMx723pvxGS9jumd2pmOODnFuU2qaFpQaeq8mNf0G//EXrZC0o6ER1ZGZHWtmxxXfS7pYr9zfbZKuCt5fJembjY+wfszMJH1B0mPu/ulJ+xL7m0u6X9JiMzs1aE24QoXfWmre37yauJP0jM943wl+xsPce1Ke8eme2aLbJP1uYTCPnSvp+aDrK8y5zSnqKt1avFT4j2+PpF9K2ifpO8H210m6o+S4rKT/lnT8pPP/VdIjkh5W4Ydtj/qeanXfKlR2PxS8dkq6ruT810r6nqQngj9fE/U91fje36JCU+fDkrYHr0uT/psHny9VYUTDj5Pwm1eKuwWe8RnvO8HPeJh7T8wzXu6ZldQjqSd4b5L6gv2PqGQUX6XnvdlfzCQLAABipxW7eAAAQMyRoAAAgNghQQEAALFDggIAAGKHBAUAAMQOCQoASZKZvTjL879hZqdNs//9Zva62VwjZBwZM/uBmaWqOOcPzezqks8bzOzt9YkQQBgkKABmzcw6JaXcfWSaw96vwhwW9fZ7kgbdfayKc74o6Y9KPn9WhRV0AUSEBAXAYYKZKq8P1kR5xMzeE2yfY2YbzWynmd1uZneY2buD01YrmOnTzFJm9uWS8z8SHNcl6WYz2x60cvyVmd0fHDcQzAoqM3uzFRZ4u7cYR8n3Xh+c87CZXVPhFkpjeWvQmvI1M/uRmf1vM1ttZvcFsZ0uSe5+UNJTZnZO8Hm3pNea2fx6/G8MYGYkKAAmW6nCAmtLJF0k6fpg2vCVkhZJ+jVJH5B0Xsk550vaFrxfKulkd3+Du/+apC+5+zckDUta7e5L3X1U0ufc/c3u/gZJGUmXBed/SYXZM8+TVNoK8vsqTO/9ZklvlvQHZnZqaeDBVN+nuftTJZuXSPpwEPf7JJ3p7udI+mdJHyo5bljSb5R8fiC4LwARIEEBMNlbJH3V3ce8sODaD1RICN4i6evuPu7ueyV9v+Scdkn7g/cjkk4zs8+a2TskvaDy3mZmW83sEUlvl9RpZidIOs7d/ys45islx1+swlok2yVtVWEq9MWTvvNEST+ftO1+d8+5+y9VmAr8zmD7IyokXEXP6PAuqMmfATRQ06xmDKBhyi3fPt12SRqVlJYkd3/OzJZI+h+SPijpchXqQl75IrO0pI0qrCfyUzP7RHD+dNcwSR9y9++EiaPEL0vej5d8Htfhfwemg/MrfQbQQLSgAJjsbknvCWo+5km6QNJ9krZI+l9BLcqvSHpryTmPSTpDkszsRElz3P0WSR+TtCw45oCk44L3xSTiZ2b2KknvlgrJjaQDwWqtUmFl1qLvSFpjZkcF1zkzWL13QnB+KkiAqnWmDl/xdvJnAA1ECwqAyW5Vob7kIRVWiv1zd99rZrdIulCFf7R/pEI3y/PBOd9WIWH5rqSTJX3JzIr/B+ijwZ9flvR5MxsNvv+fVOhmeUqFJeOLfl/SP5nZLyT9Z8k1/lmFLpkHgoLa/ZLeVSb+O1Xojvpulfd9vqRPSlKQBJ2hQl0KgAiwmjGA0MzsVe7+opm9VoVWlfOD5CWjQk3K+VUO7614jeD9Oknt7v7hKs4/W9KfuPv7jvQcM1shaZm7f6y66AHUCi0oAKpxe1DIerSk3qBYVu4+amYfV6H15CezvMZvmdlHVfj7abcK86eE5u4Pmtn3zSxVRbJ0ogrdUUVzJX2qmusCqC1aUAAAQOxQJAsAAGKHBAUAAMQOCQoAAIgdEhQAABA7JCgAACB2SFAAAEDs/H/2rg3mPPpM4AAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "#And plot the log transform of both variables\n", + "loght=np.log(df.h1.astype('float'))\n", + "logQ=np.log(df.Qobs1.astype('float'))\n", + "\n", + "plt.figure(figsize=(9,6))\n", + "\n", + "plt.xlabel('log(stage (m))')\n", + "plt.ylabel('log(discharge (cms))')\n", + "plt.plot(loght,logQ,'k*')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Here, to further simplify, we will not worry about the break point, and we will look only at the flows above h1=0.54.\n", + "### Inspect the data\n", + "You can see above that even with a log transform, the data do not fall on a straight line. If I squint at the graph above, I can imagine two lines through this data, with a break about midway through the data.\n", + "\n", + "At this point, it's a good idea to look at your field site and see what physical evidence there is that the rating curve might have a break point. \n", + "\n", + "![Lyell Fork Tuolumne Plot](LyellFork_Tuolumne_flowcontrol.png)\n", + "\n", + "You can see from the photo above (and from the lecture notes), that at certain flow levels, this river stretch switches from being channel controlled to having a bedrock control, like a weir. From Open Channel Flow class (and from the lecture notes), we know that this changes the exponent in the stage-discharge equation, which would change the slope of our line here.\n", "\n", - "Identify the rows in our data frame that correspond to flows with h1 > h11, these correspond to data above the change in channel control." + "For simplicity, we will consider only the upper end of the rating curve, focusing on flows that are high enough to not be influenced by the pictured bedrock control." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "### PART 1: Linear Regression of Transformed Variables: \n", + "As illustrated above, due to the nature of our channel, we would need to separate our variables and fit two lines to it. For simplicity, we will focus here on the upper portion of the rating curve, looking only at higher flows. Work by CEE M.S. student Gwyn Perry, the transition between the two slopes is about 0.54. I also want to ignore two outlier measurements right around this transition, so I choose to look at data above a stage height of 0.59. You can change this cut-off value and see how it changes the results. Ideally, we would have a survey of the location that identifies the exact stage when the bedrock control becomes dominant." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "h11 = 0.59" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we identify the rows in our data frame that correspond to flows with h1 > h11, these correspond to data above the change in channel control." + ] + }, + { + "cell_type": "code", + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -606,6 +708,232 @@ "h_now = df.h1[df.h1 > h11]" ] }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAFzCAYAAAD7bpkSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnwklEQVR4nO3df5hddX3g8ffHAJ07wpRWRhJ+NVqgFa2AHRQX3YC6iui2aG0Xk9pibYdZqMVd24r1sdRku0+7VNZak4GpVu3WSro6imX9Uawo5qmhTpASMfXHJtbm4SaE9RfqjArz2T/umXSYzI8zkzlzz537fj3Pfeaec77n3E/OMwkfzv18P9/ITCRJkurmMe0OQJIkaTYmKZIkqZZMUiRJUi2ZpEiSpFoySZEkSbVkkiJJkmrpmHYHsFgnnXRSrl+/vt1hSJKkZbBr164HM7N/tmMdl6SsX7+esbGxdochSZKWQUT8y1zH/LpHkiTVkkmKJEmqJZMUSZJUSyYpkiSplkxSJElSLZmkSJKkWjJJkSRJtWSSIkmSaskkRZIk1ZJJiiRJXaTZbLJhwwYOHDjwqPd11HFt8SVJ0tJdd9113HnnnVx33XX09vayY8cONm/ezLZt29od2hEiM9sdw6IMDAyka/dIkrQ4jUaDiYmJecf09PQwPj5+xP5ms8kVV1zB9u3bWbt27bLGFRG7MnNgtmN+3SNJUheY76FEb28vmzZtYt++fbMe37Jly+EnLivJJEWSpC6wb98+zjzzzCP2/8iP/AgTExP09fUd8ZSk0WgQEQwPDzM5Ocnw8DARQaPRWJGYTVIkSeoC69at4+GHHwbguOOOA+CEE07grrvuYmhoaNbi2b1797Jx40Z6e3uBhZ+4LDcLZyVJ6hLnn38+l112GYODg4yMjNBsNjn33HPZunXrrOPXrVtHX18fExMT9PT0zPnEpSomKZIkdYnR0dHD7+dKTGY6ePAgQ0NDj0psVoqzeyRJUts4u0eSpA62XE3X6t68bSaTFEmSam65pgC3ayrxUvl1jyRJNTVXA7a5mq5VfZ0q+HWPJEkdaLmmALd7KvFSmaRIklRTyzUFuN1TiZfKJEWSpBqbmgK8c+fOOZuureR1VlJlNSkR0QPcCfwIrX4s78vM62eMCeBPgcuA7wFXZubd813XmhRJklaP+WpSqmzm9n3gOZn5nYg4FtgRER/JzJ3TxrwQOKt4PQMYLn5KkqQuV9nXPdnynWLz2OI187HNzwN/WYzdCZwYEeuqikmSJHWOSmtSImJNRNwDPADcnpl3zRhyKvCv07b3F/tmXmcwIsYiYuzQoUOVxStJkuqj0iQlMx/JzPOA04CnR8RTZgyJ2U6b5TojmTmQmQP9/f0VRCpJkupmRWb3ZOY3gU8Cl844tB84fdr2acD9KxGTJEmqt8qSlIjoj4gTi/cN4HnAP88Y9iHgV6LlQuBbmblyyytKkqTaqnJ2zzrg3RGxhlYy9DeZeVtEDAFk5k3Ah2lNP/4KrSnIr6wwHkmS1EEqS1Iy817g/Fn23zTtfQLXVBWDJEnqXHaclSRJtWSSIkmSaskkRZIk1ZJJiiRJqiWTFEmSVEsmKZIkqZZMUiRJUi2ZpEiSpFoySZEkSbVkkiJJkmrJJEWSJNWSSYokSaolkxRJklRLJimSJKmWTFIkSVItmaRIkqRaMkmRJEm1ZJIiSZJqySRFkiTVkkmKJEmqJZMUSVJtNZtNNmzYwIEDB9oditrAJEWSVFtbtmxhx44dbN68ud2hqA0iM9sdw6IMDAzk2NhYu8OQJFWo0WgwMTFxxP6enh7Gx8fbEJGqEhG7MnNgtmM+SZEk1c7evXvZuHEjvb29APT29rJp0yb27dvX5si0kkxSJEkrbqFak3Xr1tHX18fExAQ9PT1MTEzQ19dHZlqj0kVMUiRJK65MrcnBgwcZGhpi586dDA0NceDAAWtUuow1KZKkFbPUWhNrVFYva1IkSbWw1FoTa1S6U2VJSkScHhF3RMSeiLgvIq6dZcyPRsTfRsQ/FWNeWVU8kqT2m6vWZO3atZWcp85W5ZOUh4HXZuaTgAuBayLinBljrgG+kJnnAhcDb46I4yqMSZLUZrPVmlR5njrXitWkRMStwNsy8/Zp+14PnE4rWVkP3A6cnZmTc13HmhRJklaPttekRMR64HzgrhmH3gY8Cbgf2A1cO1uCEhGDETEWEWOHDh2qOlxJklQDlScpEXE88H7gNZn57RmHXwDcA5wCnAe8LSL6Zl4jM0cycyAzB/r7+yuOWJIk1UGlSUpEHEsrQXlPZo7OMuSVwGi2fAXYB/x0lTFJkqTOUOXsngDeAezJzBvnGPY14LnF+JOBnwL2VhWTJEnqHFU+SbkIeAXwnIi4p3hdFhFDETFUjNkC/LuI2A38PfC6zHywwpgkSUuwUBv7uc658MILeeYzn+lMHC1JZUlKZu7IzMjMp2bmecXrw5l5U2beVIy5PzOfn5k/k5lPycy/qioeSdLSLaUd/ZYtW7jrrrvYuXPnos5bSkKk1cm2+JKkOS2lHf1c5yx03pSrr76am2++mauuuopt27YtPmh1lLZPQZYkdaa9e/dy1llnHd4u045+7969vOQlL2HNmjWH961Zs4aXvvSl857XaDSICIaHh5mcnGR4eJiIoNFoLM8fRh3HJEWSNKtGo8Epp5zCl7/85cP7vve973HLLbfM245+3bp1nHzyyTzyyCOH9z3yyCOcfPLJ857n+jya6Zh2ByBJqqe9e/fy27/929xyyy1MTk7S09PD6aefzplnnrnguQcPHuQJT3gCF1xwAQCf/exnF6wxcX0ezWSSIkma1VTSAK1akh/84Ac873nPK1UnMjo6W2ushU2tzzM4OMjIyAjNZnNJ19HqYJIiSZrTSicN05ObrVu3VvpZqj9n90iSpLZxdo8kSeo4JimS1KXKNk2zuZraxSRFkrpU2S6yS+k2Ky0Ha1IkqcuU7SK7lG6z0mJZkyJJOqxs0zSbq6ndTFIkqcuUbZpmczW1m0mKJHWhqf4nO3fuZGhoaM6i2LLjpCpYkyJJktrGmhRJktRxTFIkSVItmaRIkqRaMkmRJEm1ZJIiSZJqySRFkiTVkkmKJEmqJZMUSZJUSyYpkiSplkxSJElSLZmkSJKkWjJJkSRJtWSSIkmSaqmyJCUiTo+IOyJiT0TcFxHXzjHu4oi4pxjzqarikSRJneWYCq/9MPDazLw7Ik4AdkXE7Zn5hakBEXEisA24NDO/FhGPrzAeSZLUQSp7kpKZzcy8u3j/ELAHOHXGsI3AaGZ+rRj3QFXxSJKkzrIiNSkRsR44H7hrxqGzgR+LiE9GxK6I+JU5zh+MiLGIGDt06FDF0UqSpDqoPEmJiOOB9wOvycxvzzh8DPCzwIuAFwBvjIizZ14jM0cycyAzB/r7+6sOWZIk1UCVNSlExLG0EpT3ZOboLEP2Aw9m5neB70bEncC5wJeqjEuSJNVflbN7AngHsCczb5xj2K3AsyPimIjoBZ5Bq3ZFkiR1uSqfpFwEvALYHRH3FPt+DzgDIDNvysw9EfFR4F5gEnh7Zn6+wpgkSVKHqCxJycwdQJQYdwNwQ1VxSJKkzmTHWUmSVEsmKZIkqZZMUiRJUi2ZpEiSpFpasHC2WE/nIuAUYBz4PDCWmZMVxyZJkrrYnElKRFwCXAf8OPA54AGgB7gc+MmIeB/w5lm6yEqSJB21+Z6kXAb8xtTif9NFxDHAi4H/QKujrCRJ0rKaM0nJzN+Z59jDwAerCEiSJAnK1aQMAM/m0TUpH8/Mr1ccmyRJ6mJzzu6JiCsj4m7g9UAD+CKtupRnAbdHxLsj4oyVCVOSJHWb+Z6kPBa4KDPHZzsYEecBZwFH1KxIkiQdrflqUrbOd2Jm3rPs0UiSJBUWbOYWEf8jIvoi4tiI+PuIeDAifnklgpMkSd2rTMfZ5xe9UF4M7AfOBuac+SNJkrQcyiQpxxY/LwPe66weSZK0Ehacggz8bUT8M63px1dHRD8wUW1YkiSp2y34JCUzrwOeCQxk5g+B7wI/X3VgkiSpu5Vp5raGVjO39UU7/Ck3VhaVJEnqeqW+7qH19c5uwJWPJUnSiiiTpJyWmU+tPBJJkqRpyszu+UhEPL/ySCRJkqYp8yRlJ/CBiHgM8EMggMzMvkojkyRJXa1MkvJmWrN7dmdmVhyPJEkSUO7rni8DnzdBkSRJK6nMk5Qm8MmI+Ajw/amdmekUZEmSVJkyScq+4nVc8ZIkSarcgklKZr5pJQKRJEmabsGalIi4PSJOnLb9YxHxsUqjkiRJXa9M4Wx/Zn5zaiMzvwE8fqGTIuL0iLgjIvZExH0Rce08Yy+IiEci4mWlopakRWg2m2zYsIEDBw60OxRJi1AmSXkkIs6Y2oiInwDKzPR5GHhtZj4JuBC4JiLOmTmoWBvojwGfzkiqxJYtW9ixYwebN29udyiSFqFM4ewbgB0R8ali+98DgwudlJlNWjODyMyHImIPcCrwhRlDXw28H7igbNCSVEaj0WBiYuLw9vDwMMPDw/T09DA+Pt7GyCSVseCTlMz8KPA0YDvwN8DPZuainnpExHrgfOCuGftPBV4C3LTA+YMRMRYRY4cOHVrMR0vqYp/5zGfo7++n0WgA0Nvby6ZNm9i3b1+bI5NUxpxJSpFYAJCZD2bmbZn5t5n5YHE8IuK0hT4gIo6n9aTkNZn57RmH3wK8LjMfme8amTmSmQOZOdDf37/QR0oSACMjIxw6dIjx8XF6enqYmJigr6+PtWvXtjs0SSXM93XPDcV6PbcCu4BDQA9wJnAJ8FzgemD/XBeIiGNpJSjvyczRWYYMALdEBMBJwGUR8XBmfnDxfxRJapn5NQ/AxMQEa9assXhW6iBzPknJzF8E3gj8FLAV+DSthOXXgS8Cz8nM2+c6P1qZxzuAPXN1p83MJ2Tm+sxcD7wPuNoERdLR2rt3Lxs3bqS3txf4t6959u/fz+jobP+/JKmO5i2czcwv0CqcXYqLgFcAuyPinmLf7wFnFNeetw5FkpZq3bp19PX1MTEx4dc8UgcrM7tnSTJzBxCLGH9lVbFI6j4HDx5kaGiIwcFBRkZGaDab7Q5J0iJFpy1uPDAwkGNjY+0OQ5IkLYOI2JWZA7MdK9PMTZIkacWVWbsnIuKXI+L3i+0zIuLp1YcmSZK6WZknKduAZwIvL7YfojXbR5KWzPV0JC2kTJLyjMy8BpiAwwsMHldpVJJWPdfTkbSQMknKD4tFABMgIvqByUqjkrRqNRoNIoLh4WEmJycZHh4mIg63rpekKWWSlLcCHwAeHxF/COwA/nulUUlateZqtOZ6OpJmWrBPSma+JyJ20WqDH8Dlmbmn8sgkrUo2WpNUVpnZPT8OPAC8F/hr4GCxJo8kLclUo7WdO3cyNDRk8aykWS3YzC0ivgqcDnyD1pOUE4EmrcTlNzJzV7UhPprN3KTu1Ww2ueKKK9i+fbtPXqRV4mibuX0UuCwzT8rMxwEvBP4GuJrW9GRJWhHOCJK6S5knKWMzM5ypfRFxT2aeV2WAM/kkReo+jUaDiYmJI/b39PQwPj7ehogkLZejfZLy9Yh4XUT8RPH6XeAbxbRkpyJLqpwzgqTuVCZJ2QicBnyweJ1e7FsD/FJVgUnSFGcESd1p3iSleFrylsx8dWaeX7xenZmHMvMHmfmVFYpTUodarvb3zgiSuk+ZmpSPAf8xM3+wMiHNz5oUqbNcffXV3HzzzVx11VVs22atvaRHm68mpUyScjPwNOBDwHen9mfmjcsZZFkmKVJnsNhVUhlHWzh7P3BbMfaEaS9JmpPFrpKOVpm2+G9aiUAkrS4Wu0o6WgsmKcWqx78LPBnomdqfmc+pMC5Jq8BUsevg4CAjIyM0m812hySpgyyYpADvAbYDLwaGgF8FDlUZlKTVYXR09PD7rVu3tjESSZ2oTE3K4zLzHcAPM/NTmflrwIUVxyVJkrpcmScpPyx+NiPiRbQKaU+rLiRJkqRyScp/i4gfBV4L/BnQB/yXSqOSJEldr8zsntuKt98CLqk2HEmSpJays3t+A1g/fXxRmyJJklSJMl/33Ap8Gvg48Ei14UiSJLWUSVJ6M/N1lUciaUU0m02uuOIKtm/fbmM1SbVWZgrybRFx2WIvHBGnR8QdEbEnIu6LiGtnGbMpIu4tXv8QEecu9nMkLc6WLVvYsWMHmzdvbncokjSvORcYjIiHgAQCeCzwfVrTkQPIzOyb98IR64B1mXl3RJwA7AIuz8wvTBvz74A9mfmNiHgh8AeZ+Yz5rusCg9LSuOCfpDpa0gKDmXlCZvYVPx+TmY1p2/MmKMX5zcy8u3j/ELAHOHXGmH/IzG8Umzux/4pUGRf8k9RpFvy6JyJeUvRJmdo+MSIuX8yHRMR64HzgrnmGvQr4yBznD0bEWESMHTpkR35pKWYu+Dc+Ps4nPvGJdoclSXMqU5NyfWZ+a2ojM78JXF/2AyLieOD9wGsy89tzjLmEVpIya4FuZo5k5kBmDvT395f9aEkzTC34t3PnTs455xyazaa1KZJqa86alMMDIu7NzKfO2Lc7M39mwYtHHAvcBnwsM2+cY8xTgQ8AL8zMLy10TWtSpKNjbYqkOllSTco0YxFxY0T8ZEQ8MSL+J60i2IU+NIB30CqMnStBOQMYBV5RJkGRdPSsTZHUKcr0SXk18EZgO62ZPX8HXFPivIuAVwC7I+KeYt/vAWcAZOZNwO8DjwO2tXIaHp4rm5K0PGbWpkxMTNDX12fPFEm1U2btnu8C1wFExBrgscW+hc7bQSupmW/MrwO/Xi5USctlqjZlcHCQkZERms1mu0OSpCOUqUn5a2CIVkv8XcCPAjdm5g3Vh3cka1IkSVo9jrYm5ZxiVs7lwIdpfV3ziuULT5Ik6UhlkpRji1k6lwO3ZuYPaXWilSRJqkyZJOVm4Ku0WuPfGRE/Acza70RSdZrNJhs2bODAgQPtDkWSVsSCSUpmvjUzT83My7LlX4BLViA2SdO4MKCkbjPfAoO/nJl/FRH/dbbjc/U+qZqFs+o2Nl+TtJottXD2scXPE+Z4SVoBNl+T1K3m7JOSmTcXP9+0cuFImsnma5K61ZxJSkS8db4TM/O3lj8cSbOx+ZqkbjRfx9mp9XkuAs6h1RYf4BcpsXaPpOUzOjp6+P3WrVvbGIkkrZz5vu55N0BEXAlcUvRHISJuorV+jyRJUmXK9Ek5hUcXyh5f7JMkSapMmVWQ/wj4XETcUWxvAP6gsogkSZIotwryOyPiI8Azil3XZaYtLyVJUqXmm92zPjO/ClAkJbfOOB7AqZm5v9IIJUlSV5rvScoNEfEYWsnJLuAQ0AOcSast/nOB6wGTFEmStOzmm93zixFxDrAJ+DVgHTAO7AH+D/CHmXlkr25JkqRlMG9NSmZ+AXjDCsUiSZJ02IKFsxHx0ll2fwvYnZkPLH9IkiRJ5aYgvwp4JjA1BfliYCdwdkRszsz/VVFskiSpi5VJUiaBJ2XmQYCIOBkYpjUl+U7AJEWSJC27Mh1n108lKIUHgLMz8+vAD6sJS5IkdbsyT1I+HRG3Af+72H4ZcGdEPBb4ZlWBSZKk7lbmSco1wDuB84DzgXcD12TmdzPzkgpjkzpas9lkw4YNHDhgg2ZJWooFk5TMTGAH8Ang48CdxT5J89iyZQs7duxg8+bN7Q5FkjpSLJRvRMQvATcAnwQCeDbwO5n5vsqjm8XAwECOjY2146OlUhqNBhMTR/Y57OnpYXx8vA0RSVJ9RcSuzByY7ViZr3veAFyQmb+amb8CPB1443IGKK0me/fuZePGjfT29gLQ29vLpk2b2LdvX5sjk6TOUiZJecyMpm3/r+R5Uldat24dfX19TExM0NPTw8TEBH19faxdu7bdoUlSRymTbHw0Ij4WEVdGxJW01u358EInRcTpEXFHROyJiPsi4tpZxkREvDUivhIR90bE0xb/R9BqslqKTQ8ePMjQ0BA7d+5kaGio4/88ktQOC9akAETELwAX0apJuTMzP1DinHXAusy8OyJOoLWS8uXFekBTYy4DXg1cRqs53J9m5jPmu641Kavb1Vdfzc0338xVV13Ftm3b2h3OvJrNJldccQXbt2/3KYkkLdF8NSmlkpRlCuJW4G2Zefu0fTcDn8zM9xbbXwQuzszmXNcxSVmdOrHYtJMSKkmqqyUVzkbEQxHx7VleD0XEtxcZwHpaPVbumnHoVOBfp23vL/apy3RSsWmj0SAiGB4eZnJykuHhYSKCRqPR7tAkaVWZM0nJzBMys2+W1wmZ2Vf2AyLieOD9wGsyc2ZyE7N99CzXGIyIsYgYO3ToUNmPVgfppGLTTkqoJKmTVTpLJyKOpZWgvCczR2cZsh84fdr2acD9Mwdl5khmDmTmQH9/fzXBqu06pdi0kxIqSepkZdbuWZKICOAdwJ7MvHGOYR8CfjMibqFVOPut+epRtLqNjo52TDHqVEI1ODjIyMgIzaa/tpK03CornI2IZwGfBnYDk8Xu3wPOAMjMm4pE5m3ApcD3gFdm5rxVsRbOrm4Wo0pSd6nF7J7lYpKyOnXi7B5J0tE72rb4UuUsRpUkzWSSolpYbcWoq6VzriS1k0mKaqNTZveUsWXLFnbs2MHmzZvbHYokdSxrUqRlZG2NJC2ONSnSCrG2RpKWj0mKtIxWW22NJLWTSYq0zFZTbY0ktZM1KZIkqW2sSZEkSR3HJEXCviaSVEcmKVp1lpJw2NdEkurHJEWrzmISjkajQUQwPDzM5OQkw8PDRASNRmMFIpUkzcckRavGUhIO+5pIUn2ZpGjVWErCYV8TSaovkxStGktNOOxrIkn1dEy7A5CW01TCMTg4yMjICM1mc8FzRkdHD7/funVrleFJkhbBZm6SJKltbOYmSZI6jkmKJEmqJZMUSZJUSyYpkiSplkxSJElSLZmkqGO5KKAkrW4mKepYLgooSaubfVLUcRqNBhMTE0fs7+npYXx8vA0RSZKWyj4pWlVcFFCSuoNJijqOiwJKUncwSVFHclFASVr9KqtJiYi/AF4MPJCZT5nl+I8CfwWcQWuhwz/JzHcudF1rUjpXs9nkiiuuYPv27T71kCQB7atJeRdw6TzHrwG+kJnnAhcDb46I4yqMR23mbBxJ0mJUlqRk5p3A1+cbApwQEQEcX4x9uKp41D6NRoOIYHh4mMnJSYaHh4kIGo1Gu0OTJNVYO2tS3gY8Cbgf2A1cm5mTbYxHFXE2jiRpKdqZpLwAuAc4BTgPeFtE9M02MCIGI2IsIsYOHTq0chFqWTgbR5K0FO1MUl4JjGbLV4B9wE/PNjAzRzJzIDMH+vv7VzRILQ9n40iSFuuYNn7214DnAp+OiJOBnwL2tjEeVWh0dPTw+61bt7YxEklSp6gsSYmI99KatXNSROwHrgeOBcjMm4AtwLsiYjcQwOsy88Gq4pEkSZ2lsiQlM1++wPH7gedX9fmqlj1PJElVs+OslsSeJ5KkqrkKshbFFYglScvJVZC1bOx5IklaKSYpKm2qDuWYY46x54kkqXImKSptqg7l05/+tD1PJEmVsyalQ7RzNs1S6lCc/SNJKsOalFWgnbNpllKH4uwfSdLRMkmpuTqsILyYtXfqEK8kaXUwSam5usymKbv2Tl3ilSR1vnau3aMS6rKCcNm1d+oSrySp8/kkpQN02grCnRavJKmenN0jSZLaxtk9kiSp45ikSJKkWjJJ0ZI0m002bNhgvYkkqTImKVoSm7VJkqpm4awWZSkt8iVJmouFs1o2NmuTJK0UkxQtis3aJEkrxSRFi2azNknSSrAmRZIktY01KZIkqeOYpEiSpFoySZEkSbVkkiJJkmrJJEWSJNWSSYokSaolkxRJklRLlSUpEfEXEfFARHx+njEXR8Q9EXFfRHyqqlgkSVLnqfJJyruAS+c6GBEnAtuAn8vMJwO/WGEsC2o2m2zYsKGt3VPrEIMkSXVRWZKSmXcCX59nyEZgNDO/Vox/oKpYytiyZQs7duxg8+bNXR2DJEl1UWlb/IhYD9yWmU+Z5dhbgGOBJwMnAH+amX+50DWXuy1+o9FgYmLiiP09PT2Mj48v2+fUPQZJktqhrm3xjwF+FngR8ALgjRFx9mwDI2IwIsYiYuzQoUPLGsTevXvZuHEjvb29APT29rJp0yb27du3rJ9T9xgkSaqbdiYp+4GPZuZ3M/NB4E7g3NkGZuZIZg5k5kB/f/+yBrFu3Tr6+vqYmJigp6eHiYkJ+vr6WLt27bJ+Tt1jkCSpbtqZpNwKPDsijomIXuAZwJ52BHLw4EGGhobYuXMnQ0NDbSlcrUMMkiTVSWU1KRHxXuBi4CTgIHA9rRoUMvOmYszvAK8EJoG3Z+ZbFrructekSJKk9pmvJuWYqj40M19eYswNwA1VxSBJkjqXHWclSVItmaQsoE4N1uoUiyRJVTNJWUCdGqzVKRZJkqpWaTO3KqxU4WydGqzVKRZJkpZTXZu51VqdGqzVKRZJklaKScoc6tRgrU6xSJK0UkxS5lGnBmt1ikWSpJVgTYokSWoba1IkSVLHMUmRJEm1ZJIiSZJqySRFkiTVkkmKJEmqJZMUSZJUSyYpNeDCgZIkHckkpQZcOFCSpCPZzK2NXDhQktTtbOZWUy4cKEnS3ExS2siFAyVJmptJSpu5cKAkSbOzJkWSJLWNNSmSJKnjmKRIkqRaMkmRJEm1ZJIiSZJqySRFkiTVkkmKJEmqJZMUSZJUS5UlKRHxFxHxQER8foFxF0TEIxHxsqpikSRJnafKJynvAi6db0BErAH+GPhYhXFIkqQOVFmSkpl3Al9fYNirgfcDD1QVhyRJ6kxtq0mJiFOBlwA3tSsGSZJUX8e08bPfArwuMx+JiHkHRsQgMFhsficivrhMMZwEPLhM11rNvE/lea/K8T6V430qx/tUTl3v00/MdaDSBQYjYj1wW2Y+ZZZj+4Cp7OQk4HvAYGZ+sLKAjoxhbK5FjfRvvE/lea/K8T6V430qx/tUTifep7Y9ScnMJ0y9j4h30UpmPtiueCRJUr1UlqRExHuBi4GTImI/cD1wLEBmWociSZLmVVmSkpkvX8TYK6uKYwEjbfrcTuN9Ks97VY73qRzvUznep3I67j5VWpMiSZK0VLbFlyRJtdRVSUpE/HhE3B4RXy5+/tgc474aEbsj4p6IGFvpONut7H0qxq6JiM9FxG0rGWNdlLlXEdETEf8YEf8UEfdFxJvaEWs7lbxPp0fEHRGxp7hP17Yj1nZaxL9RpZYdWW0i4tKI+GJEfCUirpvleETEW4vj90bE09oRZ7uVuE8/HRGfiYjvR8RvtyPGsroqSQGuA/4+M88C/r7YnsslmXlep03XWiaLuU/XAntWJKp6KnOvvg88JzPPBc4DLo2IC1cuxFooc58eBl6bmU8CLgSuiYhzVjDGOij7d+9dLLDsyGpTLKOyFXghcA7w8ll+P14InFW8BoHhFQ2yBkrep68DvwX8yQqHt2jdlqT8PPDu4v27gcvbF0qtlbpPEXEa8CLg7SsTVi0teK+y5TvF5rHFq9uKwcrcp2Zm3l28f4hW8nvqSgVYE6X+7pVcdmS1eTrwlczcm5k/AG6hdb+m+3ngL4u/czuBEyNi3UoH2mYL3qfMfCAzPwv8sB0BLka3JSknZ2YTWv8gAo+fY1wCfxcRu4put92m7H16C/C7wOQKxVVHpe5V8bXYPbTWqbo9M+9auRBroezvFHC4EeT5gPdJU04F/nXa9n6OTGLLjFntVtU9aGdb/EpExMeBtbMcesMiLnNRZt4fEY8Hbo+Ify7+z2XVONr7FBEvBh7IzF0RcfEyhlY7y/E7lZmPAOdFxInAByLiKZm5quoJlunvHhFxPK2FR1+Tmd9ejtjqZLnuUxeabf2UmU8ky4xZ7VbVPVh1SUpmPm+uYxFxMCLWZWazeAQ46+rLmXl/8fOBiPgArcdnqypJWYb7dBHwcxFxGdAD9EXEX2XmL1cUctssx+/UtGt9MyI+SaueYFUlKctxnyLiWFoJynsyc7SiUNtqOX+fusx+4PRp26cB9y9hzGq3qu5Bt33d8yHgV4v3vwrcOnNARDw2Ik6Yeg88n1X2H5MSFrxPmfn6zDwtM9cDVwCfWI0JSgllfqf6iycoREQDeB7wzysVYE2UuU8BvAPYk5k3rmBsdbLgfepinwXOiognRMRxtP7d+dCMMR8CfqWY5XMh8K2pr8+6SJn71Dkys2tewONoVcx/ufj548X+U4APF++fCPxT8boPeEO7467jfZox/mJaay+1PfY63ivgqcDngHtpJby/3+64a3qfnkXrsfS9wD3F67J2x163+1Rsvxdo0ip83A+8qt2xr9D9uQz4EvB/p/5tBoaAoeJ90JrZ8n+B3cBAu2Ou6X1aW/zefBv4ZvG+r91xz/ay46wkSaqlbvu6R5IkdQiTFEmSVEsmKZIkqZZMUiRJUi2ZpEiSpFoySZH0KBHxnYVHzXv++yLiifMcvzIiTjmazygZRyMiPlUsuFb2nN+MiFdO2/6TiHhONRFKWohJiqRlExFPBtZk5t55hl1Jq+9H1X4NGM3WkgRl/QWt1WGn/BnzrwIuqUImKZJmVXTtvCEiPh8RuyPiPxX7HxMR2yLivoi4LSI+HBEvK07bRNEltVhU8V3Tzv8vxbgB4D0RcU/xtOP3I+KzxbiRovMsEXFBRNwbEZ+ZimPadW8ozrk3Iq6a448wPZaLi6cqfxMRX4qIP4qITRHxj0VsPwmQmd8DvhoRTy+2/wV4XETMttaOpIqZpEiay0uB84BzabXyv6FYT+alwHrgZ4BfB5457ZyLgF3F+/OAUzPzKZn5M8A7M/N9wBiwKTPPy8xx4G2ZeUFmPgVoAC8uzn8nrQ6ZzwSmPw15Fa125xcAFwC/ERFPmB540Q78iZn51Wm7zwWuLeJ+BXB2Zj4deDvw6mnjxoBnT9u+u/hzSVphJimS5vIs4L2Z+UhmHgQ+RSspeBbwvzNzMjMPAHdMO2cdcKh4vxd4YkT8WURcSqsF92wuiYi7ImI38BzgycVaRydk5j8UY/562vjn01qf5R7gLlqt5M+acc2TaLX7nu6zmdnMzO/Tahf+d8X+3bSSrikP8Oivo2ZuS1ohq24VZEnLZrYl3+fbDzBOa1VsMvMbEXEu8ALgGuCXaNWJ/NuFInqAbbTWWPnXiPiD4vz5PiOAV2fmx8rEMc33p72fnLY9yaP/Lewpzp9rW9IK8UmKpLncCfynogakH/j3wD8CO4BfKGpTTqa1wOSUPcCZABFxEvCYzHw/8EbgacWYh4ATivdTicSDEXE88DJoJTjAQ8VKttBayXXKx4D/HBHHFp9zdrFi+WHF+WuKJGixzubRK5/P3Ja0QnySImkuH6BVb/JPtFYn/t3MPBAR7weeS+s/3F+i9ZXLt4pz/g+tpOXjwKnAOyNi6n+GXl/8fBdwU0SMF9f/c1pfuXyV1jLzU14F/HlEfBf45LTPeDutr2fuLopsDwGXzxL/39H6aurji/xzXwS8CaBIhM6kVaciaYW5CrKkRYuI4zPzOxHxOFpPVy4qEpgGrRqVixY59XfOzyjeXwesy8xrF3H++cB/zcxXLPWciHgJ8LTMfOPiope0HHySImkpbiuKW48DthQFtGTmeERcT+spyteO8jNeFBGvp/Xv1L/Q6q9SWmZ+LiLuiIg1i0iYTqL11dSUY4A3L+ZzJS0fn6RIkqRasnBWkiTVkkmKJEmqJZMUSZJUSyYpkiSplkxSJElSLZmkSJKkWvr/syW/WuxhmY8AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "#And plot the log transform of both variables\n", + "loght=np.log(h_now.astype('float'))\n", + "logQ=np.log(Qobs_now.astype('float'))\n", + "\n", + "plt.figure(figsize=(9,6))\n", + "\n", + "plt.xlabel('log(stage (m))')\n", + "plt.ylabel('log(discharge (cms))')\n", + "plt.plot(loght,logQ,'k*')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With linear regression, we can only solve for two of the three unknowns in our rating curve equation:\n", + "$Q = a(h-b)^c$\n", + "which we have log transformed to be\n", + "$log(Q) = log(a) + c*log(h-b)$\n", + "\n", + "First, we will assume that b is 0.28 m. In practice, this is often estimated from field surveys as the maximum height of water in the measuring pool when flow stops.\n", + "\n", + "We then use the same code from basic linear regresssion (Lab 4.3), where our calculated slope will be c, and our calculated intercept will be log(a)." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "B0 : 2.8277\n", + "B1 : 2.6902\n" + ] + } + ], + "source": [ + "x=loght\n", + "y=logQ\n", + "n = len(x)\n", + "\n", + "B1 = ( n*np.sum(x*y) - np.sum(x)*np.sum(y) ) / ( n*np.sum(x**2) - np.sum(x)**2 ) # B1 parameter, slope\n", + "B0 = np.mean(y) - B1*np.mean(x) # B0 parameter, y-intercept\n", + "\n", + "print('B0 : {}'.format(np.round(B0,4)))\n", + "print('B1 : {}'.format(np.round(B1,4)))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# Now, how do we find 95% confidence intervals? We do this for our estimates of logQ\n", + "# Again, borrowing from Lab 4.3\n", + "y_predicted = B0 + B1*x\n", + "residuals = (y - y_predicted)\n", + "# sum of squared errors\n", + "sse = np.sum(residuals**2)\n", + "# standard error of regression\n", + "s = np.sqrt(sse/(n-2))\n", + "\n", + "# create an array of x values\n", + "p_x = np.linspace(x.min(),x.max(),100)\n", + "\n", + "# using our model parameters to predict y values\n", + "p_y = B0 + B1*p_x\n", + "\n", + "# calculate the standard error of the predictions\n", + "sigma_ep = np.sqrt( s**2 * (1 + 1/n + ( ( n*(p_x-x.mean())**2 ) / ( n*np.sum(x**2) - np.sum(x)**2 ) ) ) )\n", + "\n", + "# our chosen alpha\n", + "alpha = 0.05\n", + "\n", + "# compute our degrees of freedom with the length of the predicted dataset\n", + "n = len(p_x)\n", + "dof = n - 2\n", + "\n", + "# get the t-value for our alpha and degrees of freedom\n", + "t = st.t.ppf(1-alpha/2, dof)\n", + "\n", + "# compute the upper and lower limits at each of the p_x values\n", + "p_y_lower = p_y - t * sigma_ep\n", + "p_y_upper = p_y + t * sigma_ep" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAFzCAYAAAD2cOlVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABnCklEQVR4nO3dd3iUxRbH8e8QCEmA0EvoIL1IR8SCCCoiqKBeEewFYy/YCwjYxYqIYu8dFZCuICIggoLSS0IPvSWkJ3P/mI0ETMIm74YUfp/nycO+u+/OHvdy4TBz5oyx1iIiIiJSVJQo6ABEREREckPJi4iIiBQpSl5ERESkSFHyIiIiIkWKkhcREREpUpS8iIiISJFSsqADCKQqVarY+vXrF3QYIiIiEgCLFy/eba2tevTzxSp5qV+/PosWLSroMERERCQAjDEbs3pey0YiIiJSpCh5ERERkSJFyYuIiIgUKUpeREREpEhR8iIiIiJFipIXERERKVLyLXkxxoQYYxYaY5YaY5YbY4Zncc8gY8zfvp95xpg2mV7bYIz5xxizxBij/c8iIiIC5G+flyTgbGttnDGmFDDXGDPFWrsg0z3RQDdr7T5jzPnAOOCUTK93t9buzscYRUREpIjJt+TFWmuBON9lKd+PPeqeeZkuFwC18yseERERKR7ytebFGBNkjFkC7ARmWGt/z+H2G4Apma4tMN0Ys9gYMziHzxhsjFlkjFm0a9eugMQtIiIihVe+Ji/W2jRrbVvcjEpnY0yrrO4zxnTHJS8PZnr6NGtte+B84DZjzJnZfMY4a21Ha23HqlX/c/yBiIiIFDPHZbeRtXY/MBvodfRrxpiTgXeAi6y1ezK9Z5vv153Ad0Dn4xGriIiIFG75uduoqjGmgu9xKNATWHXUPXWB8cBV1to1mZ4vY4wpl/EYOBdYll+xioiISC4kJcF998GHHxbIx+fnbqMI4ENjTBAuSfrKWjvJGBMJYK19ExgKVAbeMMYApFprOwLVge98z5UEPrPWTs3HWEVERMRfwcGwcCGUzM80InvGbQoqHjp27GgXLVJLGBERkYBbvhyGDoX334fwcEhOdklMPjLGLPZNahxBHXZFRETk2A4dgrlzYeVKd53PiUtOCma+R0RERAq/996DvXtdfUvnzrBxI4SEFHRUmnkRERGRbMyaBVOnQnq6uy4EiQsoeREREZEMe/bArbe6GRaAN9+EGTOgROFKFwpXNCIiIlJwDh2CTz+FX39112XKgNv5W6io5kVERORE9ttvMHMmDBsGdevCpk1QvnxBR5UjzbyIiIicyKZOhXffhf373XUhT1xAyYuIiMiJJSUFXnrJNZkDeOQRt/25QoUCDSs3tGwkIiJyIklMhFGjICbGbX8ODS3oiHJNMy8iIiLF3ZYtrqbFWihXDhYvhuef9zzsxj2HWBlzMAAB5o6SFxERkeJuxgyXrCzznXEcEeFpF1F8ciqjpq3mnJfmMOyH5QEK0n9aNhIRESmOpk51My3nnw/XXAM9e0KdOp6GtNYyZdl2npy0gm0HEunXrhYPnd8sQAH7T8mLiIhIcZOe7gpxK1VyyUuJEp4Tl7U7Yhk2YTnz1u+heUQ4r17Rjk71KwUo4NxR8iIiIlIcJCTAmDFw222uCPf776FGDc/Dxiam8OrMtXwwbwNhwUGMvKglA0+pR1CJgmtep+RFRESkOFi4EO6/H+rXh0svdQ3nPEhPt3z311aembKKPYeSGNCpLvef15RKZQruNOkMSl5ERESKqjVr4O+/XbLSrZsryG3Z0vOwy7YeYNiE5SzeuI+2dSrw3rUdObl2Be/xBoiSFxERkaLqkUdg/nzo2xdKl/acuOw7lMyo6av5bOEmKpcJ5vlLT+bS9rUpUYBLRFlR8iIiIlJUWAtffeVmWWrUgFdfhaAgl7h4kJZu+XzhJkZNX01sYirXdq3P3T2bUD60VIACDywlLyIiIkXFpk1w1VXw0EMwYgTUquV5yMUb9zL0h+Us33aQLg0rMfzCVjStUS4AweYfNakTEREpzA4cgK+/do/r1YNff3Xdcj3aGZvIvV8t4ZKx89kTl8zI8+uz7dOHKM8hYmJi6NatG9u3b/f8OflBMy8iIiKF2fPPu58uXVyvllNO8TRcSlo6H87bwCsz15KUmsatZ53Ebd0bcevgG/h1zhweeughwsLCmDt3LiNGjOCNN94I0H9I4BhrbUHHEDAdO3a0ixYtKugwREREvPnzTwgLg2bN3MzLunXQoYPnYX9bt5snJixn7c44ujWpyrC+LWhZtyqJiYk5vi8kJISEhIT/PB8TE8OAAQP48ssvqRGAnjJHM8YsttZ2PPp5LRuJiIgUJomJrivuo4+66/LlPScuW/cncNunfzLond9JTE1j3FUd+OC6TjSsWpacJjHCwsIYNGgQ0dHRWb4+cuTIf2dojictG4mIiBS0tDSYPBn69IGQEBg/PiD9WhJT0njn1yjGzFpPurXce04TBp/ZkJBSQf/eEx0dzZlnnsm6deuOeG/p0qVJTEwkPDz8P7MqZUNCuDMpiUNAOjB27FjGjh2b7QxNoGnmRUREpKB9/jlceCHMmuWuTzsNKlTwNORPK3dw3itzGDV9Dd2aVOWnId24s0fjIxIXgIiICFJTUwEIDnbdc8uVK8fvv/9OZGRklkW7a6OiuLpqVboFubGONUMTaJp5ERERKQg7d8LWrdCuHQwYAGXLQvfunofdsPsQIyat4OdVOzmpahk+vqEzZzSumuN72rVrR+/evRk8eDDjxo0jJiaGNm3aMGbMmMM3RUW5XU5jxhBRsybPXHghY95/n5CQkGxnaPKLCnZFREQKQpcuEBsL//zjTn32KD45lTGz1vH2nGhKBRnu7tmEa7rWJ7hkgBZZFi2Cnj3dgY9nnUX//v2JiIg4IuEZP358YD7LJ7uCXSUvIiIix8u8edCxIwQHu4MUw8PdjqI8iomJ4fIBA7h55Ju8/ts2Yg4k0q9dLR4+vxnVwkNyNU6Wu4a++AK2bYN773XXsbFQ7vg1sNNuIxERkYK0eLGrZXn7bXfdubOnxAXggadeZk3N83h0chQVwoL5OvJUXr68ba4SF8hh19CUKfDNN66gGI5r4pITzbyIiIjkl+RkWL7c1bVYCx99BJdd5nq4eBBWvjIhHftTrn0f0pMT2P/rx8QtmUpI6eBc7fYJDQ09osdLeeAJYGxwMKuTktxMS1iYOz+pAGjmRURE5Hi77TY4+2zXaM4YuOYaT4lLerrlm8VbaDbkM8I7XkTiyllse/tm0lfPZtDAK3K92ycqKoqBAwcS5oupakgIg0uVYuEzz7gbypUrsMQlJ0peREREAmnjRti3zz0eMgQ++cQ1mvNo2dYDXPrmPO77ein1qpSla9xcdk8ZTbBNzvNun4iICE5OSOCB+HhCQkKISk5m6JVXUj6jxqWQyrfkxRgTYoxZaIxZaoxZbowZnsU9xhjzmjFmnTHmb2NM+0yv9TLGrPa99lB+xSkiIhIwe/dC69Yw3PdXXrNmcMEFnobcdyiZR777h76vz2XT3nheuPRkvo3sSmLMGiIjI1mwYEG2/Vj8UX/lSu4qU4Y/pk0jMjKSqP37PcV7PORbzYsxxgBlrLVxxphSwFzgLmvtgkz39AbuAHoDpwCvWmtPMcYEAWuAc4AtwB/AFdbaFTl9pmpeRESkQKxcCc2bu8cffOCWiurW9TRkWrrls4WbGDVtNXFJqVxzan3uPqcx4SGlvMWakgKvvw6nnuq2ayckQGpqoSnGzSy7mpd8a1JnXVYU57ss5fs5OlO6CPjId+8CY0wFY0wEUB9YZ62N8gX/he/eHJMXERGR4270aLeV+J9/3EzLtdd6HnLRhr0Mm7Cc5dsOcmrDygy/qCVNqgcouUhKghdfhMsvd8lLaGhgxj2O8rXDrm8GZTHQCBhjrf39qFtqAZszXW/xPZfV81meAW6MGQwMBqjrMcsVERHxy6FDbidOjRquO25qKjRs6HnYnQcTeXbKKsb/tZWI8iG8PrAdF7SOwC1meLBpE4wbByNGuE6+ixZB9eqe4y0o+Zq8WGvTgLbGmArAd8aYVtbaZZluyep/DZvD81l9xjhgHLhlI28Ri4iIHEN6OpxyCtSvD5MmQdWqcM89noZMSUvnw3kbeGXmWpJT07mt+0nc1r0RYcEB+mt61ix46SU329K6tUu6irDjcraRtXa/MWY20AvInLxsAepkuq4NbAOCs3leRESkYGzZArVru1b+jzziuaYlw2/rdjNswnLW7Yyje9OqDOvbkvpVyngb1FqYMMHF2rcvXHWVa+1fq1ZAYi5o+bnbqKpvxgVjTCjQE1h11G0TgKt9u466AAestTG4At3GxpgGxphgYIDvXhERkePvp5+gQQOYMcNdDxwIp5/uacit+xO49dPFDHrnd5JT03nn6o68f11n74kLuORl5EhXmAsuiSkmiQvk78xLBPChr+6lBPCVtXaSMSYSwFr7JjAZt9NoHRAPXOd7LdUYczswDQgC3rPWLs/HWEVERI5krTv5uXp1l6jcdx+0bet52MSUNN6eE8WY2esAGHJOE246syEhpTw2g4uNhVdfdcXDYWHwww9QrZrneAsjHQ8gIiKSlUGD3A6iP/+EkoH5t/5PK3cwfOIKNu2Np3frGjx6QQtqVQjQbp85c+Css+Dbb6Ffv8CMWcCO+1ZpERGRImffPnfSc1CQWxravt0tuXgUvfsQIyYuZ9bqXTSqVpZPbjiF0xtX8R7v0qXu7KSBA+HMM2H1amjc2Pu4hZySFxEREYDoaHfS81NPweDBnjvjAsQnpzJm1jrenhNNcMkSPHZBc67pWp9SQQEqOX3ySfjjD7j0UggOPiESF1DyIiIiJ7qDB91sS/36blfOKVm2FcsVay2T/o7h6ckriTmQSP/2tXioVzOqhYd4GzgtDd57D/r0gYgI1yCvdGmXuJxAlLyIiMiJa9Qo1/9k5Up3eOJLL3kecvX2WIZNWMaCqL20iAhn9BXt6Fi/UgCCxTWbu/122LXLbdcu4v1a8krJi4iInFhSU91PSAh07x6wupaDiSm8MmMtH87fQNnSJRl5cSsGdq5LUAmP3XG3b4epU92xAw0auO64rVp5jrcoU/IiIiInjoQE6NoVzjsPnn0WOnRwPx6kp1u+/XMLz01dxZ5DyQzoVJf7z2tKpTIBWsp59VV45RXo1cvNtLRuHZhxizAlLyIiUvwlJrqZltBQlwR06hSQYf/ZcoChE5bx16b9tKtbgfev7Uzr2uW9DzxzpktUWrWChx+G668/YZeIsqI+LyIiUrz9+CNcdx0sWBCQwxMB9h5K5oVpq/nij01ULhPMQ+c3p3+7WpTwukQE7tDH+vXd7NAnn3gfrwhTnxcRETmxJCe7XTht2kC3bq53i0dp6ZbPft/IqOlriEtK5bquDbj7nMaEh5TyNnBCAnzxhatrKVMGpk+H5s09x1tcKXkREZHixVq4+mqXvHz5pTtM8euvPQ+7aMNehv6wnBUxB+l6UmWeuLAlTaqXC0DAwGefwY03QrNmcOqp0K5dYMYtppS8iIhI8ZCW5mZXjHG1IikpLpEx3pZydh5M5Jkpq/jur63ULB/CmIHt6d26BsbjuKxaBXv3ugLia645nLjIManmRUREir4VK+Cii+CDD+C00wIyZHJqOh/Mi+bVmWtJSbMMPrMht3Y/ibDgAPy731o3uxIU5LY+e02EiinVvIiISPGTnu56tNSr53qgBCgJmLt2N8MmLGP9rkP0aFaNx/u0oH6VMt4GtdYdmnjhha4W56OP3A4iJS65puRFRESKppdfhu+/h1mzDhe5erRlXzxP/biSKcu2U69yGO9e05Eezat7jxVg7ly47DLX3v+66+DkkwMz7glIyYuIiBQdGaUOxkDlyq4YNz4eypb1NGxiShrj5kTxxux1ANx/XlNuOL0BIaU87lDaswf+/tt18j3jDJgyBc4919uYopoXEREpIvbuhcsvd7MWAwcGZEhrLTNX7mTkpBVs2hvPBa0jeOSC5tSqEBqQ8bn8cvjpJ9i82TXIk1xRzYuIiBRNGTuGypd3v6amBmTY6N2HGD5xObNX76JxtbJ8duMpdG1UxfvA8+ZB06ZuZuipp+Cxx5S4BJj3k6hERETyy48/uu3D8fFuZ860aXD11cTExNCtWze2b9+e6yEPJaXy3NRVnPfyHBZt2MdjFzRn8l1nBCZx2brVNcR7/nl33aiRziLKB0peRESk8CpXzu0o2rnTXft25owcOZK5c+cyYsQIv4ey1jJx6TZ6vPgLY2evp2+bmvx8XzduPKMhpYI8/HWYnOyKhgFq1YLx4+Hxx/M+nhyTal5ERKTwSEuDu+92ScBDD7nnMjWaCw0NJTEx8T9vCwkJISEhIdthV2+PZdiEZSyI2kvLmuGMuKglHepVCkzMQ4fC00/D2rVuu7YETHY1L5p5ERGRQiNm505+/vJL4jZvPvxkpj4oUVFRDBw4kLCwMADCwsIYNGgQCxYsyHIZ6WBiCsMnLqf3a7+yanssT/VrxYTbT/eeuGzY4H4A7roLfvhBictxpORFREQK1vLl0KMHbNvGyJEjOXf3bh7IZlUgIiKC8PBwEhMTCQkJITExkfDwcN56660jlpHS0y1fLdrM2aNm88G8DQzoVIdZQ85i0Cn1CPJ68nNSkqvDufded125MlxwgbcxJVe0bCQiIgVr3To2NGnCVdYy96iXsloO6t+/PxEREQwePJgOHTqQlpZ2xOvBNRpR+dxbCY5oQvu6FRhxUSta1SrvLUZrYcGCw2cPTZnizk+qU8fbuJIjLRuJiEjh8dZbcN997nGjRpTetIm6WSwHRUdH/+et48ePZ8yYMbRp04bNmzf/u4xUIjScahfcTcTVL1PjpBa8eFkbvons6j1xAfjiC3eAYkZh7vnnK3EpQOrzIiIix9+6dfDPP+7k51KliKhdO8vloBo1auQ4TEREBOXCwynZrDu1zrgSExxG/eRoJgyPJDyklLcYDx6Ebdvcac/9+8Pbb7suuVLgNPMiIiL5b9cuuOoqWLLEXT/9NEydCqUOJxg7duwgMjKSBQsWEBkZ6VcPlz827OWX0qdS8ZxbaFuvCt0SfqPSxlneExdwsyv/+5/bql26NNx4I5TUv/kLA9W8iIhI/tu/3x1EOGIEXHut5+F2HkzkmSmr+O6vrdQsH8JjfVpwfqsaGK8nNC9bBs2bu4Z4s2dDWBh07uw5XskbHQ8gIiLH16xZ8OWXMHYsVKgAa9ZASIinIZNT0/lgXjSvzlxLSprl9u6NuLX7SYQFB+Cvsz//hI4d4c03YfBgOOss72NKvlDyIiIi+WPlSpg5E3bsgBo1PCcuv67dxbAJy4nadYgezarxeJ8W1K9SxluMaWmu/qZpU2jXDl55xS0VSaGmZSMREQmM5GR47jlo3971PUlLcwW5HpOWzXvjefLHFUxbvoP6lcMY2rcFZzerHpiYBw+G7793CUx4eGDGlIDRVmkREckX/x6SuGOH21L800/uhaCgbBOXJUuWUKFCBf7+++9sx01MSePVmWvp+dIvzFmzm/vPa8q0e870nrhs3w6xse7xrbfC6NHuDCUpMpS8iIhI3q1bx7qePfn9118Z8cwzMH8+vPTSMd925ZVXcuDAAQYOHPif16y1TF++nXNe/oWXZ66hZ4vq/DSkG/2bleHcHmfn6STpf+3e7bY+P/mku27bFi6//IgjCKTwy7dlI2NMHeAjoAaQDoyz1r561D33A4N8lyWB5kBVa+1eY8wGIBZIA1KzmjY6mpaNRESOn5CQEM5KSuIroAewKNPz2R2SmNNuIGstUbviGD5xBb+s2UXjamUZfmFLujaqAsCtt97KW2+9xc0338wbb7yRu2Cjow+fPfT669CrFzRqlLsx5LjLbtkoP5OXCCDCWvunMaYcsBi42Fq7Ipv7+wL3WGvP9l1vADpaa3f7+5lKXkRE8pm18O23EB/PNT/9xEcffUTVoCB2paURFhZGv379GDVqVLbN5ZYsWcLFF1/Mxo0b/32ufv36fPHNd8zaGcK7c6MIKRnE3ec04epT61EqqESeT5L+12uvwQMPwIoV0LBhnv/T5fg77lulrbUxQIzvcawxZiVQC8gyeQGuAD7Pr3hERCQwpl9+OcHp6Xzku97lO1soPj7+mF1x27ZtS5kyR+4QCm7clbum7WH7wUQu7VCbB3o1pVq5w7UyUVFR3HfffXz//ffEx8cfkSRlKyEBDh2CKlXg0kshPh5q1szzf7MULsel5sUYUx9oB/yezethQC/g20xPW2C6MWaxMWZwDmMPNsYsMsYs2rVrVwCjFhERwLXJf+QR2LMHjOHkv//mnQED/j2HKCgoiN69e3PNNdf4VY+yb98+WrZsySsffEWDG18lqf1AqpQL5ttbujLqsjZHJC6Q/UnS2SZJqanQqRPccou7rlkTHnrI864nKTzyvc+LMaYsLim521p7MJvb+gK/WWv3ZnruNGvtNmNMNWCGMWaVtXbO0W+01o4DxoFbNgpw+CIismkTjBoFLVvCoEHUaNmSchUq/JtMJCcnU69ePb/rUFau38jLM9bw2oKNhNdpxgPnNePyTnUIKpF9PUzG0QGDBw9m3LhxxMTE/Pem7dtdP5mSJeGee7REVIzla58XY0wpYBIwzVqbbfm5MeY74Gtr7WfZvP4EEGetzWGOUDUvIiIB89dfbufQrbe6661boVatf1/u378/ERERRyQT48ePz3HI9HTLN4u38NzUVeyLT2bgKXUZck5TKpYJ9h7vtGnQty/8/DOcfrr38aRQKIiCXQN8COy11t6dw33lgWigjrX2kO+5MkAJX61MGWAGMMJaOzWnz1TyIiISILfe6pq3rVkDZct6Hm7p5v0MnbCcpZv306FeRYZf2JJWtcp7G9Rat5RVpYqraXn8cVeYWz1ADeykwBVEk7rTgKuAs40xS3w/vY0xkcaYyEz39QOmZyQuPtWBucaYpcBC4MdjJS4iIuJBWhqMGwerVxMTE8MFS5eyY9asYyYu/zaoy6bWZU9cEg99+zcXv/Eb2/Yn8NL/2vBN5KneExdw/Vl693anPoeFwYsvKnE5QeTnbqO5wDG7/lhrPwA+OOq5KKBNvgQmIiL/tWcP3H8/REYyMjaWqQsWMPzVV49ZxzJy5Ejmzp3LiBEjjrg3NS2dzxZuYtS01cQnp3Hj6Q24s0djyoWU8hbn/v1QvrxrKjdgAOzde8y3SPGjs41ERE5UO3e6U5/vuAOAVqVLszw5+T+3ZdVPJafeK7+s2MrQH5axansspzeqwrC+LWhcPQDt99esga5d4eWX4aqrvI8nhZ7ONhIRkSN99BEMGQJr1wIwY8MGBg4c+O8W6LCwMAYNGkR0dPR/3hoVFfWfey+96kauHvsz/3trPrGJqYwd1J6Pb+jsPXE56Nuo2qiRm21p29bbeFLk5ftWaRERKUR++QVKlXIzGHfeCRdeCI0bA7nrp3LEvWFlKdX6PP6M6EPQuv3ceXYjbjmrEaHBQd7jffJJePtt1x23TBnX2l9OeEpeREROFCkpcP310LQpTJ4MwcHQpMkRt/jVTyXTvZfc9hhbq3dla2wqZeM2MnH4tdStHOY9zvR0KF0aund3nXJ1cKJkopoXEZHiLCnJLQ9dfz0EBcHKlVC/PoSGehp28954nvxxBdOW76B+5TCG9W1J92bVvMcbGwtdurjloccf9z6eFGnH/WwjEREpBKZNg8GDISIC+vSB5s09DZeYksabv6xn7Oz1lDCG+89ryo1nNKB0SY9LRAkJLqEqVw7OPx/at/c2nhRrSl5ERIqb6GhYvx569nRdZ+fNg1NP9TSktZbpK3YwctIKtuxLoM/JETzSuzk1K3ibwQHgm29cU7y//nJdfHM6cFEEJS8iIsXP4MGwbp37CQrynLis3xXH8IkrmLNmF02rl+Ozm06h60lVvMVoLSQnu7qW9u1doqW6FvGTkhcRkaLOWpg4Ec46C8LD4Y033BJMkLelnLikVEb/vJb35kYTUjKIoX1acNWp9SgV5LHLhrXQrx9UrAjvv+8OUPwsy6PtRLKk5EVEpKhbtQouugiefRYefPDfrc95Za1lwtJtPD15JTsOJnFZh9o80KsZVcuV9hZnSorbpm0MdOzojh6wVjMukmtKXkREiqK4OFfLcu65rgh3+nQ38+LRypiDDJuwnIXRe2ldqzxjr+xA+7oVvce7eLGbbfn+e7dM9Nhj3seUE5aSFxGRoujhh+Gdd2DTJqhaFc45x9NwB+JTeHnmGj6av4HyoaV4ul9rLu9Uh6ASHmdF0tLc8lWjRtCqlWZZJCCUvIiIFBX//AMVKkCdOvDIIzBokEtcPEhPt3y9eDPPTV3N/vhkBp1SjyHnNqFCWLD3eIcOdbNDM2a4wxQnT/Y+pghKXkREioaDB+G006B/f/jgA9e3JSLC05BLNu9n2A/LWLrlAB3rVWT4RZ1pWbO8tzjT093sijEuyWrW7PCuIpEAOWbyYoypBpwG1AQSgGXAImttej7HJiJyYktPhzlzDu8i+uor6NTJ87B74pJ4fupqvly0marlSvPy5W24uG0tjNclnW3bXOHwo4/CxRfDTTd5jlUkK9kmL8aY7sBDQCXgL2AnEAJcDJxkjPkGeNFae/A4xCkicuJ56y3XvO2PP9zunF69PA2XmpbOp79v4sXpq4lPTmPwmQ254+xGlAsp5S3OjB1D1apB5cqqa5F8l9PMS2/gJmvtpqNfMMaUBPoA5wDf5lNsIiInnt27Yc8ed3jiNde4XigdOnge9veoPQybsJxV22M5vVEVnriwBY2qlfMe76efwiuvwNy5bmlo6lTvY4ocQ7bJi7X2/hxeSwW+z4+AREROWNa6U5TLlnWFrmFh7oBCD7YfSOSZKSv5Yck2alUI5c0r23Neyxrel4gyZluqVnWzLfv3Q/Xq3sYU8ZM/NS8dgTM4suZlprV2bz7HJiJyYli8GNq1gxIl3CxGRITnpZfk1HTe+y2a0T+tJSXdcmePxtzS7SRCgz0eoJiY6GaETjkF7r3X9Zk591xvY4rkUk41L9cCdwLRwGJgNa7m5XTgQWPMMuDxrJaVRETET7/84gpyP/0UBg6EHj28D7lmF8MnLCdq9yF6Nq/O0D4tqFs5zNugGTMtpUu73i3p2rMhBSenmZcywGnW2oSsXjTGtAUaA0peRERyIzkZoqLcNuIzzoDXX3e7dDzavDeekZNWMH3FDhpUKcP713Wie9Nq3uOdNw/uuQd+/BGqVIGvv1ZRrhSonGpexuT0RmvtkoBHIyJyIrj6alfgunatO0Dxtts8DZeYksbY2et585f1BJUwPNirGdefXp/SJT0uEWWoUAESEtxW6CpVlLhIgfOn5uV54ElcvctUoA1wt7X2k3yOTUSk+Mho4x8aCkOGwJVXusceWGuZtnwHIyetYOv+BPq2qckjvZsRUd7buFgLTzwB8fHwwgvQogUsXaqkRQoNf841P9fXy6UPsAVoAmS7E0lERI6yZYtbInrhBXfdqRP06eNpyHU747j6vYVEfrKYsqVL8vlNXRh9RTvviQu4JGXPHrdt29rDz4kUEv4cD5DRvag38Lm1dq/nLXYiIieC6Gho0ABq14Znn3VdZz2KS0pl9E9reXduNKGlghjWtwVXdalHySB//i2agw0b4Oab3W6n5s3htdfc7ieRQsif5GWiMWYVbtnoVmNMVSAxf8MSESninnkGnnwSVq92ycudd3oazlrLhKXbeOrHleyMTeKyDrV5oFczqpYL0JlBZcrA+vXup3lzJS5SqB0zebHWPmSMeQ44aK1NM8YcAryXxYuIFDeHDkFSElSq5JrLlS4dkMZtK2MOMmzCchZG7+Xk2uV566oOtKtb0Xu8n3/uTnx+7z1Xj7N6NQQFqMhXJB/5U7AbhGtSV993LECGl/ItKhGRoiYpCdq0gdNPd6c+N2jgmrh5cCA+hZdmrObjBRupEBbMM/1bc3nHOpQoEaCl+61bYdUqiI2FcuWUuEiR4deyEW6Z6B9AXYlERDLbvh1q1HCzLPfdB61bex4yPd3y1aLNPD9tNfvjk7mqSz3uPacp5cM8HqB48KDb6fS//8E558Ddd7sES0tEUsT4k7zUttaenO+RiIgUNT/84BKBefPc4YmRkZ6HXLJ5P8N+WMbSLQfoVL8iwy88hRY1wwMQLC7BmjvXbX0+5xwo6c9fASKFjz+/c6cYY8611k7P92hERAq79HR3CGGlSq6t/+23Q716nofdHZfE81NX8dWiLVQrV5qXL2/DxW1reT9Acf58GD0aPvrIJS9LlrhfRYowf5KXBcB3xpgSQApgAGutDdA/BUREipB+/VzyMns2lC8PL77oabjUtHQ+XrCRl2asISE5jZvPbMgdPRpTtnSAZkViYuC339xW6EaNlLhIseDP/zteBE4F/rE2o1vRsRlj6gAfATVwtTLjrLWvHnXPWcAPuMMfAcZba0f4XusFvAoEAe9Ya5/197NFRALqwAEID3eN2i6/3J1NFAALovbwxITlrNoeyxmNqzCsb0saVSvrbdC0NJdQVa/uTn/u1w/OP99zN1+RwsSf5GUtsCw3iYtPKjDEWvunMaYcsNgYM8Nau+Ko+3611h7RatK3w2kMcA6uq+8fxpgJWbxXRCR/LVsG3brB229D//7u5GePYg4k8PTkVUxcuo1aFUJ588oOnNeyuvclInDFt5MmuVmWa65xCZcSFylm/EleYoDZxpgpQFLGk9baHLdKW2tjfO/FWhtrjFkJ1AL8SUA6A+ustVEAxpgvcL1llLyIyPERFwdly7q2/pdeCk2beh4yKTWN9+ZuYPTPa0lNt9zVozGR3U4iNNjjFuWNG91ZRK+84paypkxxTedEiil/9sdFAz8BwUC5TD9+M8bUB9oBv2fx8qnGmKXGmCnGmJa+52oBmzPds8X3XFZjDzbGLDLGLNq1a1duwhIRydqwYdC+PSQmuh05b70FLVse+305mL16J+e/8ivPTV1F15OqMPOebtxzThPviQu4M4i+/RYWLXLXSlykmPOnw+5wLx9gjCkLfIs7ifrgUS//CdSz1sYZY3oD3wONcUXB/wklm/jGAeMAOnbsmNulLRERJyXF/VqqFJxxhqtryfVq+X9t3hvPiEkrmLFiBw2qlOH96zrRvWk1z+Py3XduxuXuu9027S1bXF2OyAngmDMvxpgZxpgKma4rGmOm+TO4MaYULnH51Fo7/ujXrbUHrbVxvseTgVLGmCq4mZY6mW6tDWzz5zNFRHJt3z5o29YdRgjQs6c7m8hDrUhCchovzVhDz5d+4bd1u3mwVzOm3n1GYBIXgO+/h08/hdRUd63ERU4g/tS8VLXW7s+4sNbuM8Yc8/99xlWevQuszK4+xhhTA9hhrbXGmM64ZGoPsB9obIxpAGwFBgDeq+RERDJLSnJbhytWdEW5zZt7HtJay7TlOxg5aQVb9yfQt01NHu3dnBrlQ7wNHBvrDnqMjHRHD4weDWFhajQnJyR/ftenGWPqWms3ARhj6pHNEs5RTgOuAv4xxizxPfcIUBfAWvsmcClwizEmFXdq9QDfrqZUY8ztwDTcVun3rLXL/f/PEhE5hi+/dK3y//rLHUr4xhueh1y3M47hE5fz69rdNKtRji8Gd6FLw8oBCBa3XfuNN1ziEhmpmRY5ofmTvDwKzDXG/OK7PhMYfKw3WWvnknXtSuZ7Xgdez+a1ycBkP+ITEfFfSoqrazn5ZOja9fCyiwdxSam89tNa3psbTVhwEE/0bcGVXepRMsjjmUFLlrhtz489BrVrQ1SUS7RETnDGn/YtvjqULrhkZL61dnd+B5YXHTt2tIsyqu1FRDJLT3d9WmrXhtez/DdTrllr+X7JVp6ZvIqdsUlc3rEO9/dqSpWyAepiO3w4jBkDy5craZETkjFmsbW249HPZzvzYoypb63dAOBLViYd9boBallrtwQ4VhGRwElLg6Ag17ytRYuAJQHLtx3giQnL+WPDPtrULs+4qzvStk4F77G+8w60aQNdusADD8Bdd0EFj+OKFDM5LRu94DvP6AdgMbALCAEaAd2BHsAw3M4gEZHCZ+FCuOIKmDjRJS5PP+15yP3xybw0Yw2fLNhIhbBgnrukNZd1qEOJEgHojpuYCCNHwoUXuuQlNFTdcUWykG3yYq29zBjTAhgEXA9EAPHASlwtylPW2sTjEqWISG6kp7uZlgYNoG5dt6vIo7R0y1eLNvPCtNXsj0/mqi71uPecppQPK+Vt4O3b3dEDjz7qmsvNn++WtkQkW37VvBQVqnkREYYNc4Wu33/vzvUJgL827WPYhOX8veUAnetX4okLW9KiZoB2+3z0Edx4o0taOnQIzJgixUSua15ERIoMaw8nKhUrQo0abldRcLCnYXfFJvH81FV8vXgL1cNL8+qAtlzYpqb3AxR//tktEfXuDVde6Tr6NmjgbUyRE4hmXkSkaNu6Ff73Pzfjcu65ARkyNS2dj+Zv5OUZa0hMTeP60xtwx9mNKVs6AP/esxY6d4aQEPj1V+/jiRRj2c28eGxCICJSQDL+4VWlipt1iY/P9RAxMTF069aN7du3//vc/PV7uOC1uYyYtIK2dSsw5a4zefj85t4Sl8REeOklF6Mx8M03MGNG3scTOcH5c7aRMcZcaYwZ6ruu62vlLyJSML74Anr0cA3mSpd2MxgXX5zrYUaOHMncuXMZMWIEMQcSuP2zP7ni7QXEJaXy5pUd+Oj6zjSqVtZ7vIsXu26+P/zgruvVczMvIpIn/vxT4g0gHTgbGAHE4g5b7JSPcYmIZK90abejaO9eqFYt14W5oaGhJCb6NksGleTTP3cxafhkTAnD3ee14pazTiKkVJC3GNeudYXDl10Gp50G//wDrVp5G1NEAP+WjU6x1t4GJII7mBHwVgUnIpIbCQlw7bUwbpy7vvhimDXLJS55EBUVRb9+/Qg7qSM1rx9DxbOupUraLr67qQP3nNPEe+IC8PjjcOedbskIlLiIBJA/yUuKMSYI32GMxpiquJkYEZHjIyQEtm1zMy3gZlo87PhJCS7P0vAuVL30CQB2fT2MLil/075pvbzHaK077HHbNnf90kvw559aHhLJB8fcbWSMGQRcDrQHPsSdBP2Ytfbr/A8vd7TbSKQY+f13ePhh+O47KF/+cOM5DxKS06hz3g2U6XARNj2NA799wcFFP0B6KiVKlCAtLS3vg2/eDI0awf33w5NPeopTRJw87zay1n4KPAA8A8QAFxfGxEVEipmSJV1CsHGju/aQuFhrmboshp4v/ULZUy6jatI29n1yNwcXfktYSDCDBg1i69atuR94/3432wJQp44rHB4+PM9xioh//NltVAnYCXwOfAbsMMZ47IctInIUa+GZZ2DECHfdoQOsWgUnn+xp2HU7Y7nq3YVEfvIn5UJK8uXgLnRKXUb8nhhCQkJITEwkPDycGjVq5H7wF1+EQYNg0yZ33bmzOwRSRPKVP7uN/gTqAPsAA1QAYowxO4GbrLWL8y88ETlhGAOrV0Ny8uGOuR4SgdjEFEb/vI735kYTFhzE8AtbMuiUupQMKsFzO3YQGRnJ4MGDGTduHDExMf4PvGiROyyxZUu3RNSvnzs/SUSOG39qXt4EvrPWTvNdnwv0Ar4CXrXWnpLvUfpJNS8iRcymTXD33TBqFDRs6Fr6l/I2sWut5fslW3l68ip2xyXxvw51uL9XU6qULe093sREqF8funaF8eO9jyciOfJytlFHa21kxoW1drox5mlr7b3GmAD8aSAiJ6wSJWDBAli2zCUvHhOX5dsO8MSE5fyxYR9tapfn7as70rZOBW8xpqa65nL9+7udQ99/Dy1aeBtTRDzxpwJurzHmQWNMPd/PA8A+3/ZpbZkWkdz57js32wJQuzZER8OFFx5xS1Zt+3OyPz6Zx79fRt/Rc4nadYjnLmnNd7ee5j1xAdfN99JLXV8ZgC5dIDxAJ0qLSJ74M/MyEBgGfO+7nut7Lgj4X/6EJSLF1j//wC+/QFwclC3ruuUeJXPb/jfeeCPbodLSLV/8sYlR01ZzICGFq0+tzz09m1A+zOOegm3b3IGPnTrBFVdApUrQvbu3MUUkYHKsefHNrnxorb3y+IWUd6p5ESmE4uJg6FC45BLXJj852RXiZlGMe0Tb/kxCQkJISEg44rnFG/cxbMIylm09SOcGlRh+YUuaRwRoRqRLFzh4EJYv99QMT0S8yVOfF2ttGlDVGKPjAEQkb4yBb7+F335z18HB2e4iioqKYuDAgYSFhQEQFhbGoEGDiI6O/veeXbFJ3Pf1Ui4ZO49dsUm8OqAtXw7u4j1x+flnl1gBjBkDEycqcREppPxZNtoA/GaMmQAcynjSWvtSfgUlIkXcX3+5c4jGjIEyZdwMRtljn84cERFBeHg4iYmJ/+nBkpKWzsfzN/LyjDUkpqYR2e0kbj+7EWVLuz/GYmJiGDBgAF9++WXue7YsWuROqR4zBm691fWYEZFCy5+C3W3AJN+95TL9iIhkbdkyt5U4Kspd+5G4ZNjh68GyYMECIiMj2b59O/PX76HPa3MZMWkF7epVZOrdZ/LQ+c3+TVzgyDoZvyQkuCMIADp2dJ1yb7jB7zhFpOAcs89LUaKaF5ECkpYGb78NVau62hZrITbW866cmAMJPPXjSib9HUPtiqEM7dOCc1pUx2RazslNncwRrrvO7XzatEm7h0QKqTyfbWSMqWqMecEYM9kY83PGT/6EKSJF1rvvutoWcLUiHhKCpNQ0xsxax9mjfmHGih3c07MJM+/txrktaxyRuIB/dTL/Wr0a9uxxjx9+2PVsUeIiUuT4s2z0KbAKaAAMx9XA/JGPMYlIUbB9O9xzD8THuwLcadPg0089Dztr1U7Oe3kOL0xbzZlNqjDz3m7c1bMxIaWyLvLNqU7mCLt3Q7t2MHKku27SBM46y3O8InL8+ZO8VLbWvgukWGt/sdZeD3TJ57hEpLBbuxbGjj28i6hSJU+7czbuOcSNH/7BdR/8QQlj+PD6zrx1VUfqVAo75nuzqpMB3PLVn3+6x1WqwHvvuRkXESnS/DnbaIG1tosxZhrwGq6A9xtr7UnHI8DcUM2LSD6bORM2bjxc2LpjB1SvnqehMnYHffjJ54xfGcubc6IoVcJwZ4/GXHdaA4JL+vNvq2MYNQoeesjtdmra1Pt4InJceTnb6EljTHlgCDAaCAfuCXB8IlIUvPWWqxu59lq3VJTHxAVgxMiRLN6ZRq/Rv5FYIoyL29bk4d7NqR4e4i3Gffvg0CF39MB110HFitC4sbcxRaRQ0W4jEcleYiK89JJLVmrWdMWuZcq4AwrzKDQ0lLQyVanY82ZC67cleWc0e2e8idm9PufdQf5IS4NmzVyyMnmyt7FEpMDleebFGFMVuAmon/l+X+1LTu+rA3wE1MAd4DjOWvvqUfcMAh70XcYBt1hrl/pe2wDEAmlAalbBi0g+27YNRoyAcuXgjjugcmVPw8UmpnD3h3P4YvF20pLi2TN9LOlrfuHSiy9i1KiZeR949Wq3LBQUBM8+C40aeYpTRAo3f5aNfgB+BWbiEgl/pQJDrLV/GmPKAYuNMTOstSsy3RMNdLPW7jPGnA+MA07J9Hp3a+3uXHymiHi1dq3bOXT77dCwIaxaBfXrexoyPd3y3V9beXbqKnbHJVE7ZTML3nmUkmmJJCcnZ707yF+TJ8MFF8DUqXDeea7PjIgUa/5UxIVZax+01n5lrf024+dYb7LWxlhr//Q9jgVWArWOumeetXaf73IBUDuX8YtIoL3zDjz2mNtaDJ4Tl2VbD3DZW/MZ8vVSalYI5ftbT6PqhpkMvmbgf3cH+Ss11TWXA+jZE55+Grp29RSniBQd/uw2ehKYZ63N8wKyMaY+MAdoZa09mM099wHNrLU3+q6jgX2ABd6y1o471ueo5kUkD6yFr75ytSJt2rjOuIcOQV5nQnz2HUpm1PTVfLZwE5XCgnnw/GZc2r42JUoE4LDDiy6Cdetg6VIo6c8EsogURbmueTHGxOISBwM8YoxJAlJ819Za61dbSmNMWeBb4O4cEpfuwA3A6ZmePs1au80YUw2YYYxZZa2dk8V7BwODAerWretPSCKSWWysq2e5+GJ3mGK5cu4nj9LSLV/8sYkXpq0mNjGVa06tzz3nNKF8aClvcW7b5nY3BQXBbbe5BCub06lFpHjL191GxphSuEMdp2V3CrUx5mTgO+B8a+2abO55Aoiz1o7K6fM08yLip/374eOPXV2LMa6upXFjz8nA4o37GDZhGcu2HuSUBpUYflFLmtUIQPv9NWugfXt44QW45Rbv44lIkeDlbKN+vj4vGdcVjDEX+/E+A7wLrMwhcakLjAeuypy4GGPK+Ip8McaUAc4Flh3rM0XET19/DXfddbj7bLNmnhKXnbGJDPlqKZeMncfu2GReu6IdXwzu4j1xyaiFadwY7rvPFeSKyAnPn5qXJdbatkc995e1tt0x3nc6bpfSP7it0gCPAHUBrLVvGmPeAS4BNvpeT7XWdjTGNMTNxoBb2vrMWvvUsf5jNPMikoNFi9xSS7durh/KihXQuvW/L2d0vP3yyy/93vmTkpbOh/M28OrMtSSmpnHD6Q254+xGlCkdgDqUxx8/3BSvYkXv44lIkeOlw25WszPHfJ+1di6uPiane24Ebszi+SigjR+xiYg/rIXrr4eyZWHePDfLkilxARg5ciRz585lxIgRvPHGG8ccct663QybsJy1O+Po1qQqw/q2oGHVst7ijI93sZYpA/37u3jDjn22kYicWPyZeXkP2A+MwRXw3gFUtNZem9/B5ZZmXkQySU11dS0DB0Lp0q6upWZNCD9yKSc0NJTExMT/vD0kJCTLjrfb9ifw1I8r+fGfGOpUCmVon5b0bF4N4+FQRgDi4lxCddll8Pzz3sYSkWIhzzUvuGQlGfgS+BpIBG4LbHgiEnBz57rZlm++cdfNmv0ncQGIiopi4MCBhPlmOMLCwhg0aBDR0dFH3JeUmsaYWevo8eIvzFy5g3t6NmHGPd04p0V1b4nLPl+rp7Jl4cYbXcM5EZEc+LP8cwh4CMAYEwSU8T0nIoXNtm3wzz+usPWss+DXX+G003J8S0REBOHh4SQmJhISEkJiYiJBQUFcfvnl/9a/zFq1k+ETl7NhTzy9Wtbg0QuaU6dSAJZzvvgCbroJ/vrLtfR/9FHvY4pIsefP2UafAZG4owEWA+WNMS9Za1/I7+BEJJfuuMPNuGzc6A5PPP30Y78H2LFjB5GRkQwePJhx48YxZcoUNm7cyIMjR1Gy0+X8tGonJ1Utw8c3dOaMxlW9xWitKxwuWxbOPBOuugrKlz/2+0REfPzebeQ7RLED7iDFxdbak49HgLmhmhc5IU2f7nqgVKkC0dGQng4nnZSnoTLqX0yp0pTvchnhnS/BpqVw6PeviZnzBcEl/VlpzoG1cP75rgne1197G0tEij0vu41K+ZrNXQy8bq1NMcbkX2c7EfHfxo3Quzc89BA8+SQ0aOBpuPXr13Pd0NdYUbo5QeWqkLhqDmeE72X0+Ne8JS7x8W7XkDGupiUszCUyXot8ReSE5M+fRm8BG4AywBxjTD0gyzb/InIcJCTAlCnucb167jTlxx/3POzaHbHc9+MmVlc5g7SEWPZ+/Ti7Jo6iatlSeT/xGVx/mXr1YI7vdI877oAbblDiIiJ55k/B7mvAa5me2ug7i0hECsLIkW4rcXQ01KnjTlX2IDYxhVdnruWDeRsICw6izo55tAo7QOT49xk3bhwxMTF5Gzgx0dXdtGjhYqxc2VOcIiIZsq15McZcaa39xBhzb1avZ9fyvyCp5kWKrdWr3enJJ50Ee/fC33+73UQepKdbvvtrK89OXcXuuCQGdKrDfec2pXLZ0t7jfeghmDUL5s+HEh7rZETkhJWXmpcyvl/zfrysiHiXmAhnnOF25nzzDVSq5DlxWbb1AMMmLGfxxn20rVOBd6/pyMm1K3iLMzXVde41Btq1czUtyclu9kVEJIDy9VTp400zL1JsWOtmLs4+213PnAknnwzVqmV5u7/nEu07lMyo6av5bOEmKpcJ5sFezbikfW1KlPBYfxITA+eeCw884LY+i4gEQK5nXowxr2X3GoC19s5ABCYiWfjkE7j6apfAnHXWMetajnUuUVq65fOFmxg1fTWxialc27U+d/dsQvnQUt7iTE11y1nVq7valipVvI0nIuKHnGpervE9PA1ogTseAOAyXJ+Xe/I/vNzRzIsUaXv3uhmMli3dcsvXX8MVV+RYM+LPuUSLN+5l6A/LWb7tIF0aVmL4ha1oWiMAq8EffwxPP+12E5Upc+z7RURyKddnG1lrP7TWfgg0Brpba0dba0cDPYC2+RapyImqd2+XrFgLwcEwaNAxi11zOpdo58FE7v1yCZeMnc+euGRGX9GOz2/q4j1xSU93v550kjsv6ZBOCxGR48ufJnU1cUW7e33XZX3PiYhXixdDmzZu6eWFF9zBibnof5LVuURlw8szaW08r8xcSlJqGrecdRK3d29EmdL+/N89B8nJ7sTnDh1g6FDo2hW++87bmCIieeDPn2bPAn8ZY2b5rrsBT+RbRCInij/+gM6dYdw4dzjhGWfkaZjM5xI9/e63zA5qwtQfV3JW06oM7dOChlXLeoszoxNucDBUrapziESkwPm128gYUwM4xXf5u7V2e75GlUeqeZFCLyXF9Wxp1colBePGueWhst4SjK37E3j6x5X8+E8MdSuFMbRPC3o0r4bx2sX2t9/g5pth2jSoVcvbWCIiuZSX3Ub1rbUbAHzJyg9HvW6AWtbaLQGOVaT4uvlmmDgRoqLc4YQ33+xpuMSUNN75NYrXZ60D4N5zmjD4zIaElAryFmfGbEtEhCvG3btXyYuIFBo5LRu9YIwpgUtaFgO7gBCgEdAdV7g7DFDyIpKTzZvdUkt4ONx9N1x4oeeZFoCfVu5g+MQVbNobz/mtavDoBc2pXTHMe7yPPQY7d7pZoYYNYcECnUMkIoVKtsmLtfYyY0wLYBBwPRABJAArgR+Bp6y1/92jKSKH7drl+p/ccos7j+jkk92PBxt2H2LEpBX8vGonJ1Utw8c3dOaMxlUDFDBu1sVat6uoRAklLiJS6KjDrkh+WLsWGjd2j996C847D+rX9zRkfHIqY2at4+050ZQKMtzVszHXdm1AcEmPZwdt2ADXXgsvvuh2EmUsGYmIFLC8nG2U8cb+WTx9APjHWrszEMGJFCsvv+wOJlyxwvVC8VjXYq1l8j/befLHFcQcSKR/u1o8dH4zqoUH6MygihVhzx7XIA+UuIhIoefPVukbgFOBjK3SZwELgCbGmBHW2o/zKTaRouPQIYiPd1uJBwxwsxd16ngeds2OWIb9sJz5UXtoERHO6Cva0bF+Je/xfvqp69Hy9deuHufvv5W0iEiR4U/ykg40t9buADDGVAfG4rZOzwGUvMiJLTXVLbe0bAnffut26Nx7r6chDyam8MqMtXw4fwNlS5dk5MWtGNi5LkFeD1DMEBfninIPHIAKFZS4iEiR4k/yUj8jcfHZCTSx1u41xqTkU1wihd/WrW77cMmS8PDD0KiR5yHT0y3j/9rKs1NWsudQMld0rst95zalUplgbwPHxrqdTn36QL9+rine4MFKWkSkSPInefnVGDMJ+Np3fSkwxxhTBtifX4GJFGpTp0LfvvDTT3DmmXDNNcd+zzEs23qAoT8s489N+2lXtwLvX9uZ1rUD1M02NBSWLHE7n+CYZyaJiBRm/iQvtwH9gdMBA3wIfGvdNqXu+RibSOGSng67d0O1ai5hGTLkcDLgwb5DybwwfTWfL9xE5TLBjLqsDf3b1aKE1yWi+fPd9uwvvoDSpeH3390skYhIEXfMP8mstdYYMxdIBiyw0Ban/dUi/rrsMti40SUBYWHw7LOehktLt3y+cBOjpq8mNjGV609rwF09GxMeUiow8cbGwl9/QXS0O/1ZiYuIFBP+bJX+H/ACMBs38zLaGHO/tfabfI5NpODt2eO2EpcoAVddBQcPeq4TiYmJ4eKbhlC++42s2ZVAl4aVGHFRK5pUL+ct1vT0wydT33ILnHsurFnjDlQUESlG/Pmn2KNAp4yeLsaYqsBMQMmLFG+rV0OXLvDKK66m5eKLPQ+582Ail74wgR2tBnFgx15ev7IrF7SO8H6AIrikas4cqFzZJS+gxEVEiiV/qvZKHNWMbo+f7xMpmvbvd782buw6z3bq5HnIlLR0Kne9jI7DJrKlRHUOzPuSNa9cSd82tQgL83Ae0aZNbkZo1y6XvHzzDXz0ked4RUQKM3+SkKnGmGnGmGuNMdfizjWanL9hiRSQESOgdWvXB6VECdct12NR7ty1uzn/1V8pd+a1VE7fx77PhrD/148JLRXEoEGDiI6OzvvgcXEwaRIsXuyuQ0M9xSoiUhQcM3mx1t4PjANOBtoA46y1Dx7rfcaYOsaYWcaYlcaY5caYu7K4xxhjXjPGrDPG/G2MaZ/ptV7GmNW+1x7K3X+WFEcxMTF069aN7du3B3bg5GRISnKPe/Z0sy0B2Eq8dX8Ct3yymCvf/Z3k1HTevaYjXZL/4tCODYSEhJCYmEh4eDg1atTI3cA//ABPPeUet2jhTq3u1ctzvCIiRYVff0Jba7+11t5rrb3HWvudn2OnAkOstc2BLsBtvlOqMzsfaOz7GYzr3IsxJggY43u9BXBFFu+VE8zIkSOZO3cuI0aMCNygsbHulOeMnUNdu8LIkW43UR5t2LyVNgMf4uxRs5i1eidDzmnC9HvOpEfz6uzYsYPIyEgWLFhAZGRk3hKxGTPc8lBysrsuWzbPsYqIFEXZniptjInFbY3+z0u4HdThufogY34AXrfWzsj03FvAbGvt577r1bizk+oDT1hrz/M9/zDuQ5/J6TN0qnTxFBoaSmJi4n+eDwkJISEhIW+DxscfTlAefBC6d/c8e2Gt5aeVO7nzg1+IL1GGGikxfPv4IGpV8LiUExcHTz7pioabN3fnKJUura3PIlLsZXeqdLYzL9bactba8Cx+yuUhcakPtAN+P+qlWsDmTNdbfM9l97ycgKKiohg4cOC/ha1hYWHeakW++grq1oUtW9z1c895Tlyidx+i5oCR3PjRIvbv3cOOLx7l95duonbFMEK91qEkJMC4cTBtmrsuU0aJi4ic0PJ915AxpizwLXC3tfbg0S9n8Rabw/NZjT/YGLPIGLNo165d3oKVQikiIoLw8HASExPzXitiLWTM3nTuDL17Q1CQ59jik1N5fuoqznt5DhWbdqL+/j858OWDJG5c6i3J+vtveOwx97hqVVi71p1NJCIi+Zu8GGNK4RKXT62147O4ZQtQJ9N1bWBbDs//h7V2nLW2o7W2Y9WqVQMTuBQ6nmpFrIULL4Sbb3bX9eu77cQREXmOx1rLxKXb6PHiL7wxez192kQw+/7uNDfbSIw/5K0gF1xdy5tvusMfwfVuERERwL8mdXliXNetd4GV1tqXsrltAnC7MeYL4BTggLU2xhizC2hsjGkAbAUGAAPzK1Yp/MaPH09MTAwDBgzgyy+/9C8hSE52TdqMcYW45Tx2sPVZvT2WYROWsSBqLy0iwhl9RTs61q8EHE6yBg8ezLhx44iJifFv0PR0+OADaNDA1d/ceSdcdx1UqhSQmEVEipNsC3Y9D2zM6cCvwD9Auu/pR4C6ANbaN30JzutALyAeuM5au8j3/t7AK0AQ8J619qljfaYKdou3W2+9lbfeeoubb76ZN954I+eb//gD+vWD77+Hjv+p9cqTAwkpvDJzDR/N30jZ0iW577ymDOxclyCvByiC26rdsiWccQa8/7738UREioHsCnbzLXkpCEpeiqdc7TZKTXXFrAcOwJVXul06bdp4+vz0dMu3f27huamr2HMomSs61+W+c5tSqYzH1vu7dsHrr8Pjj7uYt26FmjU9n50kIlJc5Hq3kUhh4fduo8ceg3POcTUu5cvDxImeE5d/thzgkjfncf83f1O3UhgTbz+dp/u1znPickSjvV9/haefhgUL3Iu1ailxERHxg/ZbSqGX426jtDTXDdcYaNjQNZ1LTnZ9UDzYeyiZF6at5os/NlG5TGlGXdaG/u1qUcLjEtFnt9xClV9/ZcSIEbwxZow79blBA09jioicaLRsJEVC//79iYiIOKIQdvzo0dC3r5tx6d8/IJ+Tlm757PeNjJq+hrikVK7tWp+7ejYmPKSUp3Ezlr5mAeWAjDlQT432RESKueyWjTTzIkXC+PGHd9qPGT3azbakprrtzqW8JRYZ/tiwl2E/LGdFzEG6nlSZJy5sSZPqHncoJSfD2LFEL1nCkBEjuGn8eLYmJhIWFka/fv0YNWpUQGIXETmRKHmRouWTT+Cll2DePAgJgR9/9DzkzoOJPDNlFd/9tZWa5UN4Y1B7zm9VAxOI+pMVK+Cee6gRFkZ4eDhRyckEh4SQ5KUHjIjICU7JixQN1rq6lurV3WzLwYMuefEgOTWdD+ZF8+rMtaSkWW7v3ohbu59EWLDH/1ts3Ai//QYDB0Lbtq5bbqtW7OjfP289YERE5AiqeZHCLTHRbXk+9VQYMiRgw/66dhdPTFjO+l2HOLtZNYb2aUH9KmUCM/hNN7lTnzdtClhjPBGRE5G2SkvRkpFUh4S4+pYSgfmtumVfPLd8spir3l1Iarrl3Ws68t61nbJNXI7Y2pxTrBMmQMbW7aeecrMtSlxERPKFkhcpfH77zXXFzTho88sv4Z57PA2ZmJLGaz+tpceLs5m6dDORXWsy7e4z6dG8eo7vGzlyJHPnzmXEiBHZ37R7N1xxBbz8sruuVg3q1Mn+fhER8UTJixQ+lSq5s34yZjs8FM5aa5mxYgfnvPwLL81YQ8WErWx9+2Y2TX2bkFLZnyodGhqKMYaxY8eSnp7O2LFjMcYQGhrqbjh0CD7/3D2uWhVmz4YXX8xznCIi4j8lL1LwrIVHHz08u9K8Ofz5J7Ru7WnY6N2HuO6DP7jpo0VEr13Djs8f4feXBpNyYOd/k5GjHLOr75gxriB35Up33alTwLZsi4hIzpS8SMEzBuLiXHfcjFoXD7Mth5JSeW7qKs59+RcWbdjHYxc0Z/6wC+l/eqtjHzHgk1VX36YpKdTImA264w63Xbt58zzHKSIieaOt0lIwoqLcrpzRo6FFC3jlFc/n+lhrmfh3DE//uJLtBxO5pH1tHjy/KdXKuS3V2R4xkI0dO3b8u7X57TffZPBHH8GePTBzJoSGuh1QIiJy3Cl5kYIRHg5btsCGDS558Zi4rNp+kCcmLGdB1F5a1gxnzKB2dKhX6Yh7Micj/vRZGf/NN+5wx9ateX3sWLjuOnd+koiIFCj1eZHj5+OPYcYM+PBDl6ykp3veAn0gIYWXZ6zh4wUbKRdSkvvPa8qATnUJ8niAIuC2P190EYwfD/36eR9PRERyRWcbScHbvdvNtMTFuR4oHhKX9HTLN4u38NzUVeyNT2bQKXUZck5TKpYJ9hbjnj2wbh2ccgr06QPff+8OfxQRkUJDMy+Sfw4cgLvvhgED4LzzIC3NJSwel4iWbt7PsAnLWbJ5P+3rVmDERa1oVat8YGI+7zy3gygqCkoqtxcRKUiaeZHjLyQE/vgD2rVzSUFQ9n1V/LH3UDIvTFvFF39spnKZ0oy6rA3929WihNclogUL4OSTISwMXnjBJVhKXERECi1tlZbAmjsX/vc/SE2F0qXhr7/gzjs9DZmWbvl4/ga6j5rN14u2cMNpDfj5vm5c2qF2lomLXy39M6xZA127wquvuuuTT4ZWrTzFKyIi+UvJiwTWrl2waJE7WRk8N277Y8Ne+oyey+M/LKdVrXCm3HUGj/VpQXhI9uMes6V/SgosXOgeN2kCn33m+raIiEiRoJoX8SY1FZ57DmrWdFuJrYXkZDfr4sGOg4k8M3kl3y/ZRs3yITx6QQt6t66ByaFeJjQ0lMTExP88HxISQkJCwuEn7r0X3nzTHaRYPeezjUREpODoVGnJH0FBbvvzggXu2hhPiUtyajpv/bKes0fNZvKy7dxxdiNmDunGBSdH5Ji4wDFa+m/ZAjt2uBvvusudS1StWp7jFBGRgqOqRMm9qCgYOhRefx0qVIApU1zHWY/mrNnFExOXE7XrED2aVWNo3xbUq1zG7/dn1dI/PDycGmXLQr16rmfLe++5x/XqeY5XREQKhmZeJPf274cff4QlS9y1x8Rl8954bv54EVe/t5D0dMv713bi3Ws75SpxyZDRRXfBggUMv+wyV7Rbtqw7SHHoUE9xiohI4aCaF/HPl1/Cpk1w//3uOi7OJQUeJKak8eYv63lj1jpSkpOIPKM+d/c+mdIlvW2pBtwMy403uq3aHTp4H09ERI479XkRb6ZOhVWrXLFrUJCnxMVay4wVOxgxaQVb9iVQI2Ubf777GBsPXULpC9/Ie4zx8a5Dbp06cNll7pTq1q3zPp6IiBRKmnmRrO3fD0884Xq0NGzoZlpCQz03movaFcfwiSv4Zc0uUndvYs+MN0nc9PcR9/xnd5A/rHUt/UuXhjlzPHfxFRGRgqfdRpI7hw7BBx/Azz+767JlPSUuh5JSeXbKKs57ZQ5/btzH431aMP+Ji+h/equsdwf5a8MGl7gYA48+Ck8+qcRFRKSY07KRHLZwoSvEHT4catVyu4oqVfI0pLWWCUu38czkVWw/mMgl7Wvz0PnNqFrObafOcndQjRr+DT5/Ppx5pmsyd9llbjeRiIgUe5p5kcOmToW333Z1I+A5cVm1/SCXj1vAXV8soUq5YL69pSsv/q/Nv4kLHLk7KDIy8tgt/a2Fbdvc486d4ZFHXAIjIiInDNW8nMhSUuCNN6BTJ3e+T2Kie65cOU/DHkhI4eUZa/ho/gbCQ0vxwHnNuLxTHYK8HqAIMHgwzJwJy5cHpLeMiIgUXtptJP+VnAyjRsEll7jkJSTE/eRRerrlnZ+X8eyUlRBchoGn1GXIOU2pWCbYW5z797sTn4ODYdAgl2wFexxTRESKrHxbNjLGvGeM2WmMWZbN6/cbY5b4fpYZY9KMMZV8r20wxvzje01TKYG0ebMrbE1PhzJlXJ3Lyy97Hnbp5v30GzuPp2duImHnJk6N/ZUnL27tPXHZvt0dnvjKK+66Wze46SbPu55ERKToyrdlI2PMmUAc8JG1ttUx7u0L3GOtPdt3vQHoaK3dnZvP1LKRHz75xP3lP38+tG3rebg9cUm8MG01ny/cSNqh/eyf/T6Hls/69/U8bXsG2Lv3cM3N0KHQrx+0a+c5XhERKTqO+1Zpa+0cYK+ft18BfJ5fsZzQrIVJk2DCBHc9aBCsW+c5cUlNS+fDeRvoPmo23yzewsXNwkn+7lHSo9wBjXna9pxhzBho1OjwQYojRihxERGRfxV4zYsxJgzoBdye6WkLTDfGWOAta+24AgmuOLDW/eVfrhxceKHrgVKrlqchF0bvZegPy1i1PZbTGlVm+IUteemJB9m1bTNA3rY9p6W5guEyZaBnT7dNWwW5IiKShcKwVbov8Ju1NvMszWnW2vbA+cBtviWoLBljBhtjFhljFu3atSu/Yy2UYmJi6Nat2+FtxnFx8NRTrtFciRIwfrw7+dmjHQcTueuLv/jfW/M5mJDC2EHtGX9HdxpXD2fs2LH/3peYmIgxJsttz/+JFdwOpy5d4L773HXTpvDiixAe7jlmEREpfgpD8jKAo5aMrLXbfL/uBL4DOmf3ZmvtOGttR2ttx6pVq+ZroIXVyJEjmTt3LiNGjHBPLF0Kjz9+OGGpXdvT7pzk1HTe+mU9Z4+azZRl27nz7Eb8NOQszm8dQVRUFAMHDvxPl9wtW7Ywfvz4nGM9dMg9WaqUq2np0SPPMYqIyImjQJeNjDHlgW7AlZmeKwOUsNbG+h6fC4wooBALtdDQUBITEwFoBRwYOxYzdqwrkl29Gho39vwZc9bs4omJy4nadYiezasztE8L6lYO+/f1iIgIv7rkZo4VYM3YsewdO5Zzg4NZlJTkms2JiIj4IT+3Sn8OzAeaGmO2GGNuMMZEGmMiM93WD5hurT2U6bnqwFxjzFJgIfCjtXZqfsVZlGWe9XgMeNYYrhkwwBXJekxcNu+NZ/BHi7j6vYWkp1vev7YT71zT8YjEJYM/XXIzYq3oq2NZHRrK6nr1mDJnjqc4RUTkxJNvMy/W2iv8uOcD4IOjnosC2uRPVMVIejoRkydTJyiIxMREHixdmkPJyVxWsaL/RbJZSExJ481f1jN29npKGMMDvZpyw+kNKF0y+74qmZeHxowZk+U9ERER3PLnn1ydkMDFpUuzLSmJj3v35o1TTslzrCIicmIq8N1GkkebN8Ntt9GmYUMiIyMZPHgw48aNIyYmJk/DWWuZvmIHIyetYMu+BPq2qckjvZsRUd7jjp+UFFfTAiwLCaFW584sGDuWce++m+dYRUTkxKazjYqSnTth8mS49lp3/fff0Lq12/7swfpdcQyfuII5a3bRtHo5nriwJaeeVNl7vGvXwvnnw2uvQe/e3scTEZETis42Kg5efdWdRXTOOa5Xy8knexouLimV0T+v5b250YSUDGJonxZcdWo9SgV5LIVKS3Pt++vVg5YtoWxZb+OJiIhkopmXwm7WLKhc2SUqsbGwdSs0a+ZpSGstE5Zu46kfV7IzNolLO9TmwV7NqFqutPd4R48m5e23Oa98eT77+mtP9TciInJiO+7HA0gAxMfDgAHw9NPuulw5z4nLypiDXD5uAXd9sYTq4SGMv7Uroy5r4y1xsdYd9AhQrx5/JSay5LffDvedERERCSDNvBQ2iYnw5Zdw9dWuluWvv1zC4rFV/oH4FF6euYaP5m+gfGgpHujVjP91rENQCW/1Mhw8CJdcAv36ETpkyBG9XDLk+XBGERE5oWnmpaj44gtXkDt3rrtu185T4pKebvnyj02c/eJsPpq/gSu71GPWfWdxRee63hKXjKS3XDkoXx5CQrLttpunwxlFRESyoYLdwmDtWreT6LTT4Kqr4KST4IwzPA+7ZPN+hv2wjKVbDtCxXkU+uqgzLWuW9x7v9OnwwAOuHqdiRfjmGwAiwK9uuyIiIl4oeSlo1sIVV7h+KEuWuF06HhOXPXFJPD91NV8u2kzVcqV56X9t6NeuFsbjlmqsdUtZ1aq52aA9e1zykklGt12vfWdERESyo5qXgmAtfPcdXHABlC7t+rVUrQoREZ6GTU1L59PfN/Hi9NXEJ6dx/ekNuOPsRpQLKeU93rvugrAwePbZw895TYZERERyoD4vhcn8+a7I9e234cYbPfdrAfg9ag/DJixn1fZYzmhchWF9W9CoWrkABItLUhIToUSJw0mLEhcRESkgSl6Ol3373LJQ9+7QtavrlHvuuZ6H3X4gkWemrOSHJduoVSGUN69sz3kta3hfIlq1CgYPhvfeg0aN4K23lLCIiEihoOTleLntNpg61Z1JVKaMa5vvQXJqOu/9Fs1rP60lNd1yZ4/G3NLtJEKDsz9AMVfCw2H7dti0ySUvSlxERKSQUM1Lfvr9d7dzqEoVWL/edcht29bzsLNX72TExBVE7T7EOS2q8/gFLahbOcx7vO+842J++213nZ7ulopEREQKgGpejreYGLdr6K674IUXXBLj0ea98YyYtIIZK3bQoEoZ3r+uE92bVgtAsD7bt0NUFCQkuN1ESlxERKQQ0sxLIKWkwG+/wVlnueuJE93jct4KZxOS0xj7y3re/GU9JUsY7ji7MdefXp/SJT0uEe3dC3ff7YqGzzwTUlPdVm0tEYmISCGgDrt5FBMTQ7du3di+ffuxb376aejRwy0RAfTt6ylxsdYyddl2ur/wM6/9tJYzG5bnpyHduOWsk7wnLgAhIbBgAaxY4a5LllTiIiIihZ6Sl2MYOXIkc+fOzf6QwU2bIKP9/e23u/4tDRt6/tx1O+O4+r2FRH6ymLj9u9n5xSMEL/qEiPLezjhizhy48kpIS3N9W5Ytg8hIz/GKiIgcL1o2ykZoaOixDxlMToYGDaBDB5gwISCfG5eUyuif1vLu3GiSE+LYP+djYv+aDDY96xhy6/PP4bHH4OefoV69gMQsIiKSH7RslEvZHjIYFeWazAEEB7sdOqNHe/48ay3f/7WVs0fN5q05UfRvX4spt3ehb7NwwkJDjowhNwcdpqS45azPP3fXAwbA8uVKXEREpMjSbqNsREREZH3I4K+/wuWXw4wZ0LOn534tACu2HeSJCctZuGEvJ9cuz1tXdaBdXXdmkOeDDoOC3KxQu3buDCVjXK2LiIhIEaWZlxxkHDK48OefGXb55a5o9+KLYdw46NbN8/j745MZ+sMy+oz+lbU7Y3m2f2u+v/W0fxOXzDEsWLCAyMhI/wqHo6Ph+ushLs5td545E8aO9RyviIhIYaCaF3906wa7d8M//wSk90l6uuWrRZt5ftpq9scnc2WXetx7ThMqhAUHIFhg3jzo1QsmTXJboEVERIogNanzYsQId/pzABKXJZv3M+yHZSzdcoBO9Ssy/MJTaFEz3HuM48e7JnO33urOTtq8GcqX9z6uiIhIIaPkxR8BWCLaHZfE81NX8dWiLVQrV5pXB7TlwjY1vR+gmOGrr2DDBrftuUQJJS4iIlJsqeYlj/xtXpeals77v0XTfdRsxv+5lZvPbMjP953FRW1reUtcYmPh4YdhyxZiYmK4YPNmtn/zjVr6i4hIsae/6fLomM3rgAVRe+gzei7DJ66gbZ0KTL37TB7u3ZyypQMw4bV7N7z6KkyezMiRI5m6YAEjnn7a+7giIiKFnAp2c8mf5nXbDyTy1OSVTFy6jVoVQnm8TwvOa1nd+xLR33/D9Olw330A1A8JYWNSUo6xiIiIFFVqUhcg2Tavi44mKTWNN2av4+wXZzNt+Xbu6tGYmfd2o1erGoGpbfn8c3juOXegIjA/OjrbWERERIorFezmUnbN61YeKMGIT34levchzm1Rncf7tKBOpTBvH5aeDh9+CG3aQPv2rq3//fdDpUo5xpKrJnYiIiJFjGZe8iBz47irb7mH30q04rr3/8AAH1zXiXFXd/SeuIBrMvfQQ/Duu+66TJl/E5esYvG7iZ2IiEgRppqXPEpITmPsL+t585f1lCxhuLNHY64/rQHBJT3mg7t3wwcfwJAhrpV/VJQ7/DFQW6pFRESKiONe82KMec8Ys9MYsyyb188yxhwwxizx/QzN9FovY8xqY8w6Y8xD+RVjXlhrmboshp4v/cJrP62lV8sa/DzkLCK7neQ9cQH47js327Jkibtu2FCJi4iISCb5WfPyAfA68FEO9/xqre2T+QljTBAwBjgH2AL8YYyZYK1dkV+B+mvdzjiGT1zOr2t307R6OT6/qQunnlTZ+8C//QaJidCjhzuT6MwzoWlT7+OKiIgUQ/mWvFhr5xhj6ufhrZ2BddbaKABjzBfARUCBJS+xiSmM/nkd782NJjQ4iCf6tuDKLvUoGRSAmZb0dNfSv1w5l7wEBSlxERERyUFB7zY61RizFNgG3GetXQ7UAjZnumcLcEp2AxhjBgODAerWrRvwAH9du4shXy1lZ2wS/+tYmwd6NaNK2dLeBk1Jgfffh6uvhpAQ+OYbiIgITMAiIiLFXEHuNvoTqGetbQOMBr73PZ9VgUe2VcXW2nHW2o7W2o5Vq1YNeJBVypamVsVQvru1K89f2sZ74gIwfz7cfDN8+627btwYypbN8lZ/jyEQERE5URRY8mKtPWitjfM9ngyUMsZUwc201Ml0a23czEyBaB4RzvhbutKubkVvA23ZAj/+6B6feSYsWAADBx7zbf4cQyAiInIiydet0r6al0nW2lZZvFYD2GGttcaYzsA3QD0gCFgD9AC2An8AA31LSjk6nlulc+2yy2DOHNi40S0VHYM/xxCIiIgUZwWxVfpzYD7Q1BizxRhzgzEm0hgT6bvlUmCZr+blNWCAdVKB24FpwErgK38Sl0Jp2jTYtcs9HjXKzbb4kbhAzscQiIiInMjyc7fRFcd4/XXcVuqsXpsMTM6PuI6bzZuhTx/Xzv/pp6FevVy9Xa3/RUREsqbjAQIpIQEm+3KuOnXczMuwYXkeTq3/RURE/kvHAwTSY4/Bs8/C+vW5nmkRERGRIx33mpcTxpo1sG6dezxkCMyYocRFREQkHyl58SIpyW17vu8+d12xInTvXrAxiYiIFHNKXnLLWpg92/1aujR88gm8+WZBRyUiInLCUPKSW19+6WZXfvrJXffsCdoBJCIictwoefHHgQOwwncu5CWXwIcfanlIRESkgBT0wYxFQ9++sHs3LFsGpUq5AxVFRESkQCh58cdTT0FoKJTQRJWIiEhBU/LijzPOKOgIRERExEdTCSIiIlKkKHkRERGRIkXJi4iIiBQpSl5ERESkSFHyIiIiIkWKkhcREREpUpS8iIiISJGi5EVERESKFCUvIiIiUqQoeREREZEiRcmLiIiIFClKXkRERKRIUfIiIiIiRYqx1hZ0DAFjjNkFbAzQcFWA3QEaqzjT9+QffU/+0ffkH31P/tN35Z/C+j3Vs9ZWPfrJYpW8BJIxZpG1tmNBx1HY6Xvyj74n/+h78o++J//pu/JPUfuetGwkIiIiRYqSFxERESlSlLxkb1xBB1BE6Hvyj74n/+h78o++J//pu/JPkfqeVPMiIiIiRYpmXkRERKRIUfLiY4ypZIyZYYxZ6/u1Yjb3bTDG/GOMWWKMWXS84yxo/n5PvnuDjDF/GWMmHc8YCwN/vidjTIgxZqExZqkxZrkxZnhBxFqQ/Pye6hhjZhljVvq+p7sKItaClIs/n94zxuw0xiw73jEWJGNML2PMamPMOmPMQ1m8bowxr/le/9sY074g4ixofnxPzYwx840xScaY+woiRn8peTnsIeAna21j4CffdXa6W2vbFqVtZQGUm+/pLmDlcYmq8PHne0oCzrbWtgHaAr2MMV2OX4iFgj/fUyowxFrbHOgC3GaMaXEcYywM/P3/3QdAr+MVVGFgjAkCxgDnAy2AK7L4/XE+0Nj3MxgYe1yDLAT8/J72AncCo45zeLmm5OWwi4APfY8/BC4uuFAKNb++J2NMbeAC4J3jE1ahc8zvyTpxvstSvp8TrQjNn+8pxlr7p+9xLC4hrnW8Aiwk/Pr/nbV2Du4voBNJZ2CdtTbKWpsMfIH7vjK7CPjI9/+5BUAFY0zE8Q60gB3ze7LW7rTW/gGkFESAuaHk5bDq1toYcH9YAtWyuc8C040xi40xg49bdIWHv9/TK8ADQPpxiquw8et78i2tLQF2AjOstb8fvxALBX9/PwFgjKkPtAP0PUmGWsDmTNdb+G9y6889xV2x+g5KFnQAx5MxZiZQI4uXHs3FMKdZa7cZY6oBM4wxq3z/2ik2vH5Pxpg+wE5r7WJjzFkBDK1QCcTvJ2ttGtDWGFMB+M4Y08paW6zqFQL0/zuMMWWBb4G7rbUHAxFbYRKo7+kEZLJ47ugZTH/uKe6K1XdwQiUv1tqe2b1mjNlhjImw1sb4phN3ZjPGNt+vO40x3+Gm4opV8hKA7+k04EJjTG8gBAg3xnxirb0yn0IuEIH4/ZRprP3GmNm4eoVilbwE4nsyxpTCJS6fWmvH51OoBSqQv59OMFuAOpmuawPb8nBPcVesvgMtGx02AbjG9/ga4IejbzDGlDHGlMt4DJxLMfuLxg/H/J6stQ9ba2tba+sDA4Cfi1vi4gd/fj9V9c24YIwJBXoCq45XgIWEP9+TAd4FVlprXzqOsRUmx/yeTmB/AI2NMQ2MMcG4P3MmHHXPBOBq366jLsCBjGW4E4g/31PRYa3Vj2vUVxlXxb/W92sl3/M1gcm+xw2Bpb6f5cCjBR13Yfyejrr/LGBSQcddGL8n4GTgL+BvXBI8tKDjLqTf0+m46e2/gSW+n94FHXth+558158DMbiCyy3ADQUd+3H6fnoDa4D1GX8uA5FApO+xwe20WQ/8A3Qs6JgL6fdUw/f75iCw3/c4vKDjzupHHXZFRESkSNGykYiIiBQpSl5ERESkSFHyIiIiIkWKkhcREREpUpS8iIiISJGi5EVE/GKMiTv2XTm+/xtjTMMcXr/WGFPTy2f4GUeoMeYX30F1/r7ndmPMdZmuRxljzs6fCEXkWJS8iEi+M8a0BIKstVE53HYtrm9JfrseGG/d0Qz+eg932m6G0eR8orqI5CMlLyKSK74upS8YY5YZY/4xxlzue76EMeYNY8xyY8wkY8xkY8ylvrcNwtcV1ncY5QeZ3n+P776OwKfGmCW+2ZGhxpg/fPeN83XaxRjTyRjztzFmfkYcmcZ9wfeev40xN2fzn5A5lrN8szBfGWPWGGOeNcYMMsYs9MV2EoC1Nh7YYIzp7LveCFQ2xmR1FpGI5DMlLyKSW/2BtkAb3JEGL/jO2+kP1AdaAzcCp2Z6z2nAYt/jtkAta20ra21r4H1r7TfAImCQtbattTYBeN1a28la2woIBfr43v8+riPoqUDm2ZMbcG3fOwGdgJuMMQ0yB+5ri97QWrsh09NtgLt8cV8FNLHWdgbeAe7IdN8i4IxM13/6/rtE5DhT8iIiuXU68Lm1Ns1auwP4BZcsnA58ba1Nt9ZuB2Zlek8EsMv3OApoaIwZbYzphWtFnpXuxpjfjTH/AGcDLX1nQZWz1s7z3fNZpvvPxZ1fswT4HddSv/FRY1bBtT3P7A9rbYy1NgnXNn267/l/cMlYhp0cuax19LWIHCcn1KnSIhIQJpfPAyTgThjHWrvPGNMGOA+4Dfgfrg7l8EDGhABv4M6g2WyMecL3/pw+wwB3WGun+RNHJkmZHqdnuk7nyD8jQ3zvz+5aRI4TzbyISG7NAS731ZhUBc4EFgJzgUt8tS/VcYdyZlgJNAIwxlQBSlhrvwUeB9r77okFyvkeZyQYu40xZYFLwSU+QKzvZGBwJ+NmmAbcYowp5fucJr7T3//le3+QLznKrSYceYr80dcicpxo5kVEcus7XD3LUtxpzw9Ya7cbY74FeuD+Ql+DW7o54HvPj7hkZiZQC3jfGJPxj6eHfb9+ALxpjEnwjf82bulmA/BHps+/AXjbGHMImJ3pM97BLfP86Svu3QVcnEX803FLXDNz+d99GjAcwJcgNcLVwYjIcaZTpUUkYIwxZa21ccaYyrjZmNN8iU0orgbmtFxuUc72M3yPHwIirLV35eL97YB7rbVX5fU9xph+QHtr7eO5i15EAkEzLyISSJN8RbXBwEhf4S7W2gRjzDDcrMsmj59xgTHmYdyfXxtx/WH8Zq39yxgzyxgTlItEqgpuiStDSeDF3HyuiASOZl5ERESkSFHBroiIiBQpSl5ERESkSFHyIiIiIkWKkhcREREpUpS8iIiISJGi5EVERESKlP8DuPINAdXsydgAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# First let's make a plot in the log-transformed space\n", + "plt.figure(figsize=(9,6))\n", + "\n", + "plt.xlabel('log(stage (m))')\n", + "plt.ylabel('log(discharge (cms))')\n", + "plt.plot(loght,logQ,'k*')\n", + "plt.plot(p_x,p_y)\n", + "plt.plot(p_x,p_y_lower,':r')\n", + "plt.plot(p_x,p_y_upper,':r')" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAFzCAYAAAAKZcKfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABddElEQVR4nO3dd3hVxdbH8e8kkAaEXkLvvUhHEFFEQawgWEBULIhcRa9i7+Br5doVwd5FBRErdgHpvTdDJ6GXQHrOvH9MQlFITpKzU3+f58lD9imzJ+d6yWJmzVrGWouIiIhIQRCU3xMQERERyaDARERERAoMBSYiIiJSYCgwERERkQJDgYmIiIgUGApMREREpMAokd8T8EelSpVs3bp183saIiIiEgALFy7cY62tfLLnCkVgUrduXRYsWJDf0xAREZEAMMZsPtVz2soRERGRAkOBiYiIiBQYCkxERESkwFBgIiIiIgWGAhMREREpMBSYiIiISIGhwEREREQKDAUmIiIiUmAoMBEREZECQ4GJiIiIFBgKTERERKTAUGAiIiIiBYYCExERETnG58vX2yswEREREefAAejeHT79NN+m4FlgYoypZYz53Riz2hiz0hhze/rjjxljthtjlqR/9fVqDiIiIpINISFQqhREROTbFEp4OHYqcJe1dpExpgyw0Bjzc/pzL1hrx3p4bxEREfHXzp1QrpwLSKZNA2PybSqerZhYa2OstYvSv48DVgM1vLqfiIiI5EBCgtu+ufFGd52PQQnkUY6JMaYu0BaYm/7QrcaYZcaYd4wx5U/xnmHGmAXGmAW7d+/Oi2mKiIgUP+HhcNdd8J//5PdMADDWWm9vYExp4E/g/6y1k40xVYE9gAXGAFHW2uszG6NDhw52wYIFns5TRESkWFmxAtLSoE2bPL+1MWahtbbDyZ7zMscEY0xJYBLwsbV2MoC1dudxz78JfOvlHEREROQfrIUhQ9y2zcKF+b59czzPAhNjjAHeBlZba58/7vEoa21M+mU/YIVXcxAREZGTMAYmToQSJQpUUALerph0A4YAy40xS9IfewC4yhhzGm4rZxNws4dzEBERkQw//gjLlsE990Djxvk9m5PyLDCx1s4EThaGfe/VPUVERCQTX34JixbByJEQFpbfszkpT3NMREREpABISYGSJeGNNyA+vsAGJaCS9CIiIkXb669Dt25w6JDLKYmMzO8ZZUqBiYiISFFWpw7UrevKzRcCCkxERESKGmth5Ur3/QUXwOefF+jtm+MpMBERESlqnn4aOnSAdevyeybZpuRXERGRouamm6B0aWjUKL9nkm1aMRERESkKjhyBF14Anw8qVYLbbitwxdP8ocBERESkKPjiCxg1CmbPzu+Z5Iq2ckRERIqCa6+Ftm3zpSlfIGnFREREpLDauBHOPhu2bHHbNoU8KAEFJiIiIoXXoUMuKNm1K79nEjDayhERESlsdu2CKlXcCsmaNa7cfBGhFRMREZHCZOFCaNDAFU2DIhWUgAITERGRwqVlS7juOjjzzPyeiScUmIiIiBQGX3zhOgOHhsIrr0C1avk9I08oMBERESnoVq+GK6+El17Kk9tZa1kdcyhP7vVPCkxEREQKumbN4Jdf4O67Pb9VcqqP/05cwiWv/UX07sOe3++fFJiIiIgURElJcMMNMG+euz77bCjh7WHaQ4kpXPfuPKYs2cHt5zSiXqVSnt7vZBSYiIiIFEQHD8IffxwLTDwWczCBy9+YzbyN+3j+8jb85+yGmHzotaM6JiIiIgXJgQNQtqyrU7JsGZTyftVidcwhhr47n8NJqbw3tBNnNKrk+T1PRSsmIiIiBcXevdCxIzz+uLvOg6Dkrw17uPwN1/jvi+Gn52tQAloxERERKTgqVIBLLoHevfPkdpMXbePeScuoX6k0713fkaiy4Xly38woMBEREclvM2ZAw4YQFQVjx3p+O2str/y2ged/XkfXBhV5Y0h7IsMKRgVZbeWIiIjkp8OHoV8/uOOOPLldSpqPeyct4/mf19G/XQ3eG9qpwAQloBUTERGR/FW6NHz9tatV4rG4xBRGfLyIGev3MPKcRvy3V6N8OXmTGQUmIiIiec1aeOghOO00GDgQunXz/JYxBxMY+u58Nuw6zLMDWnN5h1qe3zMnFJiIiIjkteRkV6Pk0CEXmHhs1Y5DXP9ewTgOnBUFJiIiInnlyBFXvTU0FH76CSIiPL/l72t3cevHi4gML8kXw0+nWVSk5/fMDSW/ioiI5IWUFOjVC4YOddelSoHH+R0fz93Mje8voG6lUnw1oluBD0pAKyYiIiJ5o2RJuOIKqFfP81v5fJZnflzD+OnR9GxahVeuakup0MLxK79wzFJERKSwmj/fBSWnnZYnR4ITU9L478Ql/LAilmtOr8MjFzanRHDh2SBRYCIiIuKVtDS45hpX0XXmTM+3bvYcTuLG9xewdNsBHrqgGTecUa/AHQfOigITERERrwQHw+TJUL6850HJhl1xDH1vPrvjkhg3uD19Wlbz9H5eUWAiIiISSNbCmDFu++b++/OkcNrM9Xu45eOFhJYI5rNhp3NarXKe39MrCkxEREQCbf16dyzYWs9XSj6dt4WHpqygYeXSvH1dB2qW9/4IspcUmIiIiARCXBwkJUGlSvDOOy4w8TAoOf7kTY/GlXl1UFvKFKCeNzmlwERERCS3rIU+fVwgMmOG28bxUEJyGndMXMy0lTu5ukttHruoRaE6eZOZovFTiIiI5Cdj4N574cEHPd+62XUokSsmzOanVTt5+MLmjLmk5b+CkpiYGHr06EFsbOwJ3xcGWjERERHJqT//hMOH4YIL4OKLPb/dyh0HufH9BRxMSOHNIR3o1bzqSV933333MX36dO677z4iIiKYOXMmo0eP5vXXX/d8jrllrLX5PYcsdejQwS5YsCC/pyEiInKMtXDmmZCQAPPmQZC3mxA/r9rJ7Z8tpmx4Sd6+tiPNq/+7vHx4eDiJiYmZjhMWFkZCQkKmr4mJieHKK69k4sSJVKsW+GPHxpiF1toOJ3tOWzkiIiLZYS2kprotm0mT4OefPQ1KrLVMmP43wz5cQKMqpfn6P91OGpRkvPZUIiIiGDx4MBs3bszynmPGjDm6ypLXFJiIiIj4y1oYNgyuv959X6WKK57mkeRUH/dNWs6T36+hb8soJt58OlUiw075+o0bN9KwYcN/PR4aGkpiYiKRkZGZroCEh4djjOHQuHG09vkYN24cxhjCw8MD8vP4Q4GJiIiIv4yB2rWhTh3Pb7X/SDLXvjOPiQu2clvPhrxyVVvCSgZn+p6oqChSU1MBCAkJAaBMmTLMnTuX4cOHZ5kAGx0dzaBBg/gpPJzWZG+VJVCU/CoiIpKVbdtcnZJmzeDhhz2/3YZdh7nx/fnsOJDI85e3oX+7mn6/t23btvTt25dhw4YxYcIEYmJiaNOmDa+99lrmb1y8mKgXXqBcqVLsTUri87Awkv1YZQk0BSYiIiKZsRb69YPERFi61PMk1xnrdzPi40WElgji02GdaV+nQrbeP3ny5KPfZxmMHG/FCvj9d9KaN2f48OEnBDZ5SadyREREsrJokSua1qqVp7f5YPYmHv9mFY2qlOata/OgvHxamiuf37Spuz58GEqX9vaeZH4qRysmIiIiJ/Paa5CSAnfcAe3aeXqr1DQfj3+zig/nbKZXsyq8eGVbSofmwa/oUaPgvfdg7VqXyJsHQUlWFJiIiIj8k7Xwxx/uWPDtt3tazfVgfAr/+WQRMzfs4eYz63NPn6YEB3lbPfaokSOhRQsXlBQQ2soRERHJcOSIa8RXoYLLKSlZEoIzPwmTGxt2HeamDxawbX88/9evFZd3qOXZvY6aPBnmzoVnnvH+XqegrRwREZGs+Hxw3nkQEgK//QZhp64XEgh/rtvNrZ8sIiQ4iE9v6kKHutlLcs2x2bNdo8H4eIjwOIclB1THREREBNxpm9tvd3kXHm3dxMTEcGaPHrz4/VKGvjuPGuXC+frWbgELSk7ZsO/AAcioRfLUU67HTwEMSkCBiYiIFHfffgvTprnvL7/cNeTzyONjnmBVqdN4cfo2ejWryqRbugb05M1JS8lbCxddBJde6laFSpSA0NCA3TPQlGMiIiLFV1oadOoE5crBL794tlISHh5OclAolS+9n7BaLTkw6zMOzviYsLDQLBvq+Tv+yZr3HW3YN326+9m6d8/1vQJBTfxERESOl5rqjgIHB7sVk6lTPT158+PcFTS8ZQKhUY3ZPfVZUhZOZvDgQQEr9Z5RSj4ifXumdHg4X7doQWzGysmZZxaYoCQrSn4VEZHiJTXVbdfUrQvjx0NUlKe3+25ZDKO+WI8JCmLnR/di9m8lOTk5oKXeo6KiiIyMJDExkbCwMBITE2lw6BBlt2wJyPh5SSsmIiJSvJQoAZ07uy0cD/l8lud/Wst/PllE8+qRNNr4FTf0O5c5c+b41VAvu3bu3MnoAQOY99NPDLvlFh5v2xZeeSWg98gLyjEREZHiYcYMV0isSRPPb3U4KZU7Jy7hp1U7ubxDTcZc2pLQEt7VQwHcyZs6dWDAAHj7bW/vlUuqYyIiIsVbUhIMGgRt2ricEg9t3nuEYR8sZMPuwzx6UXOu61oX42H+Cta6/Jhy5Vx5+a5dvbtXHtBWjoiIFF0+n/szNBS++QY+/NDT281cv4eLX/2L2EOJvD+0E0O71fM2KNm5E3r2dKtB4LogV63q3f3ygAITEREpmhISoH9/eP11d33aaVC+vCe3stby1oxornlnLtUiw5h6azfOaFTJk3udICwMDh6EPXu8v1ce8WwrxxhTC/gAqAb4gAnW2peMMRWAiUBdYBNwubV2v1fzEBGRYiokxG1xZKyaeCQxJY0HJi9n8uLt9GlRjf9d3oZSXnYGthYmTXKrI2XLwoIFrmptEeHlT5IK3GWtbQZ0Af5jjGkO3Af8aq1tBPyafi0iIhIY8+a5RNDgYNew7tZbPbtVzMEELh8/m8mLt3PnuY15fXA7b4MSgN9/h4ED4eOP3XURCkrAw8DEWhtjrV2U/n0csBqoAVwCvJ/+sveBS72ag4iIFDP79kGvXnD33e7aw/yOeRv3cdErM/l712EmDGnPyHMaERTkYT5JSor7s2dPl8B79dXe3Ssf5UmYZYypC7QF5gJVrbUx4IIXoEpezEFERIqBChXg00/h6ac9u4W1lg9nb2LQm3MoE1aSKf/pxnktAlMo7ZSmTXPHnDMKpl1wQZFbKcng+U9ljCkNTALusNYeysb7hhljFhhjFuzevdu7CYqISOEWFweXXea2OMD90q5Y0ZNbJaakce+kZTz89Uq6N6rElP90o1HVMp7c6wT16kHDhp6uABUUngYmxpiSuKDkY2vt5PSHdxpjotKfjwJ2ney91toJ1toO1toOlStX9nKaIiJSmPl8sH49/P23p7eJPZjIFRPm8PmCbdzWsyFvX9uRsuElPbxhLLzxhvu+cWP46SeoVcu7+xUQngUmxh3cfhtYba19/rinpgLXpn9/LfC1V3MQEZEibNEi1/embFlYuBBuvNGzW83ftI8LX5nJ+p1xvHF1O+46r4m3+SQAr74Kd90FW7d6e58CxssVk27AEKCnMWZJ+ldf4GngXGPMeuDc9GsRERH/rV3r+t2MHeuuS3qzcmGt5f1Zm7hqwhxKhwbz1Yhu9GnpYdM/a2HvXvf9I4+4gKsYrJIcz7MzTdbamcCpwslzvLqviIgUA02auMJpl1/u2S0SU9J44KvlTF60nZ5Nq/DCFad5u3UDcNNNMH8+zJ3riqc1bert/QqgopnSKyIiRc/u3S7JNTraXd90k9vG8cC2/fEMeGMWkxdt5/ZzGvHWNR28D0rANeC77jpXQr+YUhM/EREpHA4ehDlzYOVKqF/fs9v8tWEPt36yiNQ0y1vXdKBXc497z7zzDpQoAddcA336uK9iTIGJiIgUbMuXQ6tW7rjshg0QHu7Jbay1jJ8ezbM/rqFB5dKMH9Ke+pVLe3Kvo3w+V8G1VCkYMqRYHAfOigITEREpuH74Afr2ha+/hosv9iwoOZyUyt1fLOWHFbH0bVWNZwe0obSXpeXXrYNq1SAy0pXNL11aQUk6BSYiIlJw9eoFzz0HvXt7dosNuw5z84cL2LjnCA/0bcpN3etjvAwSDhxwJ4oGDIA33/QsT6awUvKriIgULOvWuSZ1cXHuGPCoUZ4lg/6wPIZLXp3JgfgUPrqxM8PObOBdUGKt+7NcOXei6NFHvblPIafARERECpYtW2DGDJdP4pHUNB9P/bCaWz5eRMOqZfh25Bl0bVDJs/uxdSt07eqSdwGuugpq1vTufoWYtnJERCT/WeuSXFu3dts30dEQEeHJrXbHJTHy08XMjt7LoM61efSi5oSWCPbkXkeVLg1JSW4bRzKlFRMREcl/zz4LHTu6iq7gWVCycPM+LnxlBou27GfswDY82a+Vd0FJcrLrdePzQfnysGBBsT8K7A+tmIiISP678UYoU8Y1q/NARmn5J75bTY3y4Xw1ohPNq0d6cq+jpkyBW25xNVfOOw+CtBbgDwUmIiKSP6ZNg08/dQXGKlaEESM8uc2RpFTun7ycqUt30KtZFf53ucel5ffudT/PwIFQuzZ06eLdvYoghW8iIpI/1q51HYL37fPsFht2xXHpa3/x7bId3N27CROGeFxa/n//g+bNITbW1SVRUJJtWjEREZG8ExcHmzdDy5Zw220wbJhrVueBb5bu4N5JywgvGcwH13fmjEYenrrJcP75EBPjckokRxSYiIhI3rn6ardKsn69C0g8CEqSU308+f1q3pu1ifZ1yvPaoHZUK+tN8APAl1/CqlXwyCNutWTsWO/uVQwoMBERkbzzxBOwa5dnqyQ7DiQw4uNFLNl6gBvOqMd95zelZLDHWQu//gpLlsB990FIiLf3KgYUmIiIiHeshccec98//rhrxueRP9ft5o7PFpOSZnl9cDv6tory7F6sWeOCkPr14YUX3IkbBSUBocBERES8tW2b+9NaTxrVpfksL/2yjld+30DjKmUYd3U7b7sCJyfDuee6PJkffvBs9ae4UmAiIiKBt3SpSwCtXRvGj4fgYE+Ckt1xSdz+2WJm/b2Xge1rMvqSloSHeFQwLSnJ9ewJCYEPP4RGjby5TzGn48IiIhJYCQmuG/Dtt7vrEiU8CUrmRu/lgpdnsHDzfp4d0JrnBrbxLijZts2Vy//0U3d91llQo4Y39yrmtGIiIiKB4fO5XIvwcPjsM2jWzKPbWN6Y/jdjp62lbsVSvH99J5pFeVzFtVo1aNMGojzMWxFAKyYiIhIIe/bAGWfA55+767POgqpVA36bfUeSueH9+Tz741rObxXF17d28y4o2bvX1Vo5csSt+nz+ufu5xFMKTEREJPciI91XCe8W4hds2scFL8/grw17GXNJC169qi1lwjys4rpiBbz1Fsya5d095F+0lSMiIjnj87lf3Fdf7boB//CDJ7kkPp/lzRnRPDttLTXKhTPplq60qlk24PcBIDUVFi6Ezp2hRw9XpbZKFW/uJSelFRMREcmZhQth+HD44AN37UFQciA+mZs+WMBTP6zhvOZV+XbkGd4FJeCqt555pgtIQEFJPtCKiYiIZE9cHJQpAx07wpw57k8PLNy8j9s+Wczuw0k8fnELrjm9DsaD4AdwtUlCQuC//3VJrnXqeHMfyZJWTERExH8//AB167oS7ACdOgV8pcTns4z7428uHz+HEsFBTLqlK9d2retdUPKf/8CAAa4AXOXKcMUV3txH/KIVExER8V/79tCnj2fHZvceTuKuL5byx9rd9G1Vjacva02klwmuAE2aQLlyLmcm2KM6KOI3Y63N7zlkqUOHDnbBggX5PQ0RkeJp/Xp45x148klP8kgyzI3ey8jPFrM/PoWHL2zO1Z1re7NK4vO5/jbt2+v4bz4xxiy01nY42XPayhERkcx9+y1MmACbNnkyvOt1s56r3pxDREgJvhrRlSFdPMwnSUhwZfK/+MKb8SVXtGIiIiL/lpLiTqY0bOhyL3budNVPA2znoURu/2wxc6L3celp1XmiXytKh3qUZfD779C9u6u1sns3VKrk6QqQnFpmKyZZ/q9vjKkCdAOqAwnACmCBtdYX0FmKiEjBcd11MH06rFkDpUp5EpT8vnYXd32+lITkNJ4b0JoB7Wt6t0oyfz707AmvvQYjRrgkVymQTrmVY4w52xgzDfgOOB+IApoDDwHLjTGPG2M8bk4gIiL54q67YOxYF5QcJyYmhh49ehAbG5vjoZNTfTz5/WqGvjufKmVC+ea2bgzsUMuboOTIEfdnx47wySdwww2Bv4cEVGYrJn2Bm6y1W/75hDGmBHAhcC4wyaO5iYhIXvH54OGHXQO+hx6Cdu3c1z+MGTOGmTNnMnr0aF5//fVs32bz3iOM/HQxS7cd5OoutXnoguaElfToJMynn7q6JAsWQM2acNVV3txHAuqUgYm19u5MnksFpngxIRERyQfGwJYtrrS8tf/KvQgPDycxMfHo9bhx4xg3bhxhYWEkJCT4dYuvl2znwa9WEGRg3OB2nN/K4069HTtCr14u2JJCI8tTOcaY240xkcZ52xizyBhzXl5MTkREPPbzz7B9uwtE3n3XnVY5yZZKdHQ0gwYNIiIiAoCIiAgGDx7Mxo0bs7zFkaRURn2xlNs/W0LTamX4/vbu3gUlX3wB99/vvm/YED76CCpW9OZe4gl/jgtfb609BJwHVAaGAk97OisREfHe/v1w2WXw6KPuOpPOwFFRUURGRpKYmEhYWBiJiYlERkZSLYuk2BXbD3LRKzOZtGgbI3s25LNhXahZPiKQP8WJ5s6F336D41Z3pHDxJzDJCJ37Au9aa5ce95iIiBQ2hw+7P8uXZ89HH3HemjV+JbPu3LmT4cOHM2fOHIYPH05sbOwpk2Gttbw9cyP9X5/FkeRUPrmxC3ee14QSwR6Uz5oxA1audN8/+STMnAlhYYG/j+SJLOuYGGPeBWoA9YA2QDDwh7W2vffTc1THREQkQFaudMdmJ0yASy5hxIgRjB8/nptvvjlHyawne//uuCTu/tKVle/VrCrPDmhNhVIhgf5JnIQEqF8fTj8dJk/25h4ScJnVMfEnMAkCTgOirbUHjDEVgRrW2mUBn+kpKDAREQmQpCS48Ubafv45S5KT//W0v8ms/0yGzRDZuDONhjxBXGIKD13QjKu9quC6ZQvUquXyYRYtgsaNoXTpwN9HPJGrkvTphdRSgTONMf2BHkDDwE5RREQ8s2UL3Hyzy7sIDYUPP+T7TZtynMwKJ0mGLR3J6SPGUr7fw1QsFcLUW89gyOkedQRevNg13vvgA3fdrp2CkiLEn8qv7wCtgZVARrVXC2jNTESkMFi+HD77zBUX69QJyHkya4bj318qqgGRvW9nR5n6XHN6HR7o28yb2iQ+HwQFQZs2MGoU9O4d+HtIvvOnIUEXa21zz2ciIiKBk5QES5ZA585wwQWwcSNUqHDCSzKSWYcNG8aECROIiYnJ1i1id+7k/NueZF2pVpCSSM3tPzP6khcD9zMc75tv4JFH4I8/oGxZGDPGm/tIvvMnMJltjGlurV3l+WxERCQw7roL3nvPBSSVK/8rKAGYfFyy6GuvvZat4fcdSabSpQ/w86qddK9fkf8NbEOVyP65nfWpVa7stmsOHXKBiRRZ/gQm7+OCk1ggCXdU2FprW3s6MxERyb6UFChZEh54AM47z5NmdTPW7+auz5dyIN4luF7frR5BQR7kkkybBhs2wH/+A126uKaC6gZc5PkTmLwDDAGWcyzHREREChJr4cYbXYLrRx9B9epw8cUBvUViShrPTVvL2zM30rBKad4d2pEW1T1cvXj/fVi9GoYNc8GWgpJiwZ/AZIu1dqrnMxERkZwzxpVgT0g4aa+b3FobG8ftny1mTWwcQ7q4BNfwEA8SXGfNgjp1oEYNeOMNCAlxQYkUG/4EJmuMMZ8A3+C2cgCw1upUjohIfrIWXnoJunZ1p20yesQEkM9neW/WJp7+cQ2RYSV497qOnN20SsDvA8CBA9CnDwwcCG+/DZGR3txHCjR/ApNwXEByfOM+HRcWEclvcXHwwgvw999HjwEH0s5DiYz6Yikz1u/hnKZVeGZAayqVDg34fdi2DWrWhHLl4OuvocNJ625JMZFlYGKtHZoXExERET/Nnw/t27sVhdmzISrwnXq/Xx7DA18tJzEljScubcngzrW9KZb288/uOPN338G558LZZwf+HlKoZFn51RjzvjGm3HHX5dOLromISF6bO9etjrz9truuXj2g+SRxiSnc9flSRny8iNoVIvhuZHdvysr70s9SdO8Od9zhAi0R/Osu3NpaeyDjwlq7H2jr2YxERORfYrZscV18a9eGceNg8GD/3hcTQ5cuXTj99NOz7CA8b+M++rw4g68Wb2Nkz4ZMuqUrDSp7UOp9/Hg44wxITnZdgJ999qR1VqR48icwCTLGlM+4MMZUwL/cFBERCYQpUwhu1Yp1M2YweswYGD4c0nvUZGXMmDHMnTuXOXPmMHr06JO+Jik1jad/WMMVE2ZTItjw5S1dufO8JpQMDiImJsYFRFkENdlStar7io8P3JhSZPjTXfga4H7gS1zS6+XA/1lrP/R+eo66C4tIcRUWFkajpCSeAW4AYo97PLMuwKfq/vvP966JPcQdny1hTWwcV3asxcMXNqdU6LF/e44YMYLx48dz88038/rrr+fsh/D54OWXoWJFGDIkZ2NIkZLb7sIfAJcBO4HdQP+8DEpERIqlP/6A557jiiuuYAVwSYkSxOJ/F+Do6Gj69etHcPCxWiPBwcH079+fjRs3kuazTJj+Nxe/8hd7Difx9rUdePqy1keDkvDwcIwxjBs3Dp/Px7hx4zDGEB4enrOf56uvXCVXkSycckvGGFPaWnsYIL1Pzr965Rz/GhERCZw3e/Wie1oaX6Rfp6amAhAfH+9XF+CoqCiqVq1KWlra0cfS0tKoWrUqKSGRXPXmHOZt3EfvFlV5sl8rKv7jGHB0dDSjRo1iypQpxMfHExERQb9+/Rg7dqx/P4DPB2+9BVde6U4Pffut63UjkoXMckW+NsYsAb4GFlprjwAYY+oDZ+O2dN7EbfGIiEhurVvnqpzWq8eFa9fy4P33Y777DuLjCQ4Opnfv3lSuXNnvfI+dO3dSr149OnbsCMC8+fNZEV+G81+aAcDYgW24rF2Nk564iYqKIjIyksTERMLCwkhMTPQrIDpqxQq45RZXIn/kSChTxr/3SbF3ysDEWnuOMaYvcDPQLT3pNQVYC3wHXGutDWA2lIhIMZacDD17Qps28N13RDVoQGjFikcDg+TkZOrUqZOtPI/juwfvikvk/knL+XXNLjpXj2TswDbUqpB5Au3OnTsZPnw4w4YNY8KECcTExGR+Q58PFixwx5lbt4Z586BdO7/nKwJ+JL8WBEp+FZEi69ChY6XXf/sNmjRxfWKA/v37ExUVdUJgcHyw4a/vlsXw0JTlxCencU+fpgztWtebbsBPPAGPP+4a7zVsGPjxpcjILPlVgYmISH5ZuxZ69HD9bq64IuDDH4hP5pGvVzJ16Q5a1yzL85e3oWGVAG+ppKW50vjlysHevfD993D11eoELJnKLDBRPRIRkfzSoAH07QvNmwd86N/X7OLeScvYdySZ//ZqzIizG1Ay2J/SVdlgLZx/PgQHu4BEx4ElABSYiIjkpYULYfRo+PRTVyTtncB2+IhLTOGJb1czccFWmlQtwzvXdaRljbIBvQc+HwQFuVWRK6+EUA8a+0mx5Vf4bIw5wxgzNP37ysaYen685x1jzC5jzIrjHnvMGLPdGLMk/atvzqcuIlIIHTwIS5ZAFnVIcuKvDXvo8+IMvli4lVvOasDU27oFPijZuhU6dz5Wk+T66115fG3dSID408TvUeBeXPVXgJLAR36M/R7Q5ySPv2CtPS3963t/JyoiUmht2ABfpFck6dmTmD//pMeIEX4f+82qLHx8ciqPfr2CwW/NJbREEF/e0pV7+zQltETwSV+fK5UrQ6lSkJIS+LFF8G/FpB9wMXAEwFq7A8gye8paOx3Yl6vZiYgUBQ8+CLffDull4Mc8+ywzZ848Ze+afxozZswpXz83ei99XpzBB3M2c323enw3sjvtapc/ySi5MH8+XHWVC0bCwuD33+HCCwN7D5F0/vTKmWet7WSMWWStbWeMKQXMtta2znJwY+oC31prW6ZfPwZcBxwCFgB3pXcrzpRO5YhIoRMT44qlVaoEu3ZBcjLhjRqdtH/NqfrenKrfTVhYGHsPxvHsj2t5b9Ym6lSM4LkBbehUz6MOvVOmwG23wa+/QuPG3txDipVc9coBPjfGjAfKGWNuAn7BVXzNiXFAA+A0IAb436leaIwZZoxZYIxZsHv37hzeTkQkHyQkQPv2ruIpQJUqULMm0dHRDBo0iIj0zsBZ9b051eu/mrmcvi/N4L1Zm7iua11+uL174IOSX3+Fzz933196qatKq6BE8kCWp3KstWONMefiVjmaAI9Ya3/Oyc2stTszvjfGvAl8m8lrJwATwK2Y5OR+IiJ5KjkZQkIgPBzGjnXByXGyW+b9n69PSrNsq9yFW75cS41y4Xx6UxdOb1Ax8D+HtfB//weHD8PAgS6xNafN+0Syya/jwumBSI6CkeMZY6KstRk1jfsBKzJ7vYhIobFsmcu7+Phj6N4dBg066cuyW+Y94/XdLr2WJ37exKagUlzTpQ739ml6tBNwwEyZAmed5YqlffIJlC2r0zaS5/zJMYkD/vmigxzLEYk+xfs+Bc4CKgE7gUfTr09LH28TcPNxgcopKcdERAq8w4ddgujo0dC2bcCGjU9OPZpLUqtCOM9e1sabVZL166FpU3jsMXj44cCPL3Kc3FZ+fR7YAXwCGOBKoBqumd87uGDjX6y1V53k4bf9uJ+ISOEwcyaMGwcffAClS8M33wR0+Fl/7+G+ScvZsi+e67rW5Z4+TYgICeAqic/nCr517AiNGsEvv8AZZwRufJEc8Cf5tY+1dry1Ns5aeyg996OvtXYiEOAzaSIihcimTTBnDmzfHtBh4xJTePCr5Qx6cy5BBiYO68JjF7cIbFAC8PTT0LWrS2wFOPtsd5JIJB/581+5zxhzOfBl+vWA455TUqqIFC9z57pmdX37uoqnl10W0MTQP9bu4oHJy4k9lMhN3etx57lNCA8JYKG05GTX0bhSJbj5ZtfJuFGjwI0vkkv+5JjUB14CTscFInOA/wLbgfbW2pleT1I5JiJSIFgL3bpBYqLbAglgYujB+BTGfLeKLxduo2GV0jw7oHXgC6VZ61ZIypVzTfeU2Cr5JMc5JsaYYOAWa+1Fp3iJ50GJiEi+W7LE1fCIiHDN98qVC+gv9R+Wx/Dw1yvZH5/MrWc35LZzGga2nHxCglvVMQZuucV1AVZQIgVUpjkm1to0oH1mrxERKdK2bIFOneDJJ911nTruGG0A7IpL5JaPFnLLx4uoUiaUr//TjVG9mwQ2KFmxAho0gB9/dNfXXAMXXBC48UUCzJ8ck8XGmKnAF6T3ywGw1k72bFYiIvlt3z6oUAFq14Z33gnoL3NrLV8u3MYT360mISWNe/o04abu9SkZ7FfDd//4fBAU5PJHund31WdFCgF/ApMKwF6g53GPWUCBiYgUTRMnwk03wbx5rrbH1VcHbOit++J54KvlzFi/h451y/P0Za1pULl0wMYHXCD11lvw558QGup+HpFCwp+S9EPzYiIiIvnOWpd70aOHO3FTtWrAhk7zWd79ayP/+2kdQQZGX9KCqzvXISjIg1yPihXdV1ycW/URKUT8OZUTBtwAtADCMh631l7v7dSO0akcEfHc6NGwapVLbg1wYujqmEPcN2kZS7cdpGfTKjxxaUuqlwtg75mkJBg1Clq3dis9cCzIEimAclv59UNgDdAbGA0MBlYHbnoiIgVAaKg7dZOS4hrxBUBiShqv/raBN/78m7LhJXn5qrZc1DoKE+iAISQE1qyByMhjjykokULKnxWTxdbatsaYZdba1saYksA0a23PTN8YQFoxEZGAS0iARx6Biy92yaEBXmGYE72XByYvJ3rPEfq3q8HDFzSnfKnABDwAbNzoetq8+qo7vpyaCiUCXBlWxCOZrZj4kwKekv7nAWNMS6AsUDdAcxMRyR8+H3z1Ffz+u7sOUFByMD6Fe79cxpUT5pDi8/HhDZ14/vLTAhuUAOzfD999B4sXu2sFJVJE+PNf8gRjTHngYWAqUBp4xNNZiYh4IS7ONd276y4oVcr9Ui9TJiBDW2v5bnkMj01dxf74ZG4+sz539Goc2HLyv/3m6pKMHAnt2sHWra55oEgR4s+pnLfSv/0TqO/tdEREPDRtGtx3nyuYdtZZAQtKtu2P5+EpK/h97W5a1ojkvaEdaVkjMEXYTvDRRzB7Ngwf7vJKFJRIEZRlYGKMCQUuw23fHH29tXa0d9MSEQmQvXtdYmi3bq7h3ooV0Lx5QIZOTfPx7l+beP7ndRgDD13QjOu61qVEoAqlpaXBG2/A+edD/frwwgsuIAlQcq5IQeTPVs7XwEFgIZDk7XRERAJs6FBYsMAli4aGBiwoWb7tIPdNXsbKHYc4p2kVHr+kBTXLRwRk7KN27XIrPLt2weOPB6wUvkhB5k9gUtNa28fzmYiIBMr27e6kSqlS8OyzkJzsgpIAiEtM4X8/reOD2ZuoVDqU1we34/yW1QJ3BHjvXpg0CYYNg6goWLQIGjYMzNgihYA/642zjDGtPJ+JiEgg7N0LLVq4o8DgSsq3bp3rYa21/LgihnOfn877szdxdZc6/HJXD/q2CnBdkjffhBEjYMMGd92okWqSSLFyyhUTY8xyXE+cEsBQY0w0bivHANZam/v/p4uIBMqhQ67AWMWKrhNw794BG3r7gQQe/XoFv6zeRbOoSN4Y0p7TapUL2PgsWuTqqLRvD3fcARdeqFUSKbZOWWDNGFMnszdaazd7MqOTUIE1EcnUV1/Bdde5pntNmgRs2JQ0H+/M3MiLv6wH4M5zGzO0WwCTW8EVRmvUyH399FPgxhUpwHJUkj4j8DDGdAFWWmvj0q/LAM2BPAtMREROKi0NgoPh9NNhwACXVxIgCzbt48GvVrB2ZxznNq/Koxc1D1xya0Zxt379XGG0yZOhXr3AjC1SyPlVkh5oZ9NfaIwJAhZYa9vlwfwArZiIyEmMGuVO2kyaFNBh9x9J5pkf1/DZ/K3UKBfOYxe34NzmgesyDMA337hS+JMmQf/+gR1bpBDIbRM/Y4+LXqy1PmOMah+LSP6qVs2tmASoR4zPZ/ly4Tae+mE1cYmp3NyjPref04iIkAD9dbd7N6xfD127uhySb7+Fvn0DM7ZIEeLP/+OijTEjgXHp1yOAaO+mJCJyErt3w003ueTQs85yKyYBsjrmEA9NWcHCzfvpWLc8Yy5tSdNqkVm/MTuuvhpWrYLoaChZEi64ILDjixQR/gQmw4GXgYdwp3R+BYZ5OSkRkX8pVQr+/hs2By697XBSKi/8vI73Zm2ibHhJnhvQmgHtawbu+O/MmXDaaa50/NixLh+mZMnAjC1SRPnTK2cXcGUezEVE5ESLFsGrr7raHhERsGSJ++WeS9ZavlkWwxPfrmL34SSu6lSbe3o3oVxEAEu9r1sH3bu7o8v33w+tVA5KxB9ZnnkzxjxrjIk0xpQ0xvxqjNljjLk6LyYnIsXc+vXwww9upQQCEpRs2BXH4LfmMvLTxVSJDGXyLV15sl+rwAQliYkwY4b7vnFj+PJLuP323I8rUoz4cypnibX2NGNMP+BS4L/A79baNnkwP0CnckSKDWuPnbIZMMBdHz4ckC7AR5JSefm39bw9YyMRIcHc3acpgzrVJjgogFVVR450qzubN0OVKoEbV6SIye2pnIwN0b7Ap9bafQEtvywiksFaeP55l5MxYIArxZ7LoMRay3fLY/i/71YTczCRAe1rct/5TalUOjC9c1i1yjXXq1HDJeRefLGCEpFc8Ccw+cYYswZIAEYYYyoDid5OS0SKjeRkeO01d+KmdGlXeKxixYAMvWFXHI9OXclfG/bSPCqSV65qS4e6FQIyNgAHD0LnznDFFfDWW1C7tvsSkRzzJ/n1PmPMM8Aha22aMeYIcIn3UxORYmHpUrjzTihf3pWVr5r7YmaHk1J5+df1vDPTbduMvqQFgzvXCcy2TWoq/PEH9OrlVko++cRVnhWRgMisiV9Pa+1vxpj+xz12/EsmezkxESnCdu6E2bPh0kuhY0dYvhxatsz1sNZavl6ygye/X82uuCSu6FCLe/o0oWKgtm0AXnrJbdlkzPmiiwI3tohkumLSA/gNONn/6ywKTEQkpx58ED7/HLZudasOAQhKVu44yGNTVzJ/037a1CzL+CHtaVu7fAAmC2zY4LacmjeHYcOgQQNo0SIwY4vICbI8lVMQ6FSOSBEwbx5Urw41a0JsrMvPCEAn4APxyfzvp3V8PHcz5SJCuLdPEwa2r0VQoE7bpKa6QKRpU5g2LTBjihRzOTqVY4y5M7NBrbXP53ZiIlJMHDgAPXvCoEEwYYLrc1OtWq6GTE3z8em8Lfzv53UcSkjhmtPr8t9ejSkbEYDKqqmp8PXXrsFeiRLw4YfQqFHuxxWRLGW2lZNxRq8J0BGYmn59ETDdy0mJSBHg87kk0Z49oVw594u+U6eADD0nei+PTV3Jmtg4Tq9fkUcvbh7Y3jZffglXXQU//QTnngtnnhm4sUUkU6cMTKy1jwMYY34C2llr49KvHwO+yJPZiUjhNW4c3HorLFwI7drBOefkesjtBxJ48vvVfLcshhrlwnl9cDvOb1ktML1tVq92jQLPPBMGDnS5L7165X5cEckWf+qY1AaSj7tOBup6MhsRKdz27nVfjRvD0KFQuTK0bZvrYROS03jjz795409Xmv72cxoxvEcDwkNyX6IecIXdhgyBtDTXnyc4GM4/PzBji0i2+BOYfAjMM8Z8hTuN0w9439NZiUjhYy2cdZYrkjZrlmu6d/nluRzSNdt7+vvV7DiYyIWto7i/bzNqlAvP/XyTklxRtBtugLAweP99V7FVla1F8pU/Bdb+zxjzA9A9/aGh1trF3k5LRAqNJUugTRv3C/2FF1xSawB+uS/fdpDR37rjv82jInnxyrZ0qhfAqq1z5ritpkqVXOVWHf8VKRD8WTHBWrsIWOTxXESksJk+HXr0cNVPr7oqIDkZuw4l8ty0tXy5aBsVIkJ4sl8rruhYKzBVW2fNgk2b3OmgHj2O5b+ISIHhV2AiInJUcjJER7u6Hmec4SqhBqD6aWJKGm/P3Mjrv28gOc3HTd3rc2vPhkSGBeD4b4ZnnoF169wKSXCwghKRAkiBiYhkz5AhbuVh/XqXmzFyZK6Gy+j++/QPa9i2P4Fzm1flgb7NqFepVO7nuncvPPUU3H2368HzxhsQGemCEhEpkBSYiEjW/v7bVW0ND3cN9665xgUlubR06wHGfLuKBZv307RaGT6+sTPdGlYKwITT7dkDr74K7du7raaoqMCNLSKeUGAiIpnbvNn1iHnkEdfjpnPnXA8ZczCBZ39cy1eLt1OpdAhP9W/F5R0ClEfyySduu+axx1zJ+y1b3GkbESkUFJiIyL/5fLBqlWuuV6cOPP889OuX62GPJKXyxp9/M2F6NBa45awGjDirAWUCmUcyZ47ry/Pgg1CypIISkUJGTfxE5N/uu89tgWzYkOueNgBpPssXC7byv5/XsTsuiYvaVOee3k2oVSEi93PdtAmGD4fnnoNWrSAhAUJDISgo92OLiCdy1MRPRIqZbdsgJMStMAwb5rZvsrnaEBMTw5VXXsnEiROplh7QzFi/m//7bjVrYuNoV7sc44e0p13t8rmfr8/ngo/ISBdAbdzoApPwABRfE5F8o8BERODwYWjdGi67DN58E+rXd1/ZNGbMGGbOnMno0aO5/ZFnePL71fy5bjc1y4fz6qC2XNAqKjB9bR57zNUg+eYbqFAB1q7VSRuRIkKBiUhxZS3MnQtdurgy8q+8Al275mio8PBwEhMTAQguVZ6JG0vw3fO/Y5MTeOSyjgw5vQ6hJXIZOMTHu9UQY6BiRXfCJjnZrfIoKBEpMrQJK1JcvfiiC0RWrHDXgwdDvXo5Gmr27NlUjqpJpbOupfqwNynd8hxqHFnPz3d048bu9XMflCxf7ub2ww/u+rbbYMIEF5SISJGiFROR4mTzZte8LqP7b8WKLpckF1LTfDz47o+E9n+K4NLlSVz3F/v+fJ8LrriEJnVr5nxgnw927ICaNd2x3/POC0girogUbApMRIqL1FRXQr5VK/j+eyhXzhVKyyFrLeVb9qBU10GUrNiKlK0r2TV5DMkx6wgODiY2NjZ3873ySreas3y5Wxn58MPcjScihYICE5GiLDXVJYheeimUKAHvvutWH3Jp4eb9PPX9aspddC/hKQeJ/fZZDqycTkREBAMHD2bs2LFHT+Vky+zZ0KGDqz9yww2wb19AOhWLSOGhHBORouyTT6B/f/jzT3fdqxfUqpXj4f7efZibP1zAZeNmsXlfPP/XryVnxs/k0OqZhIWFkZiYSGRkZM6CknnzXM7L+++76969XRl51SMRKVa0YiJS1CxdCnFxbtvmqqugUiXo0SNXQ+46lMiLv65n4vythJUI4s5zG3Nj93pEhJRg0jOxDB8+nGHDhjFhwgRiYmL8H3jtWteHp29f6NgR3nsPBg7M1VxFpHBT5VeRosRaV4+kdGm3LZJLhxJTGP/n37w9cyNpPsugTrW57ZxGVCodGoDJAn36uOBkwwYd+RUpRjKr/Ko1UpHC7tAhePZZV9PDGPj0U5fcmi4mJoYePXpkKxk1MSWNt2ZEc+azv/Pa73/Tu0U1fr3zLB6/pGXugpKdO+Guu1zuCMBrr7laKgpKRCSdAhORwm7WLLj3Xvj1V3fdsiWUP1by/fhqrFlJTfPx+fytnD32D574bjWtapTl29vO4KUr21K7YgD62uza5Qq5/fGHu27QQE32ROQE2soRKWyshUmTXLO6IUPc9Zo10KzZCS87vhrr8cLCwkhISPjHkJZpK3cy9qe1bNh1mDa1ynFv7yZ0bVgp9/N94glITHR/ggtOFIyIFGvayhEpSoxx/WzeessFJcb8KygBiI6OZtCgQUREuJWOiIgIBg8ezMaNG0943cz1e7j0tb8Y/tFCrLW8cXV7pozomrugJCXl2PebN7sGexn/CFJQIiKZ8CwwMca8Y4zZZYxZcdxjFYwxPxtj1qf/GYAWoyLFwIYNruDY3r3u+uOP4bffMq3xERUVRWRkJImJiSc9yrt4y34GvTmHq9+ey57DyTx7WWum3XEmfVpWIzY2Ntt5KUf9/jvUru3mDPDGG26+qkciIn7wcsXkPaDPPx67D/jVWtsI+DX9WkSykpgIP/0ES5a460qV/EoY3blzJ8OHD2fOnDkMHz6c2NhY1sbGMeyDBfR7fRZrY+N45MLm/DaqB5d3rEWJYPdXQnbyUgC3rZQRxDRr5o7+pqa6ayW2ikg2eJpjYoypC3xrrW2Zfr0WOMtaG2OMiQL+sNZmWYZSOSZS7FgLDz7ofrk/+6x7LD4eInKegLppzxFe/GUdXy/dQemQEgw7sz7Xn1GPUqHHyhllJy/lKJ/PJdw2bQqTJ+d4fiJSfBSkHJOq1toYgPQ/tdkscryMVQZj4MABd6w24x8POQxKdhxI4P7Jyznn+T+ZtnInw3s0YMa9Z3PbOY1OCErA/7wUkpKOBSFBQXD//XD77Tman4jI8Qps5VdjzDBgGEDt2rXzeTYieeCvv1yl1p9/dv1sXnstW3kZMTExXHnllUycOJFq1aqxKy6R13//m0/mbgFgSJc6jDi7AVXKhJ1yjKzyUo566y249VZYsADat3eng0REAiCvA5Odxpio47Zydp3qhdbaCcAEcFs5eTVBkTzl88HBg67uSKNGLj8jOdk9l81k0Yy8kIfGPEWDC4fz/qxNpKRZBravyW3nNKJGuXC/xsnISzmhxHxKigtGGjeGc86BoUNd8NSuXXZ/YhGRTOV1jslzwF5r7dPGmPuACtbae7IaRzkmUiRZCz17QpkyMHVqjofJyAsJCi1FmU79iGx/MSYkjMQ1M5j79qPUrVQq93NNSXGByHnnuVM2IiK5kFmOiWcrJsaYT4GzgErGmG3Ao8DTwOfGmBuALYC6dUnxs2QJnHaaWxEZPNjljmTUI8mB5avXM/SZD9kc1oCg0FIkrZ9NlzL7GDf+EarlJij5/HNXL+XHH6FkSVdhtmrVnI8nIuIHzwITa+1Vp3jqHK/uKVLgffml6577229w9tlw4405HupwUirvz9rEmzOiOVC2NYnr55Aw7wvid6wn6uab/50X4o/4eChRAkJCXKCUlAR79riAJCfjiYhkkyq/inhtxQqYN899f+GFrldMly45Hu5IUiqv/7GB7s/8xnPT1tK2VjmabP6agdX2MfPbiUfrlWTb1q1Qrx68+667HjAApk/XKomI5Cn1yhHxks/n6ntERcGff+ZqqPjkVD6cvZnx06PZdySZs5pU5o5ejTmtVrmcD3rwIKxcCV27uu2ku+6CK66Azp1zNVcRkczkS46JSLG1ebNLEB0zxm2LfPYZ1KmT4+Hik1P5aM5mxv8Zzd4jyXRvVIn/ntuYdrUD0NHhxhtdwLR1K4SGwvPP535MEZFcUGAiEmiLFrlf8P36QadOOT5SeyTJBSQTpruA5MzGlbn9nEa0r5OLgGTXLnjhBRg1CipWhEcfdceTQ0NzPqaISAApMBHJreRkuOceaNECbroJLr3UddOtXv1fRc/8cTgplQ9mb+KtGRvZF6iAJOPUz+7dMHasC5YGDnSl5EVEChAFJiI5lZLijtGWLAnLlkF4egEzY6B6deDEZnivv/56psPFJabw/qxNvDVzIwfiUzirSWVGntMod1s21sItt7gjyc8/74KnrVt1wkZECiwFJiI58cEH8PDDsHw5REa6zr8lTt0Mb9y4cYwbN+6kzfAOxqfwzl8befevjRxKTOWcplW47ZxGuUtq3bIFatd2QVLJkifMTUGJiBRkCkxE/LVnDwQHu/LxzZpBt26u7kdk5Im/+HHN8EaNGsWUKVOIj48nIiKCPn36sGPHDmJjY6lWrRr7jiTz1oxoPpi9mcNJqZzXvCq39WxEq5plczfPt9+Gm2+GNWugYUN3PFlEpJBQYCLij/37oUEDty3y9NPQsSN88skpX36yZnhr165l9erVPDDmGeqffxMfz91CYmoafVtGcWvPhjSLiszZ3KyF775zJ39atYILLoDHH4fKlXP4w4qI5B/VMRE5ld27YeZMd7oGXLffnj3daokf+vfvT1RUFG+99RbJyckER1ambOfLKN36PAgKJnHtTGZOeJhGVcvkbp5xcW7b5vLLYfz43I0lIpIHMqtjosBE5FRuvdVti+zY4bZvcmjOqo2MfP0bdobXASyJq//kjIoJvPbMYzkrGw/w8ccwbZrLdQHXf6dFC5dPIiJSwGUWmKgkvUiGvXvhzjth9Wp3/cADsHhxjoOSVTsO8Z9PFnHVh6vYU6ouh5f8wJ73bmXPDy9TtVRQ9oOSfftcJVlwqzl//+1WS8A1BVRQIiJFgAITkYxVQ5/PrZDMnOmuq1d35eSzacGmfVz/3nz6vjyDP9fuZniPBrSI/owrGwcz65fvctbLZuFCqFXL5ZIA3HYb/PUXlMnlNpCISAGjrRwp3u67D9avh0mT3PWhQ+6UTTp/C6RZa5m+fg+v/b6BeRv3UT6iJEO71ePa0+tSNiKHKxlLl7qk27POcjVT7r7bJd82aZKz8URECgj1yhE53saNULeuq/FRsaLbDklLc0eBI088GZNVgbQ0n+WHFTGM++NvVu44RFTZMB65sDlXdqpFREgu/u9lLVx3nduemTfP/fniizkfT0SkkNCKiRQvP/4IffvCb7+5lYhT+GeBtAwZBdISU9KYtGgbE6ZHs3lvPPUrl2L4mQ24tG0NQkrkcIf0t9/gmWdg6lTXu2blSredlIvEWxGRgkjJr1K8zZ59LG/krLPgscfcCZZMREdHM2jQICIiIgCIiIhg8ODBLFu9ntf/2ED3Z3/nwa9WUC68JG9c3Z5f/tuDyzvWyn5QcvDgsQRWn8+Vi9+82V23aKGgRESKHW3lSNHm88G117riYz//DGFh8MgjWb7tnwXSUkpEsL1yJ/q9u5LDSal0b1SJl644jdMbVMQYk7O57d7tKrPedx/cfz+ccw6sWAFB+veCiBRfCkyk6PnjD3j1VfjsM1cqfvJkl1OSTTt37mTwiFGUbNmH36IPsxG4qGkVbj6zPi1r5LBs/IIF7jjykCGuMuv990OfPu45Y9yXiEgxpsBEiobkZJcwGhoKBw64gmNbtkD9+tCyZbaGstYyf9N+yl98P7+u2UXYtkSGdK3HjWfUp3bFiOzPzdpjAcfLL8Ovv8KVV7qE1vvuy/54IiJFmJJfpfDbtQvat4d77nH1PXw+FwwEB2drmDSfZdrKWCZMj2bJ1gOUjyjJtV3rcs3pdalQKiRnc/vrL7jhBpd0W7cuxMRAqVL/Ov0jIlKc6LiwFD27d8Py5a53TZUq0L+/a2AH2c7RiE9O5cuF23hrxka27IunTsUIRl/SgoHtaxEekr3gBnDJqz4f1KvnethUruzqkdStC1FR2R9PRKQY0YqJFE6XX+6O127f7rZvcmDXoUTen72Jj+du4UB8CqfVKsfNZ9bnvBbVCA7KYa5HUpILPvr2hY8+ytkYIiJFnFZMpPD7+293zHfsWKhaFUaPhscfz1FQsjY2jrdmRPP1kh2k+Hyc17wqN3avT4c65XN2wuabb9yJn5dfdvN5/31o0yb744iIiAITKcCshYQEiIiA1FQXAAwe7E6xZLOHTUbJ+LdmRDNj/R7CSwZzZadaXN+tHnUrlcr+3HbvdlVjg4Jg1Sr46adj5ewvuij744mICKCtHCmofD7o2tWtPIwf7x5LSIDwcL/entHj5v2PP2XOjlTenrmR9bsOU6VMKNd2rcugTrUpn9OE1lmz4OyzXYXW3r3diaASJVR/RETET9rKkcLh4EH45Re47DL3S/6SS1xH3Qx+BiUADz3xLMtsbc57eQ7JQaE0i4rk+cvbcGHr6tmvzurzwZdfQunSLnekQwcYORIaNXLPh+QwwBERkX/RiokUHI8/7r42bXKnWXKgTO3mhLbqQ6nmZ0JQMAkb5nFo/teY3etJSEjI3mDJyS7osBZat4YGDWDKlBzNS0REjlGvHAmImJgYevToQWxsbKAGhKuucrU+AEaMcJVR/QxKMuazbUcMPyyP4fI3ZlNx0HNEtjiTxBW/sOPN4Rz58XkGnNmajRs3Zm9uzz8PTZpASoorjvb99zBpUjZ/QBERyS4FJuK3MWPGMHPmTEaPHp3zQdLSXEACLlF07lyIjnbXlStDu3Z+D/XImKdYmlSZc1+cyS0fL2LHwQQeuqAZPeN+Zc9P4yiRsI/ExEQiIyOpVq1a5oMlJrrTNPv3u+uWLeHCC+HIEXddq1a2C7aJiEj2aStHshQeHk5iYuK/Hg8LC8v+9sh557lf9hmrJGlp2f6FX7pGI0Jb9qZUi7MJCgkjcctyDi2cit26lIT4I/Tv35+oqCiGDRvGhAkTiImJYfLkyScfLKNc/JIl0LYtTJgAN92UvZ9JRESyJbOtHAUmkqWYmBhGjRrFV199RUJCAuHh4fTv35+xY8dmvRKxZYtbiXjgAReATJrkgoHLLstWw7o0n+XX1Tt5b9YmZv29lyBfKkdW/8n+uVMoeWQn/fr1828+GXw+GDDA5Y0895x7bN486NhRjfRERDymUzmSK1FRUURGRh5dHUlISMh8e8Ra94s/ONht1Tz+uKs90rGjC0iy4UB8MhPnb+XDOZvZtj+B6mXDuLt3ExZ98TLvff8aISEhJCYn+7dds2cPzJwJl17qTv3UrOmKtWXo1ClbcxMRkcDTiolkKVtbOfv2wbnnuu2Q4cNdYbSYmBOP/fphxfaDfDh7M1OWbCcp1UfnehW4rmtdzm1elRLBQf5v1xzf2ffuu+HFF918KlXK1nxERCRwtJUjuZKxlTNlyhTi4+OJiIg4ceskJgZWr3YN9ax11Vn793dbJdmQnOrjhxUxfDB7Mws37ye8ZDD92tXgmtPr0LSa/914M4qrTXrsMSrdcQe8+65Lqt2+3SW3tmyZzU9AREQCSVs5kisZWzmJiYmEhYX9+6TLiBEwezZs2+YqoH7ySbbG33EggU/mbuGz+VvYcziZuhUjeOiCZgxsX4uyESWzN9nFi3n3ySeZOXMmT330Ef8rXRoOH3bP1ajhvkREpMDSion45fitk59Gj6bvn3/SYs0atyWyahWULHmsEqoffD7LX3/v4cPZm/ll9U4scE7TKgw5vS7dG1YiKDvdfdO3a0qFhbE2KYmFwKXHPZ2j00MiIuIZrZhIrk1++WWXMFq9Om1Gj3bbNJs3u8CkeXO/xzkQn8yXC7fx8dwtbNxzhAqlQri5RwMGdapNrQoR2Z/YCy+4kz4zZrBh40bGXXcdb8+YAQkJJ2w5iYhI4aACa5K1w4ddFdSnnnLXLVq4VZL27f16+44dO+hywZWMeH8OnZ/8lSe+W035iJK8eMVpzL6/J/f2aXo0KMmyuuyRI/Dee64gGrjAqF49OHKEqKgodtevT2xS0sm3nEREpMBTYCIn9+GHMGqU+750aXj7bbjzzmPP+1Hr43BSKh/P3UyfF/4gttUQpq3ayYD2Nfl+ZHcmj+jGpW1rEFrixOJqp6wum5rq/pw1C4YOhR9/dNdDhri5li4NwM6dOxk+fDhz5sxh+PDhgSufLyIieUI5JnLMunUuT8QYuP9++P13mDHD5Y9kw4rtB/lk3hY+mrGWoJBwkndtJG7x9xxZ9Qc2OeGkOR+nOpJcPjSUfS1awMCBcN99rj7K3LnQpYsKoYmIFFJq4idZmzTJbdfMn++uR4+GOXP8DkqOJKXy2bwtXPLqTC58ZSaTFm7jwjY1aL1zGgcn3svhJT8QXsIwePDgkzbUi46OZtCgQURERNAVuCkkhMGDB7Nq0yZXmK1uXffCoCA4/XQFJSIiRZSSX4urI0fgiSfgzDPh/PNdUbTnnoOGDd3zfgYkK7Yf5NN5W/h6yQ4OJ6XSqEppHr2oOf3b1qRsRElumfPBqY8ZHycqLOzokeRbgoI4MzmZZ0qXdq99441A/uQiIlKAKTApTlJTYetWlywaFgYTJ7o/zz/fdfrNyCnJQlxiClOX7uCzeVtZvv0goSWCuKBVFIM616Z9nfKY41YzMnI+jq/Q+i8ffwxDh5LWsyfDhw+nbf/+vDBxIjG7dgXqJxcRkUJCOSbFSf/+7jTN6tVuKyQhAcLD/XqrtZZFWw4wcf4WvlkaQ0JKGk2rleGqTrW59LQaRwuhZVRdnThx4qlPw+zd6074XH6560+zaZNbFbnjDtAJGhGRIk91TIqrhQth7Fh45x0XgIwY4Y7+ZvSP8SMo2XckmcmLtvH5gq2s23mYiJBgLm5TnSs61aJtrXInrI7AiadqXn/99WNPHDgAu3ZB48YQGurKxNer5wKTunXh6acD+7OLiEihpBWTombrVoiIgIoV4Y8/4Mor4YcfoG1bv4dI81lmbtjD5/O38tOqWFLSLG1qleOqjrW4sE11Sof+O57NstFfy5auk++vv7on4uPdPEVEpNjRiklxsXOnW4V47DF46CHo0cMFKn4msi5Ys4kbnniTsGY92BmXQrmIkgzpUpcrOtaiSbUymb43Ojr6hEZ/w0JC+G/ZspRbssS9YOxYqFz52BsUlIiIyEkoMCns/u//3AmbJ590KxJvvAG9ernnjMkyKElITuPHlTF8Pn8bs6P3Ymt0ofL+7bx2TR96Na/yrwJopxIVHs6FW7fybYKrU3I4KYmUiAiqhYW5F/Tpk5ufUkREigkFJoVNYqIrenbuue560yYXmGS48cYsh8hIZP1y4Va+XRpDXFIqKftjOLLiVw4v/4UtcXu48H9+NL9LSXHzKVMGVq7kqhkzONK7Nx2feYYJEybwaEwMkytUyN3PKyIixYpyTAqLjITVJ5+EBx+EjRtd0mjG45nIOCnzytsfMnOba6IXvfsI4SWD6dsqip71wvn4xcf5On0b5vjmd6c8WXPkCDRoADffDI8/7uaxfDm0aqXiZyIikinlmBRma9e6BNaXXnLF0IYOdeXYa9d2z2cRBCQkpzH8qbdZU60XF7+5DIyhU90KDO/RgL6too4msn6bXtws00Jor7wCsbFu+6hUKbj1VujW7dg8WrcO9E8vIiLFjAKTgsbng2nT3FHes86CWrVc8bOUFPd8VJT7yoS1lvmb9nPhraMJa9yVoIi2lKywiwOzJ3Jk+a98k7ifz/+xRXPSQmj797t+Of37uxetXAmbNx9bpXnoIQ8+ABERKc60lZMPTlqE7OBBKFvW/dJv0gSaNoWpU7M17ua9R5i8aDuTF29j674EgnwpHFr5Jwmr/yB+41IiIsJP2KI56TwSE13CbHAwPPss3HsvREe70z5pae5xERGRXNBWTgHzryJk//0vTJkCf//tmtR9+y3UqePXWAfik/l2WQxfLd7Ows37MQYSNy8lbtkvxK+bhU1JOvra+Pj4E7Zo/jWPRYugZ0/44guXXHvddXDOOcca6CkoERERjykwyUMZRchaAe8DI8aNY9y4cVxUsiRTn3zSbdeEhrrqqJlISk3j9zW7+WrxNn5fs5vkNB+Nq5bm3j5NueS06piEdowatZIpG4OJT4Hg4GB69+5N5cqViY2NPTqPEOBZYMG4cZhx4ygbGsqBIUPcsWOAKlXcl4iISB5RYJJXli1j06xZ3Dl2LDsmTeL8pCTah4ZSY8AAxo4dm2WPGJ/PMn/TPqYs2cF3y3ZwKDGVSqVDubpLHfq3q0GL6pHHysOXCz/aqTcsLIzk5GTq1KnjVkXWrmXvwIGM/PZbpnz1FT0TEkgtUYLgK67wax4iIiJeUmDiJZ/Pbc3s2gVt21L1/vuJjIzks+Rk6oWGciQlhZtPdvrlOGtj45iyZDtTl+xg+4EEwksG07tFVfq1q0m3BhUpERx00vedkMz61FNM/PxzHnnkEardfTcVlywhsm9fEpOS6B4aSoIf8xAREckLSn71yqBBUKIEfPCBu54yBc44g/7DhhEVFXXC6ZfJkyef8NZt++OZunQHU5fsYE1sHMFBhu6NKtGvbQ3ObV6ViJBsxJMvvUTKXXdRzefjiuHDef3226FMGfrfemuW8xAREfFCZsmvCkwCZdYsd8z38cfd9ejRLjB54AG/3r7ncBLfL49h6pIdLNi8H4B2tctxadsa9G0VRaXSof7NY9UqV1/kxRcJ79yZRomJXAC8ARxIf0mWFV1FREQ8pFM5XrAWFixwXXtLlHCBySuvwMiRrrPvI49kOcShxBSmrYjlm2Ux/LVhD2k+S+Oqpbm7dxMual2d2hX9aHSXlAQffeSOGJ9xBlSq5LaOdu8+2ljv5ZNUdBURESmI8iUwMcZsAuKANCD1VFFTgZSRN/L993DhhfDjj9C7N9xyC9x2mztVk4n45FR+Xb2Lb5bu4I+17kRNrQrhDDuzPpecVp2m1SKznsO+fbB9uyv/HhwM998PV1zhApMqVWDFCgCi4IQk2FNWdBURESkg8nPF5Gxr7Z58vP/JC4ydyoEDrqbHDTfAiBHu+3feceXhwZVoP8X4H3z8KasPBPHtsh38unoXCSlpVCkTyuAutbm4TXVOq1Xu2ImaU0lMhIxOvRdd5K4XLnSrNYsXQ/XqJ33bSSu6ioiIFFD5kmOSvmLSwd/AxKsckxEjRjB+/Hhuvvlmd5T2n955B5KTYfhwt3UzaJArzz5wYJZjJ6WmMejOMUzfeJiyzc8g1ZSkQqkQzm9ZjYvaVKdj3QoEB/nZ7O6ZZ+CFF2DrVleVdfp0iIiADoVnoUlERCRDgUt+NcZsBPYDFhhvrZ2Q2esDHZhkFBj7pyqhoez89lvo1cs9cOmlrovuzz/7NW5yqo+ZG3Zz1T3PEVq/I0FhpUlLiCN+3Szi18yEXetIOHI464Hmz4f77nO5I1FR8Ntvbg733+/65oiIiBRiBTH5tZu1docxpgrwszFmjbV2+vEvMMYMA4YB1M7opBsgGUmhU6ZMgfh4CA+nX//+jC9TBvr2hZ07oXx5FxicZIvmeEmpacxcv4fvlsfw86qdxCWmUqVtL8L3rWP91Fc4sHYuEWEhXNavH2PHfn/yQQ4dgrffhrPPhtNOc/fcvh22bHGBSc+e7ktERKSIO3l1Lo9Za3ek/7kL+ArodJLXTLDWdrDWdqhcuXJA7x8VFUVkZCRnJiSwG2iUnhR65LrrGN6iBbEZqymlS7suuv+QmJLGTytj+e/EJXQY8ws3vL+AX1bt5Lzm1Xjnug4sfPg82qeu5tDa2YSFlPhX0mlMTAzXdOzInl9+cQMa41ZDpk1z182bw+rV0LlzQH9uERGRgi7PV0yMMaWAIGttXPr35wGj83oeO3fupPV11xEfH88lJUqwPDaWx95/nzeXLSNozJh/5ZzEJ6fy+5rd/LAiht/X7OJIchplw0vSu2U1LmgVRbeGlQgpEXTC+Cckne7YATt2QPXqjBk9mscXLCD2hhuotHkzlCkDmzcf61EDJw2IREREiro8zzExxtTHrZKAC4w+sdb+X2bv8brA2qlyTsIjK/DxH0v5cUUsf67bTWKKj4qlQjivRVXObxnF6Q0qUvIUJeGBY0eLAQYNYv1nn9E4/fPuBmwEdqCCZyIiUrwUuOTX7PI6MImJiTmac5JICOVadKd214uIi6hBqs9SNTKU3i2qcX7LKDrWLX/K/jQn+OADV/V13Tp3gmbaNA6sXcuts2fz1dSp/yp4ptoiIiJSXBTE5NcCJSWkLHsrtSGyXxsq1WiGMUEk+Y5w/Rn16NOyGqfVLEdQVkd7V6yAe+6Bl16CRo3YExnJ/NRU2m/YQJXWraF3b8r17k2Z1atV8ExEROQU8iX5tSBISE7jhZ/X0efF6Zz53O+sCW9BpahaXNWqLN2PzKDJhs94oG8z2tUuf/KgJC4OnngCZs5012XKuNWRbdsAeOSnn7hw924ee+ONE96WkXsyZ84chg8fTmxsrNc/qoiISKFRbLdy0nyW05/6lboVS3Fei6qc17xa5r1prIWffnKVVs85xxVeq1LFrZJkNOqzlvCIiJPmqyiPRERExNFWzkkEBxmm33M2YSWDT/2iXbtgwwbo2tWdkrn7bqhWzQUmISGuEmuZMsdeb8wJNVLUOE9ERCR7im1gAvw7KLHWBSKNGrnrW291WzXbt7vA5MsvoVatY68/PihJl1EjRXkkIiIi2Vdsc0yOOnzYHesFeOopaNYM9u931w8+CN99d+y1jRtDeHiWQyqPREREJGeKbY4J4HJGLroIZs2C9u1dtdW5c12TvixK0YuIiEjOKMfkVNq2hZEjXV8ccKslzZrl75xERESKseIdmFSuDM89l9+zEBERkXTKMREREZECQ4GJiIiIFBgKTERERKTAUGAiIiIiBYYCExERESkwFJiIiIhIgaHARERERAoMBSYiIiJSYCgwERERkQJDgYmIiIgUGApMREREpMBQYCIiIiIFhgITERERKTCMtTa/55AlY8xuYHMO3loJ2BPg6cgx+ny9pc/XW/p8vafP2FuF+fOtY62tfLInCkVgklPGmAXW2g75PY+iSp+vt/T5ekufr/f0GXurqH6+2soRERGRAkOBiYiIiBQYRT0wmZDfEyji9Pl6S5+vt/T5ek+fsbeK5OdbpHNMREREpHAp6ismIiIiUogUicDEGNPHGLPWGLPBGHPfKV5zljFmiTFmpTHmz7yeY2GW1edrjLk7/bNdYoxZYYxJM8ZUyI+5FkZ+fL5ljTHfGGOWpv/3OzQ/5llY+fH5ljfGfGWMWWaMmWeMaZkf8yysjDHvGGN2GWNWnOJ5Y4x5Of3zX2aMaZfXcyzM/Ph8mxpjZhtjkowxo/J6fl4o9Fs5xphgYB1wLrANmA9cZa1dddxrygGzgD7W2i3GmCrW2l35Md/Cxp/P9x+vvwj4r7W2Z97NsvDy87/fB4Cy1tp7jTGVgbVANWttcn7MuTDx8/N9DjhsrX3cGNMUeM1ae06+TLgQMsacCRwGPrDW/iuoM8b0BW4D+gKdgZestZ3zdpaFlx+fbxWgDnApsN9aOzZvZxh4RWHFpBOwwVobnf4X9WfAJf94zSBgsrV2C4CCkmzx5/M93lXAp3kys6LBn8/XAmWMMQYoDewDUvN2moWWP59vc+BXAGvtGqCuMaZq3k6z8LLWTsf9N3kql+B+qVpr7RygnDEmKm9mV/hl9flaa3dZa+cDKXk3K28VhcCkBrD1uOtt6Y8drzFQ3hjzhzFmoTHmmjybXeHnz+cLgDEmAugDTMqDeRUV/ny+rwLNgB3AcuB2a60vb6ZX6Pnz+S4F+gMYYzrh/vVZM09mVzz4/XeICBSNwMSc5LF/7k+VANoDFwC9gYeNMY29nlgR4c/nm+Ei4C9rbWb/epIT+fP59gaWANWB04BXjTGR3k6ryPDn830a9w+XJbgth8VoRSqQsvN3iAgl8nsCAbANqHXcdU3cvyz/+Zo91tojwBFjzHSgDW7vWTLnz+eb4Uq0jZNd/ny+Q4GnrUsI22CM2Qg0BeblzRQLtSw/X2vtIdxnTPp22cb0LwmM7PwdIlIkVkzmA42MMfWMMSG4X45T//Gar4HuxpgS6dsNnYHVeTzPwsqfzxdjTFmgB+6zFv/58/luAc4BSM99aAJE5+ksC68sP19jTLn05wBuBKanBysSGFOBa9JP53QBDlprY/J7UlJwFfoVE2ttqjHmVmAaEAy8Y61daYwZnv78G9ba1caYH4FlgA94y1p70qNXciJ/Pt/0l/YDfkpflRI/+fn5jgHeM8Ysxy2L32utLawdRfOUn59vM+ADY0wasAq4Id8mXAgZYz4FzgIqGWO2AY8CJeHo5/s97kTOBiCe9NUp8U9Wn68xphqwAIgEfMaYO4DmhTm4LvTHhUVERKToKApbOSIiIlJEKDARERGRAkOBiYiIiBQYCkxERESkwFBgIiIiIgWGAhMRyTVjzB3pNYK8vk+UMebbbL5nrDFGTSVFCgkFJiISCHcAngcmwJ3Am9l8zyvAfR7MRUQ8oDomIuI3Y0wp4HNcWfFgXPG3qsBYYC2u9cPZxphxQEcgHPjSWvto+vv7As8De4BFQH1r7YXp474CtMIVfnzMWvuvKsLGmGigmbU2yRhzHa7VezDQEvgfEAIMAZKAvhl9m4wxC4ELrLWxAf9QRCSgtGIiItnRB9hhrW1jrW0J/GitfRnX++Rsa+3Z6a970FrbAWgN9DDGtDbGhAHjgfOttWcAlY8b90HgN2ttR+Bs4Ln0YOUoY0w9YL+1Num4h1sCg4BOwP8B8dbatsBs4Pgu4ouAboH4AETEWwpMRCQ7lgO9jDHPGGO6W2sPnuJ1lxtjFuE69bYAmuMaD0ZbazMa5B3f8PE84L70Dr9/AGFA7X+MGQXs/sdjv1tr46y1u4GDwDfHzbPuca/bhevOLCIFXKHvlSMiecdau84Y0x7X++QpY8xP1trRx78mfWVjFNDRWrvfGPMeLtAwmQxtgMustWszeU1C+jjHO371xHfctY8T/34LS3+/iBRwWjEREb8ZY6rjtks+wuWVtEt/Kg4ok/59JHAEOJjeDfn89MfXAPWNMXXTr684buhpwG3GGJN+n7Ynuf06TlwFyY7GgBp3ihQCWjERkexohcv/8AEpwC3pj08AfjDGxKQnvy4GVgLRwF8A1toEY8wI4EdjzB5g3nHjjgFeBJalByebgAuPv7G19ogx5m9jTENr7QZ/J2yMKQk0xHVgFZECTqdyRCTPGGNKW2sPpwcfrwHrrbUvZOP9/YD21tqHsvmedtbah7M/YxHJa9rKEZG8dFN6gutKoCzulI7frLVf4VZTsqME7iixiBQCWjERERGRAkMrJiIiIlJgKDARERGRAkOBiYiIiBQYCkxERESkwFBgIiIiIgWGAhMREREpMP4fBqf8ZsZWYQQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Now we transform each piece back into the original form\n", + "Q_predict=np.exp(p_y)\n", + "Q_predict_upper=np.exp(p_y_upper)\n", + "Q_predict_lower=np.exp(p_y_lower)\n", + "x_topredict=np.exp(p_x)\n", + "# Plot the original data and then the prediction lines\n", + "plt.figure(figsize=(9,6))\n", + "\n", + "plt.xlabel('stage (m)')\n", + "plt.ylabel('discharge (cms)')\n", + "plt.plot(h_now,Qobs_now,'k*')\n", + "plt.plot(x_topredict,Q_predict)\n", + "plt.plot(x_topredict,Q_predict_lower,':r')\n", + "plt.plot(x_topredict,Q_predict_upper,':r')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Work on your own\n", + "Now, repeat the above but with different estimates of what b must be. Create one plot with multiple estimates of b and see how uncertain that component of the equation is in relation to parameters a and c." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Part 2: Brute force parameter estimation\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -1264,7 +1592,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.4" + "version": "3.8.12" } }, "nbformat": 4, diff --git a/modules/module7.md b/modules/module7.md index e0aefec4..78672e77 100644 --- a/modules/module7.md +++ b/modules/module7.md @@ -18,9 +18,25 @@ Download the lab and data files to your computer. Then, upload them to your Jupy ## Homework 7: -### Problem 1: Application of Bayes Theorem with MCMC +### Problem 1: ENSO Phases +Following Lab 7-1 and Lab 7-2, +A) Use the time series of the phase of the El Niño Southern Oscillation (ENSO) from 1900-2021 to create a lag-1 Markov model of the ENSO phase. +where the observed Phases of ENSO are as follows: + +1: warm (El Niño) +2: neutral (ENSO neutral) +3: cool, (La Niña) + +B) Using this Markov model and a random number generator, simulate 5,000 years of ENSO data. + +C) Using this randomly generated data, answer the following questions. + + - According to the model, what is the probability that three warm ENSO years would occur in a row? + - What is the large-sample probability that three cool ENSO years would happen in a row? (Try refreshing the numbers several times to increase the sample size if the condition never happens.) + +### Probelm 2: Rating Curves and Application of Bayes Theorem with MCMC -Following the Lab 7-3, explore how the rating curve and the 95% confidence intervals for the Lyell Fork streamflow site change depending on the method you use: +Following the class discussion and Lab 7-3, explore how the rating curve and the 95% confidence intervals for the Lyell Fork streamflow site change depending on the method you use to determine the rating curve: - Least squares linear regression fitting (with transformed variables) using h0 = 28 cm - Make 95% confidence intervals around this regression fit @@ -29,15 +45,14 @@ Following the Lab 7-3, explore how the rating curve and the 95% confidence inter - Direct monte carlo parameter estimation - Bayesian MCMC fitting -Create plots and discuss the differences in the results from these three methods. - -### Problem 2: Air Temperature Observations in Complex Terrain +Using the code in Lab 7-3, create plots and discuss the differences in the results from these three methods. -[See problem 2 in Module 8](/data-analysis/modules/module8.html) +### Problem 2 grads: Work on your term projects (CEWA 565) -### Problem 3: Statistics Synthesis (CEE 465) +### Problem 2 undergrads: Statistics Synthesis (CEE 465) +(Your final exam questions will look similar to this.) You are given the below dataset of annual peak flows on the Sauk River: ![Sauk River Plot](lab7/sauk-river-plot.png)