From 2214ce48d8264f835dd91bfcd7048a7b746c2b8d Mon Sep 17 00:00:00 2001 From: shuanghuab Date: Thu, 12 Jan 2023 23:21:42 +0000 Subject: [PATCH] code for dose response analysis blog --- local_optimization/NLDF/dose_resp.dat | 42 +++ local_optimization/NLDF/dose_resp.ipynb | 371 ++++++++++++++++++++++++ 2 files changed, 413 insertions(+) create mode 100644 local_optimization/NLDF/dose_resp.dat create mode 100644 local_optimization/NLDF/dose_resp.ipynb diff --git a/local_optimization/NLDF/dose_resp.dat b/local_optimization/NLDF/dose_resp.dat new file mode 100644 index 0000000..44658fc --- /dev/null +++ b/local_optimization/NLDF/dose_resp.dat @@ -0,0 +1,42 @@ +0.0005874,7.327 +0.0005874,-6.043 +0.0005874,-0.428 +0.002937,3.110 +0.002937,-2.340 +0.002937,3.103 +0.01468,-11.981 +0.01468,-4.885 +0.01468,-0.566 +0.03283,-4.475 +0.03283,-5.747 +0.03283,-1.660 +0.07341,-4.922 +0.07341,-4.183 +0.07341,-0.981 +0.1642,6.090 +0.1642,-10.287 +0.1642,-1.919 +0.3670,-11.767 +0.3670,-5.345 +0.3670,-7.134 +0.8207,-3.393 +0.8207,-2.099 +0.8207,-4.490 +1.835,-6.812 +1.835,-5.808 +1.835,-9.666 +4.103,-0.600 +4.103,-5.707 +4.103,-4.203 +9.175,4.313 +9.175,-8.110 +9.175,0.388 +20.52,-10.838 +20.52,-22.657 +20.52,3.820 +45.87,-34.218 +45.87,-30.895 +45.87,-23.929 +91.74,-30.911 +91.74,-44.233 +91.74,-36.776 \ No newline at end of file diff --git a/local_optimization/NLDF/dose_resp.ipynb b/local_optimization/NLDF/dose_resp.ipynb new file mode 100644 index 0000000..b57f9bf --- /dev/null +++ b/local_optimization/NLDF/dose_resp.ipynb @@ -0,0 +1,371 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bf21964c", + "metadata": {}, + "source": [ + "# Modelling dose–response relationships using data fitting\n", + "\n", + "The following is the Python code used to fit a nonlinear regression model to dose–response data of a chemical taken from the US National Toxicology Program library." + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "793e1544", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import warnings\n", + "from naginterfaces.library import opt\n", + "from naginterfaces.base import utils" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "id": "bacb5e9b", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAixUlEQVR4nO3de7hcVX3/8ffH3ICEmwSCQC6QBCkBRDxBUPsTlRZQKtoikmNIqlVEuVnxWChaTRAL0gdUrLSpIFK5JSkqArUgreANyQn3RIKRW4KA4RauhoR8f3+sNWQ4cy5zkjOz55z9eT3PPLP32rfvzJkz31lr7b22IgIzM7Nqrys6ADMzaz1ODmZmVsPJwczMajg5mJlZDScHMzOr4eRgZmY1nBzMSkrSzyR9vOg4rDU5OVhTSXpQ0kuSnpP0jKRfSTpOkj+LLSz/3Q4uOg5rHv9DWhH+KiK2BCYCZwH/AFxYbEiJpGFFx2DWCpwcrDARsToirgY+DMyWtBeApK0lXSJplaSHJH2hUrOQNEXSTZJWS3pC0pWV/UnaQ9INkp6StEzSUX3FIOliSRdIuk7SC8C7JO0k6b/y8R+QdFLV+vtL6pT0rKTHJZ2byydJCknHSvqDpEclfa5qu1GSvp6X/SFPj8rLDpK0UtIpkv6Yt/1o1bbvlbQ017Ye6bLfwyXdUVUL26eX1/oXku7N7923AFUtmyzpfyU9md/XSyVtk5f9JzAB+LGk5yV9PpcvkPRY3t/Nkqb19X7bIBIRfvjRtAfwIHBwN+UPA5/K05cAPwK2BCYB9wF/l5ddDpxO+mGzGfCOXD4aWAF8FBgOvBl4Atizj3guBlYDb8/73AJYDPwTMBLYDbgfOCSv/2vgmDw9BjggT08CIsc3GtgbWFV5rcBc4BZgB2B74FfAGXnZQcC6vM4I4L3Ai8C2efmjwJ/n6W2B/fL0m4E/Am8FhgGz8/s7qpvXORZ4DjgyH+Pv8zE/npdPAf4CGJXjuxn4em9/N+Bj+W80Cvg6cEfRny8/Bu5ReAB+lOvRS3K4JX/pDwNerv5SBz4J/CxPXwLMA3bpsv2HgZ93Kft34Et9xHMxcEnV/FuBh7uscxrw3Tx9MzAHGNtlnUpy2KOq7GvAhXn698B7q5YdAjyYpw8CXgKGVy3/Y1XieTi/B1t1OeYFlQRTVbYMeGc3r3MWcEvVvICVleTQzfofAG7v6+9WtXyb/Pq3Lvoz5sfAPNysZK1iZ+Ap0i/cEcBDVcseyssBPk/6YrtV0hJJH8vlE4G35uaVZyQ9A3wE2FHS7rk5pPrxn1X7X1E1PRHYqct+/hEYl5f/HbA7cK+kRZIO7/I6qvf1ELBTnt6pm9e0U9X8kxGxrmr+RVLNBOBvSLWJh3KT2oFVsZ7SJdbxXfZbsVN1bJG+0V+dlzRO0hW52epZ4Pukv0W3JA2TdJak3+f1H8yLetzGBpfhRQdgJmk66cv/F6SmoLWkL76leZUJwCMAEfEY8Im83TuAn0q6mfRFd1NE/EUPhxnTQzmkX7wVK4AHImJqtytG/A6YkftA/hpYKGm7qlXGA/dWxf2HPP2H/JqWdLOsVxGxCDhC0gjgBGB+Ps4K4MyIOLOO3TyatwFAkqrnga+S3oe9I+IpSR8AvlUdRpf9tQNHAAeTEsPWwNNU9WPY4OaagxVG0lb5l/cVwPcj4u6IeIX05XempC0lTQQ+S/oli6QPSdol7+Jp0pfWeuAaYHdJx0gakR/TJf1ZP8O6FXhO0j9I2jz/Qt4rJzAkzZS0fUSsB57J26yv2v6LkrbInbMfBSod5pcDX5C0vaSxpD6N79fxHo2U9BFJW0fEWuDZquP9B3CcpLcqGS3pfZK27GZX1wLTJP21pOHAScCOVcu3BJ4HVkvaGejosv3jpP6X6vXXAE+S+mm+2tdrscHFycGK8GNJz5F++Z4OnEv6Iq04EXiB1BH8C+Ay4KK8bDrwG0nPA1cDJ0fE/RHxHPCXwNGkX+SPAWeTOkvrlpPT4cC+wAOkmsx3SL+MAQ4FluTjfwM4OiJeqtrFTcBy4EbgXyLi+lz+FaATuAu4G7gtl9XjGODB3HxzHKm5jIjoJNWivkVKlMuBv+3hdT0BfIh06vCTwFTgl1WrzAH2I3XOXwtc1WUX/0xKbs/ks6UuITWNPUKq4d1S52uxQUKp6dHMNoWkSaRkMqJL34HZoOSag5mZ1XByMDOzGm5WMjOzGq45mJlZjSFxncPYsWNj0qRJRYdhZjaoLF68+ImI2L67ZUMiOUyaNInOzs6iwzAzG1QkPdTTMjcrmZlZDScHMzOr4eRgZmY1nBzMzKyGk4OZmdUodXJYvRqmTUvPZma2QamTw7XXwtKlcN11RUdiZtZaSpkc2tthzBiYPTvNz5qV5tvbi43LzKxVlDI5zJ0LEybAiBFpfsQImDgRzjij2LjMzFpFKZPDlCkpQaxdC6NHp+c5c2Dy5KIjMzNrDaVMDgDz56fEcOqpsH49XHpp0RGZmbWO0iaHjg5Ytgx22y0lh7a2oiMyM2sdpU0O552XmpEqndJf/rI7pc3MKkqbHNwpbWbWs9ImB3dKm5n1rLTJATZ0Ss+Zk54XLCg6IjOz1jAkbvazsTo64PzzYdw4mDkTVqwoOiIzs9ZQ6prD9OkpMUB6buYZSx7XycxaWamTQ5E8rpOZtTInhyZr5XGdXJsxswonhyZr5VNoXZsxs4pSJ4eHH4ZRo9Jzs7TiKbStXJsxs2KUOjmcfTa8/DKcc05zj9tqp9C2cm3GzIqhiCg6hk3W1tYWnZ2dda8/aRI89FBt+cSJ8OCDAxZWjxYtSl/G48bB44+nU2iLHttp4UKYMSPVpNasgcsvhyOPLDYmM2ssSYsjottvn1LWHC68EEaOfG3ZyJFw0UXNOX6Rp9D2pNVqM2ZWrFImh/e8B0444bVlJ5wA7353MfG0gsootaeckp47OoqOyMyKVMrkAOmXMsDhh792vqxasTZjZsUp7fAZX/kKvOUtsNdecM89cNttRUdkZtY6SpscKqdtQkoQe+1VXCxmZq2mtM1KZmbWMycHMzOr4eRgZmY1nBwK4kHuzKyVtWxykHSopGWSlks6teh4BpoHuTOzVtaSyUHSMOBfgcOAPYEZkvYsNqqB4UHuzGwwaMnkAOwPLI+I+yPiZeAK4IiBPkgRTTse5M7MBoNWTQ47A9V3dF6Zy14l6VhJnZI6V61atVEHWbAgNe0sXNi/7TYlqbTikN1mZl21anLoU0TMi4i2iGjbfvvt+7VtpWnnk59M8x//ePqirrdpZ1P7CzzInZm1ulZNDo8A46vmd8llA+K55+DFF2H9+g1lL76YynszUP0FHuTOzFpdqyaHRcBUSbtKGgkcDVw9UDs/7zzYbLPa8htv7P2LfqD6CzzInZm1upZMDhGxDjgB+B/gt8D8iFgyUPufMgXe+MbXlkmw6669f9G7v8DMyqIlkwNARFwXEbtHxOSIOHOg97/ddq+tPUTAs8/C2LG9b+f+AjMrg5ZNDo02bFi6f3S1lSvhsMN63879BWZWBqVNDiNHvrZDuuKWW3rvZHZ/gZmVQSmTQ3t76nzuzqhRvijNzKyUyWHuXJg0KXVCV3vd62DdOncym5mVMjlUzjqKeG358OGwxRbuZDYzK2VygA1nHX3pS+l5771TB/XZZ7uT2cystMmhowN+/3u4775Ug1i6NJWfeCIcdJBHSTWzcittcqicddTRkfoZRo5M5fVc9ewb9ZjZUFfa5FCxeHFqTlqzpv6rnn2jHjMb6kqbHCqD6B13XJpfvx5eeCFN99Qh7Rv1mFlZlDY5VEZmfeWV15bvsUfPHdK+UY+ZlUVpk8N558E229SW33tvzx3SHnjPzMqilMmhvT3VEJ5+unZZRO+1AQ+8Z2ZlUMrkULlCujtS77WBjg5YtAguuig9+5oIMxuKSpkcpkyBs86qLZfSaa2nnNLz6arTp6eksHQpdHZ64D0zG5pKmRwAPvvZ2rKIdPbRN77R/emqPlvJzMqitMnh619PN/zp6o1vhJkzu08APlvJzMqilMmhvT196XftkB42DHbYoecE4LOVzKwsSpkcKjWAimHD0nDd73xnWjZ3brpiGtLV09UJwGcrmVkZKLqOWz0ItbW1RWdnZ7+2WbgQjj463dxnzZpUY7jsMpg3LyWA6ovjhg2Do45KyxctSoll3Dh4/HFYscKd0mY2OElaHBHdfoOVsuYAKQGMGZNqCSNHwqOPpsfcuTB+PGy2WVpvs81SMqj0K/g2oWZWBsOLDqAoHR2pz+BLX0rPkPohRo2CffeFlStTs9GaNfC1r7lfwczKpbQ1h+nT4Zxzuu983mqrxvcreNhvM2tlpa05wIazj2bM2FBLmDMnJYhKv8LMmalfYaBVX0cxY8bA79/MbFOUtuZQ0d3ZR43sV/CFdGY2GJT2bKWKZp99tHw5vP/98OCD8NJLsPnmsOuucPXV7tcws+by2Uq9aPbZR76QzswGg9InhyL4Qjoza3V9JgdJ3YxANPQ18myijg5YtiyN/rpsmYf9NrPWU0/N4RZJCyS9V5IaHlGL6G5U1oHiC+nMrNXVkxx2B+YBxwC/k/RVSbtvykElfUjSEknrJbV1WXaapOWSlkk6ZFOOszF8NpGZWR3JIZIbImIG8AlgNnCrpJskHbiRx70H+Gvg5upCSXsCRwPTgEOBb0satpHHqEvX5iMPy21mVmefg6STJXUCnwNOBMYCpwCXbcxBI+K3EbGsm0VHAFdExJqIeABYDuy/MceoV9fmI59NZGZWX7PSr4GtgA9ExPsi4qqIWBcRncC/DXA8OwPV1yOvzGU1JB0rqVNS56pVq/p9oN6aj8p6NpGH9DCzinqGz/hCRMyvLpD0oYhYEBFn97SRpJ8CO3az6PSI+FE/46wREfNIfSG0tbX1+0q+uXPhjjvSxWjr1r22+eipp+D88xs7fEYr8pAeZlZRT83h1G7KTutro4g4OCL26ubRW2J4BBhfNb9LLhtwvTUfle1sInfCm1lXPSYHSYdJOh/YWdI3qx4XA+saFM/VwNGSRknaFZgK3NqgY73afHTqqbB+PVx6aaOO1NrcCW9mXfVWc/gD0An8CVhc9bga2KRTTCV9UNJK4EDgWkn/AxARS4D5wFLgJ8DxEfFKz3vaNJ/6VKoZ7LhjSg5DvYbQE3fCm1lXfQ68J2l4RDSqpjAgNnbgvbe9DX7963Qb0FdegeHD081+3v/+dEvQMjnqKLj+evjiF1ON4ZBD4Mori47KzBqpt4H3euyQljQ/Io4CbpdUk0EiYp8BjLGp2tvTKKh/+lOar9wvOqK8zSkdHeXshDez7vV2ttLJ+fnwZgTSTNVnKr30UiqT0qOszSnTp2+YHjduQ4e8mZVTj30OEfFonvwbYG1EPFT9aE54jVHdxj4sX399zDHluqbBzKw39VznsCVwg6SngCuBBRHxeGPDarzKmUqzZ8N3v5uamJYtc3OKmRn0405wkvYBPkyqSayMiIMbGVh/bEyHdLPvAGdm1mo2qkO6G38EHgOeBHYYiMCK5DZ2M7Oe1TPw3qcl/Qy4EdgO+MRgPlPJzMz6Vk/NYTzwmYi4o8GxmJlZi+jtOoetIuJZ4Jw8//rq5RHxVINjMzOzgvRWc7iMdI3DYiCA6luEBrBbA+MyM7MC9ZgcIuLw/Lxr88IxM7NWUE+H9I31lJmZ2dDRW5/DZsAWwFhJ27KhWWkrerg7m5mZDQ299Tl8EvgMsBOp36GSHJ4FvtXYsMzMrEi99Tl8A/iGpBMj4vwmxmRmZgWr5zah6yVtU5mRtK2kTzcuJDMzK1o9yeETEfFMZSYingY+0bCIzMyscPUkh2GSXr3GQdIwYGTjQjIzs6LVM3zGT4ArJf17nv9kLjMzsyGqnuTwD6SE8Kk8fwPwnYZFZGZmheszOUTEeuCC/DAzsxLoMzlImgr8M7AnsFmlPCI8tpKZ2RBVT4f0d0m1hnXAu4BLgO83MigzMytWPclh84i4kXRL0Yci4svA+xoblpmZFameDuk1kl4H/E7SCcAjwJjGhmVmZkWqp+ZwMmkAvpOAtwAzgdmNDMrMzIpVz9lKi/Lk88BHGxuOmZm1gnpqDmZmVjJODmZmVqOeO8G9vZ6y/pB0jqR7Jd0l6QddRn09TdJyScskHbIpxzEzs41TT82hu3s5bOr9HW4A9oqIfYD7gNMAJO0JHA1MAw4Fvp0H+jMzsybq7TahBwJvA7aX9NmqRVsBm/SFHRHXV83eAhyZp48AroiINcADkpYD+wO/3pTjmZlZ//RWcxhJup5hOLBl1eNZNnyZD4SPAf+dp3cGVlQtW0kP96uWdKykTkmdq1atGsBwzMyst9uE3gTcJOniiHiovzuW9FNgx24WnR4RP8rrnE4aluPS/u4/IuYB8wDa2tqiv9ubmVnP6rlCepSkecCk6vUj4t29bRQRB/e2XNLfAocD74mIypf7I8D4qtV2yWVmZtZE9SSHBcC/ke7h8MpAHFTSocDngXdGxItVi64GLpN0LrATMBW4dSCOaWZm9asnOayLiIG+l8O3gFHADfkOpLdExHERsUTSfGApqbnp+IgYkIRkZmb1qyc5/FjSp4EfAGsqhRHx1MYeNCKm9LLsTODMjd23mZltunqSQ2WQvY6qsgB8sx8zsyGqnoH3dm1GIFa81avhbW+DX/0Ktt666GjMrEj1DJ+xhaQv5DOWkDRV0uGND82a7dprYelSuO66oiMxs6LVe5vQl0lXS0M6tfQrDYvImq69HcaMgdm5AXHWrDTf3l5sXGZWnHqSw+SI+BqwFiCfeqqGRmVNNXcuTJgAI0ak+REjYOJEOOOMYuMys+LUkxxelrQ5qRMaSZOpOmvJBr8pU1KCWLsWRo9Oz3PmwOTJRUdmZkWpJzl8CfgJMF7SpcCNpAvYbAiZPz8lhjlz0vOCBUVHNDSsXg3TpqVns8GknrOVbpB0G3AAqTnp5Ih4ouGRWVN1dMD558O4cTBzJqxY0fc21rfqTv4ZM4qOxqx+9Zyt9EHSVdLXRsQ1wDpJH2h4ZNZU06enxADpua2t2HgGO3fy22BXV7NSRLxaKY6IZ0hNTWbWA3fy22BXT3Lobp16rqw2Ky138ttgV09y6JR0rqTJ+XEusLjRgZkNdu7kt8GsnhrAicAXgStJp7PeABzfyKDMhgJ38ttg1mtykDQMuCYi3tWkeMyGjOnTN0yPG7ehw99sMOi1WSnfS2G9JA/DZmZWIvU0Kz0P3C3pBuCFSmFEnNSwqMzMrFD1JIer8sPMzEqiniukv5fHVpoQEcuaEJOZmRWsniuk/wq4gzS+EpL2lXR1g+MyM7MC1XOdw5eB/YFnACLiDnyLUDOzIa2e5LC2eviMbH0jgjEzs9ZQT4f0EkntwDBJU4GTgF81NiwzMytSPTWHE4FppBv8XAasBj7TwJjMzKxgPdYcJG0GHAdMAe4GDoyIdc0KzMzMitNbzeF7QBspMRwG/EtTIjIzs8L11uewZ0TsDSDpQuDW5oRkZmZF663msLYy4eYkM7Ny6a3m8CZJz+ZpAZvneQEREVs1PDozMytEj8khIoY1MxAzM2sd9ZzKOuAknSHpLkl3SLpe0k65XJK+KWl5Xr5fEfGZmZVdIckBOCci9omIfYFrgH/K5YcBU/PjWOCCYsIzMyu3QpJDRDxbNTuadPtRgCOASyK5BdhG0huaHqCZWcnVM3xGQ0g6E5hFuuK6chvSnYHqO+2uzGWPdrP9saTaBRMmTGhorGZmZdOwmoOkn0q6p5vHEQARcXpEjAcuBU7o7/4jYl5EtEVE2/bbbz/Q4ZsNiNWrYdq09Gw20Br5+WpYcoiIgyNir24eP+qy6qXA3+TpR4DxVct2yWVmg9K118LSpXDddUVHYkNRIz9fRZ2tNLVq9gjg3jx9NTArn7V0ALA6ImqalMxaXXs7jBkDs2en+Vmz0nx7e7Fx2dDQjM9XUWcrnZWbmO4C/hI4OZdfB9wPLAf+A/h0QfGZbZK5c2HCBBgxIs2PGAETJ8IZZxQblw0Nzfh8KSL6XqvFtbW1RWdnZ9FhmL3GwoUwYwaMGgVr1sDll8ORRxYdlQ0VA/H5krQ4Itq6W1ZUzcFsyJs/H0aPhjlz0vOCBUVHZENJoz9frjmYNciiRanqP24cPP44rFgBbd3+RjPrv4H4fPVWcyjsOgezoW769A3T48alh9lAafTny81KZmZWw8nBzMxqODmYmVkNJwczM6vh5GBmZjWcHMwayAPv2WDl5GDWQB54zwYrJwezBvDAezbYOTmYNYAH3rPBzsnBrAGmTEkJYu3aNO7N2rVpDJzJk4uOzKw+Tg5mDeKB92ww89hKZg3S0QHnn5/GvJk5Mw2MZjZYODmYNYgH3rPBzM1KZmZWw8nBzGyQauRFlk4OZmaDVCMvsnRyMDMbZJpxkaWTg5nZINOMiyydHMzMBplmXGTp5GBmNgg1+iJLX+dgZjYINfoiSycHM7NBqNEXWbpZyczMajg5mJlZDScHMzOr4eRgZmY1nBzMzKxGoclB0imSQtLYPC9J35S0XNJdkvYrMj4zs7IqLDlIGg/8JfBwVfFhwNT8OBa4oIDQzMxKr8iaw3nA54GoKjsCuCSSW4BtJL2hkOjMzEqskOQg6QjgkYi4s8uinYHq6/xW5rLu9nGspE5JnatWrWpQpGZm5dSwK6Ql/RTYsZtFpwP/SGpS2mgRMQ+YB9DW1hZ9rG5mZv3QsOQQEQd3Vy5pb2BX4E5JALsAt0naH3gEGF+1+i65zMzMmqjpzUoRcXdE7BARkyJiEqnpaL+IeAy4GpiVz1o6AFgdEY82O0Yzs7JrtYH3rgPeCywHXgQ+Wmw4ZmblVHhyyLWHynQAxxcXjZmZga+QNjOzbjg5mJlZDScHMzOr4eRgZmY1nBzMzKyGk4OZmdVwcjAzsxpODmZmg9Tq1TBtWnoeaE4OZmaD1LXXwtKlcN11A79vJwczs0GmvR3GjIHZs9P8rFlpvr194I7h5GBmNsjMnQsTJsCIEWl+xAiYOBHOOGPgjuHkYGY2yEyZkhLE2rUwenR6njMHJk8euGM4OZiZDULz56fEMGdOel6wYGD3X/iorGZm1n8dHXD++TBuHMycCStW9L1Nfzg5mJkNQtOnb5geNy49BpKblczMrIaTg5mZ1XByMDOzGk4OZmZWw8nBzMxqKCKKjmGTSVoFPLSRm48FnhjAcAYzvxeJ34fE78MGQ/W9mBgR23e3YEgkh00hqTMi2oqOoxX4vUj8PiR+HzYo43vhZiUzM6vh5GBmZjWcHGBe0QG0EL8Xid+HxO/DBqV7L0rf52BmZrVcczAzsxpODmZmVqPUyUHSoZKWSVou6dSi42kWSeMl/Z+kpZKWSDo5l79e0g2Sfpefty061maQNEzS7ZKuyfO7SvpN/lxcKWlk0TE2g6RtJC2UdK+k30o6sIyfCUl/n/8v7pF0uaTNyviZKG1ykDQM+FfgMGBPYIakPYuNqmnWAadExJ7AAcDx+bWfCtwYEVOBG/N8GZwM/LZq/mzgvIiYAjwN/F0hUTXfN4CfRMQewJtI70mpPhOSdgZOAtoiYi9gGHA0JfxMlDY5APsDyyPi/oh4GbgCOKLgmJoiIh6NiNvy9HOkL4GdSa//e3m17wEfKCTAJpK0C/A+4Dt5XsC7gYV5lbK8D1sD/w+4ECAiXo6IZyjhZ4J0n5vNJQ0HtgAepYSfiTInh52B6nsnrcxlpSJpEvBm4DfAuIh4NC96DBjg24e0pK8DnwfW5/ntgGciYl2eL8vnYldgFfDd3MT2HUmjKdlnIiIeAf4FeJiUFFYDiynhZ6LMyaH0JI0B/gv4TEQ8W70s0jnOQ/o8Z0mHA3+MiMVFx9IChgP7ARdExJuBF+jShFSSz8S2pNrSrsBOwGjg0EKDKkiZk8MjwPiq+V1yWSlIGkFKDJdGxFW5+HFJb8jL3wD8saj4muTtwPslPUhqVnw3qd19m9ykAOX5XKwEVkbEb/L8QlKyKNtn4mDggYhYFRFrgatIn5PSfSbKnBwWAVPzWQgjSZ1OVxccU1PkdvULgd9GxLlVi64GZufp2cCPmh1bM0XEaRGxS0RMIv39/zciPgL8H3BkXm3Ivw8AEfEYsELSG3PRe4CllOwzQWpOOkDSFvn/pPI+lO4zUeorpCW9l9TmPAy4KCLOLDai5pD0DuDnwN1saGv/R1K/w3xgAmkI9KMi4qlCgmwySQcBn4uIwyXtRqpJvB64HZgZEWsKDK8pJO1L6pgfCdwPfJT0A7JUnwlJc4APk87qux34OKmPoVSfiVInBzMz616Zm5XMzKwHTg5mZlbDycHMzGo4OZiZWQ0nBzMzq+HkYGZmNZwczMyshpODbRRJO0q6QtLvJS2WdJ2k3YuOq78kfVnS54qOoy/5XgufHojtJP1qgGLaXNJNefj7vtb9pKTIFxtWyo7PZYdIurlqeAprAU4O1m95WIEfAD+LiMkR8RbgNIb4iJ0F2wboNjko6el/uWa7iHjbAMX0MeCqiHiljnX3Bu4E9gCQtAXpyuNVwB2ke0V8eIDisgHg5GAb413A2oj4t0pBRNwZET+X9Nl8B617JH2mslzSpHyHsYsl3SfpUkkHS/plvsvY/lXrXJrvRLYwf4lU9lGz77zNPVXrfC7XBiblffxHvqvX9ZI2z+ucnmP4BVAZS6iGpFmS7pJ0p6T/rCOOno5Xsx9JMyXdKukOSf9e+fXdy37OAibn9c/J6y2TdAlwDzBe0g9zLW6JpGNzuK/ZLh/j+Y19LV18hKoxhiT9TFLly3+76r8LsA9p+Ik98vxJwAJgfUQ8Dvww789aRUT44Ue/HqR/7PO6KX8Labym0cAYYAnw5rxsEmmsmr1JP0oWAxcBIg2R/MO8TgBvz9tcRBrvqMd9523uqYrhc8CXq463by6fD8ys2s8WwFbA8soxuryWacB9wNg8//o64ujueDX7Af4M+DEwIpd9G5jV5X3qup+ur3MSaVysA6rKKjFuTkoY23XdLi9/fmNeS5d9jAQe61K2Enhdnn4XcHnVslXAbsBPSLWZ24GDgJ/m5cOAVUV/tv3Y8HDNwQbSO4AfRMQLEfE8abjjP69a/kBE3B0R60lfRDdG+ma4m/SFBLAiIn6Zp7+f91nPvrvzQETckacX52P8ed7Pi5HuYdHTSLzvBhZExBMAsWGwud7i6O543e3nPaQv5kWS7sjzu/URd3ceiohbquZPknQncAtpOPqpPWxX0d/XUm0s8ExlRtJE4JH8t4VUU7grLxsPPBkR9wM7AB3A+cDupL89kZqmXpa0ZR8xW5O4A8g2xhI2DF/cH9WjWK6vml/Phs9i15Eg+xoZch2vbR7drIfjvUL6Rd1I9R5PwPci4rRN3M8Lr+4wdfQeDBwYES9K+hmvfS/6q68YXuqy/zeRk0H2FuDKPL03OQkAz5FunrM/aUTk26q2GQX8aRNitgHkmoNtjP8FRlW1ayNpH1LH4geUxsIfDXyQNDR4f0yQdGCebgd+kad/3sO+Hwd2yG3co4DD+9j/zXk/m+dfqX/Vy2v8kKTt8ut7fR9x9KS7/dwIHClph0pZ/uXdm+eA3n5Vbw08nRPDHsABdWzX39fyqoh4GhgmqZIg9iUnC0lTSU2FlYSwT9X0OcAJuabwatLI788TkW6wYy3AycH6LTcFfRA4WOlU1iXAPwN/AC4GbiXdG+I7EXF7P3e/DDhe0m+BbYEL8jFv627f+ctkbi6/Abi3j9hvI/2ivRP4b9JNn7pbbwlwJnBTbqo5t7c4ejlezX4iYinwBeB6SXfluN/QR9xPAr/MHcfndLPKT4Dh+X07i9S01Ot2/X0t3bieDc1+bwJel1/jP5FukFO5SdDepD4QIuKaiPh1Lt+TVAuF1EdxbT+ObQ3m+zlYy5A0CbgmIvYqOhbrm6T9gL+PiGMk/Q7YLyKe28h9XQWcGhH3DWiQttFcczCzjZJrHv8naes0u9GJYSTwQyeG1uKag5mZ1XDNwczMajg5mJlZDScHMzOr4eRgZmY1nBzMzKyGk4OZmdVwcjAzsxr/H7/sPieDOW4AAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# read in the data\n", + "df = pd.read_csv('dose_resp.dat', header=None) \n", + "data = df.values\n", + "t, y = data[:, 0], data[:, 1]\n", + "\n", + "# plot input vs output\n", + "plt.scatter(t, y, marker='*', color='b')\n", + "plt.xlabel(r\"Compound concentration ($\\mu M$)\")\n", + "plt.ylabel(\"Percent activity\")\n", + "plt.title(\"Dose–response data\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "82510304", + "metadata": {}, + "source": [ + "## Fit the model with the good starting point:" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "id": "061e8555", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " E04GN, Nonlinear Data-Fitting\n", + " Status: converged, an optimal solution found\n", + " Final objective value 1.219791E+03\n", + "\n", + " Primal variables:\n", + " idx Lower bound Value Upper bound\n", + " 1 -inf -3.81844E+01 inf\n", + " 2 -inf -3.36193E+00 inf\n", + " 3 -inf 3.25418E+01 inf\n", + " 4 -inf -3.39512E+00 inf\n" + ] + } + ], + "source": [ + "# create a handle for the model\n", + "nvar = 4\n", + "handle = opt.handle_init(nvar=nvar)\n", + "\n", + "# register residual structure\n", + "nres = len(t)\n", + "opt.handle_set_nlnls(handle, nres)\n", + "\n", + "# define the residual callback function\n", + "def lsqfun(x, nres, inform, data):\n", + " rx = np.zeros(nres, dtype=float)\n", + " t = data[\"t\"]\n", + " y = data[\"y\"]\n", + " \n", + " # fit the hill function to the data\n", + " for i in range(nres):\n", + " rx[i] = (y[i] - (x[0] + ((x[1] - x[0]) / (1 + (x[2] / t[i])**x[3]))))\n", + " \n", + " return rx, inform\n", + "\n", + "# define the residual gradient\n", + "def lsqgrd(x, nres, rdx, inform, data):\n", + " t = data[\"t\"]\n", + " nvar = len(x)\n", + "\n", + " for i in range(nres):\n", + " rdx[i*nvar] = -1 + (1 / (1 + (x[2] / t[i])**x[3]))\n", + " rdx[i*nvar + 1] = -(1 / (1 + (x[2] / t[i])**x[3]))\n", + " rdx[i*nvar + 2] = (x[1] - x[0]) * ((x[3] * x[2]**(x[3] - 1)) / (t[i]**x[3] * (1 + (x[2]**x[3] / t[i]**x[3]))**2))\n", + " rdx[i*nvar + 3] = (x[1] - x[0]) * np.log((x[2] / t[i])) * (x[2] / t[i])**x[3] * (1 / (1 + (x[2] / t[i])**x[3])**2)\n", + " \n", + " return inform\n", + "\n", + "# create the data structure to be passed to the solver\n", + "data = {}\n", + "data[\"t\"] = t\n", + "data[\"y\"] = y\n", + "\n", + "# set loss function to l2-norm and printing options\n", + "for option in [\n", + " 'NLDF Loss Function Type = L2',\n", + " 'Print Level = 1',\n", + " 'Print Options = No',\n", + " 'Print solution = Yes',\n", + " ]:\n", + " opt.handle_opt_set (handle, option)\n", + "\n", + "# mute warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "# use an explicit I/O manager for abbreviated iteration output:\n", + "iom = utils.FileObjManager(locus_in_output=False)\n", + "\n", + "# set initial guess and solve\n", + "x = [-40., -5., 30., -5.]\n", + "\n", + "sol = opt.handle_solve_nldf(handle, lsqfun, lsqgrd, x, nres,data=data, io_manager=iom)" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "92bc7d4d", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABUv0lEQVR4nO3dd3hUZfbA8e+ZmfQGCb1DAKWHpqLwA7uIDbFgW3RdEctasbC6igV7byirK9agIGBnFWygooB06U16LwnpM+f3x70JExJIgCSTkPN5nvvM7ffcyWTOvO9773tFVTHGGGOCeUIdgDHGmMrHkoMxxpgiLDkYY4wpwpKDMcaYIiw5GGOMKcKSgzHGmCIsORhjjCnCkoMpFyKyWkQyRSRNRHaJyC8iMkRE7DNnTBVg/6imPJ2rqnFAU+AJ4B7grdCG5BARb6hjqMxExBfqGExoWXIw5U5Vd6vqZ8ClwCARaQ8gIgki8q6IbBWRNSJyf37JQkRaisiPIrJbRLaJyEf5+xORY0XkWxHZISJLROSSkmIQkdEiMlJEvhKRvcDJItJARD5xj79KRG4JWv84EZkpIntEZLOIPOfObyYiKiKDRWSDiGwUkaFB20WIyAvusg3ueIS7rI+IrBORO0Vki7vtNUHbni0if7qlrfX77fccEZkTVArreJBzbRf0/mwWkX8FvQePBq3XR0TWBU2vFpF7RGQesNcdH7ffvl8UkZeC/n5vueexXkQezU+6B/v7mSpCVW2wocwHYDVwWjHz/wJucMffBT4F4oBmwFLgWndZKnAfzg+YSKCnOz8GWAtcA/iAzsA2oG0J8YwGdgMnufuMBmYBDwDhQAtgJXCmu/6vwFXueCxwgjveDFA3vhigA7A1/1yBh4HpQB2gNvAL8Ii7rA+Q564TBpwNZAA13eUbgV7ueE2gizveGdgCHA94gUHu+xtRzHnGufu5033f4oDjg96DR4PW7QOs2+9vNgdoDEThlPgygDh3udfdd/57MQF4w30f6gC/A9cf7O9nQ9UZrORgKtoGINH9hTkQGKaqaaq6GngWuMpdLxfny6mBqmap6jR3/jnAalV9W1XzVHU28AlwcSmO/amq/qyqAZwv9dqq+rCq5qjqSuA/bkz5x28pIrVUNV1Vp++3r4dUda+qzgfeBi5z518BPKyqW1R1K/BQ0Dnl7/dhVc1V1a+AdOCYoGVtRSReVXeq6h/u/MHAG6r6m6r6VfUdIBs4oZhzPAfYpKrPuu9bmqr+Vor3Jt9LqrpWVTNVdQ3wB9DfXXYKkKGq00WkLk5yu819H7YAz+/3/hX39zNVhCUHU9EaAjuAWji/ntcELVvjLge4GxDgdxFZKCJ/d+c3BY53q1d2icgunC/keiLSWkTS9xveC9r/2qDxpkCD/fbzL6Cuu/xaoDWwWERmiMg5+51H8L7WAA3c8QbFnFODoOntqpoXNJ2BUzIBGIDzhbvGrZLpERTrnfvF2ni//eZrDKwoZn5prd1v+kP2Jb7L3en8mMKAjUExvYFTgoAD//1MFWGNTqbCiEh3nC//aThVQfm/Lv90V2kCrAdQ1U3Ade52PYHJIvITzpfXj6p6+gEOE3uA+eBUB+VbC6xS1VbFrqi6DLjMbQO5EBgnIklBqzQGFgfFvcEd3+Ce08Jilh2Uqs4AzheRMOBm4GP3OGuBEao6ohS7Wcu+X+/724tTnZavXnFh7Dc9FnhWRBrhlCDyE9ZanNJLrf2SXf65FPv3U9XlpTgHUwlYycGUOxGJd395jwHeV9X5qurH+fIbISJxItIUuAN4393mYvcLCWAnzpdWAPgCaC0iV4lImDt0F5E2hxjW70Ca2+gaJSJeEWnvJjBE5EoRqe1WQe1ytwkEbf9vEYkWkXY47R/5Da6pwP0iUltEauG0abxfivcoXESuEJEEVc0F9gQd7z/AEBE5XhwxItJPROKK2dUXQH0Ruc1tHI8TkePdZXOAs0UkUUTqAbeVFJdbNfYDTtXZKlVd5M7fCHyDkzjiRcQjIski0ts9nwP9/UwVYcnBlKfPRSQN51fmfcBzOF+k+f6J82t2JU5p4kPgv+6y7sBvIpIOfAbcqqorVTUNOAPn1/EGYBPwJBBxKIG5yekcIAVYhVOSeRNIcFc5C1joHv9FYKCqZgbt4kdgOTAFeEZVv3HnPwrMBOYB83Hq7B+ldK4CVovIHmAITnUZqjoT51f4KzhftMuBqw9wXmnA6cC5OO/NMuBkd/F7wFychudv2JfQSvIhcBr7qpTy/Q2nMf9PN65xQH13WbF/v1Iez1QComoP+zGmtESkGU4yCSuuOsWYo4WVHIwxxhRhycEYY0wRVq1kjDGmCCs5GGOMKeKouM+hVq1a2qxZs1CHYYwxVcqsWbO2qWrt4pYdFcmhWbNmzJw5M9RhGGNMlSIiaw60zKqVjDHGFGHJwRhjTBGWHIwxxhRxVLQ5GGOqntzcXNatW0dWVlaoQznqRUZG0qhRI8LCwkq9jSUHY0xIrFu3jri4OJo1a4aIhDqco5aqsn37dtatW0fz5s1LvV21rlbauBF694ZNm0IdiTHVT1ZWFklJSZYYypmIkJSUdMgltGqdHB55BKZNg4cfDnUkxlRPlhgqxuG8z9UyOURFgQiMHAmBgPMq4sw3xhhTTZPDypVw+eUQ7T4TKzoarrgCVq0KbVzGmIoVG1v0wYHPPfccbdu2pWPHjpx66qmsWXPA+8SOatUyOdSvD/HxkJUFkZHOa3w81CvuoYnGmGqlc+fOzJw5k3nz5nHRRRdx9913hzqkkKiWyQFg82YYMgQ+/xzq1oXVq0MdkTGmMjj55JOJdqsVTjjhBNatWxfiiEKj2l7KOn6883rjjU6isH77jAmdhz5fyJ8b9pTpPts2iOfBc9sd0T7eeust+vbtW0YRVS3VNjlERTnVSflGjnSGyEjIzDzwdsaY6uH9999n5syZ/Pjjj6EOJSSqbXJYuRKGDoWJEyEjw2mU7t8fnnkm1JEZU/0c6S/8sjZ58mRGjBjBjz/+SERERKjDCYlqmxysUdoYU5zZs2dz/fXXM2nSJOrUqRPqcEKm2iYH2NcoPXgwjBrl3DFtjKk+MjIyaNSoUcH0HXfcwVdffUV6ejoXX3wxAE2aNOGzzz4LVYghU62TQ36jNMCrr4YuDmNMaAQCgSLz7rjjjhBEUvlU20tZQ836dTLGVGaWHELE+nUyxlRmlhwqWGXu18lKM8aYfJYcKlhl7tfJSjPGmHzVOjnMmQM1asC8eRV3zMp4CW1lLs0YY0KjWieHK6+E3budX/IVKf8S2unTnddQV+NU5tKMMSY0quWlrPs/92Lhwn3zVMv/+JXtEtrKWJoxpiLExsaSnp5eaN7w4cOJjY1l6NChXH311Xz77besXLmSiIgItm3bRrdu3Vi9ejWrV6+mTZs2HHPMMQXb/v7774SHh1f0aZSLallymD0bmjYtPK9ZM5g7NyThVAqVrTRjTGXh9Xr573//W+yy5ORk5syZUzAcLYkBqmnJISUFYmIKz4uJgY4dQxJOpVDZSjPGVBa33XYbzz//PNddd12oQ6lQ1TI5AOzcCe3awQMPOFfn7NgR6oiMqca+vhc2zS/bfdbrAH2fOOLdNGnShJ49e/Lee+9x7rnnFlq2YsUKUlJSADjppJN49Sj6ZVVtk8OGDfvGL7kkdHEYYyq/YcOGcf7559OvX79C8/OrlY5G1TY5GGMqkTL4hV+eWrVqRUpKCh9//HGoQ6kwlhyMMaYU7rvvviIlh6OZJQdjTLVVXJfdB9KuXTu6dOnCH3/8URGhhZwlB2NMtVVcl93BRo8eXWh6fNBlfc2aNWPBggXlEValUC3vc6gMrJM7Y0xlVmmTg4icJSJLRGS5iNwb6njKmnVyZ4ypzCplchARL/Aq0BdoC1wmIm1DG1XZsE7ujDFVQaVMDsBxwHJVXamqOcAY4PyyPkgoqnaskztjTFVQWZNDQ2Bt0PQ6d14BERksIjNFZObWrVsP6yD33gs//eS8HoojSSrWyZ0xpiqorMmhRKo6SlW7qWq32rVrH9K2+VU7777rTL/zzqFV7Rxpe4F1cmeMqewqa3JYDzQOmm7kzisTB+qWu6TuusuqvWD8eKdzu06dnNfgTu+MMRUnNja20PTo0aO5+eabAbj66qsZN25cseuvXr2aqKgoUlJS6NSpEyeeeCJLliwB4IcffiAhIYGUlJSCYfLkyYDTw2tKSgrt27fn3HPPZdeuXeV8hoevsiaHGUArEWkuIuHAQOCzstr5qlVFn+kAkJ198C96ay8wxuTL71dp7ty5DBo0iMcee6xgWa9evQp15X3aaacBEBUVxZw5c1iwYAGJiYmVuqO+SpkcVDUPuBn4H7AI+FhVF5bV/uvXh4YNC88TKfmL3toLjDHF2bNnDzVr1jykbXr06MH69WVWIVLmKu0d0qr6FfBVee2/e3dISHCeAuccD777ruTt8tsLBg+GUaOcxmljzJF58vcnWbxjcZnu89jEY7nnuHsOuk5mZmZBl9sAO3bs4LzzzivV/vO7605LSyMjI4PffvutYNnUqVML7feTTz4hOTm5YNrv9zNlyhSuvfba0p1MCFTa5FDevv7a+eUfbONGaNQI8vIOvJ09FMeYo0d+NU++0aNHM3PmTACkmLrn4HnB3XV/9NFHDB48mEmTJgFOtdIXX3xRZPv8ZLR+/XratGnD6aefXoZnU7aqbXI4UOOz3+9UMUVGQmZmxcZkTHVV0i/8UEhKSmLnzp0F0zt27KBWrVrFrnveeedxzTXXlLjP/GSUkZHBmWeeyauvvsott9xSZjGXpUrZ5lDeoqKcxucDLbNGZmNMnz59+Oijj8jJyQGcUsXJJ59c7LrTpk0rVG1UkujoaF566SWeffZZ8g5WVRFC1bLksHIlDB0KH35YdFl2tjUyG2PgnHPOYdasWXTt2hWv10tycjKvv/56wfL8NgdVJTw8nDfffLNg2f5tDvfffz8XXXRRof137tyZjh07kpqaylVXXVXu53OoREu6uL8K6Natm+bXE5bWDTfA66+D1+tUJSUkOPMvuAD27LF7D4wpb4sWLaJNmzahDqPaKO79FpFZqtqtuPWrZbUSOFcd3XgjzJrlvCYmQlqac++CJQZjTHVXLauVYF8CiIoqfNXSyJHOYA3SxpjqrNqWHPL9+ivUrr3vzujS3PVsD+oxxhztqn1yeP552LrVKSWU9q5ne1CPMeZoV20bpKOi4ORzJnBCm8WAogiK+16I4PFooQ6YVCAnBxTnHokclOwAZPo97M3z0KKNj7CocOLj4mmQWI/m9RrStkFz4vM7YjLGFGIN0hXrUBukq22bgyp0SPCTRIQz7aYGRZ0b5PxaaB6AiqLiTPtR8rwB8sKda6BZ57yksY0lrGQJMAnwqAcVISdMyYn1ElOnJu2aH8NZKT2Ii7TEYYypnKptcli1Cvqd1IHjG29DEMCDR8QZFw8+r3DB+R7AmScioMLcucJffwlRYR7CxUdyYy/NmkGuX/HnBsjN9eNXPwH1k0cuuZ5c0iSTPYEM0nKyYOcWFi7Zwp9fTyPgEdJi/NRu3oRLep5Nct36IX5XjKl+Jk6cSP/+/Vm0aBGDBg0iOzubHTt2kJmZSUO3h86JEyfSqFEj6tevz7XXXssTTzxRsH2fPn1IT08v6HZj5syZDB06lB9++CEUp1NmqmVyyL9CqWPDrdSv+SOCIhJAxKlcCgtTwnwBZk9XVBUNBAiufqsb1Pnilh3OUJywyChiExJJiq5P88h6+KhBXo6HXPWT7slguzeNsLQ9ZM9fy3vz3sDvgbQkD32O78MF3f+vnN8FYwxAamoqPXv2JDU1taDzvPw+ll555ZWC9b7++mtat27N2LFjefzxxwv1s7Rlyxa+/vpr+vbtW+Hxl5dqmRxWroSbboLMyWlcvHhOqbfToFcVIeAR8jyC3+Mhz+vB7/WS5/WQ4/OSEx5OTkQWObvS2BO2nk0eIVP2JRif+GgQn0zrmNYEwuJJD8tmg3cHvq17mPPld/zx1Xek1fJy4ann0evYTmV6/sYYR3p6OtOmTeP777/n3HPP5aGHHjrguqmpqdx6662MHDmSX3/9lRNPPLFg2V133cWIESMsOVR19etD3brwVVZbHt18v1NyyC89oHgEHnhA+eC9AFdeocTGqpMRAgFA+exTZeZM5bhuSr+z1e3Fzylh5LdYa24ugcwMNDOLQGYmgcxMcjMy2JOZzp696ezKzWb37t0sjPqTnDAvXgmjQXhT2iV0ICMukjW+rWzauospqROYGPExjTu146YzLyHMVy3/ZOYot+mxx8heVLZddke0OZZ6//rXQdf59NNPOeuss2jdujVJSUkF3WXsLysri8mTJ/PGG2+wa9cuUlNTCyWHHj16MGHCBL7//nvi4uLK9DxCpdp+04waBc0T8qjZcq+TGhQC6kFEeOJJ4YsvhL/Ux9S1Hq66SgDhH4OFnBwnlQTqe/hmnfC/UUJYGLz7rtNWAeJe5eS+enwQFgW+KOfVHdQbSSBHyd2Zxs6li/lrzizWrl7B/K2fkrtVqB/ZgrY1UtgeE2AJG9g1YymP/v4QHFOXewb8g+iIyBC/g8ZUffmlAYCBAweSmppabHL44osvOPnkk4mKimLAgAE88sgjvPDCC3i93oJ17r//fh599FGefPLJCou/PFXb5LBuHTxx1Vxe7HlX0YU/wRXxcMW57vSnzsub/Q6yw3EHWVYMAbyANyyG+tFJ1E+syfGNktCoOmzMSWTxqkwWLfucwNZwOsel4K3Zgj/D1rNj6Taee2QE0jCGu68baiUJc1Qo6Rd+edixYwffffcd8+fPR0Tw+/2ICE8//XSRdVNTU5k2bRrNmjUDYPv27Xz33XeFnsdwyimncP/99zN9+vSKOoVyVS2/WfIbpMO9/fjw1xVulZJTtXTaqQFqJCj/m6RkZStRkcrZfQPc9y+ldi3l3w8oH41RwsPBnxdg4KXKgw/sq1pCA0HjCoE8yMuC3AzIzYLcTMjLdF6z0yBzJ2Rsh4wdkLEd2b6cBrvX00D9nNwU1mfGM2/XSpaurk+LuBQ6JB7LgvD1bN+cxrP33U/D5nW58vrbin0wiTHmwMaNG8dVV13FG2+8UTCvd+/eTJ06tdB6e/bsYerUqaxdu5aICOfS97fffpvU1NQiD+u5//77GTJkCC1atCj/Eyhn1TI55HfZPWZMBNsyIoiMhMaNoWVLeP9Tp8fW5VvdlfdATjTUbu1MLtwAp1+67zGhc9cDtcs4QH8epG1Adq6h0a41NNq6mIy/5rPgzy+Ys+43mkcfT9ukY/gjYjUrNu3m5Vtu4cSTe9G1/8WWJIwppdTUVO65p/BDhgYMGEBqairHH398wbwJEyZwyimnFCQGgPPPP5+7776b7P0eDHP22WdTu3ZZfyGERrW9Q/qGG5wv9/Bw5xkO9erBH39A8+ZFHx8KlaQjPlXydqxlwRfvM+fHBTSOPI69NWsy17cGCfhptXo+5wz5BzEnnR3iQI0pmd0hXbGsy+5S2rwZhgyB6dOhbVunM72HH3ZKFZdf7nTAB6XriK/CiOBLakLKoH9x1ah3STy/Jbu2/MYp6c2pTU0WJXfk7THfMvvqXgSmj3aqsYwx5jBU25IDFO2uO5/X6zQXhIc7/Sldfz289loZBFoOstLTmf7JGNJ+Xk9EUktmhK9E/Xl0mzuN3m3XEnXeTdDt7xBxdFxeZ44eVnKoWFZyOAQHKiWccca+UsWQIeXTNXdZdfsdGRtLn0H/4MQH/s5uXUSftKbU9CQwvWsvPt54EutefhZ9rj1Me95KEsaYUqvWyaF+fad77qyswt11f/UVvPoqdOrkvJbHk+HKutvvpEaNGfDYQ0SeUoOW2/x0yG3CqhZNGdvgcmbPakJg0kPwSneYP869msoYYw6sWicHKNz2UF6lhGBRUc69cSNHOjdcjxzpTOc/bOhIeDxeup5zHicMvwpf9kb6ZLYmIy6GScecwnd/9iI3LxY+uRbevxB2rjnyAxpjjlrVPjmMH1/+pYRgFdHgnVCnHuc8cQ/SIpfe6U2J8sbwc4eufPp7Q7JT/gVrf4fXesBvo6wUYYwpVrVPDhXtQFVZB3vy3OHw+nz0vn4Q9Qe2JWVXPA0DSSxo156Pxy1g78mjoWkP+PouGHO5cwOeMdXUunXrOP/882nVqhUtWrTg5ptvJjs7m9GjR3PzzTcXWrdPnz4EX/wyZ84cRIRJkyYVWk9EuPPOOwumn3nmGYYPH86IESNISUkhJSUFr9dbMP7SSy8xfPhwnnnmGQCuvvpqxo0r3O3C6tWriYqKKtgmJSWFd999t6zfjgKWHEKgIquykrt3o/N9F9B4Tx7H5jVgRctkPnxrHNta3QVnPQHLvoXXezmlCWOqGVXlwgsv5IILLmDZsmUsW7aMzMxM7r777lJtH9zdd7CIiAjGjx/Ptm3bCs2/7777mDNnDnPmzCEqKqpg/JZbbinV8ZKTkwu2mTNnDn/7299Kd6KHocTkICJJ5Xb0SqysriYqTkVXZdWoW49TnrqBRM2jc25T1jeqz5j3PuIvT2e49hvw+mB0P5g7pnwDMaaS+e6774iMjOSaa64BwOv18vzzz/Puu++Snp5+0G1VlbFjxzJ69Gi+/fZbsoKui/f5fAwePJjnn3++XOMvT6XpPmO6iMwB3ga+1qPhxohSCL6aqLLe43AowiIiOf3RG5j2whh6bElmetIKxn00gfOvuIjk676Hj/8GE66HrUvglH+DxwqVpuJM/Xgp29Ye/Mv4UNVqHEuvS1ofdJ2FCxcW6YU1Pj6eZs2akZeXd9Btf/nlF5o3b05ycjJ9+vThyy+/ZMCAAQXLb7rpJjp27FjqUkhprFixgpSUlILpl19+mV69epXZ/oOV5hugNTAKuApYJiKPicjB3/ESiMjFIrJQRAIi0m2/ZcNEZLmILBGRM4/kOIejPK8mCjURodftl5HYriY9s1qSHhvNxA/Hs3TlOrhqAnS9GqY9BxMGgz831OEaE1K7d+8udn5+/2WpqakMHDgQ2Nfdd7D4+Hj+9re/8dJLL5VZTPtXK5VXYoBSlBzcksK3wLcicjLwPnCjiMwF7lXVXw/juAuAC4E3gmeKSFtgINAOaABMFpHWquo/jGOUysaNMHAgfPSR0yic3ynfxImQkeFcTdS/P7jtREeF7leew8Ivp9LzJ2FazHImvjeGc6/7G23OeQFqNIUpD0F2Olz8tvP8CWPKWUm/8MtL27ZtizT87tmzh02bNtGlSxfGjClc1bpjxw5q1aqF3+/nk08+4dNPP2XEiBGoKtu3byctLa3Qw35uu+02unTpUlBtVZWUqs1BRG4VkZnAUOCfQC3gTuDDwzmoqi5S1SXFLDofGKOq2aq6ClgOHHc4xyit/W9Gq6iriUKtXb9eNDm7HT0zk8mKjODTt95nwV8roNcd0O85WDoJPrjYSRLGHKVOPfVUMjIyCq768fv93Hnnndx8880cd9xx/Pzzz2xyGx5nzpxJdnY2jRs3ZsqUKXTs2JG1a9eyevVq1qxZw4ABA5gwYUKh/ScmJnLJJZfw1ltvVfi5HanSVCv9CsQDF6hqP1Udr6p5qjoTeL2M42kIrA2aXufOK0JEBovITBGZuXXr1uJWOaiDVR9V9I1xodKydwrJA7rQM6MFueHhfPbGO/Tut55Nja+FC/8Da36GMZc5z54w5igkIkyYMIFx48bRqlUrkpKS8Hg83HfffdStW5cXX3yRs88+m5SUFG677TZSU1PxeDykpqbSv3//QvvK7+57f3feeWeRq5ZK4/rrr6dRo0Y0atSIHj16APvaHIIvgS0vJXa8JyKXqOrH+827WFXHlrDdZKC439v3qeqn7jo/AEPdRIOIvAJMV9X33em3cBrBD/qctcPpeG/jxgNXHx1tpYSSrPtjCUtTf2da9ErIy+WvPbfy5su1Ye5HTiN1y9Ng4Afgiyh5Z8aUUmXseO+XX37hsssuY8KECXTp0iXU4ZSpQ+14rzRXK90LfLzfvGHAQZODqp5Win3vbz3QOGi6kTuvzFWX6qPSaHXSMXSqE+CRgfBL5ArqJ7xIWMzd+AKXkjk1A764Dcb9HS5+x7ns1Zij1IknnsiaNda1DBykWklE+orIy0BDEXkpaBgNHPwar8P3GTBQRCJEpDnQCii3u7Pyq48+/xzq1oXVq8vrSJXbypWQ3LMNL3/ehuOzm+P1+rhz6FMsXJwB3a5xbpZb/IVzR3X1uJLZmGrvYG0OG4CZQBYwK2j4DDiiS0xFpL+IrAN6AF+KyP8AVHUhTinlT2AScFN5Xqn06quwYAG8/76TKNxnh1c7+aWoL5d0Y8KkJLpnNyNSfIx673Fy8/LghBvgpNtg5n/h5xdDHa4xpgKUps3Bp6rlVVIoE4f7sB+fD/zFpJ5K8UjQCnbhhU6SGDwYFjz6KrGtopgd/hcZUfDUPcOdVvtProWF4+Git6H9haEO2VRxlbHN4WhWZg/7EZH8dobZIjJv/6HsQq54+VcqFZcYKs0jQStYcJcel48ZQp2Ny2mT24DoTBj28uPOHdMXjIQmPWDCEFg7I9QhG2PK0cGqlW51X88Bzi1mqLL27zYbnEeDilTfRulg4vXS7dGbqL9kCc38tYnYns2j746EsEgY+CHE1YOPr4K0zaEO1RhTTg6YHFR1ozs6AMhV1TXBQ8WEVz6Cr1TK70LowgvhhhuO3nsaDlVYw4Z0uvoUmizdQL1ADfwrNvPalx9DdKJzWWvmLhg7CPJyQh2qMUdk06ZNDBw4kOTkZLp27crZZ5/N0qVLAXjhhReIjIws1JVGSV15p6enc/311xfsr0+fPvz2228AxMbGFtouf1/ffvstPXr0IL+a3+/307lzZ3755ZdyO++SlOYmuDicrjOmisjNIlK3vIOqCPlXKv3xB9x4I+TlVUwPqVVJwjn9OKZJLi3XZpOgMWz/fRGfTP8O6nWA81+Bv36Fb+4LdZjGHDZVpX///vTp04cVK1Ywa9YsHn/8cTZvdkrFqampdO/enfGH8MXwj3/8g8TERJYtW8asWbN4++23S7wJ7vTTT6dp06YFd1K//PLLdOvWjRNPPPHwT+4IlaZvpYeAh0SkI3Ap8KOIrDvM+xgqjeC/9auvhi6Oyq7eA/8m49zz0KiLmV07j0VfTeWPuo3o0uEi2DAbfn0FGnSBlMtCHaoxh+z7778nLCyMIUOGFMzr1KkT4NyNnJ6ezmuvvcaIESNK1T/SihUr+O233/jggw/wuNUSzZs3p3nz5iVu+/zzz9OzZ0969OjBK6+8wu+/h/YZK4dyR9MWYBOwHahTPuGYysYbF0f9h4aTc/0QAmffyYy4jXz131QaDb2FOqc9BBvnwpd3QMOuUDs0naeZqu/70aPYsmZlme6zTtMWnHz14IOus2DBgiJdducbM2YMAwcOpFevXixZsoTNmzdTt+7BK04WLlxY8JS34mRmZhbqcnvHjh2cd955ANSvX5/bbruNHj168NJLL5GYmHjQY5W30nS8d6PbzcUUIAm4TlU7lndgpvKI7d2bmuefT4NvX6NLRkPyPMrI514gV3H6YAqLcu6gzs0qcV/GVBX5XXJ7PB4GDBjA2LFOpxD5XXbv70DzgwU//W3OnDk8nN/jp+umm27C7/dz9dVXH3H8R6o0JYfGwG2qOqecYzGVWN1h95L+y880W/Uj2S17MTdyHY88+RgP3/eAc4nrh5fAt/+Gs58OdaimCirpF355adeuXZEuuwHmz5/PsmXLOP300wHIycmhefPm3HzzzSQlJbFz585C6+d35V2jRg3mzp2L3+8/YOnhYDweT6mSTEU42H0O8e7o08BfIpIYPFRMeKay8NaoQf0HH0QW/M6x8X5a5tbBkxvgkdeeh9Znwgk3we+jYPGXoQ7VmFI75ZRTyM7OZtSoUQXz5s2bxy233MLw4cNZvXo1q1evZsOGDWzYsIE1a9bQvXv3A3blnZycTLdu3XjwwQcLrjxavXo1X35Z9f4vDlatlP+shlk43WgEd6Fx6Lcjmyov7rTTiD+7L6S+SvtWydT31yCweTevTxwDpz0I9VNg4o2wu1z6SjSmzOV32T158mSSk5Np164dw4YN44cffijSJXf//v0ZM2bMQbvyBnjzzTfZvHkzLVu2pH379lx99dXUqVP1mmlL7D6jKjjc7jPMocvbsYOVfc8mvFVL/jrhMmZtWchesuh+wWmc2bQuvN4LGh/nPHa0khSPTeVk3WdUrDLrPiNo4ymlmWeqB19iInXuGkrmzFl0rJ/FsZ6GeMXLHxOmsCwnAs54BFZ+DzOr3pOvjDH7HKzNIdJtW6glIjWD2huacYCns5nqIeHCC4nq2pUtTz3NqTefSbus+uR6Aox79b+ktR8IyafAN/+G7StCHaox5jAdrORwPU77wrEUbm/4FHil/EMzlZV4PNQf/iD+vXvZ9syz9H3oSjpkNSDbF+DFJ57Ff86L4Alz2h8C5dbjujGmHB2sb6UXVbU5zmM8W6hqc3fopKqWHKq5iFatSLrmGnZPnEjO3DmcdvulHJNdlzwJ8Ph/Up1LWtdOh19eDnWoxpjDUJq+lQIiUiN/wq1iurH8QjJVRa0bbyCsUSM2DR9OTGIkPS8+i0Z5NfHvzeSl2enQ5lz4fgRsWRTqUI0xh6g0yeE6Vd2VP6GqO4Hryi0iU2V4oqKo98C/yVm1ih1vvUXjLs3p3Lk7NTWWPavWM6HmAAiPhc/+adVLxlQxpUkOXgm6ZU9EvEB4+YVkqpLY//s/4s46i20jXydnzRq6DjiR1oktCMfHkqkL+bPrXbBuBvz+n1CHakwR+3ehDfDTTz/RpUsXfD5fsXdPVxelSQ6TgI9E5FQRORVIdecZA0DdYcOQsDA2PfIoqspZt11Iy4BzBdOXP+wkrfmpMOVh2PVXqEM1pkRNmjRh9OjRXH755aEOJaRKkxzuAb4HbnCHKcDd5RmUqVrC6tah9q23snfaNNL+9z8ALnhoEMdm1WGvL5fXl3XBr8Dnt8FRcNOlObo1a9aMjh07FtzxXF2V5nkOAWCkOxhTrJqXX8auiRPYPOIxYnr2xBsby/kPXUPu8DdZGrmVF3Ov444VL8C8j6DTwFCHayqZXZ+vIGfD3jLdZ3iDGGqcm1ym+6xOSnOHdCsRGScif4rIyvyhIoIzVYf4fNQfPpy8bdvY+tJLAIRHRtD3xktpkluTPSK8LZfBpHshfWuIozXGlKQ0XXa/DTwIPA+cDFxD6aqjTDUT1bEjNQZeys73P6DGBRcQ2bYtNZvU4sSzTyVj0iTWeRrwbUYyp0+6By76b6jDNZWI/cKvfErzJR+lqlNwOulbo6rDgX7lG5apqurcfjvemjXZOPwh1O9cvnrsSe1pe0xHIjWMWXoyi+b/DEu/CXGkxpiDKU1yyBYRD7BMRG4Wkf5A0eu/jAG88fHUvedusubNY5f75CyAU648gxYxjfGL8pV/ENs/uxdyMkIYqTGQkZFBo0aNCobnnnuOGTNm0KhRI8aOHcv1119Pu3btQh1mSJTYZbeIdAcWATWAR4B44GlVnV7u0ZWSddlduagqf119DVmLFpH81Zf4atUqWPbBfSNZ5ttMzdwwbuqdi+/M4aEL1ISUddldscq8y25VnaGq6aq6TlWvUdUBlSkxmMpHRKj34AMEMjPZ8nThx4YOfHgwydlJ7AzPZdQPAetaw5hKyhqWTbmIaNGCpGv/zu5PP2Pv9N8K5nu9Xgbc93ca5ySwJcLDWy+PhkAgdIEaY4plycGUm1pDhhDWuDGbHn4YzckpmB8dF8MZ11xM7bxY1nuiSX35iRBGaYwpTmnuczipNPMOhYg8LSKLRWSeiEzYr9fXYSKyXESWiMiZR3IcE1qeyEjq/ft+clauZPt/3y60rHGrRnTrdSIxGsHq7QEmT/o6RFEaY4pTmpJDcR3yH2kn/d8C7VW1I7AUGAYgIm2BgUA74CzgNbejP1NFxf7f/xF3xhlsGzmSnLVrCy07/swTSa5dFxWY+/Mclq1cE6IojTH7O9hjQnuIyJ1AbRG5I2gYDhzRF7aqfqOqee7kdKCRO34+MEZVs1V1FbAcOO5IjmVCr+6/hiFeL5sedTrmC3bBP6+hmV9I9+Qw6b8fsWPH7hBFaYwJdrCSQzjO/Qw+IC5o2ANcVIYx/B3Ir1NoCAT/vFzHAZ5XLSKDRWSmiMzcutW6Y6jMwurVo9Yt/2Tvjz+R9u23RZZf/sDtNM+A7b4Mxjz7H7KD2ieMKW8jRoygXbt2dOzYkZSUFH777Tf69OlDkyZNCv2YueCCCwp18b1w4UJOOeUUjjnmGFq1asUjjzxSaP2JEyfSsWNH2rRpQ4cOHZg4cSIAN910EykpKbRt25aoqChSUlJISUlh3LhxXH311UW6Cc8/5urVq2nfvn2R+MeOHUu7du3weDyU6SX9qnrQAWha0joH2G4ysKCY4fygde4DJrDvfotXgCuDlr8FXFTSsbp27aqmcgvk5uqK8y/Qpb37aF5aetHliyfpW3eN0AcffFBHDns2BBGaivbnn3+GOgT95Zdf9IQTTtCsrCxVVd26dauuX79ee/furR06dNCpU6eqqurOnTv1uOOO05iYGFVVzcjI0BYtWuj//vc/VVXdu3evnnXWWfrKK6+oquqcOXM0OTlZV65cqaqqK1eu1OTkZJ07d27BsVetWqXt2rUrFM+gQYN07NixheblH7O49VWd93Hx4sXau3dvnTFjxgHPtbj3G5ipB/heLU2bQ4SIjBKRb0Tku/yhFEnnNFVtX8zwKYCIXA2cA1zhBgmwHmgctJtG7jxTxTkd8z1I3ubNbHul6CPI5ZgzuaLjQhpmRbEpfA9v3P9iCKI01c3GjRupVasWERERANSqVYsGDRoAMHDgQMaMGQPA+PHjufDCCwu2+/DDDznppJM444wzAIiOjuaVV17hiSecK++eeeYZ/vWvf9G8eXMAmjdvzrBhw3h6v/t+ykKbNm045phjyny/pel4byzwOvAmUCbPehSRs3CeCdFbVYP7UPgM+FBEngMaAK2A38vimCb0olJSqHHxxex47z0S+l9A5H4f6Ijznmbgih68n30TG8N38p+HRnLdgzeEKFpTkb7++ms2bdpUpvusV68effv2Peg6Z5xxBg8//DCtW7fmtNNO49JLL6V3794AnHrqqVx33XX4/X7GjBnDqFGjeOSRRwCnSqlr166F9pWcnEx6ejp79uxh4cKFDB06tNDybt268eqrr5YY91133cWjjz56KKdaLkpTcshT1ZGq+ruqzsofjvC4r+C0X3wrInNE5HUAVV0IfAz8ifO0uZtU1R4+fBSpc8fteBMS2PTgcHT/m9/iGxB35r0M8LxMYl4UGwNbeOfp0SGJ01QPsbGxzJo1i1GjRlG7dm0uvfRSRo8eDTg3bPbs2ZMxY8aQmZlJs2bNKiSmp59+mjlz5hQMoVKaksPnInIjTttAdv5MVd1xuAdV1ZYHWTYCGHG4+zaVm7dGDercdRcbhw1j1yefUPPiiwuvcNxg6swdw1lbPueLnHNYm76W1Fc+4rKbLw1NwKZClPQLvzx5vV769OlDnz596NChA++8807BsoEDB9K/f3+GDx9eaJu2bdvy008/FZq3cuVKYmNjiY+Pp23btsyaNYtOnToVLJ81a1aV6sSvNCWHQcBdwC/ALHewXu7MYUu44Hyiu3dnyzPPkrtlS+GFHi+c+yKtAwvpUXcjEepl1dZljPvPxJDEao5uS5YsYdmyZQXTc+bMoWnTpgXTvXr1YtiwYVx22WWFtrviiiuYNm0akydPBiAzM5NbbrmFu+92nqA8dOhQHn/8cVavXg04Vxo99thj3HnnneV8RmWnNB3vNS9maFERwZmKtXEj9O4NZVz1W4SIUO+hh9CsLDYNf6jIvQ80SIHjb6DHttF0bF0brwrL1y3k0/e/Kt/ATLWTnp7OoEGDaNu2LR07duTPP/8sVEoQEYYOHUqtoJ6FAaKiovj000959NFHOeaYY+jQoQPdu3fn5ptvBiAlJYUnn3ySc889l2OPPZZzzz2Xp556ipSUlCOKd8mSJYW6GB87diwTJkygUaNG/Prrr/Tr148zzyybjiVK02V3NHAH0ERVB4tIK+AYVf2iTCIoA9Zld9m48UZ44w24/np47bXyP972/77NlqeeosHTT5Nw7jmFF2anw6vHQ2Q842vewZJFS/DioVPHEzjzolPKPzhT7qzL7opV5l124zwmNAc40Z1eD4S+Kd2UmagoEIGRI50OUkeOdKajosr3uImD/kZUp05sfvRR8va/kTEiFvo9A1v+5MJGa2nWojl54mf+vN/48Ytp5RuYMaZUySFZVZ8CcgHcS0+lXKMyFWrlSrj8coiOdqajo+GKK2DVqvI9rni91H/8MQKZmWx8qJjqpWP6Qptz4ccnuey8HjRs0JAsyWXW7z/zy9e/Fb9TY0yZKE1yyBGRKEABRCSZoKuWTNVXvz7Ex0NWFkRGOq/x8VCvXvkfO6JFC2rfegvpk6ew58ti2hT6PgWeMPjyTgYNvoZ6SXXYK9lM//UHpn7+a/kHaMpVSdXapmwczvtcmuTwIM49B41F5ANgCs4NbOYosnkzDBkC06c7r+XdKB0s8eqriezUkc2PPFK0eim+AZz6AKz4DhZ8wj9uGUK9pFpkSDYzZvzEjxMqdxVTRTXyV0WRkZFs377dEkQ5U1W2b99OZGTkIW1XYoM0gIgkASfgVCdNV9VthxVlObEG6aove8UKVl04gOju3Wk86g3EE/S7JeCHN0+D3Wvhpt8hOpH/vPo6m7dsIUrD6dy+B6dc0jt0wR9ERTfyVyW5ubmsW7eOrKysUIdy1IuMjKRRo0aEhYUVmn+wBunSXK3UH/hOVXe70zWAPqo6sSyCLguWHI4OO1NT2fTQw9Qddi+JgwYVXrhxHozqA52vhPNeAuCt1//Dpo0bidBwOrbuzhlXnlrxQR9AVJRTPbe/yEjIzKz4eIwpzpFerfRgfmIAUNVdOFVNxpSpGgMHEnvyyWx55lmyFi8uvLB+R+hxI/zxDqxx2hquHXIdDZo0Iltymbv0d75+Y1IIoi5eqBr5jSkrpUkOxa1Tmm43jDkkIkL9EY/iqZHA+qFDCez/07vPMEhoDJ/fCnnONRHXXPt3GrVoQi55zN0wi4nPfBKCyIsKZSO/MWWhNMlhpog8JyLJ7vAcThcaxpQ5X2IiDR5/gpzlK9jy1H7dG4fHwDnPw7Yl8OOTBbMHDRrEsV3aIgFlYdqfjH34QwL7d+oXAqFs5DfmSJWmzSEG+DdwGs7lrN8CI1R1b/mHVzrW5nD02fz4E+x45x0avvQi8W6f+QUm3ADzPoLrpkCDzgWzv5v6A7P+N41Mr59j/c0YMPwqvL7S/P4xpno67AZpEfECk1X15PIKrixYcjj6BHJyWHPlVeSsWEGzcWOJcB+aAkDmTnj1BIhOgsE/gC+8YNGcBXOYkvoFab482uY24bx7riAyIaLiT8CYKuCwG6TdZykERCShXCIz5gA84eE0euF5JCyM9bfcSiAj6JlQUTWd6qUtC2Hqs4W2S2mfwsU3/o3YHOHP8L8Y9+TbbF5mzxg35lCVpsydDswXkbdE5KX8obwDMyasQQMaPPMM2cuXs3H48MI3Sx17NnS4GKY+A5vmF9quScMmXDfsNmKy/CwP38Rn73zIgsmLKjh6Y6q20iSH8ThtDj+x73kO1iBtKkRsz5Oo9c+b2fPZ5+xMTS28sO9TTili4o3gzy20KCE+gaGPPUS0P8AG706m/PQF3735fQVGbkzVVprnObyD8+jO6ar6Tv5Q/qEZ46g1ZAgxvf+PzY89zt7pQR3uRSdCv2dh0zz4+YUi24nHw92PPEzNGrGkSSYz1/7K5w+Nw59jT541piQlJgcROReYg9O/EiKSIiKflXNcxhQQj4eGzzxDeLOmrLv1VnLcp2sB0PZ8aHsB/PgUbFpQ7Pa33D6U5A6tCATymB1YyMQH32X3+j0VErsxVVVpqpWGA8cBuwBUdQ5gT4IzFcobF0fjkSMRj4e1Q27Av3v3voX9noXIGjB+MOQW30/PZRddRt+BFxCVE2B+xBomvPYei6csrZjgjamCSpMccoO7z3CF/g4jU+2EN25Mo5dfImf9etbddhuak+MsiKkF57/qXL303SMH3L5Tu47c/O+7iczOYbVvK9/8+CmTn51EIM8+zsbsrzTJYaGIXA54RaSViLwM/FLOcRlTrOhu3aj/8MNk/DqdDf+6D82/E7r1GdD9H/DrK7DyhwNuHxUdwz2PjSAhIZI0Mvltzww+f+BD0jdUmns6jakUSpMc/gm0w3nAz4fAbuC2cozJmIOq0f8Cat9+O3u++ILNjz+x7xLX0x+BpFbOHdSZOw+4vYhw+x33ktKrC2G5fmaHL2f8a+8wZ8K8CjoDYyq/AyYHEYkUkduAp4C/gB6q2l1V71dV64DdhFTS4OtIHDSIne+9x/bXX3dmhkfDgP/A3i3wxe1QQtcw/U4/h+vu+icRWdms9G7hu9lf8+XwCeSk51TAGRhTuR2s5PAO0A2YD/QFnqmQiIwpBRGhzj13k3D+eWx98SV2vP+Bs6BBZ6f31oUTYG7qwXcC1KyZxL2PP0ad+jXI1ixm6jw+few9lk1eXs5nYEzldsC+lURkvqp2cMd9wO+q2qUigyst61up+tLcXNbddjvpU6ZQ9/77SbzyCufJce+eD+tnOX0v1T6mVPuat3AuX737EVlR4SQGYmntbcUpd/QlPC685I2NqYIOt2+lgltOVTWvzKMypgxIWBiNnn+O2FNPZfOjj7LjvffB44UL/wNh0fDxIMjJKHlHQMd2nbj7sUeIifKQrnv5LTCHiY+/y+KvlpTzWRhT+Rys5OAH8i/hECAKyHDHVVXjKyTCUrCSg9GcHNbfeSdp307e95jR5VPg/QHQ+QrnUtdDMHveH3zzwTgyo8JJCERzTKA5PQefRXyTuHI6A2Mq3hE9Q7oqsORgwKliWn/nUNK++YakIddT+9Zbke9HwE9PwwWvQ8plh7a/QIDnXnycnO25ZPugeV4dmtduy4k39MQXbg9DNFVfpUsOIvIIcD7OzXRbgKtVdYOICPAicDZOKeVqVf2jpP1ZcjD5NC+PTQ89zK6xY0kYcCH1H/g38uGFsGE2XPc91Dn2kPe5evUKPnz9DfIiovHipU1OY9r3PImWfZPxeOxhQqbqqozJIV5V97jjtwBtVXWIiJyNc1/F2cDxwIuqenxJ+7PkYIKpKttefoVtr71GbJ8+NBx+F573zoCIeLjuO4iqcVj7nfDZWJb+MpfMSB8JgWha5TWh64A+1O9qD4Y2VdNhP+ynvOQnBlcMzuNHwSlNvKuO6UANEalf4QGaKk1EqH3LP6k3/EHSf/qJ1YNvJ+f/noNdf8En/3CuZjoM/c+7mNuG309MtJDtz2Bm+GImTEzl6wfGs2eNdeRnji4ha3MQkRHA33DuuD5ZVbeKyBfAE6o6zV1nCnCPqhYpFojIYGAwQJMmTbquWbOm4oI3VUb6zz+z/o47EREa3nAGMategJ53wGkPHtF+t2/dzOsvPw8aRa5XaeKvRePwFpxwdS/iGlujtakaQlJyEJHJIrKgmOF8AFW9T1UbAx8ANx/q/lV1lKp2U9VutWvXLuvwzVEi9qSTaP7xR3hrJfHXU5+wPb0POvU5WDD+iPabVLsu9z38BBdfdQnhWRmsl+38nPc7n7zxLt88+gVp69MA2LgReveGTZvK4myMKaw8P18hv1pJRJoAX6lqexF5A/hBVVPdZUuAPqq68WD7sDYHUxJ/+l42DruXtG8nE9MikgZdN+O7cRLU71gm+58+fSpTJnxGICIWv0dp6q9Fg7Cm/LypFyPeqMH118Nrr5XJoYwpcOON8MYbHPbnqzI2SLdS1WXu+D+B3qp6kYj0wylF5DdIv6Sqx5W0P0sOpjRUlZ0ffsiWJ5/E68ulQR8l5sHJkNCozI4xY8YvfDtuIoGIGPI8SkN/Igl76/PWpBP5ZlFDIiMhM7PMDmeqqagoyCqmh7tD/XxVxuTwCXAMzqWsa4AhqrrevZT1FeAsnEtZrymuvWF/lhzMochavJj1t9xEzl8bqNkhgtqvf4U3qUGZHmPKlN/48ctxRMTEkutVagZiaJxdj8YdutLlkmPwhnvL9Himetm4EYYOhYkTISMDoqOhf3945hmodwgXz1XGq5UGqGp7Ve2oqueq6np3vqrqTaqarKodSpMYjDlUkcceS/OJn5N4wSnsnJ/Fyr5nkv7jD2V6jFNPPZ6tmU/z8htXsGNjFln+vcyLWsGUZRMY/9Bovn9mMmlr08r0mKb6qF8f4uOd0kNkpPMaH39oiaEkIW9zKAtWcjCHK+OjJ9n47H/I2RNG/Dn9qDN0KGFl9B924YXOP/HgwTDq9d3AKzSNzyQj0rm7umEgkdqB2rTr1Z3k01rg8dkNdab0Cn2+RjmlifGHeJ1FpatWKmuWHMyRCHzzKNvfeJ3tSxIgLIJa1w8m8Zpr8ERElP2x/H6++noC86f+jobHkONVojScZnl1SIhpQLdzj6NWh1plflxjimPJwZiDUYWv7yHnuzfZsv440mavIaxRI2r/82bizzkH8ZZP+8DWjet5551RsDPA3kgvKlAzEEuDvERq12tB1wFdiWto90yY8mPJwZiSBALw6U0w90P2NrmJzZ8vJnvRIsJbtKD2zTcRd9ZZSDn1oxTw+5n647f8Mvk7fESz1318RD1/DWr5a1CrfjO6nNeF+KaVpiNkc5Sw5GBMafjzYNw1sOgztN+LpO1owLZXXiZ72XIiWrcm6dq/E9+3LxJefg//yc3O4rPPxrFs5nzEF01mmPP/WSsQR+28GtRIakT3ft2p0bqGdfpnjpglB2NKKy8bxlwOyydDv2fRLtew5+tJbHt9JDnLV+CrW5eaV15BzUsvxRtfvr/ks9LTmfjZR6yeuxSvN5q94c7/anwgirr+GkSH1aB15w60Pq0lYTH2tDpz6Cw5GHMocrNg7CBYOgnOegJOuAENBNg7bRrb336bjF+nI9HRxPc9i5oXX0xkp044t+iUY0jZWUz+9kvm/jKTsEAEeyOEgChe9VA3kEC8P5aEWg3pfHIKdTrUsSufTKlYcjDmUOXlwCd/h0Wfw+kPw0m3FizKWrSIHe+/z56vvkYzM4lo1ZIaF11EfL9++GqV/5VGgYCfZQvn8dXXX+DflkUgPIIMXwCASA2jjj+e6EAMsbXr0bVPF+p2qmNVUKZYlhyMORz+XBg/GBaOh153win/hqASgj89nT1ffcWucZ+QNW8eeDxEH38c8Wf1Je6M0/HVrMnGjTBwIHz0UdneoBQsO2Mv30z+ij9/m01Yno+ccB9ZXidZhKuP2v54ogPRRMYlckznNrQ+KZnwOKuGMpYcjDl8/jz48g744x3odDmc9xJ4w4qslrV0KXu+/pq0r74mZ80a8HqJOf44vt7ah8e/7c3Zf29aYR3v7d21k2+/n8Sy2QvxZfvICfOS6XOeYSEKNTSGBH80Pokivl49uvZMoV77enitKqraseRgzJFQhR+fgh8eg+RT4ZJ3ISL2AKsq2YsW8e9TJ9EncjItIlYBsDqnKT+l/x+/5/bim7Vd8MbGVFj4Wenp/PLLD8yZ8QeyOw+8PjLCINfjlC686qFmIIaYQAQ+iSSyRk1ap7ShzQktCY+1EsbRzJKDMWVh1jvwxe1Qrz0MTIWEhgdcNb9jtBmf/0V371ROjv+R46J/J0yzweslsm1bort3J7pbN6K7dsGbkFBhp6Gq7Ni4nik/TWH9wlV4soWAz0uGL1CQMABiA5HEBSKJCITjCY8irm4d2qQcS8vOTfG5XYCYqs2SgzFlZen/YNy1EBYFl74HTU444Ko33OD0eRMeDjk5cNM/Mnl80Gz2zphBxowZZM2dh+bmAhDerBmR7dsT1aE9kR06ENmmDZ6oqIo6KzQQYMemjfw642dWzV9GIC0PES/ZPsjw5qFBF2PFBCKICUQQEQjD4w0nPD6eeslN6NStDYmN7P6LqsSSgzFlactiGHMZ7FoL/Z6FroOKXa2kjtEC2dlkzp1L5h+zyVwwn6z5C8jbvNlZ6PEQkdyC8JYtiWjVigj3NbxJk3LrzqM4uVlZrFq2iN9+n8Guv7ZCVgARD7leIcOXR67sK2mICjEaQVQgjHD1IfjwREYQU682rdq3pH1KK8KjrJqqMrHkYExZy9wJ4/4OK76Drlc790OEHfkv/dwtW8hasJCsBfPJWryE7GXLyF23zmn3ACQ8nPAWLQhv0oTwJo0Ja9LEGW/cGF+9ehWWOAIBPzs3bWLW7Jmsmr+CvN2ZiB8CXiHPC5meANmevELb+NRLlIYREfDhUy+CF/H58MZGE1+/Ni3aNKd9u5ZERFoCqSiWHIwpDwE/fPcITHse6rSDi0dD7dZlf5iMDLJXrCR7+XKyly0je8Vycv9aS+66dQXVUgASFkZYw4aENWiAr149wurVxVfXfa1XD1+dOnhr1Cj3G/Y0EGDv7l0sX7KYhbP/JH3zTgJZeaAQ8DjJI9vjJ1NyC1VXAXhUiNIwwgM+fOpB1OP0aeXzIlGRRCXFU6tJfVq3aUGLBvXxVmAp6mhkycGY8rRsMky4HnIz4OxnIOXyQvdDlBf1+8nbtImctWvJ+esvcteuJeevteRu2kjeps3kbd3qdCgYRCIi8NWujTcpEV/NROc1MRFvzUR8SYl4E53BV7MmnvgEPDHRZZ5MVJXsjL1s3biBP+ctYvOqDWTvyiCQ60eAgECeB3I8AXI8frLIhWJC8KgQoT7C1Ys34MGDgAp4BLxeiPDhi4kiOimBWo3r0KRZQ1rUa0h0RGSZnk9VZsnBmPK2ZyOMvw5WT4U250K/5yG2dkhD0rw88rZtI2/TJnI3bSZvy2ZyN24ib9s2/Dt2kLdjR8ErQSWQQrxevHFxeOLj3dc4vHHxeBPi8cTF442PwxMbhyc62hli3NeoKDzR0Uj0vulDrfLSQIDM9DS2bdjM0oXL2LZuC5m70/Fn5YI/gLq5IM+j5ImS6wmQI35yKFoiCRamXsLUg1c9eFTwqCDuBuoB9XjA50EiwvDGRBAeH0N8nZrUrl+bhrXq0jCpFnGR0Yd0LuXlSG+ytORgTEUI+OGXl+D7xyAizilFtL8w1FGVSFUJpKfj376dvB078e/YTt7OnQT2pOHfs4dA2h78e9Lwp+0hsHsP/rQ0AnucVy3uKfcHIJGRBUnDEx2FREQiERF4IsKR8AgkIgKJCMcTEeEuc8f3XxbujIsvDAkLQ8J8iM+HeoTMnFzSdmayftMOduzcw969WeRm56J5AQLud13AAwFR/O6QJwHyJECu+Mkhj4CU/J3oUcGrHrxIQYLxQEGSKXhv3QSmPkF9XiTMgycyDF9sFBFxkUQnxBIbG09cTBwJsfHUiImhZnQcCdGxhPlKvlz4xhvhjTfg+us5rJssLTkYU5G2LIaJN8CGP6Dt+dD3KYgrp74zQiyQk0MgLY1ARgaBjEwCGXsJZGSgmZnuvPz57nhmBpo/np2NZueg2dkEcvaNO9P7ximP7ygRxOeDsDAnsfh85EbEkh1Zg+zIeDIjY8iJiCbXF0FemI+ALwy/10vA6yHgETfBuENQsskjgD8/2eAnDz9+CZQcTzHyk44XQTR/2qlhk4Bztzuq7N7l46FRwwDnedKZmYfyNhw4OdidLMaUtTrHwrXfOqWIHx6H5d/Byf+C4waD9+j6l/OEh+NJSoKkpHLZv6pCbm6hZBFwXzU3D83NRfNyIS8PzR9yct3xoPkF84Lm5+Y6+3DnaV4e+AMQ8KN5fjSQCf69aMAPeQE0x79vuT8Afj8aKPqq/gC5hJEjUeR5I8mVCLJ9EeSGhZHrDcfvC8fv8xLwegl4PKjXS8AjqEdQEXfcSTwqQUlIdN/gdRMRAZrGeIiOhv794Zlnyu69P7o+qcZUFl4f9LrDKTl8fQ/8bxjMfs+pamp2UqijqzJEBMLD8YaHQ2zxXZYcrVTVuaAgEADVQtN+v59778jjy/G5ZEsEWVkQH1+2nTvarYzGlKekZLhiLFz6AWSnweizYcwVsHVpqCMzlZyIIF6v064S7rS3eKKi8MTEEBYfz6odiZx6WV0mfleDIUNg06YyPr61ORhTQXIy4NdX4OeXIHcvdL4S+gyD+AahjsxUUwdrc7CSgzEVJTwaet8Nt86B466HOanwUmeYNAz2bAh1dMYUYsnBmIoWUwv6PgH/nAXtB8Bvb8CLneDzW2HHylBHZwxgycGY0KnZFC54DW6ZDZ2vckoSL3eFsVfDml/L5xJOY0rJkoMxoVazKZzzHNw2D3rc7HTm9/ZZ8EYv+ONdp63CmApmDdLGVDY5e2Hex/D7KNjyJ0TWgA4XQ8pl0KBLhfTbZKoHu0PamKpIFdb8DDPegsVfgj8bah8LnS6DjpdCfP1QR2iquEp7tZKI3CkiKiK13GkRkZdEZLmIzBORLqGMz5iQEoFmPeHit2HoUjjnBYhMgMkPwnNt4L994dfXnIcOGVPGQnaHtIg0Bs4A/gqa3Rdo5Q7HAyPdV2Oqt6ga0O0aZ9i+AuaPhT8/c+68/t8wp7qp7XnQ6kyo08aqnswRC1m1koiMAx4BPgW6qeo2EXkD+EFVU911lgB9VHXjwfZl1Uqm2tq2HBZ95gwbZjvz4hpAy1Oh5WnQoo+TWIwpRqXreE9EzgfWq+rc/R4k0hAILiOvc+cVSQ4iMhgYDNCkSZPyC9aYyqxWS6cPp153wO71sGIKLJ/slCpmvwfihYZdoEkPaHoSNDkeomqGOmpTBZRbchCRyUBx3UDdB/wLp0rpsKnqKGAUOCWHI9mXMUeFhIbQ5W/O4M+D9TOdRLFqKkwf6fQSi0DddtD0RGh0HDToDIktwGNXtZvCyi05qOppxc0XkQ5AcyC/1NAI+ENEjgPWA42DVm/kzjPGHAqvD5qc4AwAuZmwfhas+cW5Amr2B86lsgARCdCgk5MoGnSB+p2gRlNLGNVchVcrqep8oE7+tIisZl+bw2fAzSIyBqchendJ7Q3GmFIIi3KufGrW05n258HWxU47xYY/nNdfX4OA+7jQsBjnuRR12kCdtvteY+taY3c1Udme5/AVcDawHMgArgltOMYcpbw+qNfeGbpc5czLy4bNC2HTPOdpdlv+hKXfwOz3920XmQCJyU5VVJL7mtjCmRedaInjKBLy5KCqzYLGFbgpdNEYU435IpzG64b73V60dxtsWeQMWxc7nQOumwELx4MGPQIzIgFqNoGExhDfEBIauUNj5zWuHni8FXtO5rCFPDkYYyq5mFrQvJczBMvLhl1/Ocli+wrYscK5IW/XX067RtbuwuuLF+LqQ1xdp3oqprbzGlsn6NUdD4+puPMzxbLkYIw5PL4IqNXKGYqTtQf2rIfd62D3WudS2z3rIX2zk0DWzYS9W4FiLjYMi3YuuY1KdO7TiE50x2u64zX3TUfVhMh4iIhz2kqsIb1MWHIwxpSPyHhnqNPmwOv48yBju5Mw9m6B9C3u+DbI3AkZO5zXLYv2jav/IAcViHATRXFDZILzGh7rPHwpLHiICpoX5SSasCjwRVbLhGPJwRgTOl6fU80UV7d066tC9p59iSJzB2Tucp7Pnb3HfU0rPJ21yym55M/PST/0OAsSRvS+hOGLcF694c64N9ydHw7eCHd5hDt+gHm+SPCGgSfMeQ0e94Q570/BtC9ofni5JyxLDsaYqkPE+fUfmYBzu9RhCPidBJGb6XSPnpvpDvnjGc4zNPLH84dC8zKdXnLzciB3l/Pqz3baYfKy9y3zZ0MgryzfgX3E4ySKk26BU+4v891bcjDGVC8eb1CCqQABf9GEkZcDeVnOuD/XHXKcROLPde438ee5rznuvP2W+XOc8cbl0zepJQdjjClPHq/TlkF0qCM5JNWvlcUYY44SGzdC796waVPZ79uSgzHGVFGPPALTpsHDD5f9vi05GGNMFRMV5bTNjxwJgYDzKuLMLyuWHIwxpopZuRIuvxyi3WaM6Gi44gpYtarsjmHJwRhjqpj69SE+HrKyIDLSeY2Ph3rFPUHnMFlyMMaYKmjzZhgyBKZPd17LulHaLmU1xpgqaPz4feOvvlr2+7eSgzHGmCIsORhjjCnCkoMxxpgiLDkYY4wpwpKDMcaYIiw5GGOMKUJUi3lEXxUjIluBNYe5eS1gWxmGU5XZe+Gw98Fh78M+R+t70VRVaxe34KhIDkdCRGaqardQx1EZ2HvhsPfBYe/DPtXxvbBqJWOMMUVYcjDGGFOEJQcYFeoAKhF7Lxz2Pjjsfdin2r0X1b7NwRhjTFFWcjDGGFOEJQdjjDFFVOvkICJnicgSEVkuIveGOp6KIiKNReR7EflTRBaKyK3u/EQR+VZElrmvNUMda0UQEa+IzBaRL9zp5iLym/u5+EhEwkMdY0UQkRoiMk5EFovIIhHpUR0/EyJyu/t/sUBEUkUksjp+JqptchARL/Aq0BdoC1wmIm1DG1WFyQPuVNW2wAnATe653wtMUdVWwBR3ujq4FVgUNP0k8LyqtgR2AteGJKqK9yIwSVWPBTrhvCfV6jMhIg2BW4Buqtoe8AIDqYafiWqbHIDjgOWqulJVc4AxwPkhjqlCqOpGVf3DHU/D+RJoiHP+77irvQNcEJIAK5CINAL6AW+60wKcAoxzV6ku70MC8H/AWwCqmqOqu6iGnwmch6BFiYgPiAY2Ug0/E9U5OTQE1gZNr3PnVSsi0gzoDPwG1FXVje6iTUDdUMVVgV4A7gYC7nQSsEtV89zp6vK5aA5sBd52q9jeFJEYqtlnQlXXA88Af+Ekhd3ALKrhZ6I6J4dqT0RigU+A21R1T/Ayda5xPqqvcxaRc4Atqjor1LFUAj6gCzBSVTsDe9mvCqmafCZq4pSWmgMNgBjgrJAGFSLVOTmsBxoHTTdy51ULIhKGkxg+UNX8p9FuFpH67vL6wJZQxVdBTgLOE5HVONWKp+DUu9dwqxSg+nwu1gHrVPU3d3ocTrKobp+J04BVqrpVVXOB8Tifk2r3majOyWEG0Mq9CiEcp9HpsxDHVCHcevW3gEWq+lzQos+AQe74IODTio6tIqnqMFVtpKrNcP7+36nqFcD3wEXuakf9+wCgqpuAtSJyjDvrVOBPqtlnAqc66QQRiXb/T/Lfh2r3majWd0iLyNk4dc5e4L+qOiK0EVUMEekJTAXms6+u/V847Q4fA01wukC/RFV3hCTICiYifYChqnqOiLTAKUkkArOBK1U1O4ThVQgRScFpmA8HVgLX4PyArFafCRF5CLgU56q+2cA/cNoYqtVnolonB2OMMcWrztVKxhhjDsCSgzHGmCIsORhjjCnCkoMxxpgiLDkYY4wpwpKDMcaYIiw5GGOMKcKSgzksIlJPRMaIyAoRmSUiX4lI61DHdahEZLiIDA11HCVxn7VwY1lsJyK/lFFMUSLyo9v9fUnrXi8i6t5smD/vJnfemSLyU1D3FKYSsORgDpnbrcAE4AdVTVbVrsAwjvIeO0OsBlBschDHgf6Xi2ynqieWUUx/B8arqr8U63YA5gLHAohINM6dx1uBOTjPiri0jOIyZcCSgzkcJwO5qvp6/gxVnauqU0XkDvcJWgtE5Lb85SLSzH3C2GgRWSoiH4jIaSLys/uUseOC1vnAfRLZOPdLJH8fRfbtbrMgaJ2hbmmgmbuP/7hP9fpGRKLcde5zY5gG5PclVISI/E1E5onIXBF5rxRxHOh4RfYjIleKyO8iMkdE3sj/9X2Q/TwBJLvrP+2ut0RE3gUWAI1FZKJbilsoIoPdcAtt5x4j/XDPZT9XENTHkIj8ICL5X/5JwX8XoCNO9xPHutO3AGOBgKpuBia6+zOVharaYMMhDTj/2M8XM78rTn9NMUAssBDo7C5rhtNXTQecHyWzgP8CgtNF8kR3HQVOcrf5L05/Rwfct7vNgqAYhgLDg46X4s7/GLgyaD/RQDywPP8Y+51LO2ApUMudTixFHMUdr8h+gDbA50CYO+814G/7vU/772f/82yG0y/WCUHz8mOMwkkYSftv5y5PP5xz2W8f4cCm/eatAzzu+MlAatCyrUALYBJOaWY20AeY7C73AltD/dm2Yd9gJQdTlnoCE1R1r6qm43R33Cto+SpVna+qAZwvoinqfDPMx/lCAlirqj+74++7+yzNvouzSlXnuOOz3GP0cveToc4zLA7UE+8pwFhV3Qag+zqbO1gcxR2vuP2civPFPENE5rjTLUqIuzhrVHV60PQtIjIXmI7THX2rA2yX71DPJVgtYFf+hIg0Bda7f1twSgrz3GWNge2quhKoA9wFvAy0xvnbo07VVI6IxJUQs6kg1gBkDsdC9nVffCiCe7EMBE0H2PdZ3L8nyJJ6hsyjcPVo5AGO58f5RV2eSns8Ad5R1WFHuJ+9BTt0GnpPA3qoaoaI/EDh9+JQlRRD5n7774SbDFxdgY/c8Q64SQBIw3l4znE4PSL/EbRNBJB1BDGbMmQlB3M4vgMiguq1EZGOOA2LF4jTF34M0B+na/BD0UREerjjlwPT3PGpB9j3ZqCOW8cdAZxTwv5/cvcT5f5KPfcg53ixiCS555dYQhwHUtx+pgAXiUid/HnuL++DSQMO9qs6AdjpJoZjgRNKsd2hnksBVd0JeEUkP0Gk4CYLEWmFU1WYnxA6Bo0/DdzslhQKkob7/mxT5wE7phKw5GAOmVsV1B84TZxLWRcCjwMbgNHA7zjPhnhTVWcf4u6XADeJyCKgJjDSPeYfxe3b/TJ52J3/LbC4hNj/wPlFOxf4GuehT8WttxAYAfzoVtU8d7A4DnK8IvtR1T+B+4FvRGSeG3f9EuLeDvzsNhw/XcwqkwCf+749gVO1dNDtDvVcivEN+6r9OgEe9xwfwHlATv5DgjrgtIGgql+o6q/u/LY4pVBw2ii+PIRjm3Jmz3MwlYaINAO+UNX2oY7FlExEugC3q+pVIrIM6KKqaYe5r/HAvaq6tEyDNIfNSg7GmMPiljy+F5EEZ/KwE0M4MNESQ+ViJQdjjDFFWMnBGGNMEZYcjDHGFGHJwRhjTBGWHIwxxhRhycEYY0wRlhyMMcYUYcnBGGNMEf8PHslsmkZX5dwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# plot result from NLS\n", + "plt.title(\"Dose–response curves\")\n", + "plt.plot(t, y, '*b')\n", + "dose = np.linspace(0.0005874, 91.74, 1000)\n", + "resp = sol.x[0] + ((sol.x[1] - sol.x[0]) / (1 + (sol.x[2] / dose)**sol.x[3]))\n", + "plt.plot(dose, resp, label='L2')\n", + "\n", + "# fit all loss functions and plot results\n", + "# loss function types\n", + "lossfuns = {'HUBER', 'L1', 'LINF', 'CAUCHY', 'ATAN', 'SMOOTHL1', 'QUANTILE'}\n", + "\n", + "# turn off log printing\n", + "opt.handle_opt_set (handle, 'Print File = -1')\n", + "\n", + "for lossfun in lossfuns:\n", + " \n", + " # set option for the loss function\n", + " opt.handle_opt_set(handle, 'NLDF Loss Function Type = ' + lossfun)\n", + " \n", + " # call the solver\n", + " sol = opt.handle_solve_nldf(handle, lsqfun, lsqgrd, x, nres, data=data, io_manager=iom)\n", + " \n", + " # calculate response using fitted parameters\n", + " resp = sol.x[0] + ((sol.x[1] - sol.x[0]) / (1 + (sol.x[2] / dose)**sol.x[3]))\n", + " \n", + " # plot curve\n", + " plt.plot(dose, resp, label=lossfun)\n", + "\n", + "# show the plot\n", + "plt.xlabel(r\"Compound concentration ($\\mu M$)\")\n", + "plt.ylabel(\"Percent activity\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "937d3904", + "metadata": {}, + "source": [ + "## Fit the model with the naive starting point and different regularisation functions:" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "839d5752", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3nUlEQVR4nO3deXxU5fX48c+Z7AESZIcECCgiqOwouFuqgBtiXcGK9lsBpYootVrboqBWf1i1tUClVVGxLigqKG5gAVERQUAFVFYl7LKHkJDl/P64d8KQyTJJZkvmvHVed+5+ZpjMmed57n0eUVWMMcYYX55IB2CMMSb6WHIwxhjjx5KDMcYYP5YcjDHG+LHkYIwxxo8lB2OMMX4sOZhaS0TOE5HsEB7/RhFZFKrjl3G+wSKyWURyRKR7GM87VEQ+DNf5TO1gycFEjIhsEpFflloW1i/kKPMY8DtVra+qy0NxAhHJEhEVkXjvMlV9SVUvDMX5TO1lycEYwPfLMoLaAqsiHYQxYMnBRDn3V+4JPvPTROTBUtv8UUR+dksiQ32WJ4nIYyLyk4jsEJF/iUiKu+48EckWkT+IyHbguQBiOUNEvhSR/e70DJ91N4rIBhE5KCIbvXGIyAkissDd52cRebWM4yaJSA4QB6wUkfWVvXaf+O8SkZ0isk1EbvLZNkVE/iYiP7rnXuS+9oXuJvvc6qu+pUtrlbzO+SIyQUQ+dV/rhyLSpLL3ztQ+lhxMbdcCaAJkAMOAqSLS0V33CHAi0A04wd3mL6X2bYTzi314RScRkUbAu8A/gMbA48C7ItJYROq5yweqagPgDGCFu+sE4EPgOCATeKr0sVU1X1Xru7NdVfX4Krz2dPd1/R8wSUSOc9c9BvR0Y2kE3A0UA+e46xu61VefB/o6fTYbAtwENAMSgbEBxmtqEUsOJtLeEpF93gcwuRrH+LP7BbsA54vtahERnC/8Maq6R1UPAg8D1/rsVwyMc/c9XMk5LgbWquqLqlqoqi8D3wGX+hzrFBFJUdVtquqtHirAST6tVDVPVYPZnlIAjFfVAlWdA+QAHUXEA/wGGK2qW1S1SFU/U9X8AI5Z2esEeE5Vf3Dfs9dwkq+pYyw5mEi7XFUbeh/ArVXcf6+qHvKZ/xFoBTQFUoFlPonnfXe51y5VzQvwPK3cY/v6Echwz38NMBLYJiLvishJ7jZ3AwIsEZFVIvKbKry2yuxW1UKf+VygPk5JKhlYX41jlvs6fea3l3FOU8dYcjDRLhfnS96rRan1x7nVOl5tgK3Az8Bh4GSf5JPuU30DUJUuibfilAB8tQG2AKjqB6p6AdAS55f2v93l21X1ZlVtBYwAJvu2I1Sistdenp+BPKCs6qnKXnOFr9PEDksOJtqtAIaISJyIDADOLWObB0QkUUTOBi4BZqhqMc4X9BMi0gxARDJEpH8145gDnCgiQ0QkXkSuAToD74hIcxEZ5CapfJzqnWL3nFeJSKZ7jL04X87FAZ5zBZW/dj/ua38WeFxEWrn79xWRJGCXe/72VX2dAcZs6ghLDibajcap794HDAXeKrV+O86X7lbgJWCkqn7nrvsDsA5YLCIHgLlAR6pBVXfjJJ67gN041UWXqOrPOH9Hd7ox7MH5Er/F3bU38IV7NdIsnHaADQGetrLXXpGxwDfAl25MjwIeVc0FHgI+davb+lThdZoYIjbYjzHGmNKs5GCMMcaPJQdjjDF+LDkYY4zxY8nBGGOMn2jobKzGmjRpollZWZEOwxhjapVly5b9rKpNy1pXJ5JDVlYWS5cujXQYxhhTq4hI6bvhS1i1kjHGGD+WHIwxxvix5GCMMcZPnWhzMMbUHQUFBWRnZ5OXF2iHuaYyycnJZGZmkpCQEPA+lhyMMVElOzubBg0akJWVhTMsh6kJVWX37t1kZ2fTrl27gPeL6Wqlbdvg3HNh+/bKtzXGhEdeXh6NGze2xBAkIkLjxo2rXBKL6eQwYQIsWgTjx0c6EmOML0sMwVWd9zMmk0NKCojAlClQXOxMRZzlxhhjYjQ5bNgAQ4ZAqjvGVmoqDB0KGzdGNi5jjKlf33/U1ccff5zOnTvTpUsX+vXrx48/lnvvWtDEZHJo2RLS0iAvD5KTnWlaGrQIdBBGY0xMUFWKiwMduC90unfvztKlS/n666+58sorufvuu0N+zphMDgA7dsDIkTB7NjRvDps2RToiY0w02LRpEx07duSGG27glFNOYfPmzUycOJHevXvTpUsXxo0bV7LthAkT6NixI2eddRbXXXcdjz32WEhiOv/880l1qzr69OlDdnZ2SM7jK2YvZZ0505neequTKKzfPmOizwOzV7F664GgHrNzqzTGXXpyhdusXbuW559/nj59+vDhhx+ydu1alixZgqpy2WWXsXDhQlJSUnjjjTdYuXIlBQUF9OjRg549e/od66WXXmLixIl+y0844QRef/31Ksf/zDPPMHDgwCrvV1UxmxxSUpzqJK8pU5xHcjIcPhy5uIwxkde2bVv69HGG1/7www/58MMP6d69OwA5OTmsXbuWgwcPMmjQIJKTk0lOTubSSy8t81hDhw5l6NChQYlr+vTpLF26lAULFgTleBWJ2eSwYQOMHQtvvQW5uU6j9ODBEKJSoTGmGir7hR8q9erVK3muqtx7772MGDHimG2efPLJgI4VrJLD3Llzeeihh1iwYAFJSUkB71ddMdvmYI3SxphA9O/fn2effZacnBwAtmzZws6dOznzzDOZPXs2eXl55OTk8M4775S5/9ChQ1mxYoXfoyqJYfny5YwYMYJZs2bRrFmzoLyuysRsyQGONkoPHw5Tpzp3TBtjjK8LL7yQNWvW0LdvX8C51HT69On07t2byy67jC5dutC8eXNOPfVU0tPTa3y+3NxcMjMzS+bvvPNO5syZQ05ODldddRUAbdq0YdasWTU+V0VEVUN6gnDo1auX2mA/xtQNa9asoVOnTpEOIyA5OTnUr1+f3NxczjnnHKZOnUqPHj0iHVaZynpfRWSZqvYqa/uYrVaKNOvXyZjab/jw4XTr1o0ePXrwq1/9KmoTQ3XEdLVSJPn26zR5cqSjMcZUx3//+99IhxAyVnIIs2ju18lKM8YYL0sOYRbN/TpZL7XGGK+YTg4rVkDDhvD11+E7ZzReQhvNpRljTGTEdHK4/nrYv9/5JR9O3ktoFy92ppGuxonm0owxJjJiskG69LgXq1YdXRaOK3u9/ToBTJoU+vNVJhpLM8bEquzsbEaNGsXq1aspLi7mkksuYeLEiSQmJjJ//nwGDRpUMtxnkyZNmDt3Lvfffz///ve/adq0KQADBgzgkUceqVEcMVlyWL4c2rY9dllWFqxcGZFwokK0lWaMiQbh7rJbVbniiiu4/PLLWbt2LT/88AM5OTncd999JducffbZJXdZz507t2T5mDFjSpbXNDFAjJYcunUDn65TAGe+S5eIhBMVoq00Y0ykbNq0if79+3P66aezbNky5syZw2uvvcZrr71Gfn4+gwcP5oEHHgCcLrunT59O06ZNad26NT179mTs2LHVPvfHH39McnIyN910EwBxcXE88cQTtGvXruSc4RKTyQFg7144+WT4y1+cq3P27Il0RMYYP+/dA9u/Ce4xW5wKAyv+ZR2pLrtXrVrld4y0tDTatGnDunXrAPjkk0/o1q0bAFdddVVJqeKJJ55g+vTpADz66KP0798/sPejHDGbHLZuPfr86qsjF4cxJvpEa5fd4FQrldXJ35gxY2pUaiktZpODMaYWqOQXfqhEqsvuzp07+y07cOAAP/30EyeccAJLliwJ8BXUXEw2SBtjTKDC2WV3v379yM3N5YUXXgCgqKiIu+66ixtvvLFkmNBwseRgjDEVuPDCCxkyZAh9+/bl1FNP5corr+TgwYPHdNk9cODAoHTZLSK8+eabzJgxgw4dOnDiiSeSnJzMww8/HKRXU4VYrMtuY0w0sS67Q8O67K4lrJM7Y2o/67I7AkRkAPB3IA74j6pGpmUqRKzLbmNqP+uyO8xEJA6YBAwEOgPXiUjnyEYVHNbJnTGmNojK5ACcBqxT1Q2qegR4BRgU7JNEomrHOrkzxtQG0ZocMoDNPvPZ7rISIjJcRJaKyNJdu3ZV6ySj7yxk4ULlzt9Xre+UmiQV6+TOGFMbRGtyqJSqTlXVXqray9sTYaC8VTszXokHhJene6pUtVPTQXGskztjTLSL1uSwBWjtM5/pLguK8q7ereyq3mC1F8yc6XRu17WrM/Xt9M4YE9vq16/vt2zhwoX06NGD+Pj4Mm+eC4VoTQ5fAh1EpJ2IJALXArOCdfCNG/3HdADIz6/4i97aC4yJLeHusrs8bdq0Ydq0aQwJ48hkUZkcVLUQ+B3wAbAGeE1VVwXr+C1bQkYGgLoPJ1lU9kVv7QXG1H2bNm2iY8eO3HDDDZxyyils3ryZiRMn0rt3b7p06cK4ceNKtp0wYQIdO3bkrLPO4rrrruOxxx4LSUxZWVl06dIFjyd8X9lRe5+Dqs4B5oTq+L17Q3xKEZvWxgGKqvDxx5Xv520vGD4cpk51GqeNMaHx6JJH+W7Pd0E95kmNTuIPp/2hwm0i1WV3NIna5BBq770HeXnHvvxt2yAzEwoLy9/PBsUxpu6L5i67wyVmk0N5jc9FRU4VU3IyHD4c3piMMceq7Bd+qESqy+5oEpVtDqGWkuI0Ppe3zhqZjTFe4eyyO5rEZMlhwwYYOxb++18FBKdR2rl8KT/fGpmNMUddeOGFrFmzhr59+wLOpabTp08/psvu5s2bB6XLboDc3FwyMzNL5u+8807OPvtsBg8ezN69e5k9ezbjxo1j1aqgXaNTpphMDt6rjhzFgOD9N738crspzZhYlpWVxbfffnvMstGjRzN69Gi/bceOHcv9999f0mV3WQ3SVVXepbPZ2dk1PnZVxGRyAOeqo+t/U8g8/ZwuB3vww7L6/Pijc+/CtGmRjs4YUxsMHz6c1atXk5eXx7Bhw+pUl90xPdjPln2Had0kES2K81tnDdLGREZtGuynNrHBfqogwSM0//WnNGhYWHJndCB3PdtAPcaEVl340RpNqvN+xmy1EsCqPSvJXV3IwX3OjXBxiQXkHo5n5YFPeXj1DGS1ICIIx07nPn4FKz/pw8CRX9D/rrcRBI84eda7jUc8JHgSSI1PJSUhxZnGp5CakMpxScfRNLUpTVOa0jCpIVJWXx7GxKjk5GR2795N48aN7W8jCFSV3bt3k5ycXKX9YjY5pKRAXt7pxywrOpIIFLJnyxG2/vQJxXg713Cms4cvprggqWT7FW/3ZcXbffEk5HPB1NNLtvc+CkTJRcmn/KwdD7QkgfYk0p5ETpAkukkyrUl0/jBK/jjcaaDzCSmQ2gTqNYGmHSGjF6S1rNJ7ZEwkZGZmkp2dTXW74jf+kpOTj7kCKhAxmxzKK2V5RFh/65NIcab7dev2vQRse+S3/OG13/D28r4cPpJMSmIel3f/jEev/g8tCr3XvmrJCUScaRHKYYrJBQ5TzB6K+VmK2SVF7KKYLVLERsnjMzlEgXuIxip0L47nvKJEzimKJw2PT9B6zMT3nF5ScAjJ/Rk5klOyrCijN0dO+x2FJ15Uds+Dplao+7+mhaatWle+mQEgIU5IivdvN62pmG2Q3rYNTjoJDhw4eo+DL4kros3Y9/2W7/7gFHJWtIG4YijyUL/bjzTuH6zrjYvwJP5MXOom4lI2EVdvPZ6EA6h6KDrUgYK9p1GYcxLOsNqBSSafTvITp3vWcFXcAo73bGN2UR9+XzCCPJIqP4AxJqqNPPd47hl4UrX2rahBOiZLDk6VkneudGJQzhqQy5Db99LguFMoLpU8//55ExoOPsh5gw4y/+0G7NvdlNsu9R/eunTKLSsH+2+jpeaL2ZG/jnWHPue7+Pkcqv8i9eOa0LPhYE5J60+8JFT2Ul3dgMuYX3wXe356nkvWT6Jr80Q+OvVvqCcmPwK1Vh34LWeCrEtmzW+8K0tMfjNs2ACjRsGbb/qvExFObVePWwbW818J3LDAKXVce20ar77qvZO6XQijPQEYQGFxIQs2L+DFNS+yYMe/+e7wLEZ1G8WgEwaVNIYH5iFYkkWbOWP5P887cPadoQrcGFOLxeSlrC1bQvPmZa9Thaefrvhy1ZoOE1od8Z54+rXtx3P9n+PfF/6b5vWa85fP/sJN79/E+n3rq3aw026GTpfCgkdh74+hCdgYU6vFZHIAZyyGsng8sGVL2QkgWMOE1oSI0KdlH6YPnM74M8azfv96rpx9JS+seqFq1zIPeASKi+Czp0IXrDGm1orZ5JCdDR06+C8XcUoWZSWAaBomVEQY3GEwsy6fxTkZ5zBx6URu//h29ufvD+wA6ZnQ9RpY/iLk7gltsMaYWicmk0NKCrRqBWvX+q+78EInAXhLA75deEfjMKGNkhvx5PlPcs9p97Bo6yJ+/d6vyT4YYAddp42AwjxY/VZIYzTG1D4xmRy8JQDvcKzJyU4pYuBAmDPH+cL39qt0+PCxCcA7TOjixc40GrrQEBGGdhrK1AumsvvwbobOGcqqnwO4vLbFqdDkRPjmjdAHaYypVWL2PodbbnHaHRITnTEcWrSAr76Cdu18L3M9qrZ0xLdx/0ZumXsL+/P3868L/kXXpl0r3mH+I85j7A9Qv1l4gjTGRAXreK8MviWAzp2dq5PGj4+udoXqaJfejmkDpnFc8nGM+GgEK3auqHiHDhcAChsXhiM8Y0wtEbPJYeZMePZZ6NYNvAMqTZnitEW8+mp0tStUVYt6LXiu/3M0Tm7MqHmjKr7UtWU3SE6HDf8LW3zGmOgXs8kByi8lXHhh6NsVQt3td/N6zXn6gqdJjEtk5NyR7Di0o+wNPXHQ7hzYsCA0gRhjaqWYTg7lXX00Zw5MmgRduzrTmTODf+5w3EiX2SCTyf0mcyD/ALfOu5XcgtyyN2x7JuzfDAe2hS4YY0ytEtPJAcJ/9VG4b6Tr1LgTT5z3BOv2rePPn/657BvlWrlDG279KjRBGGNqnZhPDjNnhr6U4CsSDd5nZJzBHT3u4MMfP2Taqmn+G7Q4FSQOtlhyMMY4Yj45hFukbqS78eQbuaDtBTz51ZN8se2LY1cmpkKzzrBlWWiDMMbUGpYcIiASN9KJCA+e+SBt09ryx0/+yL68fcdu0KobbFtpfUIbY4AAkoOINA5HINEmlFcThbsqyys1IZVHz36UPfl7eODzB45tf2jWGQ7vgUM2NKMxJrCSw2IRmSEiF0ndH5+wRCS65Q6HTo07cXv325n701zeWvfW0RXN3JGkdn0XkbiMMdElkORwIjAV+DWwVkQeFpETa3JSEblKRFaJSLGI9Cq17l4RWSci34tI/5qcpzqioVvuUBt28jBOa3Eaf13yVzYf2OwsbNrJme605GCMCSA5qOMjVb0OuBkYBiwRkQUi0rea5/0WuAI4ps8GEekMXAucDAwAJotI8EfO9lG6+qi2d58RCI94eOish4iTuKPVSw1aQFI67FoT6fCMMVEgoDYHERktIkuBscBtQBPgLuC/1Tmpqq5R1e/LWDUIeEVV81V1I7AOOK065whU6eqjaOyWOxRa1GvBmJ5j+GL7F7y57k2neNTsJCs5GGOAwKqVPgfSgMtV9WJVnamqhaq6FPhXkOPJADb7zGe7y/yIyHARWSoiS3ftqnojakXVR9HYLXcoXHnilfRs3pPHvnyMXbm72CY9OffBB+rs6zXGBC6Q5PAnVZ2gqiUjyIjIVQCq+mh5O4nIXBH5tozHoCDEjapOVdVeqtqradOmVd6/ouqjSF1NFG4e8XB/3/vJL8rn4S8eZsLbV7NoU0/GjyuIdGjGmAiLD2Cbe4DXSi27F5hR0U6q+stqxLMFaO0zn+kuC7pYqT6qTFZ6Ft/89gu+OnL0ozBlahxTptaeMSyMMcFXbslBRAaKyFNAhoj8w+cxDSgMUTyzgGtFJElE2gEdgCUhOldJ9dHs2dC8OWzaFKozRbd165VWZ8/Hk+iMcpSaXFjnGuGNMVVTUbXSVmApkAcs83nMAmp0iamIDBaRbKAv8K6IfACgqqtwSimrgfeBUapaVJNzVWTSJPj2W5g+3UkUWVmhOlN0a5OZwGltO1JckEhCfB55+XExWYoyxhxV6TChIhKvqqEqKQRFdYYJBYiPh6IyUk8sVqdccQWsL15EYfcn6P3F7RxIPrfOtrUYYxzVGiZURLztDMtF5OvSj5BEGibeK5XKSgyxWp0ycya8N709qW1/oNE191hiMCbGVVStNNqdXgJcWsaj1ip9pRJAXJyTMGK5OqVV/Vb8NqElH5HD4m2LIx2OMSaCyk0OquodFuxXQIGq/uj7CE94oeF7pZLHfQeuuAJuuaXu3tMQqBub9iajsJDHvpxIUXHImnuMMVEukPscGgAficgnIvI7EWke6qDCwXul0ldfwa23QmFh3b6nIVBJx7Xjjj37+H7vD8xaPyvS4RhjIqTSBumSDUW6ANfglCSyq3kfQ0hUt0HalOH799GXr+H6ruezreAA7wx+h9SE1Mr3M8bUOtVqkC7DTmA7sBtoFozATBRKz0CA37f6BbsO7yp7WFFjTJ0XSMd7t4rIfGAe0Bi4WVW7hDowEyFpTldW3YrjuLDthUxbNY0dh3ZEOChjTLgFUnJoDdyhqier6v2qujrUQZkISjkOElJhfzZ39LyDwuJCnlr+VKSjMsaEWUX3OaS5TycCP4lII99HeMIzYSfilB72Z9O6QWuGdhrKrPWz+G6PdeVtTCypqOTgHathGU43Gr5daFjrb12WngEHnP4Ob+5yM+lJ6Tz25WMEevGCMab2q+g+h0vcaTtVbe9OvY/24QvRhF1aJux3kkNaYhoju47ki+1f8MmWTyIcmDEmXAJpkJ4XyDJTh6RnQM4OKDwCwNUnXk3btLb8benfKCyO6m62jDFBUlGbQ7LbttBERI7zaW/IopzR2UwdkZYBKBx0bpJPiEtgTI8xbNi/gZlrY/wuQWNiREUlhxE47QsncWx7w9vAP0MfmomYdDf3Hzg6ztIv2vyCHs16MGnFJA4VHIpQYMaYcKmozeHvqtoOGFuqzaGrqlpyqMvSMp3p/qPJQUS4q9dd7Mnbw7PfPhuhwIwx4RLIfQ7FItLQO+NWMd0aupBMxJWUHLKPWdylaRcGZg3khVUvsP1QjPdQaEwdF0hyuFlV93lnVHUvcHPIIjKRl9QAktOPKTl4je45miIt4p/LrfBoTF0WSHKIExHxzohIHJAYupBMVEjLKGmQ9pVRP4PrO11vN8YZU8cFkhzeB14VkX4i0g942V1m6rK0Vsc0SPv6bZffkpaUxmNL7cY4Y+qqQJLDH4D/Abe4j3nA3aEMykSBtFZwYGvZqxLTuKXrLXyx7QsWbVkU5sCMMeFQaXJQ1WJVnaKqV7qPp1XVhgir6xq0gpydJTfClXb1iVfTpkEbuzHOmDoqkDukO4jI6yKyWkQ2eB/hCM5EUForQCGn7KuSEuISGNNzDOv3r+fNdW+GNzZjTMgFUq30HDAFKATOB14ApocyKBMF0ryXs5ZdtQTQr00/ujfrzqTldmOcMXVNIMkhRVXn4Qwp+qOq3g9cHNqwTMSltXKm5TRKg3Nj3NheY9mdt5vnvn0uTIEZY8IhkOSQLyIeYK2I/E5EBgP1QxyXibSS5FB+yQGcG+MGZA3g+VXP24hxxtQhgSSH0UAqcDvQE7geGBbKoEwUSE6HhHpwwP9eh9JG93BvjFthN8YZU1cEcrXSl6qao6rZqnqTqv5KVReHIzgTQSIV3uvgK7NBJkNOGsLb6962G+OMqSMCKTmYWFXBvQ6l3dzlZhomNeThLx62G+OMqQMsOZjypWUEnBzSk9IZ03MMy3cuZ/aG2SEOzBgTaoHc53BmIMuqQkQmish3IvK1iLxZqtfXe0VknYh8LyL9a3IeU0NprZz+lYoDu+dx0AmD6NKkC39b+jcOHDkQ4uCMMaEUSMnhqQCXVcVHwCmq2gX4AbgXQEQ6A9cCJwMDgMluR38mEtJaghY5d0oHwCMe7utzH3vz9jJp+aQQB2eMCaX48laISF/gDKCpiNzpsyoNqNEXtqp+6DO7GLjSfT4IeEVV84GNIrIOOA34vCbnM9XkeyNcWsuAduncuDNXd7yaV75/hSs6XEHHRh1DGKAxJlQqKjkk4tzPEA808Hkc4OiXeTD8BnjPfZ4BbPZZl00541WLyHARWSoiS3ft2hXEcEyJAG6EK8tt3W8jPTGdBxc/SLEWhyAwY0yolVtyUNUFwAIRmaaqP1b1wCIyF2hRxqr7VPVtd5v7cLrleKmqx1fVqcBUgF69etnlMaEQQBcaZUlPSufOXnfy50//zIzvZ3DNSdeEIDhjTCiVmxx8JInIVCDLd3tV/UVFO6nqLytaLyI3ApcA/fTotY9bgNY+m2W6y0wkpDaGuEQ4WLXkADDo+EHM2TCHx5c9zjmZ59CyfmDVUsaY6BBIg/QMYDnwJ+D3Po9qE5EBOGNCXKaquT6rZgHXikiSiLQDOgBLanIuUwMlN8JVPTmICOPOGIeijF883u59MKaWCSQ5FLrjOSxR1WXeRw3P+0+c9ouPRGSFiPwLQFVXAa8Bq3FGmxtlY0dEWBXudSgto34Go3uMZtGWRbyz4Z0gB2aMCaVAqpVmi8itwJtAvnehqu6p7klV9YQK1j0EPFTdY5sgS2sF2V9We/drO17Lexvf49EvH6Vvq740SWkSxOCMMaESSMlhGE410mfAMvexNJRBmSjirVaqZrVQnCeO8WeM53DBYe7/7H6rXjKmlgik4712ZTzahyM4E17btsG558J238Hf0jKg6Ajk7q72cds3bM+YnmNYkL2AGT/MqHmgxpiQC6T7jFQR+ZN7xZJ32NBLQh+aCbcJE2DRIhg/3mdhA/cqoyre61DakE5DOKPVGUz8ciIb9tsos8ZEu0CHCT2Cc7c0OJeWPhiyiEzYpaQ4FyZNmQLFxc5UxFle3XsdSvOIhwfPfJDk+GTuWXgPBUUFNQ/cGBMygSSH41X1/wEFAO6lpxLSqExYbdgAQ4ZAaqozn5oKQ4fCxo1U+y7psjRNbcr9Z9zPmj1reHzZ4zU+njEmdAJJDkdEJAVQABE5Hp+rlkzt17IlpKVBXh4kJzvTtDRo0QKo3wwkLqAR4QLRr00/hnYayvQ103lv43uV72CMiYhAksM4nHsOWovIS8A8nBvYTB2yYweMHAmLFzvTkkZpT5zT7lDDaiVfd/W8i25NuzHus3Gs37c+aMeNRmU28htTC0gglxaKSGOgD0510mJV/TnUgVVFr169dOlSu7o2ZP5zASSkwLBZQTvkztydXD37ahokNuDli1+mfmL9oB07mtx6Kzz9NIwYAZMnRzoaY44lIstUtVdZ6wK5Wmkwzl3S76rqO0ChiFwe5BhNNAtwLOmqaJbajInnTmTzwc3cvfBuCosLg3r8SKuwkd+YWiCgaiVV3e+dUdV9OFVNJlakZ8L+bOdbLoh6t+jNH0//I59s+YRHljxSp26Qq7CR35haIJDuM8pKIIHsZ+qKRu2gMA9yth+9eilIru54NdkHs3lu1XO0btCaYScPC+rxI6XCRn5jaoFASg5LReRxETnefTyO04WGiRXHtXOme0Lzs/eOnndwQdsL+NvSv/HBpg9Cco5IKLeR35haIJASwG3An4FXcS5n/QgYFcqgTJRp5CaHvRsh68ygH94jHh4+62F+Pvwz9yy8h5T4FM7JPCfo5wm3mTOPPp9kQ2qbWqbCkoOIxAHvqOo9qtpLVXur6h9V9VCY4jPRIL21c69DiEoOAMnxyUzqN4mOjToy5n9j+HyrDRtuTCRVmBzcsRSKRSQ9TPGYaBSX4DRK7w1ta2qDxAY8fcHTZKVncfvHt1uCMCaCAmlzyAG+EZFnROQf3keoAzNRplE72Lsp5KdJT0pn6gVTaZ3WmlHzRjH3x7khP6cxxl8gyWEmTpvDQo6O52AN0rHmuHYhrVby1TilMc/1f47OjTtz14K7eHPtm2E5rzHmqEobpFX1ebdvpTaq+n0YYjLR6LgsOLwHDu+DlIYhP523BDFm/hj+8tlf2H5oOyO6jsAjgfyeMcbUVCB3SF8KrMDpXwkR6SYiwetHwdQOTTo4093rwnbK1IRUnvrFU1x2/GVMXjmZsQvGkluQG7bzGxPLAvkZdj9wGrAPQFVXADYSXKxpepIz3bk6rKdNjEvkwTMfZGyvscz7aR7D3h9G9sHssMZgTCwKJDkU+Haf4QpuPwom+h2XBfEpsPO7sJ9aRBh28jD++Yt/suXgFq6afZV1921MiAWSHFaJyBAgzh0i9CngsxDHZaKNJw6angi71kQshLMzz2bGZTM4vuHx3L3wbv786Z85VGC33BgTCoEkh9uAk3EG+PkvsB+4I4QxmWjVtFNESg6+MupnMG3ANIZ3Gc7b697m8rcvZ2H2wojGZExdVG5yEJFkEbkD+H/AT0Bf9w7pP6lqXrgCNFGkWSc4uNW5YimC4j3x3Nb9Nl4Y+AL1E+ozat4o7l5wN7tyd0U0LmPqkopKDs8DvYBvgIHAY2GJyESvZp2c6a7Ilh68ujXrxmuXvMaobqOY+9NcLn7zYqasmGJXNBkTBBUlh86qer2qPg1cCdT+ntBMzTQ/2Zlu+zqycfhIiEtgZNeRvD3obc7KOIvJKydzyZuXMOOHGRQUFUQ6PGNqrYqSQ8lflqrWrWG6TPWkZUD95rAl+oZkbZ3WmsfPe5wXB75Iy/otGf/5eAbOHMhLa14ir9BqQY2pqnLHkBaRIsB7KYgAKUCu+1xVNS0sEQbAxpAOo5eHONVKt38V6UjKpap8tvUzpn49la92fkWj5Eb8qsOvuOrEq2hZv2WkwzMmalQ0hnS53WeoalzoQjK1VmZP+P5dyN0DqY0iHU2ZRIQzM87kzIwzWbZjGdNWTeOZb5/hmW+f4ZzMc7jqxKvo26ovCZ6ESIdqTNSKyHCfIjIBGIRzM91O4EZV3SoiAvwduAinlHKjqkbvT9RYlNnbmW5eAh0HRDaWAPRs3pOezXuyNWcrr//wOm+sfYP5m+fTMKkh/bP6c1G7i+jWrJv12WRMKeVWK4X0pCJpqnrAfX47TuP3SBG5COe+iouA04G/q+rplR3PqpXCqCAPHm0LvX4DA/4a6WiqrKCogE+3fsq7G95l/ub55BXl0TSlKedknsPZmWfTt2VfUhNSIx2mMWFRrWqlUPImBlc9nOFHwSlNvKBOxlosIg1FpKWqbgt7kKZsCcnQpi+s/1+kI6mWhLgEzmt9Hue1Po9DBYf4+KeP+d/m//HBpg94Y+0bJHgS6Nm8J71b9KZX816c0uQUEuMSIx22MWEXkeQAICIPATfg3HF9vrs4A9jss1m2u8wvOYjIcGA4QJs2bUIaqyml/Xkwdxwc2ApprSIdTbXVS6jHpcdfyqXHX0pBUQHLdy5nQfYCPtv6GU8tfwqApLgkujTtQpcmXejcuDOdG3cmo34GTg2oMXVXyKqVRGQu0KKMVfep6ts+290LJKvqOBF5B3hEVRe56+YBf1DVCuuMrFopzHZ+B5NPh4H/D04fEeloQmJf3j6W7VzGsh3LWLp9KWv3rqXQvaI7PSmdTo060fG4jrRv2J526e1on96e9CT/0XS3bYNrr4VXX4UWZf01GFMDNf18VVStFJE2h2MCEGkDzFHVU0TkaWC+qr7srvseOK+yaiVLDhEwuS8kp8Nv3o90JGGRX5TP2r1rWb17dclj/b71HCk+UrJNo+RGtEtvR5sGbWhZvyUZ9TN4YUIfZr7YlOHDlX9NsUZvE1y33gpPPw0jRsDkyVXfP+qSg4h0UNW17vPbgHNV9UoRuRj4HUcbpP+hqqdVdjxLDhGwcCJ8/CDcvhwaxebwHkXFRWw9tJWN+zeyYd8GNh5wpltytjD/+vfQgiS/fTwJR7jtg7tpktKEximNaZzcmCYpTWiS0oRGyY1IT0qnXkI9u3rKVCglBfLKuLczORkOHw78ONGYHN4AOuJcyvojMFJVt7iXsv4TGIBzKetNlVUpgSWHiDiwDZ48BU4bXiuvWgq1Tdn5jB5TwAfvppB/OI6E5AI6nP0tXX/zIodTNrL78G725u8tc1+PeEhLTCMtMY30pHTneVJaybLUhFRS41NJTUilXkK9kuep8c4jJSGFegn1SI5LtraROmrbNhg7Ft56C3JzITUVBg+Gxx6rWvVSNF6t9KtyliswKszhmOpIawknD4avXoSz74J6TSIdUVTJykyiVZMkCvKdX3NHjiRw7gndmXxt95JtCooL2HN4Dz/n/czuw7vZfXg3B44cYH/+fg4cOcCB/APO9MgBNh/cXPK8WAMba0sQkuKSSIxL9Jses8yT5Lc+wZNAvCeeeE88cRJX8jzeE0+8HH0e53HWJUhCyfPS23jEQ5zEISJ48ODxePBwdFnJOvGUPI5ZhxDnicOD55jtY1nLlpCW5pQekpOdaVpacNu1Ina1kqkDzh4L386E+X+Fi/8W6Wiizo4dMHIkDB8OU6c6v/Z8JXgSaF6vOc3rNQ/4mKrKkeIjHCo4RG5BLrmFuc7U+7zw6PNDBYcoKCogvyif/KJ8jhQdOWaaX5TP/oL9HCk+csyy/KJ8CosLKSwupEiLgvyuBI9vAvGIB0FKps7/gohQ8p+bUEovD3hbdwqUvW2pZeVte8w5S53LV+nlvglREOZ9fSsnXrSf/7u5iI0fXOD3+aqpiDdIB4NVK0XQnN/Dkn/DDW9D+3MjHY0JMlWlUAuPJoviomPmSx7qrCsoLihJKr4JRlUpprjkeVnTYi0+ugylqNidll7ns32xFlPMsesU9Zt6X0vJf+73XqDbKorzfwXb+uxT1rbHLPdd5gRSEs8xUz123neZ18XtL+bqjldX69836qqVTB3yy/thw3yYcSMMmw0tTolwQCaYRIQESbB+qGKQXRJhaiaxHlz3CsQnw/OXwHfvRjoiY0wQWHIwNdf4eLjpXUhvDa8MgRcGwZrZUFCFa+qMMVHFqpVMcDRqD7+dB0uehk//Aa9eD3FJTjVTs86Qnul0tZGUBon1Iak+JKSCJx48cSBx4PE4U/G4yzxAJVelBHTVSgDbBOs4xoRbQrJTgg8ySw4meOIT4Yzb4PRbYNMnsH4ebFkOP3wAh3ZGOjpj6qYz74ALHgj6YS05mOCLi4fjz3ceXoVHIGc75OfAEe8jF4oLQYuguNidFvlMK7ueP4Ar7YJ1NV4duKrP1FEtu4bksJYcTHjEJ0LD2Os91zreM7WVNUgbE0ITJsCiRTB+fKQjMaZqLDkYEwIpKU4b95QpTo3ZlCnOfEpKpCMzJjCWHIwJgQ0bYMgQp0M0cKZDh8LGjZGNy5hAWXIwJgTC0TGaMaFkycGYEPF2vLd4sTPdvj3SERkTOLtayZgQmTnz6PNJkyIXhzHVYSUHY4wxfiw5GGOM8WPJwRhjaqlt2+Dcc0PTnmXJwRhjaqlQ3mRpycEYY2qZcNxkacnBGGNqmXDcZGnJwRhjaplw3GRpycEYY2qhUN9kaTfBGWNMLRTqmyyt5GCMMcaPJQdjjDF+LDkYY4zxY8nBGGOMH0sOxhhj/FhyMMYY4yeiyUFE7hIRFZEm7ryIyD9EZJ2IfC0iPSIZnzHGxKqIJQcRaQ1cCPzks3gg0MF9DAemRCA0Y4yJeZEsOTwB3A2oz7JBwAvqWAw0FJGWEYnOGGNiWESSg4gMArao6spSqzKAzT7z2e6yso4xXESWisjSXbt2hShSY4yJTSHrPkNE5gJldQN1H/BHnCqlalPVqcBUgF69emklmxtjjKmCkCUHVf1lWctF5FSgHbBSRAAyga9E5DRgC9DaZ/NMd5kxxpgwCnu1kqp+o6rNVDVLVbNwqo56qOp2YBZwg3vVUh9gv6puC3eMxhgT66KtV9Y5wEXAOiAXuCmy4RhjTGyKeHJwSw/e5wqMilw0xhhjwO6QNsYYUwZLDsYYY/xYcjDGGOPHkoMxxhg/lhyMMcb4seRgjDHGjyUHY4wxfiw5GGNMLbVtG5x7LmzfHvxjW3IwxphaasIEWLQIxo8P/rEtORhjTC2TkgIiMGUKFBc7UxFnebBYcjDGmFpmwwYYMgRSU5351FQYOhQ2bgzeOSw5GGNMLdOyJaSlQV4eJCc707Q0aFHWCDrVZMnBGGNqoR07YORIWLzYmQa7UTrivbIaY4ypupkzjz6fNCn4x7eSgzHGGD+WHIwxxvix5GCMMcaPJQdjjDF+LDkYY4zxY8nBGGOMH1HVSMdQYyKyC/ixmrs3AX4OYji1mb0XDnsfHPY+HFVX34u2qtq0rBV1IjnUhIgsVdVekY4jGth74bD3wWHvw1Gx+F5YtZIxxhg/lhyMMcb4seQAUyMdQBSx98Jh74PD3oejYu69iPk2B2OMMf6s5GCMMcaPJQdjjDF+Yjo5iMgAEfleRNaJyD2RjidcRKS1iPxPRFaLyCoRGe0ubyQiH4nIWnd6XKRjDQcRiROR5SLyjjvfTkS+cD8Xr4pIYqRjDAcRaSgir4vIdyKyRkT6xuJnQkTGuH8X34rIyyKSHIufiZhNDiISB0wCBgKdgetEpHNkowqbQuAuVe0M9AFGua/9HmCeqnYA5rnzsWA0sMZn/lHgCVU9AdgL/F9Eogq/vwPvq+pJQFec9ySmPhMikgHcDvRS1VOAOOBaYvAzEbPJATgNWKeqG1T1CPAKMCjCMYWFqm5T1a/c5wdxvgQycF7/8+5mzwOXRyTAMBKRTOBi4D/uvAC/AF53N4mV9yEdOAd4BkBVj6jqPmLwM4EzCFqKiMQDqcA2YvAzEcvJIQPY7DOf7S6LKSKSBXQHvgCaq+o2d9V2oHmk4gqjJ4G7gWJ3vjGwT1UL3flY+Vy0A3YBz7lVbP8RkXrE2GdCVbcAjwE/4SSF/cAyYvAzEcvJIeaJSH3gDeAOVT3gu06da5zr9HXOInIJsFNVl0U6ligQD/QApqhqd+AQpaqQYuQzcRxOaakd0AqoBwyIaFAREsvJYQvQ2mc+010WE0QkAScxvKSq3tFod4hIS3d9S2BnpOILkzOBy0RkE0614i9w6t0bulUKEDufi2wgW1W/cOdfx0kWsfaZ+CWwUVV3qWoBMBPncxJzn4lYTg5fAh3cqxAScRqdZkU4prBw69WfAdao6uM+q2YBw9znw4C3wx1bOKnqvaqaqapZOP/+H6vqUOB/wJXuZnX+fQBQ1e3AZhHp6C7qB6wmxj4TONVJfUQk1f078b4PMfeZiOk7pEXkIpw65zjgWVV9KLIRhYeInAV8AnzD0br2P+K0O7wGtMHpAv1qVd0TkSDDTETOA8aq6iUi0h6nJNEIWA5cr6r5EQwvLESkG07DfCKwAbgJ5wdkTH0mROQB4Bqcq/qWA7/FaWOIqc9ETCcHY4wxZYvlaiVjjDHlsORgjDHGjyUHY4wxfiw5GGOM8WPJwRhjjB9LDsYYY/xYcjDGGOPHkoOpFhFpISKviMh6EVkmInNE5MRIx1VVInK/iIyNdByVccdauDUY+4nIZ0GKKUVEFrjd31e27QgRUfdmQ++yUe6y/iKy0Kd7ChMFLDmYKnO7FXgTmK+qx6tqT+Be6niPnRHWECgzOYijvL9lv/1U9YwgxfQbYKaqFgWw7anASuAkABFJxbnzeBewAmesiGuCFJcJAksOpjrOBwpU9V/eBaq6UlU/EZE73RG0vhWRO7zrRSTLHWFsmoj8ICIvicgvReRTd5Sx03y2eckdiex190vEewy/Y7v7fOuzzVi3NJDlHuPf7qheH4pIirvNfW4MiwBvX0J+ROQGEflaRFaKyIsBxFHe+fyOIyLXi8gSEVkhIk97f31XcJxHgOPd7Se6230vIi8A3wKtReQttxS3SkSGu+Ees597jpzqvpZShuLTx5CIzBcR75d/Y99/F6ALTvcTJ7nztwMzgGJV3QG85R7PRAtVtYc9qvTA+cN+oozlPXH6a6oH1AdWAd3ddVk4fdWcivOjZBnwLCA4XSS/5W6jwJnuPs/i9HdU7rHdfb71iWEscL/P+bq5y18Drvc5TiqQBqzznqPUazkZ+AFo4s43CiCOss7ndxygEzAbSHCXTQZuKPU+lT5O6deZhdMvVh+fZd4YU3ASRuPS+7nrc6rzWkodIxHYXmpZNuBxn58PvOyzbhfQHngfpzSzHDgPmOuujwN2RfqzbY+jDys5mGA6C3hTVQ+pag5Od8dn+6zfqKrfqGoxzhfRPHW+Gb7B+UIC2Kyqn7rPp7vHDOTYZdmoqivc58vcc5ztHidXnTEsyuuJ9xfADFX9GUCPdjZXURxlna+s4/TD+WL+UkRWuPPtK4m7LD+q6mKf+dtFZCWwGKc7+g7l7OdV1dfiqwmwzzsjIm2BLe6/LTglha/dda2B3aq6AWgG/B54CjgR598edaqmjohIg0piNmFiDUCmOlZxtPviqvDtxbLYZ76Yo5/F0j1BVtYzZCHHVo8ml3O+Ipxf1KEU6PkEeF5V763hcQ6VHNBp6P0l0FdVc0VkPse+F1VVWQyHSx2/K24ycPUEXnWfn4qbBICDOIPnnIbTI/JXPvskAXk1iNkEkZUcTHV8DCT51GsjIl1wGhYvF6cv/HrAYJyuwauijYj0dZ8PARa5zz8p59g7gGZuHXcScEklx1/oHifF/ZV6aQWv8SoRaey+vkaVxFGeso4zD7hSRJp5l7m/vCtyEKjoV3U6sNdNDCcBfQLYr6qvpYSq7gXiRMSbILrhJgsR6YBTVehNCF18nk8EfueWFEqShvv+/KzOADsmClhyMFXmVgUNBn4pzqWsq4C/AluBacASnLEh/qOqy6t4+O+BUSKyBjgOmOKe86uyju1+mYx3l38EfFdJ7F/h/KJdCbyHM+hTWdutAh4CFrhVNY9XFEcF5/M7jqquBv4EfCgiX7txt6wk7t3Ap27D8cQyNnkfiHfft0dwqpYq3K+qr6UMH3K02q8r4HFf419wBsjxDhJ0Kk4bCKr6jqp+7i7vjFMKBaeN4t0qnNuEmI3nYKKGiGQB76jqKZGOxVRORHoAY1T11yKyFuihqgereayZwD2q+kNQgzTVZiUHY0y1uCWP/4lIujNb7cSQCLxliSG6WMnBGGOMHys5GGOM8WPJwRhjjB9LDsYYY/xYcjDGGOPHkoMxxhg/lhyMMcb4seRgjDHGz/8HndP3QUt3JaAAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# set loss function to huber\n", + "opt.handle_opt_set (handle, 'NLDF Loss Function Type = HUBER')\n", + "\n", + "# set initial guess and solve\n", + "x = [30., 30., 30., 30.]\n", + "\n", + "# regularisation types\n", + "regs = {'OFF', 'L1', 'L2'}\n", + "\n", + "# solve huber loss function with various regularisations\n", + "for reg in regs:\n", + " # set reg type\n", + " opt.handle_opt_set(handle, 'Reg Term Type =' + reg)\n", + " # call the solver\n", + " sol = opt.handle_solve_nldf(handle, lsqfun, lsqgrd, x, nres,data=data, io_manager=iom)\n", + " # calculate response\n", + " resp = sol.x[0] + ((sol.x[1] - sol.x[0]) / (1 + (sol.x[2] / dose)**sol.x[3]))\n", + " # plot curve\n", + " plt.plot(dose, resp, label='reg = '+reg)\n", + "\n", + "# show huber curves\n", + "plt.title(\"Huber loss function\")\n", + "plt.plot(t, y, '*b')\n", + "plt.xlabel(r\"Compound concentration ($\\mu M$)\")\n", + "plt.ylabel(\"Percent activity\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "f86a8a83", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEcCAYAAAAsv3j+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABBDUlEQVR4nO3deXhU1fnA8e872SYLCSSERQIEZBMB2URwVxDEDVGpC9atilvVqmi1WrXya6vFirVVBKtVCy641gV3i4gWEWSRHWQNhH1LyJ55f3/cGwhkYZLMkmTez/PcZ+bemXvOO5PJvHPPufccUVWMMcaY8jzhDsAYY0z9Y8nBGGNMBZYcjDHGVGDJwRhjTAWWHIwxxlRgycEYY0wFlhyMMcZUYMnBhJyIbBSRPuGOozZEZJ2IDAlRXV1FZIGI5IjI7aGos1zdS0Tk9FDWaeqX6HAHYCKLiDQDWgFLwx1LA3Av8F9V7R3MSkRkHXC9qn5Rtk1Vjw1mnab+syMHE2o9gZWqWhjqikWkof0Yag8sCXcQJjJZcjCh1gv4qSY7uE05Y0VkkYjsFZE3RMTrPnaMiMwQkT1uU8gFlez7WxFZBOwXkWh32z1ueftF5AURaSkiH7tNOF+4Rzj+xFZl/W69m9wyV4jI4Oq2H1buV8AZwD9EJFdEuoiIikincs95SUT+z5/3yX28rYi8IyLbRWSniPxDRP4NtAM+cOu5t1xZQ2rwHldZr2mYLDmYUOsJLKrFfr8AzgY64CSYa0QkBvgA+AxoAdwGTBWRroftezlwLtBUVUvcbRcDZwFdgPOBj4HfAek4/xdHbOOvrn43hl8Dx6tqE2AYsK6q7YeXrapnAt8Av1bVJFVdecR3yFHhfXJjjQI+BNYDmUAb4HVV/SWwATjfrecv/r5Gf+o1DZclBxNqvYCfRGSAiPxPRGaKyGvul1B1nlbVzaq6C+fLqjcwEEgCHlPVIlX9CucL8PJK9t2oqvnltv1dVbeq6iacL+HvVXW+qhYA7wL+dJhXV38pEAd0F5EYVV2nqj9Xsz1QKnufAAYARwH3qOp+VS1Q1Vl1fI3+1GsaKEsOJmRERIAeOEcOG4EzVfVUnF/OI9znjHabN3JF5ONyu28pdz8P5wvrKGCjqvrKPbYe51dxeRsrCWdrufv5lawn+fGSqqxfVVcDvwEeAbaJyOsiclRV2/2oy1+VvU8AbYH15Y6c/OXve1xVvaaBsuRgQqkDUKqq61U1u9wv+SLAB6CqU93mjSRVHX6E8jYDbUWk/Oe4HbDpsOcFa1z6autX1VdV9WScjmUFHq9uux/ygIRy661qEOtGoF0VnfLVvT/+vsemkbHkYEKpF7C4/AYRaQ8MxWmKqKnvcb4w7xWRGPe8/POB1+sWZt3rd/sdzhSROKAA52jEV9V2P+tbAFwhIlEicjZwWg1inQNkA4+JSKKIeEXkJPexrUDHmr7GGtRtGiBLDiaUDumMFpFk4N/ANapaXNPCVLUI54tqOLADeBa4SlWXBybcOtUfBzzmbt+C05l7fzXb/XGHW98eYDTwXg1iLXX37YTTAZ0FXOo+/GfgQfdspLE1eI2mERObCc6Eg9u88T7wV1X9MtzxGGMOZUcOJlwuB04Afu+eQ3/pkXYwxoSOHTkYY4ypwI4cjDHGVGDJwRhjTAWWHIwxxlTQ0EaprFTz5s01MzMz3GEYY0yDMm/evB2qml7ZY40iOWRmZjJ37txwh2GMMQ2KiKyv6jFrVjLGGFOBJQdjjDEVWHIwxhhTQaPoczDGNB7FxcVkZWVRUFAQ7lAaDa/XS0ZGBjExR5o25SBLDsaYeiUrK4smTZqQmZmJMwWIqQtVZefOnWRlZdGhQwe/94voZqXsbDjtNNiy5cjPNcaERkFBAWlpaZYYAkRESEtLq/GRWEQnh3HjYNYsePTRcEdijCnPEkNg1eb9jMjkEB8PIjBxIvh8zq2Is90YY0yEJoc1a+CKKyDBnXAxIQFGj4a1a8MblzHGJCVVnH77ySefpHv37vTq1YvBgwezfn2V164FTEQmh9atITkZCgrA63Vuk5OhVU1m5DXGNHqqis/n7yyuwdOnTx/mzp3LokWLuOSSS7j33nuDXmdEJgeArVvhppvggw+gZUtYty7cERlj6oN169bRtWtXrrrqKnr06MHGjRsZP348xx9/PL169eLhhx8+8Nxx48bRtWtXTj75ZC6//HKeeOKJoMR0xhlnkOA2dQwcOJCsrKyg1FNexJ7K+s47zu0ttziJwsbtM6b++cMHS1i6eV9Ay+x+VDIPn39stc9ZtWoVL7/8MgMHDuSzzz5j1apVzJkzB1XlggsuYObMmcTHx/P222+zcOFCiouL6du3L/369atQ1tSpUxk/fnyF7Z06deKtt96qcfwvvPACw4cPr/F+NRWxySE+3mlOKjNxorN4vZCfH764jDHh1759ewYOHAjAZ599xmeffUafPn0AyM3NZdWqVeTk5DBixAi8Xi9er5fzzz+/0rJGjx7N6NGjAxLXlClTmDt3Ll9//XVAyqtOxCaHNWtg7Fh47z3Iy3M6pUeOhCAdFRpjauFIv/CDJTEx8cB9VeX+++/nxhtvPOQ5Tz31lF9lBerI4YsvvuCPf/wjX3/9NXFxcX7vV1sR2+dgndLGGH8MGzaMF198kdzcXAA2bdrEtm3bOOmkk/jggw8oKCggNzeXDz/8sNL9R48ezYIFCyosNUkM8+fP58Ybb+T999+nRYsWAXldRxKxRw5wsFN6zBiYPNm5YtoYY8obOnQoy5YtY9CgQYBzqumUKVM4/vjjueCCC+jVqxctW7akZ8+epKSk1Lm+vLw8MjIyDqzfddddTJ8+ndzcXEaNGgVAu3bteP/99+tcV3VEVYNaQSj0799fbbIfYxqHZcuWccwxx4Q7DL/k5uaSlJREXl4ep556KpMnT6Zv377hDqtSlb2vIjJPVftX9vyIbVYKNxvXyZiGb8yYMfTu3Zu+ffty8cUX19vEUBsR3awUTuXHdXr22XBHY4ypjVdffTXcIQSNHTmEWH0e18mOZowxZSw5hFh9HtfJRqk1xpSJ6OSwYAE0bQqLFoWuzvp4Cm19PpoxxoRHRCeHK6+EvXudX/KhVHYK7ezZzm24m3Hq89GMMSY8IrJD+vB5L5YsObgtFGf2lo3rBPDMM8Gv70jq49GMMZEqKyuLW2+9laVLl+Lz+TjvvPMYP348sbGxzJgxgxEjRhyY7rN58+Z88cUXPPLIIzz//POkp6cDcPbZZ/PYY4/VKY6IPHKYPx/atz90W2YmLFwYlnDqhfp2NGNMfRDqIbtVlYsuuogLL7yQVatWsXLlSnJzc3nggQcOPOeUU045cJX1F198cWD7nXfeeWB7XRMDROiRQ+/eUG7oFMBZ79UrLOHUC/XtaMaYcFm3bh3Dhg3jhBNOYN68eUyfPp1p06Yxbdo0CgsLGTlyJH/4wx8AZ8juKVOmkJ6eTtu2benXrx9jx46tdd1fffUVXq+Xa6+9FoCoqCgmTJhAhw4dDtQZKhGZHAB274Zjj4WHHnLOztm1K9wRGWMq+Pg+2PJTYMts1ROGV//LOlxDdi9ZsqRCGcnJybRr147Vq1cD8M0339C7d28ARo0adeCoYsKECUyZMgWAxx9/nGHDhvn3flQhYpPD5s0H7//iF+GLwxhT/9TXIbvBaVaqbJC/O++8s05HLYeL2ORgjGkAjvALP1jCNWR39+7dK2zbt28fGzZsoFOnTsyZM8fPV1B3EdkhbYwx/grlkN2DBw8mLy+PV155BYDS0lLuvvturrnmmgPThIaKJQdjjKnG0KFDueKKKxg0aBA9e/bkkksuIScn55Ahu4cPHx6QIbtFhHfffZc333yTzp0706VLF7xeL3/6058C9GpqEIsN2W2MqU9syO7gsCG7Gwgb5M6Yhs+G7A4DETkb+BsQBfxTVcPTMxUkNmS3MQ2fDdkdYiISBTwDDAe6A5eLSPfwRhUYNsidMaYhqJfJARgArFbVNapaBLwOjAh0JeFo2rFB7owxDUF9TQ5tgI3l1rPcbQeIyBgRmSsic7dv316rSu67D2bOdG5roi5JxQa5M8Y0BPU1ORyRqk5W1f6q2r9sJEJ/lTXtuKcS8/LLNWvaqeukODbInTGmvquvyWET0Lbceoa7LSCqOnv3SGf1Bqq/4J13nMHtjjvOuS0/6J0xJrIlJSVV2DZz5kz69u1LdHR0pRfPBUN9TQ4/AJ1FpIOIxAKXAe8HqvC1ayvO6QBQWFj9F731FxgTWUI9ZHdV2rVrx0svvcQVIZyZrF4mB1UtAX4NfAosA6ap6pJAld+6NbRpc+g2kSN/0Vt/gTGN37p16+jatStXXXUVPXr0YOPGjYwfP57jjz+eXr168fDDDx947rhx4+jatSsnn3wyl19+OU888URQYsrMzKRXr154PKH7yq631zmo6nRgerDKP/54SElxZoFz6oOvvjryfmX9BWPGwOTJTue0MSY4Hp/zOMt3LQ9omd1Su/HbAb+t9jnhGrK7Pqm3ySHYPv7Y+eVfXnY2ZGRASUnV+9mkOMY0fvV5yO5QidjkUFXnc2mp08Tk9UJ+fmhjMsYc6ki/8IMlXEN21yf1ss8h2OLjnc7nqh6zTmZjTJlQDtldn0TkkcOaNTB2LFQ2LEphoXUyG2MOGjp0KMuWLWPQoEGAc6rplClTDhmyu2XLlgEZshsgLy+PjIyMA+t33XUXp5xyCiNHjmT37t188MEHPPzwwyxZErBzdCoVkcmh7KwjgKgopymp7G964YV2UZoxkSwzM5PFixcfsu2OO+7gjjvuqPDcsWPH8sgjjxwYsruyDumaqurU2aysrDqXXRMRmRzAOevootEF7O+wkp+/aUH24qbs3xXHd+u3MeiX67nmXyA4k294BEAQcbZ5xL0vIGXbRdzHnPveGA+JsdEkxkWTGBdFYlw0qQmxtEzx0jrFS3pSHNFREdmqZ0yjMWbMGJYuXUpBQQFXX321DdndGLzzDrw18y0uHTwKX0ncge2rZrRk1YyWeKILufzpCxA8oAACKiAeUEE1FsWL+ryUajylmkCJJpKvaeSWtKCgqAn5RZBXVFpp/VEeITMtgW6tkzmmVRP6Z6bSp11T4qKjQvMGGGPqrDEP2R2xyQGgRfRmBt5/GT9MeJ7SwiR8xV48sfk07/M5HS99gg1xO/ABPgQV8OHkiYI96fz83BO0vWUsMU13Vll+U4VOEkur6CYc5W1Jmrct3rhOFMf2IjevGau27WdR1h4+WuRcLBEfE8UJHVO54LijGHZsKxLjIvrPYyKYqiKVDWNgaqU2M35G9LfPqSfeS6dJ8J1zEgJeLxQVebm4z3k8e/twUF+lyy23e1my2sspP77Jnx/6ifyCXeQV7iWvYA978rayc/9WdhbsZEfRPnYW55JVtI1FRdvZn7sU56JvaKpC95hmjDq2E11an0JR1CnMW5fPl8u3cde0hcTHLGZk3zbceGpH2qclVv0ijGlkvF4vO3fuJC0tzRJEAKgqO3fuxOv11mi/iJ1DOj6+4kVwAB4PjBhR+WB4Ve1zxGsiVNG8XezauogN2fNYtW0BS/b9zNLivayO9lAiQrQqPaKSOCGtJ0elj+R/647iPwuyKfH5GNkng9+e3ZUWyTX74xrTEBUXF5OVlUVBZf9spla8Xi8ZGRnExMQcsr26OaQjNjl4vZVf6+DxOGcvVSY72zkF9r33IC/PGXhv5Eh44olanvqqSsGOFSxa/jazs77h+/0bWRwt+ETI0ChOS+5Ofuk5vLYgjZgo4a6hXbn2xEw8Hvs1ZYypu+qSQ8Q2K61dC926wb59h273+aq+QjrgA++J4E3vxoD0BxgAoMqerO/578IX+GzbXN7Yt4gS+Yk+neNoXtiPxz4azIwV2/jrqOPsKMIYE1QReS5lfDwcdVTFxFCmuiukgzpRjwhN2w5k5HnPM/G6+cwYNoX7mvWjqLSAmd7vSOvyKL6c3/OLv7/B/A27A1ixMcYcKiKblbKz4dZb4d13Kz4m4nzpP/ts9ftfdhm88UZorqTW0hIWLniRqUte4jPfPqIVjslJZWTf33PxyWcFPwBjTKNUXbNSRB45tG4NLVtW/pgqTJpU/TzRdZ0mtKYkKpre/cYw/qrv+ODUpzknpiVLk3fxp1V38vuXz2XfnvWhCcQYEzEi8sgBnGEzKrtK3eOBTZucL/5Jk+DGGw8eRdT6bKUgWLdhDuOm38kP3r2k+JSbWwzi0qFPExWbENpAjDENlh05VCIrCzp3rrhdxDmyqGye6Po0TWhmuwH8/bqvGZD3K9IL4/jzzu8ZPWUgyxa+EvpgjDGNTkQmh7IO6VWrKj42dKiTAMrmki4/hHd9myY0ITaap351G/lFT9Fl68lsEeXy+X/hialDKNgb2kG6jDGNS0Qmh7IjgLLpWL1e5yhi+HCYPt35wi9rJsrPPzQBBPVspVpIiovmxWuOZ2PxxbDtD1yQ2ImXS7Zy6VvDWD73ufAGZ4xpsCK2z+Hmm505oGNjnYvhWrWCH3+EDh3qT79CTSzK2sOo5/7HgA6p3HTCWh6a+yd24eP2uPZcPeIVPAlp4Q7RGFPPWJ9DJcofAXTv7pyd9Oij9atfoSZ6ZTTlofO7882qHSzc0Z+3f/ElZyS248miDdz86unsWft1uEM0xjQgEXvkAFWffRQV5ZzSGhsLRUWHnrFUn6kqv35tPp8s3sK0GwfSt10z3p7zJH9a9i/SS0t5stNojj31AaeH3RgT8ezIoQpVHSUMHRr8foXqrqOoLRHhzxf1pHWKl7FvLqKwxMclJ9zNK0Mm4Yv2ctWa13n3tfOhMDdwlRpjGqWITg5VnX00fTo88wwcd5xzW9kIrXUVrAvpkr0xPH5xL9bu2M+Tn68EoEfGSUwb9QV9E9vwUPF6Jvz7dHy77cI5Y0zVIjo5QOjPPoqPd1p1KruOIlBO6tScywe045/frDkwBlOzhDQmXvIRo1oO4sWYQu6ZdjYF62YFrlJjTKMS0X0O4RDwYb+rkFNQzNAJM0mKi+aj208hNtr5HaCqvDLnCf667BV6FhXxdL/7Set3TeAqNsY0GNbnUI+E6kK6Jt4Yxo3owaptubz83boD20WEq0+4hydPHMfK2Diu+fExsmf8KbCVG2MaPEsOYRCqpqwh3VtyRtd0/vblKrbtO/S0rCFdLmTy0OfZGevlqp//zbqP7qh8sCljTEQ6YrOSiKSp6s4QxVMrwWhWCvWw3MGydsd+hk2YyXm9WvPkpb0rPL58x1JunP5LKMlnUpO+dLvoXxAVU7EgY0yjU9dmpdki8qaInCMRNNt3qIflDpYOzRO54dQOvDN/E3PX7arweLfm3XnpgjeJjW3CdTnzWfDaSCi2uXuNiXT+JIcuwGTgl8AqEfmTiHSpS6UiMkpEloiIT0T6H/bY/SKyWkRWiMiwutRTG6E4myjUbj2jE61TvDz64VJ8vopHih2aduSVC9+jWXwqNxWvZeFrI6G4Ho8VYowJuiMmB3V8rqqXAzcAVwNzRORrERlUy3oXAxcBM8tvFJHuwGXAscDZwLMiElXLOvxy+MVoDXX4jOokxEZz99CuLMray0c/ZVf6nNZJrXnxgrdI9aZyU+kGfnp1BBTtD3Gkxpj64ojJQUTSROQOEZkLjAVuA5oDdwOv1qZSVV2mqisqeWgE8LqqFqrqWmA1MKA2dfjr8Oaj+jYsd6CM7NOGbq2aMP7TFRSVVN7x3DKxJS9e8CZNvanc6NvEkqkjoDAnxJEaY+oDf5qV/gckAxeq6rmq+o6qlqjqXCDQY0K3ATaWW89yt1UgImNEZK6IzN2+fXuNK6qu+ai+DcsdCFEe4bfDu7FhVx6vfl/51dHZ2XDpOa14bMBbJHtTuUGzWTrlAhtuw5gI5E9yeFBVx6nqgdljRGQUgKo+XtVOIvKFiCyuZBkRgLhR1cmq2l9V+6enp9d4/+qaj955J/jDZ4TD6V3SGdQxjae/Wk1OQXGFx8uOoib9tSUvnP86TeJTuVG2sua1i60PwpgI409yuK+SbfcfaSdVHaKqPSpZ/lPNbpuAtuXWM9xtAddYm4+qIyLcf043du0vYvLMNQe2V3YUldGkDV9f+zVRccncULqBza9fBiWFYYzeGBNKVSYHERkuIn8H2ojI0+WWl4CSIMXzPnCZiMSJSAegMzAnSHUdaD764ANo2RLWrQtWTfVHr4ymnNurNS/MWsvOXOfLvqqjqHXrPEwa/jL5sYmMKVzBzreugtJg/emNMfVJdUcOm4G5QAEwr9zyPlCnU0xFZKSIZAGDgI9E5FMAVV0CTAOWAp8At6pqaV3qqs4zz8DixTBlipMoMjODVVP9cueQzhQUlx44eqjuKKpraleeGfY8W2PjuWnfAnLevR58QfuTGGPqCX+ukI5W1Xr9c7G2V0hHR0NpJd9z9X1K0EC4840FfLw4m5n3nkGLJl4uushJEmPGONOnZmcf2tcya9MsbvviFnoV5DOp1Vl4L3jGJg0ypoGr1RXSIjLNvTtfRBYdvgQl0hApa2OvLDE09Gsa/HX74M4UlyoTZ/wMHLkT/uQ2J/OnUx9nvtfLPZs+pfTLBn7puDGmWtU1K93h3p4HnF/J0mAd3sYOztSgIo2/U7pMh+aJXNSnDVO/38CWvf4NlzG8w3DuG3A/MxIT+POyf6H/mxjkKI0x4VJlclDVsktpLwaKVXV9+SU04QVH+TZ2j/sOXHQR3Hxz47imwV+3D+6Mz6c889/Vfu9zxTFXcG33q3kjuQkvzv4j/PRWECM0xoSLP6eyNgE+F5FvROTXItIy2EGFQtmZSj/+CLfcAiUljeuaBn+0TU3gF8e35fUfNpC1O8/v/X7T/y6Gtx/KU6nN+OjT38DPXwUvSGNMWPg9E5yI9AIuxTmSyFLVIcEMrCYa0kxw9c3mPfmcPn4GF/Vtw2MX9/J7v6LSIm789Fcs2LaASTv2MuDy96BN3+AFaowJuEDNBLcN2ALsBFoEIjATfkc1jefyAW15c14W63f6P9BebFQsTw3+B+2T2/ObtBRWvT4KdvjfPGWMqd/8GXjvFhGZAXwJpAE3qKr/PzFNvXfrGZ2I9gh/+3JVjfZLiUth4tDn8cY34+ZmXrZOHQk5W4MUpTEmlPw5cmgL/EZVj1XVR1R1abCDMqHVItnLLwe25735m1i9rWaD7LVOas2zQyeTE5vALQkl5L46ygbqM6YRqO46h2T37nhgg4ikll9CE54JlZtOPxpvTFSNjx4AuqV2Y8IZf+Pn2FjGajbFb15rw2wY08BVd+RQNlfDPJxhNMoPoWG9v41M86Q4rjkxkw8XbWb5ln013v/ENify0IkP8228lz/umYt+dDf4ebKDMab+qe46h/Pc2w6q2tG9LVs6hi5EEypjTu1IUmw0Ez5fWav9L+p8ETf0vIG3myTxws/vwKwJAY7QGBMq/nRIf+nPNtPwNU2I5VendODTJVtZvGlvrcq4rc9tnJN5Dn9LbcpHs8fDojcDHKUxJhSq63Pwun0LzUWkWbn+hkyqmJ3NNHzXndyBlPgYnqzl0YOIMO7kcfRr0Zffpzdn7sd3wNpvAhylMSbYqjtyuBGnf6Ebh/Y3/Af4R/BDM+GQ7I1hzKkd+Wr5Nuat312rMmKjYvnbmU/TJrkdd7RMY81bv4RtywIcqTEmmKrrc/ibqnYAxh7W53CcqlpyaMSuOTGTtMTYWvc9gHsNxFmTiPY25Za0JHa+OgpyImjgKmMaOH+uc/CJSNOyFbeJ6ZbghWTCLTEumptPP5pZq3fw/ZqdtS4no0kG/xj8LDtj4rgtsYT8qZdAYU4AIzXGBIs/yeEGVd1TtqKqu4EbghaRqReuHNieFk3i+OtnK/F3/K3K9EzvyWOn/YXFsbHc79tC6bSr7RoIYxoAf5JDlMjBKb9EJAqIDV5Ipj7wxkRx6xmdmLNuF9+s2lGnsga3G8y9x9/Ll4nxPLF7Hnx0p10DYUw9509y+AR4Q0QGi8hg4DV3m2nkLhvQloxm8fz54+WU+ur2ZX5l9yu58pgrmZKSzNTV78A3TwQoSmNMMPiTHH4L/Be42V2+BO4NZlCmfoiLjuLes7uxLHsf787fVOfyxvYfy5ltz+TxtFS+mv1XWPDqkXcyxoTFEZODqvpUdaKqXuIuk1S1ktmXTWN0fq/WHNe2KU98uoL8orr92aM8UTx26mMcm3Ysv23Zgp8+uQtWfxGgSI0xgeTPFdKdReQtEVkqImvKllAEZ8JPRHjgnGPYsq+AF2bV/c8eHx3P3wf/g7TE1vy6VQuy3roGNs+ve6DGmIDyp1npX8BEoAQ4A3gFmBLMoEz9MqBDKkO7t2TijJ/ZnlNY5/Kaxzfn2bMmUhybyC0tmrH31VGwa20AIjXGBIo/ySFeVb/EmVJ0vao+Apwb3LBMfXPf8G4UlviY8EXtL4wrr2NKR/525tNkxUTzm5QYiqZcBPtrf02FMSaw/EkOhSLiAVaJyK9FZCSQFOS4TD3TMT2JKwe257U5G/gpq3aD8h3u+FbHM+6k/2NubDQPReeir46CIv+nKjXGBI8/yeEOIAG4HegHXAlcHcygTP1019AupCXG8fv/LMZXx1Nby5zb8Vxu73M7HyV6+Ufez/DWdXaRnDH1gD9nK/2gqrmqmqWq16rqxao6OxTBmfol2RvDA+d2Y8HGPUybuzFg5V7f83ou7nwxk5sm8+7mb+Cju+wiOWPCzJ8jB2MOuLB3GwZ0SOXxT5aze39RQMoUER4Y+AAnHnUif0hvznfL3oCv/xKQso0xtWPJwdSIiDBuRA/2FZTw+CfLA1ZujCeGv572V45u1pm7WrVmxbfjYd7LASvfGFMz/lzncJI/22pCRMaLyHIRWSQi7x426uv9IrJaRFaIyLC61GOCo2urJlx/cgde/2Ejs+o47lJ5SbFJPDP4GRLjU7k1oy1bP74Llr4fsPKNMf7z58jh735uq4nPgR6q2gtYCdwPICLdgcuAY4GzgWfdgf5MPXPnWV3o2DyR3769iNzCwHUgt0psxTNDniUnOpZbM9qx/53r4eevAla+McY/1U0TOkhE7gbSReSucssjQJ2+sFX1M1Ut+0aZDWS490cAr6tqoaquBVYDA+pSlwkOb0wU40f1YvPefP48PbCzvHVL7cZfT/8rqz3K3UcdRfHro2HjnIDWYYypXnVHDrE41zNEA03KLfuASwIYw3XAx+79NkD502CyqGK+ahEZIyJzRWTu9u3bAxiO8Ve/9qlcf3IHpn6/IaDNSwAntzmZBwc+yLdRpYxrkY5OvQS2LA5oHcaYqkVX9YCqfg18LSIvqer6mhYsIl8ArSp56AFV/Y/7nAdwhuWYWtPyVXUyMBmgf//+dt5jmNw9tCtfLd/GXdMWMP2OU2ieFBewsi/pcglb9m9h0qJJJDVtwj3/Holc9wmkHR2wOowxlfOnzyFORCaLyGci8lXZcqSdVHWIqvaoZClLDNcA5wGj9eBUY5uAtuWKyXC3mXrKGxPFP67oy578Yu6atjBgF8eVubX3rVzR7Qr+He/huQQPvHIh7LWPhDHB5k9yeBOYDzwI3FNuqTURORtnTogLVDWv3EPvA5eJSJyIdAA6A9bYXM8d0zqZh87rzsyV25n8TWAH7BURfjvgt4w4egTPJsXyiicP/n0h5G4LaD3GmENV2axUTomqTgxwvf8A4oDP3RlIZ6vqTaq6RESmAUtxmptutbkjGobRJ7Tju593MP7TFfRt14wBHVIDVrZHPDxy4iPkleQxns9J3LWdi18+H67+EJLSA1aPMeYgOdLk8e7ZSduAd4ED4zWr6q6gRlYD/fv317lz54Y7jIi3r6CYC5/5lj15xfzn1pNom5oQ0PKLS4u57b+38d2m7/jLjr2cHX8UXP0BJDYPaD3GRAoRmaeq/St7zJ9mpatxmpG+A+a5i30TmwqSvTH886r+lJT6uOGVuQG9/gEgJiqGCadPoE+LPtyf3pQv87PglRGQV29+pxjTaPgz8F6HSpaOoQjOhFZ2Npx2GmzZUvsyOqYn8czovqzalsttr/5IcakvcAHizCT3zOBn6N68B2PTU/li/wZ45QJLEMYEmD/DZySIyIMiMtld7ywi5wU/NBNq48bBrFnw6KN1K+eUzuk8OuJY/rtiO/e+tSjgZzAlxSYxacgkjm3ek3tapPJ57jqnkzp/d0DrMSaS+TtNaBFworu+Cfi/oEVkQi4+HkRg4kTw+ZxbEWd7bY0+oT1jh3bh3fmbePTDpRypb6umkmKTeG7Ic/Ro3ot7W6Txec4aeOl8yLULIo0JBH+Sw9Gq+hegGMA99VSCGpUJqTVr4IorIMHtP05IgNGjYW0dp3W+9YxO/OrkDrz03bqgJYiJQybSo3kv7klP4/O8DfDSObBvc0DrMSYS+ZMcikQkHlAAETmacmctmYavdWtIToaCAvB6ndvkZGhV2fXtNSAiPHjuMVx7Uib/+nYdv3s3cDPIlUmKTeK5s56jZ3ov7mnejI+Ld8CLZ8PudQGtx5hI409yeBj4BGgrIlOBL3EuYDONyNatcNNNMHu2c1uXTunyRISHzuvOLacfzWtzNnD76/MpKA7spSuJMYk8d9Zz9G7Zh9+mNmGaJw9eHA47VgW0ntoIRCe/MeFwxOscAEQkDRiI05w0W1UDO8paHdl1Dg3DpK9/5rFPltO7bVMm/7I/6U0CNw4TQEFJAXd/fTczs2ZyR24x1+eVwC/fhVY9A1pPTdxyC0yaBDfeCM8+G7YwjKlUddc5+HMR3EjgK1Xd6643BU5X1fcCHGetWXJoOD5ZvIU731hAamIsz13Zj54ZKQEtv9hXzIOzHmT62ulcm+/jzt37kMumQodTA1rPkcTHO81zh/N6IT8/pKEYU6W6XgT3cFliAFDVPThNTcbU2Nk9WvHmTYNQVS6a+C3//GZNQDuqYzwx/PmUP3Np10v5V7yHP6Q3p3TKxfDTWwGrwx/B6uQ3JlT8SQ6VPcefMZmMqVSPNilMv+MUTu/agv/7aBnXvfQDW/ZW8jO7ljzi4YETHuCGnjfwdkwJv2mbSd4718N3/whYHUcSrE5+Y0LFn+QwV0SeFJGj3eVJnCE0jKm1pgmxTP5lP/5wwbH8b81Oznrya/49e33AzmYSEW7vezsPnPAAMz1FXNehKzu+/D188jvnYo4QCFYnvzGh4E+fQyLwe2AIzumsnwN/VNX9wQ/PP9bn0LCt37mf3737E9+u3slxbZvywDnHBHRU1xkbZ3Dv1/fQTIWJG36mY6fhMHISxCYGrA5jGqJad0iLSBTwhaqeEazgAsGSQ8Onqrzz4yb+8ulytu4r5KzuLbl3WFc6t2wSkPKX7FzCr7/8NYWF+3hq8yYGpHSCy1+HlIwj72xMI1XrDml3LgWfiAT2lBJjDiMiXNwvgxljz+CeYV353887OWvCTG54ZS7z1td9zKRj045l6jlTadGkLTe2bsW0wmyYfAZs/CEA0RvT+PjTrPQfoA9Oc9KBpiRVvT24ofnPjhwan137i3jp27W8/L/17M0vpl/7Zlw+oB3n9mxNfGxUrcvdV7SPe2fey7ebvmVUIdy/dSsxFzwNx10WwOiNaRjqep3D1ZVtV9WXAxBbQFhyaLz2F5bwxg8bmTJ7PWt27KdJXDTn9z6K83sdxYAOqUR5aj7MV6mvlKfnP82Li1+kr8bx1w0/03zAzXDWHyAqJgivwpj6qU7JwS0gHminqisCHVwgWHJo/FSVOWt38cYPG5m+OJuCYh/Nk2I5q3srhh7bkhM6pJIQW7MzrD9Z+wm///b3JKvy1Mb19GzZF0b9C5KPCtKrMKZ+qeuRw/nAE0CsqnYQkd7Ao6p6QcAjrSVLDpFlf2EJM1ZsZ/ribP67fBt5RaXERAl92zXj5E7NGXR0Gj3apOCNOXLz0/Jdy7njqzvYlreFu3fnMLrIg1zyYsivqDYmHOqaHOYBZwIzVLWPu22xqvYIeKS1ZMkhchUUlzJn7S6+Xb2DWat3sGTzPgCiPUK31k04LqMpvds25ZjWyRydnlRpf8Xewr08OOtBZmTN4MxiD49mbyTl9AfgpDvB48+lQMY0THVNDrNVdaCIzC+XHBapaq8gxForlhxMmZ25hcxbv5uFWXtYsHEPizbuJcedy1oE2jZLoEvLJDq1aELH5olkNIunbWoCLZPjeH3FVCbMe5IWRPFE1gZ6ZpwIFz4Hya3D/KqMCY66JocXcIbpvg+4GLgdiFHVmwIdaG1ZcjBV8fmUNTv2s3JrDiu35rBqWy6rtuawdsd+iksPfvY9Aq1T4klNzWZL3D8p0p3ctCeHq/aXsm7gY3i6n0fzpDiaJsQQE2VHE6ZxqGtySAAeAIa6mz4F/k9VAzcYTh1ZcjA1VVzqY8veAjbuyiNrdz4bdzu3Wbvz2JK7m13e1/A0WUinAnhq+2a+LTiVcSVXko+XhNgoUuJjSPbGOLfxMSTHR5MSH0MTbwwJsVEkxEYRHxNFQmy0c9/d5tyPJiEmirgYD7FRHqIt2ZgwqS45VHl6h4h4gZuATsBPwCBVLQlOiMaEVkyUh7apCbRNTajiGefy7soPePyHP3Gxty1jd8zhf6Wr+azzw6yI6c6+/GL2ukvW7jxyskvYm19MbmHN/0U84sQTG+0hLtpz4H5MlJM8YqPL3UZ7iPII0R4hqvwiQnSU4BHnMY/n0NsoEaI8ngrPiRLweAQRQQCPCCJOTIJzX0ScdXEe55DnOfvJYeseTyX749xSyX5wcO5hOeTsZDlkW/mHxN146Lay5x26H9U855Bt5Z5/eBmV1UNlZRwW36HbKtZdVXz+SvbGkJIQ+FOwqzxyEJE3cOaN/gYYDqxT1d8EPIIAsCMHEyxb92/loe8e4rvN3zGoGH6/ZTNt+10Pgx+qdGwmn0/JLy4lr6iU/KJS8opLDt4vKiWvqOTA/aJSH0UlzlJc6qPQvS0q8VFUevD+4dtLSpVSn1Kq7u3hiyqlpc5tiU/x+Zxb0zjddNrR3De8W632rVWzkoj8pKo93fvRwBxV7VurCILMkoMJJlXljRVvMGHeBHwlBdyyaye/lGZEX/B36Hh6uMPzi6riUyokkBKfDwV8qqDgU1Cc56oqqjhLuW1Onilbd/Ytu4VD15WD5ZTfXzlYtlOaunGWi7lc7OXXy69oua0HylIqPP/wMrSSwiqvu5p6Div7kBIPKevQ8qtrya/qoeqa/7u2akKvjKZVF1qNWjUr4Rw1lAVWIjU91jGmkRARLut2Gae3PZ0/fv9HnmQG00vhkddHcWyPy2DIHyAhcKPIBoOI04RUmyvKTWSqrifsOBHZ5y45QK+y+yKyL1QBGlNftEpsxdNnPM2E0yewM6k5V7RpxZ/Xf8jef/SDeS+HbJ4IY0KhyuSgqlGqmuwuTVQ1utz95FAGaUx9ISIMaT+E/1z4PqO6XsrryU04r0Uyb8z4HSUvDIbN88MdojEBEZZz6ERknIgsEpEFIvKZiBzlbhcReVpEVruP18s+DmOaxDbhwYEPMu38N+ncqh//1zyVX3i2MeeVs+GDOyBna7hDNKZOwnWC9XhV7aWqvYEPgYfc7cOBzu4yBpgYnvCM8U/X1K68MOwFnjz9SfJS2vCr1i24LWs6K5/tDzMeh8LccIdoTK2EJTmoavk+i0QOdtKPAF5Rx2ygqYjY2AWmXhMRzmp/Fv8Z+QF39L2DecmpXNKqKfctmcTGZ/rC3H9BqV0iZBqWsF2aKSJ/FJGNwGgOHjm0ATaWe1qWu62y/ceIyFwRmbt9+/bgBmuMH+Ki4ri+5/V8fMlnXNvjOr5MbsYFqXGM+98f2PLcQFj8NvhKwx2mMX4JWnIQkS9EZHElywgAVX1AVdsCU4Ff17R8VZ2sqv1VtX96enqgwzem1lLiUriz351Mv/gThqZdx+PPvMWQkiQemTGWDZYkTABlZ8Npp8GWLYEvO2jJQVWHqGqPSpb/HPbUqTgD+gFsAtqWeyzD3WZMg5OekE7Ox3eSt6ofiV8/wQfJKZyfUMC9X9/DiudOgJ/esuYmUyfjxsGsWfDoo4Ev26+Z4AJeqUhnVV3l3r8NOE1VLxGRc3GOIs4BTgCeVtUBRyrPrpA29U18PBRUMjRldGwJff95Anm+Ik7Ky+fK0nhO7H8rnj5XQlxS6AM1DVJVny+vF/Lz/S+nuiukw9Xn8JjbxLQIZ7TXO9zt04E1wGrgeeCWMMVnTJ2sWQNXXAEJ7rh+CQkwejRsXB/NZ7/4iluPu4UVKenc3ARG/PQUr07qzf7PHoScILQPmEanqs/X2rWBq6Nmk+4GiKpeXMV2BW4NcTjGBFzr1pCc7Py683qd2+RkaNUKIIWbet/Mr3pez6frP2Xqgsn8OWctf9/0Lhe+9CoXtzqJTgNvg4zjaz5Ep4kI1X++AsMGkjcmSLZuhZtugtmzndvDOw1jomI4r+N5vHbR+0w9Zyqntj2d15OTGJm3gNHTr+Tt5wewf/azUJgTnhdg6rUjfb7qKix9DoFmfQ6msdhVsIsPVrzNu8um8nPhTuJ9PoblF3NBq4H0638rnrYD7GjCBEydZoJrCCw5mMZGVVm0fSHvLnyejzd/Sx6ltCgpYVhpLMMzz6ZH/5uR1Mxwh2kaOEsOxjRgecV5zFz7CR8vmcI3e1dRLJBRXMzw6DSGdL6QY3pfhySmhTtM0wBZcjCmkdhXtI8vl03j4+XT+D4/G59Aq5ISTo9O5Yz2Qzi+743EJB8V7jBNA2HJwZhGaFf+TmYumcp/f/6I7/I3UyCQ5PNxsiRxylEnMqjnVaS37hPuME09ZsnBmEauoDif2cum8d+V7zAjZw273PMQO5UKJyZ3ZFDH4fQ79nLi42wqFnOQJQdjIohPfaxc8yXfLX+D73Ys5EfNp1iEGFX6RiXTr3lP+h49nF4dhhIfkxDucE0YWXIwJkyys+Gyy+CNNwJ7gVJN5O/fzo+LXuG79V8we/8GVkUJKkK0QvfoJvRL60HfjmfTJ3MwKXEp4QnShIUlB2PC5JZbYNIkuPFGePbZcEcDqLIvewELlr3JvOzv+TE/m8UxHkrcayfae+I5NrkDPVqfQI/2p9Mt7Rjio+PDHLQJFksOxoRYoAZGCzqfj4LNP/LT8ndYkD2bxXmbWRwtbIt2RtaJAjrFNqNH6jF0zziJLum96NS0E0mxNkhgY2DJwZgQy86GsWPhvfcgL88ZGG3kSHjiifA1L/nF54Odq9m+5ksWb/yan3atYElpDovjYtkXFXXgaW2im9A5pQNdWvahS3pPujTrQrsm7YjyRFVTuKlvqksOYRl4z5jGLhQDowWFxwPpXUhP78IZJ9zMGQB5u9CNP7Bl4yxWbv2RVXvXsVK3sDJvJ9/sWEip2yQVKx7aedPJTOlIh+bHktm0A5nJmWSmZJIca2dJNTSWHIwJkrKB0caMgcmTnaOJBikhFek6jNZdh9EaOA0gdxtkL6Jw0zzWZs9h5a6VrC7Zw9qYXFbvzWJG9ncH+jEAUmOSyEzOpH3TTmQkt6VNUhvaJLUho0kGad40xMaLqnesWckYExgF+2D7Cti2lOJtS9m0bRFr965lXWku62JiWBcTzYaYGHZEHdr05PXE0CbxKNokt6dNEydptE5sTcvElrRMaEnz+OZEe+x3bDBYs5IxJvi8ydD2eGh7PDFApruQtwu2L4dty2DXGvJ3rGTznp/ZlLeVLA9siolmU85eNu38mR9josk97CDCg9A8Po2Wia1pmdCSFgktDiSOFgktSItPI82bRnJssh2BBJAlB2NMcCWkQvsTnQWIB44Gji4tgb0bYdfPsNNZdOdq9u1Zz5a8LWylhK3RUWyNimZbTg5bY7NZGxPHbA/k4qtQTYwnhlRv6oFkcfhtanwqad40msY1pWlcU2KiYkL7PjQwlhyMMeERFQ2pHZyl0xAABEgBUlTpun877NngLHs3uved2/37NrLVV8DWqCh2ll9i97Mzdgc7oqNYIbBLiymh8qbz+Oj4A4kiOS754P1Y9763KSmxKaTEpZAcm0xSbBJJMUnER8fXmyOUYF5kacnBGFP/iEBSC2fJqNgknqhKx4I9dMzZAjnZztzb+zY7t2XrOdlozhb2ibIzynMggez1eNgTHcPeOB97C0rYG72XPR4PK0TZqyXs9RXhqyKhAERJ1IFEkRSTRFJsEk1imhzY1iS2SYXH46PjSYhJcG6jEw7cj/HU7ehl3DiYNQsefTTwF1lah7QxpvHy+SBvh5Mw9m+H/Tvd2+3O9v073HX3fvF+fECOR9jniWKPx8OeKA+5Hg85sQnkxiaQGxNHTkwsuVHR5HqiyBWnmStHS8n1FZFbWkhpJc1elYnxxFRIHNWte6O8eKO9XHncxRQXVrympKYXWVqHtDEmMnk8B49A/FGUhydvByn7t5OyfydtC/ZA/m7Id28PrO+GnHL3fcUHilAgX4T9HiHH4ySW/Kho8mMTyItJID/GS15MHHnRMeRHRZMn0eT7POQVl5JfkkMee9mupeRrKXmlReT7isgvLaRESw/U0fHxiWx5fSz7fjwTLYo/5CLLQLHkYIwxZWITILYdNG3n/z6qUJx3IIFI/m4SCnNIKMwhvXAfFOZUsuyDvBwo3AlFuc624rxqqykCCjxCoQgF0du5V3bzRnEccbHFFBTEBPwiS0sOxhhTFyIQm+gsKW1qX05pCRSVSyBFeU7iKM6Dojxii3KJLc6Dov1QtJ/CD3px87BvGHN9MZO/GhLwiywtORhjTH0QFQ3xzZzFD++cdfD+MxcHPhxP4Is0xhjT0FlyMMYYU4ElB2OMMRVYcjDGGFOBJQdjjDEVWHIwxhhTQViTg4jcLSIqIs3ddRGRp0VktYgsEpG+4YzPGGMiVdiSg4i0BYYCG8ptHg50dpcxwMQwhGaMMREvnEcOE4B74ZDhD0cAr6hjNtBURFqHJTpjjIlgYUkOIjIC2KSqCw97qA2wsdx6lrutsjLGiMhcEZm7ffv2IEVqjDGRKWjDZ4jIF0Blw0A9APwOp0mp1lR1MjAZnCG761KWMcaYQwUtOajqkMq2i0hPoAOw0J1NKQP4UUQGAJuAtuWenuFuM8YYE0Ihb1ZS1Z9UtYWqZqpqJk7TUV9V3QK8D1zlnrU0ENirqgEea9AYY8yR1LdRWacD5wCrgTzg2vCGY4wxkSnsycE9eii7r8Ct4YvGGGMM2BXSxhhjKmHJwRhjTAWWHIwxxlRgycEYY0wFlhyMMcZUYMnBGGNMBZYcjDHGVGDJwRhjGqjsbDjtNNiyJfBlW3IwxpgGatw4mDULHn008GVbcjDGmAYmPh5EYOJE8PmcWxFne6BYcjDGmAZmzRq44gpISHDWExJg9GhYuzZwdVhyMMaYBqZ1a0hOhoIC8Hqd2+RkaFXZDDq1ZMnBGGMaoK1b4aabYPZs5zbQndJhH5XVGGNMzb3zzsH7zzwT+PLtyMEYY0wFlhyMMcZUYMnBGGNMBZYcjDHGVGDJwRhjTAWWHIwxxlQgqhruGOpMRLYD62u5e3NgRwDDacjsvXDY++Cw9+GgxvpetFfV9MoeaBTJoS5EZK6q9g93HPWBvRcOex8c9j4cFInvhTUrGWOMqcCSgzHGmAosOcDkcAdQj9h74bD3wWHvw0ER915EfJ+DMcaYiuzIwRhjTAWWHIwxxlQQ0clBRM4WkRUislpE7gt3PKEiIm1F5L8islRElojIHe72VBH5XERWubfNwh1rKIhIlIjMF5EP3fUOIvK9+7l4Q0Riwx1jKIhIUxF5S0SWi8gyERkUiZ8JEbnT/b9YLCKviYg3Ej8TEZscRCQKeAYYDnQHLheR7uGNKmRKgLtVtTswELjVfe33AV+qamfgS3c9EtwBLCu3/jgwQVU7AbuBX4UlqtD7G/CJqnYDjsN5TyLqMyEibYDbgf6q2gOIAi4jAj8TEZscgAHAalVdo6pFwOvAiDDHFBKqmq2qP7r3c3C+BNrgvP6X3ae9DFwYlgBDSEQygHOBf7rrApwJvOU+JVLehxTgVOAFAFUtUtU9ROBnAmcStHgRiQYSgGwi8DMRycmhDbCx3HqWuy2iiEgm0Af4HmipqtnuQ1uAluGKK4SeAu4FfO56GrBHVUvc9Uj5XHQAtgP/cpvY/ikiiUTYZ0JVNwFPABtwksJeYB4R+JmI5OQQ8UQkCXgb+I2q7iv/mDrnODfq85xF5Dxgm6rOC3cs9UA00BeYqKp9gP0c1oQUIZ+JZjhHSx2Ao4BE4OywBhUmkZwcNgFty61nuNsigojE4CSGqapaNhvtVhFp7T7eGtgWrvhC5CTgAhFZh9OseCZOu3tTt0kBIudzkQVkqer37vpbOMki0j4TQ4C1qrpdVYuBd3A+JxH3mYjk5PAD0Nk9CyEWp9Pp/TDHFBJuu/oLwDJVfbLcQ+8DV7v3rwb+E+rYQklV71fVDFXNxPn7f6Wqo4H/Ape4T2v07wOAqm4BNopIV3fTYGApEfaZwGlOGigiCe7/Sdn7EHGfiYi+QlpEzsFpc44CXlTVP4Y3otAQkZOBb4CfONjW/jucfodpQDucIdB/oaq7whJkiInI6cBYVT1PRDriHEmkAvOBK1W1MIzhhYSI9MbpmI8F1gDX4vyAjKjPhIj8AbgU56y++cD1OH0MEfWZiOjkYIwxpnKR3KxkjDGmCpYcjDHGVGDJwRhjTAWWHIwxxlRgycEYY0wFlhyMMcZUYMnBGGNMBZYcTK2ISCsReV1EfhaReSIyXUS6hDuumhKRR0RkbLjjOBJ3roVbArGfiHwXoJjiReRrd/j7Iz33RhFR92LDsm23utuGicjMcsNTmHrAkoOpMXdYgXeBGap6tKr2A+6nkY/YGWZNgUqTgziq+l+usJ+qnhigmK4D3lHVUj+e2xNYCHQDEJEEnCuPtwMLcOaKuDRAcZkAsORgauMMoFhVnyvboKoLVfUbEbnLnUFrsYj8puxxEcl0Zxh7SURWishUERkiIt+6s4wNKPecqe5MZG+5XyJlZVQo291ncbnnjHWPBjLdMp53Z/X6TETi3ec84MYwCygbS6gCEblKRBaJyEIR+bcfcVRVX4VyRORKEZkjIgtEZFLZr+9qynkMONp9/nj3eStE5BVgMdBWRN5zj+KWiMgYN9xD9nPryK3taznMaMqNMSQiM0Sk7Ms/rfzfBeiFM/xEN3f9duBNwKeqW4H33PJMfaGqtthSowXnH3tCJdv74YzXlAgkAUuAPu5jmThj1fTE+VEyD3gREJwhkt9zn6PASe4+L+KMd1Rl2e4+i8vFMBZ4pFx9vd3t04Ary5WTACQDq8vqOOy1HAusBJq766l+xFFZfRXKAY4BPgBi3G3PAlcd9j4dXs7hrzMTZ1ysgeW2lcUYj5Mw0g7fz308tzav5bAyYoEth23LAjzu/TOA18o9th3oCHyCczQzHzgd+MJ9PArYHu7Pti0HFztyMIF0MvCuqu5X1Vyc4Y5PKff4WlX9SVV9OF9EX6rzzfATzhcSwEZV/da9P8Ut05+yK7NWVRe49+e5dZzilpOnzhwWVY3EeybwpqruANCDg81VF0dl9VVWzmCcL+YfRGSBu97xCHFXZr2qzi63fruILARm4wxH37mK/crU9LWU1xzYU7YiIu2BTe7fFpwjhUXuY22Bnaq6BmgB3AP8HeiC87dHnaapIhFpcoSYTYhYB5CpjSUcHL64JsqPYukrt+7j4Gfx8JEgjzQyZAmHNo96q6ivFOcXdTD5W58AL6vq/XUsZ/+BAp2O3iHAIFXNE5EZHPpe1NSRYsg/rPzjcJOBqx/whnu/J24SAHJwJs8ZgDMi8o/l9okDCuoQswkgO3IwtfEVEFeuXRsR6YXTsXihOGPhJwIjcYYGr4l2IjLIvX8FMMu9/00VZW8FWrht3HHAeUcof6ZbTrz7K/X8al7jKBFJc19f6hHiqEpl5XwJXCIiLcq2ub+8q5MDVPerOgXY7SaGbsBAP/ar6Ws5QFV3A1EiUpYgeuMmCxHpjNNUWJYQepW7Px74tXukcCBpuO/PDnUm2DH1gCUHU2NuU9BIYIg4p7IuAf4MbAZeAubgzA3xT1WdX8PiVwC3isgyoBkw0a3zx8rKdr9MHnW3fw4sP0LsP+L8ol0IfIwz6VNlz1sC/BH42m2qebK6OKqpr0I5qroUeBD4TEQWuXG3PkLcO4Fv3Y7j8ZU85RMg2n3fHsNpWqp2v5q+lkp8xsFmv+MAj/saH8KZIKdskqCeOH0gqOqHqvo/d3t3nKNQcPooPqpB3SbIbD4HU2+ISCbwoar2CHcs5shEpC9wp6r+UkRWAX1VNaeWZb0D3KeqKwMapKk1O3IwxtSKe+TxXxFJcVZrnRhigfcsMdQvduRgjDGmAjtyMMYYU4ElB2OMMRVYcjDGGFOBJQdjjDEVWHIwxhhTgSUHY4wxFVhyMMYYU8H/A2gOPHYPMbGnAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# set loss function to NLS\n", + "opt.handle_opt_set (handle, 'NLDF Loss Function Type = L2')\n", + "\n", + "# solve huber loss function with various regularisations\n", + "for reg in regs:\n", + " # set reg type\n", + " opt.handle_opt_set(handle, 'Reg Term Type =' + reg)\n", + " # call the solver\n", + " sol = opt.handle_solve_nldf(handle, lsqfun, lsqgrd, x, nres,data=data, io_manager=iom)\n", + " # calculate response\n", + " resp = sol.x[0] + ((sol.x[1] - sol.x[0]) / (1 + (sol.x[2] / dose)**sol.x[3]))\n", + " # plot curve\n", + " plt.plot(dose, resp, label='reg = '+reg)\n", + "\n", + "# show NLS curves\n", + "plt.title(r\"$l_2$-norm loss function\")\n", + "plt.plot(t, y, '*b')\n", + "plt.xlabel(r\"Compound concentration ($\\mu M$)\")\n", + "plt.ylabel(\"Percent activity\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "id": "d8984282", + "metadata": {}, + "outputs": [], + "source": [ + "# destroy the handle\n", + "opt.handle_free(handle)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1258bec0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}