{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Determining the stability of the fits\n",
"After doing a fit, it is quite common to validate the result that has been obtained, since the stability and the type of minimum (global or local) is not reflected by the FCN or by the minimizer. In this section, several tools meant to evaluate the quality of the fits, and included in the Minkit package, are discussed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Scanning the FCN\n",
"The most common operation to do after a fit is to scan the profile of the FCN for the variables involved in it. This can be done in many different ways: fixing the rest of the parameters and evaluating the minimization function, mapping the full paramter space by varing all the parameters, or doing a fit for fixed values of one of the parameters (typically done by MINOS). In the Minkit package an easy and fast solution to the first two items is proposed using the *minkit.fcn_profile* function. Let's see how to do it on a fit to a Gaussian function."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" \n",
" FCN = 28257.803801867896 \n",
" TOTAL NCALL = 92 \n",
" NCALLS = 92 \n",
" \n",
" \n",
" EDM = 5.589309293901719e-06 \n",
" GOAL EDM = 1e-05 \n",
" \n",
" UP = 1.0 \n",
" \n",
"
\n",
"\n",
" \n",
" Valid \n",
" Valid Param \n",
" Accurate Covar \n",
" PosDef \n",
" Made PosDef \n",
" \n",
" \n",
" True \n",
" True \n",
" True \n",
" True \n",
" False \n",
" \n",
" \n",
" Hesse Fail \n",
" HasCov \n",
" Above EDM \n",
" \n",
" Reach calllim \n",
" \n",
" \n",
" False \n",
" True \n",
" False \n",
" \n",
" False \n",
" \n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" \n",
" + \n",
" Name \n",
" Value \n",
" Hesse Error \n",
" Minos Error- \n",
" Minos Error+ \n",
" Limit- \n",
" Limit+ \n",
" Fixed? \n",
" \n",
" \n",
" 0 \n",
" c \n",
" -0.00837095 \n",
" 0.00995019 \n",
" \n",
" \n",
" -1 \n",
" 1 \n",
" No \n",
" \n",
" \n",
" 1 \n",
" s \n",
" 0.99453 \n",
" 0.00706184 \n",
" \n",
" \n",
" 0.7 \n",
" 1.3 \n",
" No \n",
" \n",
"
\n",
"\n",
"\n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import minkit\n",
"\n",
"x = minkit.Parameter('x', bounds=(-4, +4))\n",
"c = minkit.Parameter('c', 0., bounds=(-1, +1))\n",
"s = minkit.Parameter('s', 1., bounds=(0.7, 1.3))\n",
"g = minkit.Gaussian('g', x, c, s)\n",
"\n",
"data = g.generate(10000)\n",
"\n",
"with minkit.minimizer('uml', g, data) as minimizer:\n",
" result = minimizer.migrad()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to scan the variables, we must define the points to be evaluated as a numpy array."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7xddX3n/9c7CQgdRUKIFoEEGLAWbUWTgaj9WaGKaH+/4kyxRaniBRkV66WXUaed8VKdqe2oP6mowyAVHFqkKJpaKFINtVoDJMpFQCSmpoBMoSGI/KxAyOf3x/6eug3nspOz99l7n/N6Ph7rcdb+ru9a+7v27Xs+63tZqSokSZIkSbO3aNgFkCRJkqT5wgBLkiRJkvrEAEuSJEmS+sQAS5IkSZL6xABLkiRJkvpkybALMNf233//OuSQQ4ZdDEnSNDZu3PjPVbV82OUYBuspSRoPU9VVCy7AOuSQQ9iwYcOwiyFJmkaSLcMuw7BYT0nSeJiqrrKLoCRpwUqyV5Krk1yX5MYk72rpb0iyKUkl2b8r/3OSfD/JtW35r13bTkhyS9vvbV3phya5qqV/Ksmec3uWkqS5ZIAlSVrIHgCOq6qnAkcBJyRZA3wVeC4w2dXJv6uqo9ryboAki4GzgBcARwIvSXJky/8+4INVdTiwDXj1QM9IkjRUBliSpAWrOu5vD/doS1XVN6rqu7twqKOBTVW1uaoeBC4ETkwS4Djg4pbvPOBF/Sm9JGkUGWBJkha0JIuTXAvcBVxRVVfNsMszWpfCy5I8uaUdCNzWlef2lrYMuLeqtu+UvnMZTk+yIcmGu+++e1bnI0kaLgMsSdKCVlUPV9VRwEHA0UmeMk32rwMrW5fCPwE+26cynF1Vq6tq9fLlC3LyREmaNwywJEkCqupeYB1wwjR57pvoUlhVlwJ7tEkw7gAO7sp6UEvbCuybZMlO6ZKkecoAS5K0YCVZnmTftr438DzgW9Pk/+k2rookR9OpR7cC1wBHtBkD9wROBtZWVdEJ2k5qhzgV+NygzkeSNHwGWJKkvtm4ZRtnrdvExi3bhl2UXh0ArEtyPZ0g6Yqq+nySNya5nU6L0/VJzmn5TwK+meQ64Ezg5DZRxnbgDcDlwM3ARVV1Y9vnrcBvJdlEZ0zWxwd5QmP4HkjSvLLgbjQsSRqMjVu2cco563lw+w72XLKIC05bw6qVS4ddrGlV1fXA0yZJP5NOALVz+oeBD09xrEuBSydJ30xnlsGBG8f3QJLmG1uwJEl9sX7zVh7cvoMdBQ9t38H6zVuHXaQFx/dAkobPAEuS1BdrDlvGnksWsTiwx5JFrDls2bCLtOD4HkjS8NlFUJLUF6tWLuWC09awfvNW1hy2zK5pQ+B7IEnDZ4AlSeqbVSuX+k/9kPkeSNJw2UVQkiRJkvrEAGsXOPWtJEmSpOnYRbBHTn0rSZIkaSa2YPXIqW8lSZIkzWRgAVaSvZJcneS6JDcmeVdL/6UkX09ybZKvJDm8pT8qyaeSbEpyVZJDuo719pZ+S5Lnd6Wf0NI2JXnboM4FnPpWkiRJ0swG2UXwAeC4qro/yR7AV5JcBnwUOLGqbk7yeuD3gVcArwa2VdXhSU4G3gf8epIjgZOBJwNPAP4myRPbc5wFPA+4HbgmydqqumkQJ+PUt5IkSZJmMrAAq6oKuL893KMt1ZZ9Wvpjge+19ROBd7b1i4EPJ0lLv7CqHgD+Ickm4OiWb1NVbQZIcmHLO5AAC5z6VpIkSdL0BjrJRZLFwEbgcOCsqroqyWnApUn+BbgPWNOyHwjcBlBV25N8H1jW0td3Hfb2lsZE/q70Y6Yox+nA6QArVqzow5lJkiRJ0iMNdJKLqnq4qo4CDgKOTvIU4C3AC6vqIOBPgQ8MsgytHGdX1eqqWr18+fJBP50kSZKkBWpOZhGsqnuBdcALgKdW1VVt06eAZ7b1O4CDAZIsodN9cGt3enNQS5sqXZIkSZKGYpCzCC5Psm9b35vOZBQ3A4/tmqRiIg1gLXBqWz8J+FIbx7UWOLnNMngocARwNXANcESSQ5PsSWcijLWDOh9JkiRJmskgx2AdAJzXxmEtAi6qqs8neQ3w6SQ7gG3Aq1r+jwOfbJNY3EMnYKKqbkxyEZ3JK7YDZ1TVwwBJ3gBcDiwGzq2qGwd4PpIkSZI0rUHOIng98LRJ0i8BLpkk/UfAi6c41nuB906Sfilw6awLK0mSJEl9MCdjsCRJkiRpITDAkiRJkqQ+McCSJPVs45ZtnLVuExu3bBt2USRJGkkDvdGwJGn+2LhlG6ecs54Ht+9gzyWLuOC0NaxauXTYxZIkaaTYgiVJ6sn6zVt5cPsOdhQ8tH0H6zdvHXaRJEkaOQZYkqSerDlsGXsuWcTiwB5LFrHmsGXDLpIkSSPHLoKSpJ6sWrmUC05bw/rNW1lz2DK7B0qSNAkDLElSz1atXGpgJUnSNOwiKElasJLsleTqJNcluTHJu1r6G5JsSlJJ9u/KnyRntm3XJ3l617ZTk9zallO70lcluaHtc2aSzO1ZSpLmkgFWnzmFsSSNlQeA46rqqcBRwAlJ1gBfBZ4LbNkp/wuAI9pyOvBRgCT7Ae8AjgGOBt6RZKKp76PAa7r2O2GQJyRJGi67CPaRUxhL0nipqgLubw/3aEtV1TcAJmlsOhE4v+23Psm+SQ4AngNcUVX3tP2uoBOsXQnsU1XrW/r5wIuAywZ5XtPZuGWb4+gkaYAMsPposimMrbwkabQlWQxsBA4Hzqqqq6bJfiBwW9fj21vadOm3T5K+cxlOp9MixooVK3b9JHrkhUBJGjy7CPaRUxhL0vipqoer6ijgIODoJE8ZQhnOrqrVVbV6+fLlA3se72UmSYNnC1YfOYWxJI2vqro3yTo6Y6S+OUW2O4CDux4f1NLuoNNNsDv9ypZ+0CT5h2LiQuBD23d4IVCSBsQAq8+cwliSxkeS5cBDLbjaG3ge8L5pdlkLvCHJhXQmtPh+Vd2Z5HLgv3VNbHE88PaquifJfW3ijKuAlwN/MrATmoEXAiVp8AywJEkL2QHAeW0c1iLgoqr6fJI3Av8J+Gng+iSXVtVpwKXAC4FNwA+BVwK0QOoPgGvacd89MeEF8HrgE8DedCa3GNoEF+CFQEkaNAMsSdKCVVXXA0+bJP1M4MxJ0gs4Y4pjnQucO0n6BmDOx3VJkobDSS4kSZIkqU8MsCRJkiSpTwywJEmSJKlPDLAkSZIkqU8MsCRJkiSpTwywJEkAbNyyjbPWbWLjlm3DLookSWPLadolSWzcso1TzlnPg9t3sOeSRVxw2hrvlSRJ0m6wBUuSxPrNW3lw+w52FDy0fQfrN28ddpEkSRpLBliSJNYctow9lyxicWCPJYtYc9iyYRdJkqSxZBdBSRKrVi7lgtPWsH7zVtYctszugZIk7SYDLEkS0AmyDKwkSZoduwhKkiRJUp8YYM0xp0GWJEmS5i+7CM4hp0GWJEmS5jdbsOaQ0yBLkiRJ89vAAqwkeyW5Osl1SW5M8q6WniTvTfLtJDcneWNX+plJNiW5PsnTu451apJb23JqV/qqJDe0fc5MkkGdTz84DbIkSZI0vw2yi+ADwHFVdX+SPYCvJLkM+FngYOBJVbUjyeNa/hcAR7TlGOCjwDFJ9gPeAawGCtiYZG1VbWt5XgNcBVwKnABcNsBzmhWnQZYkSZLmt4EFWFVVwP3t4R5tKeB1wEurakfLd1fLcyJwfttvfZJ9kxwAPAe4oqruAUhyBXBCkiuBfapqfUs/H3gRIxxggdMgS5IkSfPZQMdgJVmc5FrgLjpB0lXAvwV+PcmGJJclOaJlPxC4rWv321vadOm3T5IuSZIkSUMx0ACrqh6uqqOAg4CjkzwFeBTwo6paDfwv4NxBlgEgyektoNtw9913D/rpJEkaa95SRJJ235zMIlhV9wLr6IyRuh34TNt0CfDzbf0OOmOzJhzU0qZLP2iS9Mme/+yqWl1Vq5cvXz67k5EkaR6buKXI+79wC6ecs94gS5J20SBnEVyeZN+2vjfwPOBbwGeBY1u2XwS+3dbXAi9vswmuAb5fVXcClwPHJ1maZClwPHB523ZfkjVt9sCXA58b1PlIkrQQeEsRSZqdQc4ieABwXpLFdAK5i6rq80m+AlyQ5C10JsE4reW/FHghsAn4IfBKgKq6J8kfANe0fO+emPACeD3wCWBvOpNbjPQEF5IkjbqJW4o8tH2HtxSRpN2QzqR9C8fq1atrw4YNwy6GJGkaSTa2sboLzijUUxu3bPOWIpI0g6nqqkG2YEmSRoj/NKtX3lJEknbfnExyIUkaLicumFySvZJcneS6JDcmeVdLPzTJVUk2JflUkj1b+iuS3J3k2rac1nWsU5Pc2pZTu9JXJbmhHevMNm5YkjRPGWBJ0gLgxAVTegA4rqqeChxF50b2a4D3AR+sqsOBbcCru/b5VFUd1ZZzAJLsB7wDOAY4GnhHm5gJ4KPAa4Aj2nLCHJyXJGlIDLAkaQGYmLhgcXDigi7VcX97uEdbCjgOuLilnwe8aIZDPR+4oqruqaptwBV0grUDgH2qan11Bj2f38OxJEljzDFYkrQArFq5lAtOW+MYrEm02W43AocDZwHfAe6tqu0ty+3AgV27/GqSZ9O5zchbquq2tv22rjwT+xzY1ndO37kMpwOnA6xYsaIPZyVJGhZbsCRpgVi1cilnHHu4wdVOqurhqjqKzg3rjwaeNE32vwQOqaqfp9NKdV6fynB2Va2uqtXLly/vxyElSUNigDViNm7ZxlnrNjkAXZLmWFXdC6wDngHsm2Sil8dBwB0tz9aqeqClnwOsaut3AAd3HW5inzva+s7pkqR5ygBrhDjLlyTNrSTLk+zb1vcGngfcTCfQOqllOxX4XMtzQNfuv9LyAlwOHJ9kaZvc4njg8qq6E7gvyZo2e+DLJ44lSZqfHIM1Qiab5cuuPJI0UAcA57VxWIuAi6rq80luAi5M8h7gG8DHW/43JvkVYDtwD/AKgKq6J8kfANe0fO+uqnva+uuBTwB7A5e1RZI0TxlgjZCJWb4e2r7DWb4kaQ5U1fXA0yZJ30xnPNbO6W8H3j7Fsc4Fzp0kfQPwlFkXVpI0FgywRoizfEmSJEnjzQBrxKxaudTASpIkSRpTTnIhSZIkSX1igCVJkiRJfWKAJUmSJEl9YoAlSZIkSX1igCVJkiRJfWKAJUmSJEl9YoAlSfPAxi3bOGvdJjZu2TbsomgB8PMmSVPzPliSNOY2btnGKees58HtO9hzySIuOG2N99PTwPh5k6Tp2YIlSWNu/eatPLh9BzsKHtq+g/Wbtw67SJrH/LxJ0vQMsCRpzK05bBl7LlnE4sAeSxax5rBlwy6S5jE/b5I0PbsIStKYW7VyKRectob1m7ey5rBldtfSQPl5k6TpGWBJ0jywauVS/9HVnPHzJklTs4vgGHL2JkmSJGk09dyCleSnquqHgyyMZubsTZI0NesqSdKwzdiCleSZSW4CvtUePzXJRwZeMk3K2Zsk6ZGsqyRJo6KXLoIfBJ4PbAWoquuAZw+yUJqaszdJ0qSsqyRJI6GnLoJVdVuS7qSHB1MczcTZmyRpctZVkqRR0EuAdVuSZwKVZA/gTcDNgy2WpuPsTZL0CNZVkqSR0EsXwdcCZwAHAncAR7XHkiSNCusqSdJImLEFq6r+GThlDsoiSdJusa6SJI2KGQOsJH8K1M7pVfWqgZRIkqRdZF0lSRoVvXQR/DzwV235IrAPcP9MOyXZK8nVSa5LcmOSd+20/cwk93c9flSSTyXZlOSqJId0bXt7S78lyfO70k9oaZuSvK2Hc5EkzU+7VVdJktRvvXQR/HT34yR/Dnylh2M/ABxXVfe3AcdfSXJZVa1PshrYeZaGVwPbqurwJCcD7wN+PcmRwMnAk4EnAH+T5Iltn7OA5wG3A9ckWVtVN/VQNknSPDKLukqSpL7qpQVrZ0cAj5spU3VMXD3coy2VZDHwx8B/2mmXE4Hz2vrFwC+lM9/uicCFVfVAVf0DsAk4ui2bqmpzVT0IXNjyStK8s3HLNs5at4mNW7YNuyjjoqe6aqreFkkObb0pNrXeFXu2dHtbSJKm1csYrB/Q6dee9vf/AG/t5eAtmNoIHA6cVVVXJXkTsLaq7tzpfiUHArcBVNX2JN8HlrX09V35bm9pTOTvSj9minKcDpwOsGLFil6KLkkjY+OWbZxyznoe3L6DPZcs4oLT1nirhp3Moq6atLcF8FvAB6vqwiQfo9PL4qPY20KSNIMZW7Cq6jFVtU/X3yfu3BVjmn0frqqjgIOAo5M8G3gx8CezK/auqaqzq2p1Va1evnz5XD61JM3a+s1beXD7DnYUPLR9B+s3bx12kUbO7tZVU/W2AI6j05sCOr0rXtTW7W0hSZrWlC1YSZ4+3Y5V9fVen6Sq7k2yDjiWTmvWptZ69VNJNlXV4XTuW3IwcHuSJcBjga1d6RMOamlMky5J88aaw5ax55JFPLR9B3ssWcSaw5YNu0gjox911c69LYDvAPdW1faWpbvnxEB6W8zHnhYbt2xj/eatrDlsmS2ukhaU6boIvn+abRNX96aUZDnwUAuu9qbTPeJ9VfXTXXnub8EVwFrgVOBrwEnAl6qqkqwF/izJB+h0uzgCuJpON5AjkhxKJ7A6GXjpdGWSpHG0auVSLjhtjf+sTm5WdRV0elsARyXZF7gEeFKfytazqjobOBtg9erVj5huftzYrVXSQjZlgFVVx87y2AcA57Urg4uAi6rq89Pk/zjwySSbgHvoBExU1Y1JLgJuArYDZ7TKkCRvAC4HFgPnVtWNsyzzvOBVQ2n+WbVyqd/nSfShruo+1kRvi2cA+yZZ0lqxuntI2NuiB5N1a/XzK2mhmHGSC4AkTwGOBPaaSKuq86fbp6quB542Q55Hd63/iM74rMnyvRd47yTplwKXTvccC41XDSUtVLtTV03V2wJYR6c3xYV0eld8ru1ib4se2K1V0kLWyyyC7wCeQ6fSuhR4AZ17i0xbaWk4vGooaSGaRV01aW+LJDcBFyZ5D/ANOr0swN4WPbFbq6SFrJcWrJOApwLfqKpXJnk88L8HWyztLq8aSlqgdquumqq3RVVtpjMD4M7p9rbokd1aJS1UvQRY/1JVO5JsT7IPcBc/2Z9cI8SrhpIWKOsqSdJI6CXA2tBmVvpfdKaxvZ9O33ONKK8aSlqArKskSSNhxgCrql7fVj+W5K+BfVqXCkmSRoJ1lSRpVCyaKUOStUlemuTfVNV3rbAkSaPGukqSNCpmDLDo3MTxF4Cbklyc5KQke820kyRJc8i6SpI0EnrpIvi3wN+2KWyPA14DnAvsM+CySZLUE+sqSdKo6PVGw3sD/w/w68DTgfMGWShJknaVdZUkaRT0cqPhi+jcC+SvgQ8Df1tVOwZdMElaSDZu2ebtFWbBukqSNCp6acH6OPCSiTvSS5L6a+OWbZxyznoe3L6DPZcs4oLT1hhk7TrrKknSSJhxkouqutwKS5IGZ/3mrTy4fQc7Ch7avoP1m7cOu0hjx7pKkjQqeplFUJI0QGsOW8aeSxaxOLDHkkWsOWzZsIskSZJ2U0+TXEiSBmfVyqVccNoax2BJkjQPTBlgJVkx3Y5V9Y/9L47migPqpdGyauVSv4u7wbpqfFkPSZqvpmvB+iuggHSlFbAceByweIDl0gA5oF7SPGJdNYashyTNZ1OOwaqqn6uqn29/f47OvUW+CtwPvHmuCqj+c0C9pPnCumo8WQ9Jms9mnOQiyRFJPgFcBmwEjqyqPxl0wTQ4DqiXNN9YV40X6yFJ89l0Y7CeAvwe8GTgj4BXOwXu/OCAeknzhXXVeLIekjSfTTcG6zrgNjr9248Gjk5+3MW9qt442KJpkBxQL2mesK4aU9ZDkuar6QKsV9MZKCxJ0qiyrpIkjZTpAqwLgcdU1d3diUmWAz8YaKkkSeqNdZUkaaRMN8nFmcD/NUn6LwAfHExxJEnaJdZVkqSRMl2AtaqqPrNzYlVdAjx7cEWSpPll45ZtnLVuExu3bBt2UeYj6ypJ0kiZrovgT02zbcbp3SVJ3lB1DlhXSZJGynSVz11Jjt45Mcm/A+6eJL8kaSfeUHXgrKskSSNluhas3wUuajdu3NjSVgMvB04ecLkkaV6YuKHqQ9t3eEPVwbCukiSNlClbsKrqajr3FAnwirYEOKaqrpqLwknSuJu4oepvHf8zdg8cgNnUVUkOTrIuyU1Jbkzyppb+1CRfS3JDkr9Msk9LPyTJvyS5ti0f6zrWqpZ/U5Iz027GlWS/JFckubX99QMgSfPclC1YSVZU1T8C75jD8mhEbNyyjfWbt7LmsGX+QyjNkjdUHZxZ1lXbgd+uqq8neQywMckVwDnA71TV3yZ5FZ1Wsv/S9vlOVR01ybE+CrwGuAq4FDgBuAx4G/DFqvrDJG9rj9+6G2WVJI2J6cZgfXZiJcmn56AsGhETg/Lf/4VbOOWc9c58JmmU7XZdVVV3VtXX2/oPgJuBA4EnAl9u2a4AfnW64yQ5ANinqtZXVQHnAy9qm08Ezmvr53WlS5LmqekCrHStHzbogmh0OChf0hjpS12V5BDgaXRaoG6kExgBvBg4uCvroUm+keRvk0zcf+tA4PauPLe3NIDHV9Wdbf3/AI+f4vlPT7IhyYa773Zujgne4kDSOJpukouaYl3znIPyJY2RWddVSR4NfBp4c1Xd17oFnpnkvwBrgQdb1juBFVW1Nckq4LNJntxzQasqyaRlrKqzgbMBVq9ebZ2LtziQNL6mC7CemuQ+OlcH927rtMdVVfsMvHQaiolB+Y7BkjQGZlVXJdmDTnB1wcQNi6vqW8DxbfsTgV9u6Q8AD7T1jUm+Q6c74R3AQV2HPailAfxTkgOq6s7WlfCu2Z7wQjFZbwrrI0n9MOi5BqabRXBxVe1TVY+pqiVtfeLxjMFVkr2SXJ3kujY707ta+gVJbknyzSTntsqNdJzZZmC6PsnTu451apuB6dYkp3alTzprk2Zv1cqlnHHs4VZmkkbabOqqVmd8HLi5qj7Qlf649ncR8PvAx9rj5UkWt/XDgCOAza0L4H1J1rRjvhz4XDvcWmCi3jq1K10zmOhNsTjYm0JS38zFXAODvMv9A8BxVfVU4CjghCRrgAuAJwE/B+wNnNbyv4BOZXUEcDqdGZlIsh+d2aGOoTMV7zu6prmdmLVpYr8TBng+kqT55VnAy4DjuqZefyHwkiTfBr4FfA/405b/2cD1Sa4FLgZeW1X3tG2vpzP74CbgO3RmEAT4Q+B5SW4Fntseqwfe4kDSIMzFXAPTdRGclTaT0v3t4R5tqaq6dCJPkqv5cbeKE4Hz237rk+zbulM8B7hiohJrU+iekORK2qxNLX1i1qaJSk2SpClV1Vf4yUkyun1okvyfptOdcLJjbQCeMkn6VuCXZlHMBc1bHEjqt7mYa2BgARZA60qxETgcOKv7po+ta+DLgDe1pAOB27p2n5iFabr0qWZt2rkcp9NpFWPFihW7f0KSJEmSxtZczDUw0ACrqh4GjkqyL3BJkqdU1Tfb5o8AX66qvxtkGVo5nJ1J0sB4Y25JksbHoFvHBxpgTaiqe5OsozNG6ptJ3gEsB/5jV7Y7+Ml7jUzMwnQHnW6C3elXMv2sTZI0J5xKWpIkdRvYJBdttqV92/rewPOAbyU5DXg+8JKq2tG1y1rg5W02wTXA99vMTJcDxydZ2ia3OB64fIZZmyRpTnhjbkmSRsuwb1I+yBasA4Dz2jisRcBFVfX5JNuBLcDX2qzqn6mqdwOXAi+kMwPTD4FXAlTVPUn+ALimHffdO83a9Ak6sxFehhNcSJpj3phbkqTRMQo9SwY5i+D1wNMmSZ/0OdvsgWdMse1c4NxJ0iedtUlzw3EnkjfmliRplIzCTcrnZAyW5p9RuDogjQqnkpYkaTSMQs8SAyztllG4OiBJkiR1G4WeJQZY2i2jcHVAkiRJ2tmwe5YYYGm3jMLVAUnSwuZYYGnhGYfvvQGWdtuwrw5IkhYuxwJLC8+4fO8Hdh8sSZKkQfEedNLCMy7fewMsSZI0dibGAi8OjgWWFohx+d7bRVCSpjEOfb2lhcixwNLCMy7fewMsSZrCuPT1lhYqxwJLC884fO/tIihJUxiXvt6SJM0XG7ds46x1m9i4Zduwi7LbbMGSpCl4vzdJkubOfOk5YoClgXHsisbduPT1liRpPpis58g41r0GWBqI+XIFQhqHvt6SJM0H86XniAGWBmK+XIGQJEnS3JgvPUcMsDQQ8+UKhCRJkubOfOg5YoClgZgvVyAkSePN8cDS6Fgo30cDLA3MfLgCIUkaX44HlkbHQvo+eh8sSZI0L3kvO2l0LKTvowGWpAVtPtzQUNLkJsYDLw6OB5aGbCF9H+0iKGnBWkjdFTS5JAcD5wOPBwo4u6o+lOSpwMeARwPfBU6pqvvaPm8HXg08DLyxqi5v6ScAHwIWA+dU1R+29EOBC4FlwEbgZVX14Jyd5ALmeGBpdCyk76MtWJIWrIXUXUFT2g78dlUdCawBzkhyJHAO8Laq+jngEuB3Adq2k4EnAycAH0myOMli4CzgBcCRwEtaXoD3AR+sqsOBbXSCM82RVSuXcsaxh8/rf+akUTFTr5CF8n00wJK0YC2k7gqaXFXdWVVfb+s/AG4GDgSeCHy5ZbsC+NW2fiJwYVU9UFX/AGwCjm7Lpqra3FqnLgROTBLgOODitv95wIsGf2aSNLcmeoW8/wu3cMo56xd013u7CGqoFsp0nRpNC6m7gmaW5BDgacBVwI10gqnPAi8GDm7ZDgTWd+12e0sDuG2n9GPodAu8t6q2T5K/+7lPB04HWLFixazPRZLm2mS9QhZqvWqApaFx/ItGgbcTEECSRwOfBt5cVfcleRVwZpL/AqwFBjpmqqrOBs4GWL16dQ3yuSRpECZ6hTy0fceC7xVigKWh8UqHpFGQZA86wdUFVfUZgKr6FnB82/5E4Jdb9jv4cWsWwEEtjSnStwL7JlnSWrG680vSvGGvkB8zwNLQeKVD0rC1MVIfB26uqg90pT+uqu5Ksgj4fTozCkKnNevPknwAeAJwBHA1EOCINmPgHXQmwnhpVVWSdcBJdMZlnQp8bm7OTlA1hQkAABkySURBVJLmlr1COgywNDRe6ZA0Ap4FvAy4Icm1Le0/0wmWzmiPPwP8KUBV3ZjkIuAmOjMQnlFVDwMkeQNwOZ1p2s+tqhvb/m8FLkzyHuAbdAI6jQjHAku98bvSu1QtrK7eq1evrg0bNgy7GJKkaSTZWFWrh12OYbCemjuOBZZ643dlclPVVU7TLkmSFiTvhSf1xu/KrjHAkjRvzXTDQ0kLm/fCk3rjd2XXOAZL0rxkdwZJM3EssNQbvyu7xgBL0rzkbQAk9cJZz6SOmSax8LvSOwMsjTRnrNHu8jYAkiT1xl4f/TWwMVhJ9kpydZLrktyY5F0t/dAkVyXZlORTSfZs6Y9qjze17Yd0HevtLf2WJM/vSj+hpW1K8rZBnYuGY+LL/v4v3MIp56x3HI12yUR3ht86/mesKCRJmoaTWPTXICe5eAA4rqqeChwFnJBkDfA+4INVdTiwDXh1y/9qYFtL/2DLR5Ij6dyw8cnACcBHkixOshg4C3gBcCTwkpZX84Rfds3WqpVLOePYww2uJEmahpNY9NfAughW5wZb97eHe7SlgOOAl7b084B3Ah8FTmzrABcDH06Sln5hVT0A/EOSTcDRLd+mqtoMkOTClvemQZ2T5pZdvCRJw2ZXdS0ETmLRXwMdg9VamTYCh9NpbfoOcG9VbW9ZbgcObOsHArcBVNX2JN8HlrX09V2H7d7ntp3Sj5miHKcDpwOsWLFidielOeOXXZI0TI5L0XzRy4UCJ7Hon4EGWFX1MHBUkn2BS4AnDfL5pinH2cDZAKtXr65hlEG7xy+7JGlYnI1U84EXCubenNxouKruBdYBzwD2TTIR2B0E3NHW7wAOBmjbHwts7U7faZ+p0iUtEN5IWNIgOS5F84Fj2ufewFqwkiwHHqqqe5PsDTyPzsQV64CTgAuBU4HPtV3Wtsdfa9u/VFWVZC3wZ0k+ADwBOAK4GghwRJJD6QRWJ/PjsV2S5jmvyEkaNLuqaz5wTPvcG2QXwQOA89o4rEXARVX1+SQ3ARcmeQ/wDeDjLf/HgU+2SSzuoRMwUVU3JrmIzuQV24EzWtdDkrwBuBxYDJxbVTcO8Hw0ohyAvDDZdUfSXLCrusadFwrm3iBnEbweeNok6Zv58SyA3ek/Al48xbHeC7x3kvRLgUtnXViNLVsxFi6vyEmS1DHTxWYvFMytgU5yIQ2arRgLl1fkJEnyYvMoMsDSWLMVY2HzipykUWBXdQ2TF5tHjwGWxpqtGJKkYbL1QMPmxebRY4ClsWcrhiRpWGw90LB5sXn0GGBJGkl2uZE0Dmw90KD1Uh96sXm0GGBJGjl2uZE0Lmw90CBZH46nRcMugDRoG7ds46x1m9i4Zduwi6Ieedd5SeNk1cqlnHHs4f7jq76zPhxPtmBpXvPKz3iyy40kSdaH48oAS/Oag4/Hk11uJM0njinVdKb7fFgfjicDLM1rXvkZXw7YlTQf2JNC0+nl82F9OH4cg6V5beLKz28d/zNWaiPGsXEaBUkOTrIuyU1JbkzyppZ+VJL1Sa5NsiHJ0S39OUm+39KvTfJfu451QpJbkmxK8rau9EOTXNXSP5Vkz7k/Uw2LY2g0HT8f85MtWJr3vPIzeryiqxGyHfjtqvp6kscAG5NcAfwR8K6quizJC9vj57R9/q6q/u/ugyRZDJwFPA+4Hbgmydqqugl4H/DBqrowyceAVwMfnYuT0/DZk0LT8fMxPxlgSdg/fq45Nk6joqruBO5s6z9IcjNwIFDAPi3bY4HvzXCoo4FNVbUZIMmFwInteMcBL235zgPeiQHWguEYmoVtpv8v/HzMTwZYWvBsTZl7XrHTKEpyCPA04CrgzcDlSf4Hne70z+zK+owk19EJun6nqm6kE5Td1pXnduAYYBlwb1Vt70o/cICnoRFkT4qFqdf/L/x8zD+OwdKCZ//nuefYOI2aJI8GPg28uaruA14HvKWqDgbeAny8Zf06sLKqngr8CfDZPj3/6W2s14a77767H4fUmHFc6vzj/xcLly1YWvBsTRkOr9hpVCTZg05wdUFVfaYlnwq8qa3/BXAOQAu+aOuXJvlIkv2BO4CDuw57UEvbCuybZElrxZpI/wlVdTZwNsDq1aurj6enMWBPivE1XRdA/79YuAywtODZ/7n/HNOmcZEkdFqnbq6qD3Rt+h7wi8CVdMZQ3dry/zTwT1VVbWbBRXSCqHuBI5IcSieAOhl4acu3DjgJuJBO4Pa5uTg3jQ/HpY6nmQJj/79YuAywJGxN6SevxGrMPAt4GXBDkmtb2n8GXgN8KMkS4EfA6W3bScDrkmwH/gU4uaoK2J7kDcDlwGLg3DY2C+CtwIVJ3gN8gx93N5QAWzrGVS+Bsf9fLEwGWFIPbJHpnVdiNU6q6itApti8apL8HwY+PMWxLgUunSR9M51ZBqVJ9dLSYT00egyMNRUDLGkGtsjsGiscSdp107V0WA8Nz3SBrV0ANRUDLGkGtsjsGiscSeov66Hh6CWwtQugJmOAJc3AFplH6uXGiVY4ktQf1kODMVNdZmCr3WWAJc2g1xaZhdI/3q4qkjS3rIf6r5e6zMBWu8sAS+rBTC0yCyno8IqeJM0966H+6nUGQLu8a3csGnYBpPlgPt2tfeOWbZy1bhMbt2ybdPvEFb3FwSt6kjQi5lM91C/T1We91mWrVi7ljGMPN7jSLrEFS+qDXroRjEPXjV4H9HpFT5JGy3yph/rFmwBrmAywpD6Y6Yd6XLpu9Nr9z0ksJGm0zJd6qFf9mKDCukyDYoAl9cl0P9S9Bi5zcXVxuudwQK8kja9xqYd6MV05nKBCo84AS5oDvXbdmKnCmKni62W7XSYkaeGZq3qolzyzraucoEKjzgBLmgO9/NDPVGHMVOH0UjHaZUKSFqa5qId6ydOPuqrX1inrMw2LAZY0R2b6oZ+pwpipwukleLLLhCQtXIOuh3rJ04+6ytYpjToDLGlEzFRhzFTh9BI8WSlJkqYy23qolzz9qqtsndIoS1UN5sDJwcD5wOOBAs6uqg8lOQr4GLAXsB14fVVdnSTAh4AXAj8EXlFVX2/HOhX4/Xbo91TVeS19FfAJYG/gUuBNNcMJrV69ujZs2NDXc5Xmymz7tUvjIsnGqlo97HIMg/WURtlcjMGSxsVUddUgA6wDgAOq6utJHgNsBF4E/L/AB6vqsiQvBP5TVT2nrf8mnQDrGOBDVXVMkv2ADcBqOoHaRmBVVW1LcjXwRuAqOgHWmVV12XTlsuKSpNFngGU9JUmjbqq6atGgnrCq7pxogaqqHwA3AwfSCZL2adkeC3yvrZ8InF8d64F9W5D2fOCKqrqnqrYBVwAntG37VNX61mp1Pp0ATpIkSZKGYk7GYCU5BHganZamNwOXJ/kfdAK8Z7ZsBwK3de12e0ubLv32SdIlSZIkaSgG1oI1IcmjgU8Db66q+4DXAW+pqoOBtwAfn4MynJ5kQ5INd99996CfTpIkSdICNdAAK8kedIKrC6rqMy35VGBi/S+Ao9v6HcDBXbsf1NKmSz9okvRHqKqzq2p1Va1evnz57p+QJEmSJE1jYAFWmxXw48DNVfWBrk3fA36xrR8H3NrW1wIvT8ca4PtVdSdwOXB8kqVJlgLHA5e3bfclWdOe6+XA5wZ1PpIkSZI0k0GOwXoW8DLghiTXtrT/DLwG+FCSJcCPgNPbtkvpzCC4ic407a8EqKp7kvwBcE3L9+6quqetv54fT9N+WVskSZIkaSgGFmBV1VeATLF51ST5CzhjimOdC5w7SfoG4CmzKKYkSZIk9c3A7oM1qpLcDWyZxSH2B/65T8UZlHEoI4xHOS1j/4xDOS1jf/SjjCurakEOmu1DPQXj8TnZXZ7bePLcxpPnNr1J66oFF2DNVpINo37zy3EoI4xHOS1j/4xDOS1jf4xDGee7+fweeG7jyXMbT57b7hn4NO2SJEmStFAYYEmSJElSnxhg7bqzh12AHoxDGWE8ymkZ+2ccymkZ+2Mcyjjfzef3wHMbT57bePLcdoNjsCRJkiSpT2zBkiRJkqQ+McCSJEmSpD4xwJpEkhcnuTHJjiRTTt+Y5IQktyTZlORtXemHJrmqpX8qyZ4DKON+Sa5Icmv7u3SSPMcmubZr+VGSF7Vtn0jyD13bjup3GXstZ8v3cFdZ1nalj8preVSSr7XPxfVJfr1r28Bey6k+Y13bH9Vel03tdTqka9vbW/otSZ7frzLtRhl/K8lN7XX7YpKVXdsmfd+HUMZXJLm7qyyndW07tX02bk1y6qDK2GM5P9hVxm8nubdr28BfyyTnJrkryTen2J4kZ7byX5/k6V3b5ux1XEhm85kZdT2c24ok65J8o33eXjiMcu6OHs5tZfu9vD7JlUkOGkY5d9VsfiNGXQ/n9qT2f8IDSX5nrss3Gz2c2ynt/bohyd8neepcl3F39XBuJ7ZzuzbJhiS/0JcnriqXnRbgZ4GfAa4EVk+RZzHwHeAwYE/gOuDItu0i4OS2/jHgdQMo4x8Bb2vrbwPeN0P+/YB7gJ9qjz8BnDQHr2VP5QTunyJ9JF5L4InAEW39CcCdwL6DfC2n+4x15Xk98LG2fjLwqbZ+ZMv/KODQdpzFQyrjsV2fu9dNlHG6930IZXwF8OFJ9t0P2Nz+Lm3rS4dVzp3y/yZw7hy/ls8Gng58c4rtLwQuAwKsAa6a69dxIS2z/cyM8tLj9/bsiTqh/eZ9d9jl7uO5/QVwals/DvjksMvd47nt1m/EOCw9nNvjgH8HvBf4nWGXt8/n9syJ32zgBfPsfXs0P56T4ueBb/XjeW3BmkRV3VxVt8yQ7WhgU1VtrqoHgQuBE5OEzo/hxS3fecCLBlDME9uxe32Ok4DLquqHAyjLdHa1nP9qlF7Lqvp2Vd3a1r8H3AU84s7dfTbpZ2ynPN1lvxj4pfa6nQhcWFUPVNU/AJva8ea8jFW1rutztx6Y6yuxvbyOU3k+cEVV3VNV24ArgBNGpJwvAf58QGWZVFV9mc6FmqmcCJxfHeuBfZMcwNy+jgvJyH9mZqGXcytgn7b+WOB7c1i+2ejl3I4EvtTW102yfSTN4jdi5M10blV1V1VdAzw0d6Xqjx7O7e/bbzcMpx7fbT2c2/3Voivg39D5XZk1A6zddyBwW9fj21vaMuDeqtq+U3q/Pb6q7mzr/wd4/Az5T+aRFet7W7PoB5M8qu8l7Oi1nHu1ptn1ad0YGdHXMsnRdK46fqcreRCv5VSfsUnztNfp+3Ret172nasydns1nauXEyZ73/ut1zL+ansPL05y8C7u2w89P1c63SwP5cf/gMHcvJYzmeoc5vJ1XEhm+5kZZb2c2zuB30hyO3ApnRa6cdDLuV0H/Ie2/u+BxyRZNgdlGzR/C8bfzvX42Evy75N8C/gr4FX9OOaSfhxkHCX5G+CnJ9n0e1X1ubkuz2SmK2P3g6qqJFNG3O3q0M8Bl3clv51OMLEnnW4WbwXePcRyrqyqO5IcBnwpyQ10goW+6PNr+Uk6XTd2tOS+vZbzWZLfAFYDv9iV/Ij3vaq+M/kRBuovgT+vqgeS/Ec6rYLHDaEcvToZuLiqHu5KG5XXUqNpss/MuHsJ8Imqen+SZwCfTPKUrt/mcfY7wIeTvAL4MnAHMJ/eO42hJMfSCbD6M05pRFTVJcAlSZ4N/AHw3Nkec8EGWFU12xfvDuDgrscHtbStdJq8l7QWhYn0vpYxyT8lOaCq7mz/9N81zaF+Dbikqv612bqrxeaBJH9K58d8t/SjnFV1R/u7OcmVwNOATzNCr2WSfehc3fi91rVh4th9ey13MtVnbLI8tydZQqebzNYe952rMpLkuXSC2V+sqgcm0qd43/sdFMxYxqra2vXwHDrj8ib2fc5O+17Z5/JN2JX37GTgjO6EOXotZzLVOczl67iQzOozM+J6ObdX07qaVtXXkuwF7M/09eEo6OU36Xu0FqwkjwZ+tarGZoKSacxV3aQ+S/LzdOrHF+xUZ84bVfXlJIcl2b+q/nk2x7KL4O67BjginVnu9qRTea1t/TjX0RnzBHAqMIgWsbXt2L08xyP63U/0eW7jdV4ETDq7Sh/MWM4kSye61SXZH3gWcNMovZbtPb6ETt/xi3faNqjXctLP2DRlPwn4Unvd1gInpzPL4KHAEcDVfSrXLpUxydOA/wn8SlXd1ZU+6fs+pDJ2jwH4FeDmtn45cHwr61LgeH6yJXhOy9nK+iQ6E0V8rSttrl7LmawFXp6ONcD32wWIuXwdF5Ld/syMgV7O7R+BXwJI8rPAXsDdc1rK3dPLb9L+SSb+R3s7cO4cl3FQpvqN0AhLsgL4DPCyqvr2sMvTT0kOb/+/kc6slo+ic6F6dqabAWOhLnT6O98OPAD8E3B5S38CcGlXvhcC36Zzlfj3utIPo/PP7CY6MwE9agBlXAZ8EbgV+Btgv5a+GjinK98hdK4OLdpp/y8BN9AJBv438OgBvZYzlpPO7DQ30OlzfgPw6lF7LYHfoDNw9dqu5ahBv5aTfcbodD/8lba+V3tdNrXX6bCufX+v7XcLnStOg/q+zFTGv2nfo4nXbe1M7/sQyvjfgRtbWdYBT+ra91Xt9d0EvHJQZeylnO3xO4E/3Gm/OXkt6VyoubN9F26n04LwWuC1bXuAs1r5b6BrFta5fB0X0rK7n5lxWHr43h4JfLV97q8Fjh92mft4bifRqZO+TafVoO9134DOa7d/I0Z96eHcfrql3wfc29b3GXa5+3Ru5wDb+HE9vmHYZe7jub2VTv1/LZ2LUL/Qj+edmJZQkiRJkjRLdhGUJEmSpD4xwJIkSZKkPjHAkiRJkqQ+McCSJEmSpD4xwJIkSZKkPjHAkuaJJPsmef2wyyFJmltJ3pzkp/p4vO+2e+rt7v6vSPLhQT5Pkr+fYftP1IlJnpDk4un2kfrFAEuaP/YFdinAajd79HdAksbbm4G+BVi7KsniuX7OqnrmDFl+ok6squ9V1UmDLZXU4T9W0ohI8vIk1ye5LsknkyxP8ukk17TlWS3fO5Ocm+TKJJuTvLEd4g+Bf5vk2iR/3PL+btv3+iTvammHJLklyfl0bo588DDOV5K0a5L8myR/1eqJbyb59VYHPAFYl2Rdy/fRJBuS3Djx29/Sv5vkXUm+nuSGJE9q6cuSfKHlP4fODYEn9vlsko1t2+ld6fcneX+S64BnJHllkm8nuRp41hTln+55fiPJ1a0O+59JFid57UR91vL8a8tYkvvb30cn+WLXOZ3Ysv9Endjqvm+2ffZK8qct/zeSHNt1/M8k+esktyb5o91/t7SgDfsOyy4uLgXwZODbwP7t8X7An9HuKA6sAG5u6+8E/h54FLA/sBXYAzgE+GbXMY8HzqZTgS0CPg88u+XbAawZ9nm7uLi4uPS+AL8K/K+ux49tf787UX+0x/u1v4uBK4Gf78r3m2399cA5bf1M4L+29V8Gqrs+an/3pnNRbll7XMCvtfUDgH8ElgN7Al8FPjxJ+Sd9HuBngb8E9mjbPgK8vB1vU9f+l3XVi/e3v0uAfdr6/sCmVu/tXCf+62Pgt4Fz2/qTWtn3Al4BbAYe2x5vAQ4e9vvuMn7LEiSNguOAv6iqfwaoqnuSPBc4MvnXC3z7JHl0W/+rqnoAeCDJXcDjJznm8W35Rnv8aOAIOhXJlqpaP5hTkSQNyA3A+5O8D/h8Vf3dFPl+rbU2LaET/BwJXN+2fab93Qj8h7b+7In1qvqrJNu6jvXGJP++rR9Mpx7ZCjwMfLqlHwNcWVV3AyT5FPDESco11fP8ErAKuKbVeXsDd1XV3a2nxhrgVjrB0Fd3OmaA/5bk2XQuHh7I5HVit18A/qSV41tJtnSV94tV9f12HjcBK4HbZjie9BMMsKTRtYhOK9OPuhNb5fNAV9LDTP5dDvDfq+p/7rT/IcD/18+CSpIGr6q+neTpwAuB9yT5YlW9uztPkkOB3wH+XVVtS/IJOq0xEybqj6nqju5jPQd4LvCMqvphkiu7jvWjqnp4lqf0r08FnFdVb59k24XArwHfAi6pqtpp+yl0WrpWVdVDSb7LT57vruqlfpWm5RgsaTR8CXhxkmUASfYDvgD85kSGJEfNcIwfAI/penw58KqJVq8kByZ5XF9LLUmaM0meAPywqv438MfA09um7t//fehcRPt+kscDL+jh0F8GXtqe4wXA0pb+WGBbC66eBKyZYv+rgF9sY6z2AF68i8/zReCkiToqyX5JVrZtlwAnAi+hE2zt7LF0WrseamOpJvbbuU7s9nd0AjOSPJFON/xbpsgr7TKjcmkEVNWNSd4L/G2Sh+l063sjcFaS6+l8V78MvHaaY2xN8tU2iPeyqvrdJD8LfK21et0P/AadK3KSpPHzc8AfJ9kBPAS8rqWfDfx1ku9V1bFJvkGnxec2HtmlbjLvAv48yY10xvj+Y0v/a+C1SW6mE4BM2rW8qu5M8k7ga8C9wLW78jxVdVOS3we+kM7Mtg8BZ9Dpzr6tPf+RVXX1JMe8APjLJDcAG9p5P6JOBM7q2ucjwEfbPtuBV1TVA11d8qVZySNbWiVJkiRJu8MugpIkSZLUJwZYkiRJktQnBliSJEmS1CcGWJIkSZLUJwZYkiRJktQnBliSJEmS1CcGWJIkSZLUJ/8/vNOOO8vCO0oAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"\n",
"c_scan = np.linspace(*c.bounds, 40)\n",
"s_scan = np.linspace(*s.bounds, 40)\n",
"\n",
"with minkit.minimizer('uml', g, data) as minimizer:\n",
" c_fcn = minimizer.fcn_profile('c', c_scan)\n",
" s_fcn = minimizer.fcn_profile('s', s_scan)\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 4))\n",
"\n",
"ax0.plot(c_scan, c_fcn, '.');\n",
"ax0.set_xlabel('center');\n",
"ax0.set_ylabel('FCN value');\n",
"ax1.plot(s_scan, s_fcn, '.');\n",
"ax1.set_xlabel('standard deviation');\n",
"ax1.set_ylabel('FCN value');\n",
"\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In case we want to map the two-dimensional parameter space, we just need to specify the two parameters as input variables and define the values as a grid."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAFBCAYAAACVRMOtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7gkVXnv8e9v75m9hzsDgwYZFIgYRKOoBPHyKKLCgBfweAkcPRCDmSgaTXI0wqMRgyFKTk5IiBdClABqQCR6nCiIKBBijijjERFQZAAvg+gIw1XCMLP3e/6otaUZaq3urunu3b3793meenbvtaqqV1d39eqqeutdigjMzGy8Tcx3A8zMbP65MzAzM3cGZmbmzsDMzHBnYGZmuDMwMzP62BlIOkvSOknXZeqPkHStpGskrZb0/H61xczMytSv+wwkvQC4Hzg3Ip5aU78t8KuICElPAy6IiH360hgzMyvq25FBRFwJrC/U3x8P90TbAL77zcxsnszrNQNJr5L0A+BLwO/PZ1vMzMZZ304TAUjaA/hi3WmizeZ7AfC+iHhJpn4lsBJgYtHUs5Zs/5j69eReSuE1arbQsMxiKm2z0uYstaNB28t1XVe0OTbLVDZZpp1BHSOqxwsW11eobLJcaRmVnqu+Lhqvr/tlys+Vr4qJ7tv+wPq1d0TELoVnbOvQF20Td66f6WqZb1+74ZKIWLElzztoi+a7AVCdUpK0l6RlEXFHTf2ZwJkA2+y8ezzlZX9cu57JDfXfIpMb898uE5llACYfqu8pJjLlABMb8x8aNajTpkJvtXFTfn2bMs81W1hfbhmAmfq64o+JzDJpwXzd7IB6g8yXC9DmC7D+gFqThQPtyclCO0rL1dfFosL6Cs8Vi+t3+VhcWibfvtxys1P5ZWYKdeXl6t+Tmen8e/WtT7/zx9nKDt25foZvXfL4rpaZ3PWmZVv6vIM2b52BpCcCN6cLyM8EpoE756s9ZmZ1ApildAphYehbZyDpPOAgYJmktcBJwGKAiDgDeDVwjKSNwH8BvxtOoWpmQyeYCXcGjUXE0W3qTwVO7dfzm5n1QnVksPB/pw7FNQMzs2Hm00RmZmMuCGbG4Az2wuoMshEgI/xG9vpDWIrUKUQaZS/nlNbX64ihXp+3nS1E8RTvwKlvRxRiHFWK4ipFLvVYNiy68F4Nzd4zuM30KD5NZGY25gKYcWdgZmY+MjAzG3MBvmZgZma5q0QLizsDM7OCIHzNwMxs7AXMLPy+YAQ7A0FkMy82eMfmMVztEZochzbJaNrrc5+lcM9i2GmDF9zrBHYTpbY3CDstta/02WzwPqrwXDEx5N9cPQ6lLWZB7cX68WkiMzNDzAzNr8b+cWdgZlYQDC6r+nxyZ2Bm1sY4HBnM67CXZmY2HHxkYGZWUKWjWPhHBqPZGXT5vuSij1LtFjWlK01CEpqerGwSTdQksVxxfQ0jjXKr63EkVHH86yaRRoUVFpPYFZqR3U4DjBgqBkL1+LmicK4iuwkH8D092++QpSHg00RmZgVzRwbdTJ2Q9CNJ35N0jaTVqWwnSZdKuin9XZrKJel0SWskXZuGCp5bz7Fp/pskHdtS/qy0/jVp2WLD3BmYmRUEYoaJrqYuvCgi9ouI/dP/JwBfi4i9ga+l/wEOA/ZO00rgY1B1HlRDCj8bOAA4aa4DSfP8QctyK0oNcWdgZtbGbKiraQscAZyTHp8DHNlSfm5UrgJ2lLQrcChwaUSsj4i7gEuBFalu+4i4Ko0tf27Lumq5MzAzK+jXaaK06q9I+raklanssRFxe3r8c+Cx6fFuwE9bll2bykrla2vKs0bzArKZ2cCImdKV7XrL5q4DJGdGxJmbzfP8iLhN0mOASyX9oLUyIkJqkmOnGXcGZmYFVW6irjuDO1quA9SvN+K29HedpM9TnfP/haRdI+L2dKpnXZr9NmD3lsWXp7LbgIM2K78ilS+vmT9rJE8TheonctMCpdnITk1ERHYiZuun2chPDZ8r24bZ2WZTk9fb5HVF5KfS+hq0vbFc+3r9NFJ2akyqnwag16eJJG0jabu5x8AhwHXAKmAuIuhY4Avp8SrgmBRVdCBwTzqddAlwiKSl6cLxIcAlqe5eSQemKKJjWtZVy0cGZmYFEY1OE7XzWODzKdpzEfAvEfFlSVcDF0g6Dvgx8Lo0/0XA4cAa4AHgjVXbYr2kDwBXp/lOjoj16fHxwNnAVsDFacpyZ2Bm1sZsj08xRMQtwNNryu8EXlxTHsBbM+s6Czirpnw18NRO2+TOwMysoIomGskz6l1xZ2BmVtSX00RDx52BmVlBw2iikePOwMysjZkxSFQ3mp1B5n1plNWwUDcU73/TsL9cCGTTcMVsFtTeZiatlsuss/G26P41x0T+l6Cyr3my6+epnqzHY1n3IVS0l4r7VYN9bovCVTtZf8pNtNAt/FdoZmZtjeaRgZnZAM36ArKZ2XhzaKmZmVXXDIbiAmJ/uTMwM2vDoaVmZmMuAt90tiUknQW8HFgXEY/KjyHp9cC7qYLJ7gPeEhHf7WTd3YaQjsIRnpqEA/Y6hLBpKGhGcQD7Urhnk9dVWl8hTLSR3HYqpZ5XKbx1CL5oGn6Wer1vlcJEG4WO94R6nptoGPXzU3g25TE3bwVeGBG/DXwA2HzgBzOzeRdURwbdTKOob0cGEXGlpD0K9f+35d+reORADGZmQ8PRRINzHG1ybZuZzYdgiwe5Hwnz3hlIehFVZ/D8wjwrgZUAU9ssHVDLzMwqPjLoM0lPAz4OHJYGdaiVBpI+E2CbZbsPd+IVM1tQAt+B3FeSHg98DvgfEfHD+WqHmVlZZ+Maj7p+hpaeBxwELJO0FjgJWAwQEWcA7wN2Bj6axgHdFBH7t18xPc1aWsx4OJGpG4UfCY2yXTbIQFoKR20aPtrrgeBz6yuFnBbakMtoms9mSjEGsxSC2yjkeBj0ODNpabl+/2j3kcEWioij29S/CXhTv57fzKxXfGRgZjbmIuQjAzMzG490FAv/FZqZWVs+MjAzKwgYi9xEI9cZBD1OWLXw3+PeyETKFJPRlZSidQYUQaNBJrfrtQGOgTzQm28bjEne//ZpLE4TjVxnYGY2SFVo6cL/1ejOwMysDaejMDMbc05UZ2ZmgIe9NDMbe9Wwlz4yMDMbez5NNKwyR2y5pHNRGpu2IB/KVhqnNV/X849TgxDCYthmr0M6h32M5sKRfzHsNPcel5Ielj6DjZL9TeaXaaLU9gbLFb87G4SPFpfr8/d0dc3Ap4nMzMaeE9WZmY0532dgZmbg00RmZgbOTWRmNvYcWmpmZoCHvRxa3WYtjdxYxkBMjOgYs+0MKqyzFILZNDNp7rlKYwyXKLMjl0I6J0vxj/XLFccyzq9t7BTHHS9t9lxIed/HQB6PdBQLv7szM7O2RvLIwMxskHwB2cxszPk+AzMzA3wB2czMYjwuILszMDMrCHzNYDipQYhZwwyK+eyUDdc37EeaTbKCNg1hLT1X0xDSbteXCzmlD2GipddUfL09zk7aRCE0O/uZLi3TMGtpNiR1AN/TPjIwMxtzvoBsZmbAeHQGw37iwsxsXs3dgdzN1ClJk5K+I+mL6f+zJd0q6Zo07ZfKJel0SWskXSvpmS3rOFbSTWk6tqX8WZK+l5Y5XSqPYOTOwMysjVnU1dSFdwDf36zsXRGxX5quSWWHAXunaSXwMQBJOwEnAc8GDgBOkrQ0LfMx4A9alltRaog7AzOzkqAvRwaSlgMvAz7ewexHAOdG5SpgR0m7AocCl0bE+oi4C7gUWJHqto+Iq6KKhjgXOLL0BO4MzMwK5i4gd9kZLJO0umVaWbPqvwP+DNg8nOyUdCroNEnTqWw34Kct86xNZaXytTXlWSN5ATmXhTQ38H2jcLXC+voh146eX7ZqOKh8MctoD5cprq9h23M00TCkM5eNtXRKdnK8fncV97nCpihnGO5+fb3S4ALyHRGxf65S0suBdRHxbUkHtVSdCPwcmALOBN4NnNztkzcxXp9QM7Mu9ekC8vOAV0r6EXA+cLCkT0XE7elU0Abgn6muAwDcBuzesvzyVFYqX15TnuXOwMysjQh1NbVfX5wYEcsjYg/gKOCyiHhDOtdPivw5ErguLbIKOCZFFR0I3BMRtwOXAIdIWpouHB8CXJLq7pV0YFrXMcAXSm0aydNEZmaDNMB0FJ+WtAvVGeJrgDen8ouAw4E1wAPAGwEiYr2kDwBXp/lOjoj16fHxwNnAVsDFacpyZ2BmVhDR35vOIuIK4Ir0+ODMPAG8NVN3FnBWTflq4KmdtqNvp4kknSVpnaTrMvX7SPqGpA2S3tmvdpiZWXv9PDI4G/gwVXxrnfXA22kT+/ooItuFNYk2KEYi5MZULo7h2rBu3BSSs/U6aqiRpuMjj6hyVF2DuobjHDepG0Q0USfXAUZd3zZjRFxJ9YWfq18XEVcDG/vVBjOzLde/dBTDpKMjA0m7AU9onT992ZuZLXjjcGTQtjOQdCrwu8ANwEwqDmBgnUG6e28lwOLtlraZ28ysd5zC+mFHAr+VboKYFxFxJtXdeGz92N2H4ISymY2NaD6G0yjppDO4BVgMzFtnYGY2nzzsZeUB4BpJX6OlQ4iIt5cWknQecBBVwqa1VGlWF6dlz5D0G8BqYHtgVtIfA/tGxL1NXoiZWT8EvmYwZ1WauhIRR7ep/zmPzJ3R+bqz465mynscyjYGPxI6l0vaBo2T4mX1eGzkmM0H06mQpy6XgG8kPha9Dm/Ohl/nF2mSjK5UV1pfb4xuhFA32nYGEXGOpCngSanoxohwOKiZjQ1fMwBSetVzgB9R/Q7YXdKxDi01s3Hh00SV/w0cEhE3Akh6EnAe8Kx+NszMbBhEuDOYs3iuIwCIiB9KWtzHNpmZDRVfM6islvRx4FPp/9dTRQGZmY0FXzOovIUqdepcKOl/AB/tW4vMzIaMTxMB6c7jv03TvAu6z07aPJQtN9ZyYZliqGqhchjGnCuFbpZCSBs9VeGnVo9DSBsptqEQd5ozyEysTcJHG2TvhVLW0tIyDevmaQzkoLPRy0ZdtjOQdEFEvE7S96i+gx8hIp7W15aZmQ2JMThLVDwyeEf6+/JBNMTMzOZP9gArDagMcHxE/Lh1ohpb08xs4Uuhpd1Mo6iTs20vrSk7rNcNMTMbWtHlNIJK1wzeQnUEsJeka1uqtgP+s98NMzMbFqP6a78bpWsG/wJcDHwQOKGl/L6IyA5naWa20Iz1fQYRcQ9wD3A0gKTHAEuAbSVtGxE/GUwTNyOITGRfk9CzRqFspQHRi4OHlxZrMrD4wv+18ghN98hBbacR+MZoNoB994Pb9zqcu7Rc7vugV8YlhXXbawaSXiHpJuBW4N+pEtZd3Od2mZkNh6D6JdfNNII6uYD8l8CBwA8jYk/gxcBVfW2VmdkQiehuGkWddAYbI+JOYELSRERcDuzf53aZmQ2PcY4manG3pG2BK4FPS1oH/Kq/zTIzGxaje+9ANzo5MjgC+C/gT4AvAzcDr+hno8zMhoqPDCAiWo8CzuljW8zMhs+4D24j6esR8XxJ9/HIvk5ARMT2fW9dRrchnzGR76obhbk1zT7a94G7h8wwZB+FwhW9HmcmHaRehxw3DInO1jXM3lveH7tfpmdG9Nd+N0r3GTw//d1ucM0xMxtGC/+HXCf3GZwu6TmDaIyZ2VAag2sGnRxgfRv4c0k3S/obSQ4rNbPx4s4AIuKciDgc+B3gRuDUdEeymdnC5zuQH+WJwD7AE4Af9Kc5ZmY2H9qGlkr6a+BVVPcXnA98ICLu7nfD8g3qcaK6BmO/lpPblaIyGkSAlCKQeh5R0n1YRvTj3vtRvZ9/FOT2kabRP7kIvkIwVtP9sclz9co4fCQ7uQP5ZuA5EXFHvxtjZjaUxqAz6OSn4D8BKyS9D0DS4yUd0N9mmZkNEV8zAOAjwHNI4xoA96UyM7OxoOhuGkWdnCZ6dkQ8U9J3ACLiLklTfW6XmdlwGOFw0W500hlslDRJ2hySdqF8H7+Z2QIyuqd+utHJaaLTgc8Dj5F0CvB14K/62iozs2EyBjeddZK19NOSvk01wpmAIyPi+31vWalNXY6FWg5zKyXNyryrxTGVux8vtrRcdmzkpnqcLK/UvuI+MSxJ7Aal6XafyHzYGo9Z3H0Ic5MxvxuPgVzaVxvs3z0zol/w3ShlLd2p5d91wHmtdRGxvp8NMzMbGuPcGVDlJAqq37OPB+5Kj3cEfgLs2ffWmZnNt7l0FAtc9oAtIvaMiL2ArwKviIhlEbEz8HLgK+1WLOksSeskXZepV8qIukbStZKe2fRFmJn10ziElnZyAfnAiLho7p+IuBh4bgfLnQ2sKNQfBuydppXAxzpYp5nZ4I3BBeROOoOfSXqvpD3S9B7gZ+0WiogrgdJ1hSOAc6NyFbCjpF07a7aZ2eiTNCnpO5K+mP7fU9I30xmTz8zd0yVpOv2/JtXv0bKOE1P5jZIObSlfkcrWSDqhXVs66QyOBnahCi/9XHp8dHGJzuwG/LTl/7WpzMxsqPTxNNE7gNbozFOB0yLiiVTXaY9L5ccBd6Xy09J8SNoXOAp4CtWZmI+mDmaSKlPEYcC+wNFp3qxOQkvXpwbPG0krqU4lsWiHpcxmWj2bCTGbLYTGxWT+ncst1zRsrtHYycXxYvN1vQ9JzTRwZqa3zzPKmm7zBtliixpkum0cEt0gU3BxfxzaMZB7fwFZ0nLgZcApwJ+q2mkPBv57muUc4P1Up9CPSI8BLgQ+nOY/Ajg/IjYAt0paA8zljlsTEbek5zo/zXtDrj2D2Iw5twG7t/y/PJU9SkScGRH7R8T+k1tvM5DGmZkB3V8v6PzI4O+AP+PhjA47A3dHxKb0f+vZkl+fSUn196T5c2dYuj7zMp+dwSrgmBRVdCBwT0TcPo/tMTPrlWWSVrdMK1srJb0cWBcR356n9j1KJ7mJGpF0HnAQ1UZZC5wELAaIiDOAi4DDgTXAA8Ab+9UWM7Mt0n2E0B0RURov/nnAKyUdDiwBtgf+niqQZlH69d96tmTuTMpaSYuAHYA7KZ9h6ejMy5zSHcj/QGETRMTbSyuOiOJF5qiGyHpraR4zs2HQ63sHIuJE4EQASQcB74yI10v6LPAaqlEljwW+kBZZlf7/Rqq/LCJC0irgXyT9LfA4qlD9b1Fd6dlb0p5UncBRPHwtolbpNNFqqruQlwDPBG5K036AU1ib2fgY3H0G76a6mLyG6prAJ1L5J4CdU/mfAicARMT1wAVUF4a/DLw1ImbSkcXbgEuoopUuSPNmZY8MIuIcAElvAZ4/d1FD0hnAfzR8oWZmo6ePN5JFxBXAFenxLTwcDdQ6z4PAazPLn0IVkbR5+UVUp+M70sk1g6VU57PmbiDbNpXNDzXIXth0AO4us6NWdc1CQXOhm6WQv2L4aOa5illGe5zRtKgUThkOV/213PvVMGtp9vNUDB/t/jPTj6yl2ZDyvl35rIxyioludLIZPwR8R9LlVB+ZF/BwvKuZ2cI3Bonqip2BpAngRuDZaQJ4d0T8vN8NMzMbGuN+ZBARs5I+EhHP4OGr2mZmY2UcThN1ctPZ1yS9Wj3Pb2BmNiLGIGtpJ9cM/pAqlGmTpAeprhtERGzf15aZmQ0DX0CuRMR2g2iImdnQcmdQkbSU6s62JXNlabyCwRPEovp3JhfWWQo9K4aJZrMkNgvPLIad5iL+mmZBbRKS2GvFxhfCR3NtjBHeI0vbokmW0QbLANAkE28x62+mvPGg9w2eq7C+nhnhj16n2nYGkt5ElcJ6OXANcCDVLdEH97dpZmbDYRxOE3VyAfkdwO8AP46IFwHPAO7ua6vMzGygOjlN9GBEPCgJSdMR8QNJv9X3lpmZDYsxODLopDNYK2lH4P8Al0q6C/hxf5tlZjYkHE1UiYhXpYfvTykpdqDKjmdmZgtEaTyDnWqKv5f+bsvDiesGKshHBzWJNsiNm1wtl4lOKoybPFmM/ik9V2bB3NjD0DhZWZNlcvccRql9M7P5ukHKRlaVtm2+Lnv/ZWlbNB2vOrdcKRldISInu1zDBIv5fSTfhOb7Y2aZxQP42T7mRwbfptoEAh4P3JUe7wj8BNiz760zMxsGY9AZZH/KRMSeEbEX8FXgFRGxLCJ2Bl4OfGVQDTQzm0/i4TTWnU6jqJPQ0gPTIAkARMTFwHP71yQzsyHj3EQA/EzSe4FPpf9fD/ysf00yMxsiI/xrvxudHBkcDewCfD5Nj0llZmbjwUcGEBHrqe5CNjMbTyP6Bd+NTnITPQl4J7BH6/wRMT+5iQqJ6mYX1y8yu6gQGleoy4WQNk1UN1sI35vIjSVbCvcsxbFmQwgLsXubCsnjGoQ4FsMpC8MIZodALj3XIJPYZUJIGw/50SREuGHIcS4UtPGYxdlxwkv7QWF9hW+kbkPKe2kcThN1cs3gs8AZwMcpppo0M1ug3BkAsCkiPtb3lpiZDaMRvg7QjU46g3+TdDzVxeMNc4XpWoKZ2YLn00SVY9Pfd7WUBbBX75tjZjaE3BlUdyIPoiFmZsPKRwaJpKcC+/LIYS/P7VejzMyGijsDkHQScBBVZ3ARcBjwdWB+OgPlsxTmQkh7nSWxuL5CqOpEk8yQTTJQFtbXKEMmtBnPOLNIqX3FBXPxioUsqMXwzEwoaOn1luqaLDNR+NA0CRMthoLm15f7XDQdAzkXLl0OHy2srzReeTakvM/f1GNyAbmTPfw1wIuBn0fEG4GnU41pYGa24KnBNIo6OU30XxExK2mTpO2BdcDufW6XmdnwGIMjg046g9Vp2Mt/ohrj4H7gG31tlZmZDVQn0UTHp4dnSPoysH1EXNvfZpmZDY9xiCZqe81A0tfmHkfEjyLi2tYyM7MFb5yzlkpaAmwNLJO0lIevi2wP7DaAtpmZDYcR/YLvRuk00R8Cfww8jupawVxncC/w4U5WLmkF8PfAJPDxiPjQZvVPAM6iGi9hPfCGiFhbXmkQmVCyXMbDcrhavi63XCn7aGQynUKz8L2mIX/KhRc2HJg9245SWORMIa9hacD5ifoQ0pjtPry1qMmg91Vl1+srh+2W3uPc+1h4rgYhzKWQ6GLW31z4dYP9qrQ+yGclzn0f9My4D24TEX+f7j5+Z0TslcZE3jMinh4RbTsDSZPAR6juS9gXOFrSvpvN9jfAuRHxNOBk4IONX4mZWb+MwWmiTn5q/VzSdgCS3ivpc5Ke2cFyBwBrIuKWiHgIOB84YrN59gUuS48vr6k3M5t3uYHvc9Mo6qQz+POIuE/S84GXAJ8AOklpvRvw05b/1/Loaw3fBf5bevwqYDtJO3ewbjOzwfGRAfDwgDYvA86MiC8BUz16/ncCL5T0HeCFwG3UDKAjaaWk1ZJWz9z/qx49tZlZZ8bhyKCTm85uk/SPwEuBUyVN01knchuPvFN5eSr7tYj4GenIQNK2wKsj4u7NVxQRZwJnAkw/YfmIbmozG0kj/Gu/G518qb8OuAQ4NH1R78QjxzbIuRrYW9KekqaAo4BVrTNIWib9OgTjRKrIIjOz4TIGp4k6uQP5AeBzLf/fDtzewXKbJL2NqiOZBM6KiOslnQysjohVVNlQPygpgCuBt7ZtsYDF9aGHs4vr+7ZcSFpVVwiby4TUzWayJ0K7sNMGoaWLSuGKpfDCJiGJhbi+TZtqi8vhqKWR1At7TCaENBdy2lgxy2iD7d4wbDf7XhWeq/hZKoUjZz5P5RDmbFU27LSUtbQUWlrcVzPL9Tu0VIzuqZ9udDSeQVMRcRFV2uvWsve1PL4QuLCfbTAz22LuDMzMTKUj2QXCnYGZWckIXwfohjsDM7M2fM3AzMx8ZDCUJkBTmURmU5kEdoWIoSZRD6UkXKXIi/JyDRLVFSKNtDEXhdIwOVsu0miyEOFTOM+q2fxykWtijxPVNUpGByi3DZtsP+h5VFjpc5H9nDUY57iqy7Wh9LnNVhXHM85GDWWiC607PU4DaWa28PT6DmRJSyR9S9J3JV0v6S9S+dmSbpV0TZr2S+WSdLqkNZKubc0PJ+lYSTel6diW8mdJ+l5a5nQVf/2M4pGBmdmg9f400Qbg4Ii4X9Ji4OuSLk5170ph960OA/ZO07Op8sM9W9JOwEnA/qmV35a0KiLuSvP8AfBNqhD/FcDFZPjIwMyspMujgk6ODKJyf/p3cZpKSx5Ble4/IuIqYEdJuwKHApdGxPrUAVwKrEh120fEVRERwLnAkaU2uTMwM2unD+koJE1KugZYR/WF/s1UdUo6FXRaygUH+SzQpfK1NeVZ7gzMzArm0lF0eWSwbC7TcppWbr7eiJiJiP2okngeIOmpVDna9gF+hyoP3LsH9Tp9zcDMrJ3u70C+IyL272zVcbeky4EVEfE3qXiDpH+mSvMP+SzQt1HleGstvyKVL6+ZP2vkOgMpWLS4fmzdjYvr49yKY7iWwk4zSbOKibYajiWbq5sohQk2GB85OzYyFBOmaVP3oa/FJHaFUEtlxk7Ohpy2kW1Hk2R0kB/ruBQ+2mSc40I7iuGjhfckm1iutB80GEO8ScK5tnWZkPKJqcJY2z3S65vOJO0CbEwdwVY8PETArhFxe4r8ORK4Li2yCnibpPOpLiDfk+a7BPgrSUvTfIcAJ0bEekn3SjqQ6gLyMcA/lNo0cp2BmdlA9Scdxa7AOWms+Anggoj4oqTLUkch4BrgzWn+i4DDgTXAA8AbAdKX/geohgwAODki1qfHxwNnA1tRRRFlI4nAnYGZWVvq8X1tEXEt8Iya8oMz8weZFP8RcRY1Y8FExGrgqZ22yZ2BmVk7TkdhZmZOVGdmNu6CJtFEI8edgZlZGz4yGEJSMDVdPx7vpg31L2d2Kh+GN5PJdAowsaFJGF5hfORSiGuDbJLl7JT1YY6aLIThNRkfeaZwZW2ysAc1+KVVHG2qkAU1+7pK4a2lcM9MXTFsd1GzrKW9zmab/5xlFyl+bmcy+0J53PFCXWF/JJO1dNEAQkt9zcDMbMzN3YG80LkzMDMriRiLawbOTWRmZj4yMDNrx6eJzMzMF5DNzMxHBkNpQsHW0w/V1m2Yqo9Z2zSVj5trEgLXOGyuOBB4JkSvENY3UQx/zIUQFmIIC+GP2lQfvlcKwYxCmlEV6ii5fEUAABFjSURBVLL7XSl8tMGA88UhYUvry9U1HfS+sN1z71cpfLSYHbfB56z8uW2wTCF8NAp1mq7/DC6Z3ph/sl4IYHbh9wYj1xmYmQ3cwu8L3BmYmbXj00RmZjYW9xm4MzAza8NHBmZm464/I50NHXcGZmYFVW6ihd8bjFxnMDERbDe9obbugSVTteWbNhRCSzcUQvTqV8fsQ92H7kGzgcCjlOm0GEJY/5onNzbMWpoLfyztJA3rsq+qFApaksvwmRvYHmCiEGabzYJayhZaGsC+tN3r62JxYX0Nwk5z2UehYfh1Zt+plimFlubDhxdP1Wcr3nZJ/fdBT/V42MthNHKdgZnZoPnIwMxs3PmagZmZwXiksHZnYGbWxjiElvZ1PANJKyTdKGmNpBNq6h8v6XJJ35F0raTD+9keM7NG5ga46XQaQX07MpA0CXwEeCmwFrha0qqIuKFltvcCF0TExyTtC1wE7FFa7yLNsnT6gdq6+6ana8sfnM6HQ8xMFyKNHqyva5yorkESu2J0UiGiRBvrwx8iE2UEoE35kInIJIlTIXmcZkoRL4WPXi5qqDgGcqEuF01UihgqJQHMtb1hEsDSe5J7j2cL7SsnnWsyZnH3kUalscVnp7tPRgewZEl9QrrtpgYQTTQG+nlkcACwJiJuiYiHgPOBIzabJ4Dt0+MdgJ/1sT1mZt0L0Gx30yjq5zWD3YCftvy/Fnj2ZvO8H/iKpD8CtgFeUrciSSuBlQBbPXbbnjfUzKxoRE/9dGO+x0A+Gjg7IpYDhwOflB59J1BEnBkR+0fE/tM7bjXwRprZmIsupxHUzyOD24DdW/5fnspaHQesAIiIb0haAiwD1vWxXWZmXRmHm876eWRwNbC3pD0lTQFHAas2m+cnwIsBJD0ZWAL8so9tMjPrnqOJmouITZLeBlwCTAJnRcT1kk4GVkfEKuB/Av8k6U+oDq5+L2JEt6SZLUyBcxNtqYi4iCpctLXsfS2PbwCe1806F03MsPP0r2rr7preurY8l8AO4FcP5jfB7HT9gdNM/RDMAExMFcYsLtVl1jlRCOubKI2PnAlJjE2FcM9S2GkuhHSmsJcsKoQQbqpPOgalMZBL4aP5qlxoaXEM5FLoayZMNBqGj+aS0UE+wd3sVCnBYvehoMXQ0lLSuVwyx+I4x4VkdNP5z8V2mYR0O2VCzXtFxFicJvIdyGZm7bgzMDMzdwZmZuPO1wzMzAzGI7TUnYGZWTvuDMzMxt3o3jvQjZHrDBZrhsdN31Nbd9eS+tDS+zbUZzMF2LAkH1O36aFMWF9pTOVeh+htLIQJlkIIMyGk2ljIdFoIcWRT5jUvKmQtbbgDZcNYCxlSi3JjDBdDSwthorm6xfndKQoZTUvZZ3MhpKXMpKXxjHN1M4XP0kzhczuTyUBazEy6JJ+ZdKtMZlKAHaYfrC1/zJL7ssv0RODOwMzM8AVkMzMbjwvI85211MzMhoCPDMzM2hmDIwN3BmZmJUE5L9YC4c7AzKzIoaVDaUozLJ9aX1v3y+ntasvvWZIfHe2BDfm4uU2ZENKZJflLLbnsowAThTDRXEbTiXykXXl9mRBCTRUyk87kP/ATmbrShbXS7lPMGDqTCT2cbXiJK5O1NBtySjkDaS6ENErZRwvho6W6XNhpk8yk1XLdlbevq3+XYzoffjO1pPvMpADLltxfW77b9F3ZZXrGnYGZmbkzMDMbd2NyzcChpWZmRQEx293UhqQlkr4l6buSrpf0F6l8T0nflLRG0mfSkMFImk7/r0n1e7Ss68RUfqOkQ1vKV6SyNZJOaNcmdwZmZu30fgzkDcDBEfF0YD9ghaQDgVOB0yLiicBdwHFp/uOAu1L5aWk+JO1LNb78U4AVwEclTUqaBD4CHAbsCxyd5s1yZ2BmVjJ3mqibqd0qK3NXxBenKYCDgQtT+TnAkenxEel/Uv2LVUVhHAGcHxEbIuJWYA1wQJrWRMQtEfEQcH6aN8udgZlZO70/MiD9gr8GWAdcCtwM3B0Rc+FWa4Hd0uPdgJ9WTYlNwD3Azq3lmy2TK88auQvIU9rIHlO/rK37xfQOteXrM9lMAe59KJ/R9KEN9ZtnYylr6XShrhB2mgvfm2kaWprLWppPGMnETCH8MTKDwBc++IXg0WLYaXYA+4YX8aJBaGkxa2lmcPuYyu9Os4WQ3pliaGmTLKPd1zXJTAowu6S+bmKrfPjo1kvyO8JOW+UHt991yb215U+YuiO7TM90H020TNLqlv/PjIgzH7nKmAH2k7Qj8Hlgny1r5JYZuc7AzGywGt10dkdE7N/R2iPulnQ58BxgR0mL0q//5cBtabbbgN2BtZIWATsAd7aUz2ldJldey6eJzMxKgmosjW6mNiTtko4IkLQV8FLg+8DlwGvSbMcCX0iPV6X/SfWXRXVYvgo4KkUb7QnsDXwLuBrYO0UnTVFdZF5VapOPDMzM2un9TWe7AuekqJ8J4IKI+KKkG4DzJf0l8B3gE2n+TwCflLQGWE/15U5EXC/pAuAGYBPw1nT6CUlvAy4BJoGzIuL6UoPcGZiZtdPjziAirgWeUVN+C1Uk0OblDwKvzazrFOCUmvKLgIs6bZM7AzOzos7CRUfdyHUG05phr0X1YyD/fLo+quDOjdtk13fvxiXZulwSu00PFSJDHipE+JTGH96US1RXiAwpRRplgjkmMs8DMJOJQKoqMztDIUqmeEGqkKhOmUR10XSHbJSorpRYrv41l8cyLkSZZcY5BpiZbhBNlA+Qy9bNlpbJRAwBRGY84+lCMrodt/6vbN1jC+MZP376ztry31xcH13YMwHRwV3Fo84XkM3MbPSODMzMBs6niczMzCmszczGXURH9w6MOncGZmbt+MjAzMzCRwbDZ1qT7Ll429q6X8ysqy1fv6R+foC7NxaS2G1VH3a64aH8ZtvQNOw0U7cpH/mK8tF72RDSiZlCSGdhjGH1+IdRcQzkyUxdYYzmotzLKoWWZsJHoTAucSlJYTF8tFDXKLFc92Gnmwrho7NL8l+Ek1vVh5Zuu1V+LONdtqofyxhg+ZL8eMa/OVW/f++1qLAj9ESj3EQjZ+Q6AzOzgRqTYS/dGZiZteObzrZMuzE4JZ0m6Zo0/VDS3f1sj5lZt4Lq7vduplHUtyODljE4X0o1ys7VklZFxA1z80TEn7TM/0fUJG4yM5tXET4y2ELdjsF5NHBeH9tjZtaIjwy2TN0YnM+um1HSE4A9gcv62B4zs2bG4MhgWC4gHwVcODcow+YkrQRWpn83TO665rruVn/rFjUuYxkwgMFXh74NMBztGIY2wHC0Y97bcGuhHf+vsNy/9r4pv7WlK7iPuy75aly4rMvF5vsz0LV+dgalsTk3dxTw1tyK0kDSZwJIWt3p2KL9NAztGIY2DEs7hqENw9KOYWjDsLRjs0HpG4mIFb1oy7Dr5zWDjsbglLQPsBT4Rh/bYmZmBX3rDCJiEzA3Buf3qcb4vF7SyZJe2TLrUcD5aXBnMzObB329ZlA3BmdEvG+z/9/f5WrP3MJm9cowtGMY2gDD0Y5haAMMRzuGoQ0wHO0YhjaMBPkHuZmZedhLMzMbzs5A0mslXS9pVlI2GiGX7iJdtP5mKv9MuoDdbRt2knSppJvS36U187yoJZ3GNZIelHRkqjtb0q0tdft124ZO25Hmm2l5rlUt5YPaFvtJ+kZ6366V9LstdVu0LTpIazKdXtua9Fr3aKk7MZXfKOnQ7l55V234U0k3pNf+tXTvzFxd7XvTp3b8nqRftjzfm1rqjk3v4U2Sju1jG7JpZnq1LSSdJWmdpNowc1VOT228VtIzW+p6sh0WnIgYugl4MlV88BXA/pl5JoGbgb2AKeC7wL6p7gLgqPT4DOAtDdrw18AJ6fEJwKlt5t8JWA9snf4/G3hND7ZFR+0A7s+UD2RbAE8C9k6PHwfcDuy4pdui9D63zHM8cEZ6fBTwmfR43zT/NNVNjTcDk31qw4ta3vu3zLWh9N70qR2/B3w48/m8Jf1dmh4v7UcbNpv/j4Cz+rAtXgA8E7guU384cDEg4EDgm73cDgtxGsojg4j4fkTc2Ga22nQXkgQcDFyY5jsHOLJBM45Iy3a6jtcAF0fEAw2eq5ft+LVBbouI+GFE3JQe/wxYB+zS4Lk210lak9b2XQi8OL32I6gi1TZExK3AmrS+nrchIi5vee+vorqvpte6TfHS6lDg0ohYHxF3AZcCTeLnhyLNTERcSfXjK+cI4NyoXAXsKGlXercdFpyh7Aw6VJfuYjdgZ+DuqEJbW8u79diIuD09/jnw2DbzH8WjP/SnpEPU0yRlhhXpWTuWSFot6aq5U1XM07aQdADVr8abW4qbbovc+1w7T3qt91C99k6W7VUbWh1H9at0Tt1700Sn7Xh12tYXSpq78XPg20L1aWZ6tS3aybWzV9thwZm3dBSSvgr8Rk3VeyLiC/PdhtZ/IiKk/Fhf6RfHb1PdUzHnRKovzimq8LZ3Ayf3sR1PiIjbJO0FXCbpe1Rfih3p8bb4JHBsxK8TunS8LUadpDcA+wMvbCl+1HsTETfXr2GL/RtwXkRskPSHVEdMB/fpudqpSzMzyG1hXZi3ziAiXrKFq8ilu7iT6pBwUfqVmE2DUWqDpF9I2jUibk9fcPVj7lVeB3w+Ija2rHvul/QGSf8MvDO3cC/aERG3pb+3SLqCKh34vzLAbSFpe+BLVB36VS3r7nhb1OgkrcncPGslLQJ2oPocdJMSZUvbgKSXUHWeL4yIX4/7mHlvmnwBtm1HRNzZ8u/Hqa73zC170GbLXtGPNrR4VJqZHm6LdnLt7NV2WHBG+TRRbbqLiAjgcqpz+ADHAk2ONFalZTtZx6POi6Yvzbnz9kcCXSbX67wdkpbOnXqRtAx4HnDDILdFeg8+T3We9sLN6rZkW3SS1qS1fa8BLkuvfRVwlKpooz2BvYFvdfHcHbdB0jOAfwReGRHrWspr35sGbei0Hbu2/PtKqrv/oTpqPSS1ZylwCI88ku1ZG1I7HpVmpsfbop1VwDEpquhA4J70o6RX22Hhme8r2HUT8Cqqc3kbgF8Al6TyxwEXtcx3OPBDql8W72kp34tqp18DfBaYbtCGnYGvATcBXwV2SuX7Ax9vmW8Pql8bE5stfxnwPaovvk8B2zbcFm3bATw3Pdd309/jBr0tgDcAG4FrWqb9erEt6t5nqtNMr0yPl6TXtia91r1aln1PWu5G4LAt+Ey2a8NX02d17rWvavfe9KkdHwSuT893ObBPy7K/n7bRGuCN/WpD+v/9wIc2W65n24Lqx9ft6TO3luo6zZuBN6d6UQ2udXN6rv1blu3Jdlhok+9ANjOzkT5NZGZmPeLOwMzM3BmYmZk7AzMzw52BmZnhzsAWGEk7Sjp+vtthNmrcGdhCsyNVFtOOpRuTvC/YWPMOYENF0jEpydp3JX1S0i6S/lXS1Wl6Xprv/apy2l8h6RZJb0+r+BDwm6ry5f+vNO+70rLXSvqLVLaHqpz851LdDLd7XXvMxsW85SYy25ykpwDvBZ4bEXdI2gn4MHBaRHxd0uOpUgc8OS2yD9U4AtsBN0r6GNV4C0+NiP3SOg+hSkNxANVdqaskvQD4SSo/NlryKJmNK3cGNkwOBj4bEXcARMT6lPxt3yqtEQDbS9o2Pf5SVAnhNkhaR31q7UPS9J30/7ZUncBPgB+7IzCruDOwYTcBHBgRD7YWps5hQ0vRDPWfZwEfjIh/3Gz5PYBf9bKhZqPM1wxsmFwGvFbSzlCNvQx8hWroRFJZu/GT76M6bTTnEuD3544mJO0m6TE9bbXZAuAjAxsaEXG9pFOAf5c0Q3Vq5+3ARyRdS/V5vZIqO2VuHXdK+k9VA6VfHBHvkvRk4BvpaOJ+qgyrM7l1mI0jZy01MzOfJjIzM3cGZmaGOwMzM8OdgZmZ4c7AzMxwZ2BmZrgzMDMz3BmYmRnw/wGtp6p5HYAszAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"grid = np.array([a.flatten() for a in np.meshgrid(c_scan, s_scan)])\n",
"\n",
"with minkit.minimizer('uml', g, data) as minimizer:\n",
" fcn = minimizer.fcn_profile(('c', 's'), grid)[::-1] # need to reverse for \"imshow\"\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(6, 5))\n",
"\n",
"fcn = fcn.reshape((len(c_scan), len(s_scan)))\n",
"\n",
"m = ax.imshow(fcn, extent=(*c.bounds, *s.bounds), aspect = 'auto')\n",
"ax.set_xlabel('center');\n",
"ax.set_ylabel('standard deviation');\n",
"fig.colorbar(m, ax=ax);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Evaluating residuals and pulls\n",
"Another common operation that is performed once a fit is done is to look at the residuals or the pull plots of the binned distribution. To do this correctly, it is necessary that we evaluate the resulting PDF in the bins of the distribution, and thus a binned data set must be constructed. In the following lines we will plot the result of the fit together with the distribution of the residuals below."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAFlCAYAAACDVh3MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5yUZf3/8deHBUSOohw8ALuUeMoQbQXNTmoq+tVIUxMRFUjMwMzUFPuVlmlaec5S5CAkoOI5JRPxXKKsiqiggrqrIAqJAiIiLJ/fH9csLriwOzsze83c834+HvPYmXvumXkP7M5nruu+7usyd0dERKTQNYsdQEREJBtU0EREJBFU0EREJBFU0EREJBFU0EREJBFU0EREJBGaxw5QW6dOnbysrCx2DBERySPPP//8/9y9c3375VVBKysro6KiInYMERHJI2ZW1ZD91OUoIiKJoIImIiKJoIImIiKJoIImIiKJoIImIiKJoIImIiKJoIImIiKJkFfnoYkUq88/hxdegDlzYOlSWLsWOneG3XeHvn2hbdvYCUXynwqaSEQvvQTXXAP33APLl9e9T4sWcOSRcPrpcOihYNa0GUUKhbocRTajrKwMM9twyea0bJWVcMwx0KcP3HknHH003HVX2P7557BuHSxeDNOmwYgR8J//QP/+cMABoMl0ROpm7h47wwbl5eWuqa8kX5gZtf8+Nr3dGO4wbhyceWZoaV1wQbi+zTZbftznn8P48fC738EHH8B558Ef/gDN1cciRcDMnnf38vr2UwtNJEvqa9GtXg2nngo/+Ql885swbx785jebL2a1n2+XXco4/fTwmKFD4YorQvfj0qU5f1siBUPf70SypKqq6kstuhrLl8OAAfDEE3DRRaGQlZQ0/PlqnqtDB7j5Zvj2t8MxtX79YMYM6Nkz++9HpNCooInk2LJlcMghYQTjpElw4onpP0dpaelGBbK0tJQnnqikf3/4zndCUdtllyyGFilA6nIUaaRNuxhLS0u/tM+qVfB//wevvAL33de4YgZQWVmJu2+4VFVV0bcvPP44rFkDBx4I77yT2fsRKXQqaCKNVNMlWHOprKzcZI/mHHssPPccTJkCRxyR/Qy9e4fW2SefhFGQy5Zl/zVECoUKmkjO/JmHHoIbbwxD9HPl61+He++FN98Mr7N2be5eSySfqaCJ5MCttwL8grPOgtNOa9hjGtKFuTkHHghjxoRBJ6NGNSqySMHToBCRLJs9u6aIPc6f//y9Bj9u01GS6Ro8GJ59Fq68Mox+PO64Rj+VSEFqcAvNzFqZ2XNm9pKZvWpmv0tt72lmz5rZAjO73cxaprZvlbq9IHV/WW7egkj+WL0aBg6EbbeFbt3OoWXL3Mw0sjlXXQX77RcKqgaJSLFJp8txDXCQu+8F9AH6m9l+wBXA1e6+M/ARMCy1/zDgo9T2q1P7iSTa+efDa6/BLbfAu+8+/6WRibnWsmU4NaC6OpzEvX592J7LabxE8kWDC5oHn6RutkhdHDgIuDO1fQLww9T1AanbpO4/2GqfSCOSMA8/DNdfDz//eTjvLJavfAWuvRYeeyz8hC+PyGyK4irS1NIaFGJmJWY2G1gCTAfeBD5293WpXRYCO6Wu7wS8C5C6fzmwXR3POdzMKsysYqnm8ZECtXIlDBsWlnu5/PLYaWDIkDAzyahR8MYbsdOINI20Cpq7V7t7H6Ab0BfYLdMA7j7a3cvdvbxz586ZPp1IFL/9LSxaBGPHwtZbx04TJj6+8UZo1QrOOCN2GpGm0ahh++7+MfAYsD+wjZnVjJbsBixKXV8EdAdI3d8B+DCjtCJ56IUX4Lrr4Kc/hf33b5rXrJkKa0vHxLbfPrQWH30UYFDTBBOJKJ1Rjp3NbJvU9a2BQ4B5hMJ2bGq3U4D7UtfvT90mdf+jnk9r1YhkQXV1mCS4Sxe47LKme926psKqrWYQyBlnNAOeoVmzazSLiCReOi20HYDHzGwOMAuY7u4PAOcDvzSzBYRjZGNT+48Ftktt/yVwQfZii+SHcePCgptXX13/mmZN6YtBIOt56aX9MevExRfHTiWSW1rgU2Qz6lvg06wdXbqspFcveOqpcNxqS8rKyjZqSZWWlm40/2MmC4jWl/VnP4PRo8MkybvtVn8WkXyiBT5Fcu4CliwJJzM35ISU+roJc+nii6FNm7DSdewsIrmigibSCOHz/xwGDYK+fWOn+fIgkU3ngezSBX79a3jgAXjkkUghRXJMXY4im7GlbrxBg2Dy5NVUVW1Njx7Zf/5c+OyzcJ5c+/ZhZGbtFbNz/doimVCXo0iWfdEK2ovJk6F9+3GNLmYxtGoVhvHPmQO33RY7jUj2qaCJNFDNcacBA16iQweorBwRO1LajjsuLAp68cWwbl29u4sUFBU0kTRUVMB998EvfwkdO8ZOk75mzeCSS2DBApg4MXYakexSQRNJw29/G5aG+cUvYidpvKOOgn33hd//Hj7/PHYakexRQRNpoGeegX/9Kwx9b98+dprGMwuttKqqMPekSFKooIk00G9+A507w8iRsZNk7tBD4VvfgksvhTVrYqcRyQ4VNJEGeOYZmDEjLODZtm12nrO+c8dyySwU6EWL4NZbm+xlRXJK56GJbEbtc7MGDICnnw7ddNkqaLG5h2NpK1bA/PkluFfHjiRSJ52HJpIlr74K998PZ56ZnGIGoZU2ahTMnw+dOv203uVoRPKdCppIPa64Alq3DgUtaY4+OkxWvNNON7B+veZ2lMKmgiayBVVVMHkyDB8O220XO032NWsWjgu+9BI89FDsNCKZUUET2YK//CV86J9zTuwkuXPiidC9e9MuUCqSCypoIpvVibFj4aSToFu32Flyp2XLULCffhqefTZ2GpHGU0ET2azTWb0azj03do7cGzo0nCx+9dWxk4g0ngqaSB3ClFAjOOww2GOP2Glyr107OO00uPNOeOed2GlEGkcFTaQOU6cC7FDQczam68wzw7lpf/1r7CQijaOCJrIJd7jmGoB5HHpo7DRNp7QUfvQjGD0aoE3sOCJpU0ET2cR//xuWiYFraVZkfyFnnw3LlwOcGjmJSPqK7M9VpH7XXFOz1tk/YkdpcvvvD/36AZzF+vWx04ikRwVNJKWsrAyzUu68s5qPPrqC0tLOsSNFcfbZAL144IHYSUTSo4ImklJVVcW551ZRUlJCVdX5VFZWxo4UxY9+BPCuBodIwWlwQTOz7mb2mJnNNbNXzeys1PaLzWyRmc1OXY6o9ZhRZrbAzF43s8Ny8QZEsqc1Y8aED/QePWJniad5c4CbmD4d3ngjdhqRhkunhbYOOMfd9wD2A0aYWc0ZOle7e5/UZRpA6r4TgK8B/YG/mVlJFrOLZNlAPv44mZMQp28MLVrAjTfGziHScA0uaO6+2N1fSF1fCcwDdtrCQwYAt7n7Gnd/G1gA9M0krEiuhGXPzmDPPeGAA2Knia+0tBVr107h6qs/wqy1lpORgtCoY2hmVgbsDdTM/DbSzOaY2Tgz65jathPwbq2HLaSOAmhmw82swswqli5d2pg4IhmbNQvgG5xxRlgnrNhVVlby5JMDgY6MGfOplpORgpB2QTOztsBdwC/cfQXwd+CrQB9gMXBlOs/n7qPdvdzdyzt3Ls5RZRJHGNUYFrTs1288Zqs46aTYqfLHt74Fe+4JN9wQO4lIw6RV0MysBaGYTXL3uwHc/QN3r3b39cDNfNGtuAjoXuvh3VLbRPJCVVUV7s6HHzqtWg1h+PA2tG8fO1X+MIMRI+DFFwH6xY4jUq90RjkaMBaY5+5X1dq+Q63djgZeSV2/HzjBzLYys55AL+C5zCOLZNeECfDZZ3DGGbGT5J9Bg8LExfCz2FFE6pVOC+0AYDBw0CZD9P9kZi+b2RzgQOBsAHd/FbgDmAs8BIxw9+rsxhfJjHsYybf//rDXXrHT5J927eDkkwF+zP/+FzuNyJY1b+iO7v40UNfh8mlbeMylwKWNyCXSJB59NJxrNXFi7CT564wz4IYbtmL8eDjvvNhpRDZPM4VIUfv732G77eC442InyV9f+xpstdVz/OpXr28YRKNh/JKPVNCkiO3AvffCkCHQqlXsLPlt9Oi+wK48/rjj7hrGL3lJBU2K2KlUV8Pw4bFz5L9jj4UOHeDmm2MnEdk8FTQpKrXPPWvefDjf+x706hU7Vf5r3TqMeLzzTli2LHYakbqpoElRqTn37IknnHXryhg6NHaiwnHaabBmDdx6a+wkInVTQZOiNHZsGJIelkqRhujTB8rL1e0o+UsFTYrOihUwdSoMHBi60qThTjsNXnkFNHOI5CMVNCk6t98Oq1ej7sZGGDgQ2rQB+MmGbbWPS2pIv8SkgiZFZ9w42GMP6KvFjNLWrh2ccALACaxcGbbVHJesuWhIv8SigiZFZndmzgytMy0T0zinnQbQlvbth2NmlJaWxo4kAqigSdEZQvPmaJmYDPTtG5aVKS8fjbtTWVkZO5IIoIImRWTtWoCTOfJI6No1dprCZRZaaRUVMHt27DQiX1BBk6IxbRpAVw0GyYKTToKttoIxY2InEfmCCpoUjXHjABZz+OGxkxS+bbeFo4+GyZPDWnIi+UAFTYrC++/Dgw8CTKB5gxdNki0ZOhQ++gjuvz92EpFABU2Kwj/+AdXVAONjR0mMgw6C7t1rWr4i8amgSeK5hw/dAw4AeCN2nMQoKYFTToGHH4aFC2OnEVFBkyIwcya89ppmBsmFU08NXxi04rfkAxU0Sbxx48J0TVqVOvu++lX47ndh/PhQ2ERiUkGTRFu1Cm67DY4/PkzbJNk3ZAgsWABPPx07iRQ7FTRJtDvvhE8+UXdjLh17LLRtG1ppIjGpoEmijRsXVqQOA0IkF9q0gR//GO64I3x5EIlFBU0Sa/58ePJJTUTcFIYMCd27U6fGTiLFTAVNEmv8eGjWDE4+OXaS5PvmN2GXXdTtKHE1uKCZWXcze8zM5prZq2Z2Vmr7tmY23czmp352TG03M7vOzBaY2Rwz2ydXb0JkU+vWwYQJcPjhsOOOsdMkn1lopT31FMDOseNIkUqnhbYOOMfd9wD2A0aY2R7ABcAMd+8FzEjdBjgc6JW6DAf+nrXUIvV4+GF47z0NBmlKJ58cWsRwauQkUqwaXNDcfbG7v5C6vhKYB+wEDAAmpHabAPwwdX0AMNGDmcA2ZrZD1pKLbMG4cdCpExx5ZOwkxWPHHeGwwwBOSU0zJtK0GnUMzczKgL2BZ4Gu7r44ddf7QM1KUzsB79Z62MLUtk2fa7iZVZhZxdKlSxsTR2QjS5eGCXMHD4aWLWOnKS6hRdyNRx6JnUSKUdoFzczaAncBv3D3FbXvc3cH0povwN1Hu3u5u5d37tw53TgiXzJpUljMc+hQKCsrw8w2XEpLS2PHS7SjjgL4UBMWSxRpFTQza0EoZpPc/e7U5g9quhJTP5ekti8Cutd6eLfUNpGcqZmIeN99Yc89oaqqCnffcKmsrIwdMdG22gpgEvfeC8uWxU4jxSadUY4GjAXmuftVte66Hzgldf0U4L5a209OjXbcD1heq2tSJCeefx5efhmGDYudpJiN4/PPYcqU2Dmk2KTTQjsAGAwcZGazU5cjgMuBQ8xsPvD91G2AacBbwALgZuBn2YstUrdx46BVKzjhhNhJitlL9OmjddKk6TV47V53fxrY3HwLB9exvwMjGplLJG2rV8PkyWFuwQ4dYqcpbkOGwFlnwZw50Lt37DRSLDRTiCTGPffA8uVw660HahBIZIMGhRGmaqVJU1JBk8QIH55vUV39mAaBRLbddvCDH8Ctt8Lnn8dOI8VCBU0SobISZswAGJ+arUJiGzYMPvwQ/vnP2EmkWOhPXxLhlltqZtSfUM+e0lQOOQR22kndjtJ0VNCk4K1fH2Z5P+QQ2HhyGomppAROPRWmTavGbKcNxzXLyspiR5OEUkGTgvfoo/DOO5qIOB+deipACZdeumjDcc2qqqrIqSSpVNCk4I0dCx07woABsZPIpnbeGeBxxo0Ls7iI5JIKmhS0jz4Kw/UHDQonVEs+Gsebb9aslSaSOypoUtAmTYI1a9TdmN/uol07DQ6R3FNBk4I2dizsvXe4SL76lIEDYepUWLGi/r1FGksFTQrWiy/C7Nnw4osjNDNInhs6FD79FG6/PXYSSTIVNClYY8cCfMayZTdoZpA817cv7LGHuh0lt1TQpCCtXh2On8FddOwYO43Uxyy00mbOBNg9dhxJKBU0KUj33AMffwxhiT4pBCedBM2bAwyJHUUSSgVN8lpZWdmG42O1Z5kYOxZ69gR4PGI6qUtpaelG/2c1xzW7doUjjwQ4mbVro0aUhFJBk7xWVVW14fhYzSwTb70VZgcZMgRAZ+vmm8rKyo3+z2of1wynV3TlwQdjpZMkU0GTgjN+fDgmE6ZVkkJy+OEAizU4RHJCBU0KTDNuuQUOOwy6d4+dRdIVjqFNYNo0WLw4dhpJGhU0KTCHsnBhWGtLCtV4qqvhH/+InUOSRgVNCswwOnUKqyFLoXqDAw5AExZL1qmgScFYuhTgBwweDC1bxk4jmRg6FF5/HZ55JnYSSRIVNCkYoYuqpbobE+C446BNm5rZXkSyQwVNCoJ7zYffTL72tdhpJBOlpaW0b2+sWjWWceM+oUcP/YdKdqigSUF49lmYOxc0M0jhqzlP7ZlnhgFteffdb8WOJAnR4IJmZuPMbImZvVJr28VmtsjMZqcuR9S6b5SZLTCz183ssGwHl+IyZgy0bg2g6dqTol8/+PrXAYbHjiIJkU4L7Ragfx3br3b3PqnLNAAz2wM4Afha6jF/M7OSTMNKcVq+HKZMgRNPBFgZO45kiRmcfjrAN6ioiJ1GkqDBBc3dnwSWNXD3AcBt7r7G3d8GFgB9G5FPhEmTwlpap5+++XkCpTCddBLAp9x0U+wkkgTZOIY20szmpLokaxby2Al4t9Y+C1PbvsTMhptZhZlVLA3jskU2ctNNsM8+UF6+5XkCpfB06ABt2jzAmDGfYNZ+owmoRdKVaUH7O/BVoA+wGLgy3Sdw99HuXu7u5Z07d84wjiRPP+bMqemakiSaMeN4oC1///uKDRNQizRGRgXN3T9w92p3Xw/czBfdiouA2jPtdUttE0nTT2nbFgYOjJ1DcqVvX+jdO7TENXOIZCKjgmZmO9S6eTRQMwLyfuAEM9vKzHoCvYDnMnktKT4ffQTwYwYNgnbtYqeRXKkZHDJ7NhocIhlJZ9j+FOAZYFczW2hmw4A/mdnLZjYHOBA4G8DdXwXuAOYCDwEj3L066+kl0cLMIFuru7EIDBoUTssYPTp2Eilk5nnUxi8vL/cKfUUTQtfTnnvC3LnP4a4BssVg2DC4/XZYtao97is2bC8rK9vouFppaakGAxUZM3ve3cvr208zhUhe+s9/amYG0XjuYjF8OKxaBXDiRtvrWrVcpC4qaJKXbroJ2rcHuC12FGkiffvCXnsBnK7BIdIoKmiSd5Ytg6lTvzjpVorDFzOH7M3MmbHTSCFSQZO8M3YsrFkDP/1p7CTS1AYPBljOX/8aO4kUIhU0ySvV1fC3v8F3v1szca0Uk7ZtAW5h6lR4//3YaaTQqKBJXpk2DSorYeTI2EkknhtYuxZuvjl2Dik0KmiSV66/Hrp1gx/+MHYSiWc+hx0GN94Ia9fGziKFRAVNoiorK6s1e/5uTJ8ejp01bx47mcQ0ciS89x7ce2/sJFJIVNAkqtrnGJ155mvAGk47LXYqie3ww6FnTzQ4RNKigiZ5YeVKuOUWgDvo0iVyGImupAR+9jN48kkAjQ6ShlFBk7wwcWIoanB97CiSJ4YOhVatAEbEjiIFQgVNonMPXUv77gswK3YcyRPbbhsmLYaTUisviGyZCppE9+ij8NprGqovXxZ+J9owfnzsJFIIVNAkuquvhs6d4fjjYyeRfNOnD8BTXH89rFsXO43kOxU0iWw3HnwQRowIx0tKS0trDeM3SktLYweU6K6islJD+KV+OttHIjubVq3CiDZA61zJhi81NXr06EmLFnDllXDssRGDSd5TQZNoli4FOJmTTw5djiJQ95eaG24Ix9P++9+mzyOFQ12OEs3f/gbQirPPjp1E8t2pp0LHjqGVJrI5KmgSxerV4Vs3PMBuu8VOI/muTZswJdo99wB8JXYcyVMqaBLFrbfWdDnqK7c0zMiRNXN8nhU7iuQpFTRpcuvXw1VXwd57AzweOY0Uih13hBNPBBiqE62lTipo0uT+9a9wIvUvfxk7iRSa8DvTlhtvjJ1E8pEKmjS5yy+H7t11IrWkr3dvgIe45ppwHFakNhU0aVJPPQVPPw2/+hW0bBk7jRSmy1iyBMaNi51D8k2DC5qZjTOzJWb2Sq1t25rZdDObn/rZMbXdzOw6M1tgZnPMbJ9chJfCc9ll0KULDBsWO4kUqh493gGeZuTIKsxaUFZWFjuS5Il0Wmi3AP032XYBMMPdewEzUrcBDgd6pS7Dgb9nFlOS4IUX4KGHYMmSUbRuramtpHGqqiqZNu1bQCnjx6+lqqoqdiTJEw0uaO7+JLBsk80DgAmp6xOAH9baPtGDmcA2ZrZDpmGl8JSVlW2Yl/Eb35iK2Qo+/viPG1ap1lRX0hj9+4eJi//4R9CRE6mR6W9CV3dfnLr+PtA1dX0n4N1a+y1MbfsSMxtuZhVmVrE0nJgkCVJVVYW7M2+eY3Yco0a1p0OH2Kmk0JnBhRfCG28AHBM7juSJrH21cXcHvBGPG+3u5e5e3lkT+iXWn/4UZtM/S+fESpYccwzsuivAr/G0P3kkiTItaB/UdCWmfi5JbV8EdK+1X7fUNilCb74JEyfCaaeFASEi2VBSAhdcANCHadNip5F8kGlBux84JXX9FOC+WttPTo123A9YXqtrUorMH/4ALVrUfPiIZM+gQQBvc/HFqJUmaQ3bnwI8A+xqZgvNbBhwOXCImc0Hvp+6DTANeAtYANwM/CyrqaWA7MzEiWG9sx00LEiyrEULgEuoqIB//jN2GonNPI++1pSXl3tFRUXsGJJFZrfSuvVJvPUWdO1a//4i6TJrzs47r6NNm3BqSDMNekwcM3ve3cvr20//9ZIz8+YBnMjIkSpmkkvVXHQRvPQS3H137CwSkwqa5MzvfgewivPOi51Ekm7gQNh9d7joIqiujp1GYlFBk5x4+WW44w6A6+jUKXYaSbqSkvAFau5cuO222GkkFhU0yVjt2UBqLr17/xP35XTrdkfseFIkfvSjMBv/xRfD2rWx00gMKmiSsZrZQGoujz3mwFFcfnkH3n33pdjxpEg0axZOEVmwAG6+OXYaiUGjHCVjZkbN79H69dCvH7z/fpiWaOutI4eTxKv9++cOBx4Yuh4XLID27SOHk6zQKEeJYupUqKgI35RVzKQplJaWbujqbtbMeP31o1i6NEy3JsVFLTTJWM035DVrwkizdu3C+UAlJbGTSTEyMwYOdO65B+bPh27dYieSTKmFJk3uxhvh7bfDN2MVM4npsstC9/dvfhM7iTQlFTTJimXL4Pe/h+9/Hw49NHYaKXZlZfDzn8OECbDjjodv6JLU6tbJpoImWfGb38DHH8OVV4a1qkRiu/BC6NgRFi8+n/XrwwhcrW6dbCpokgW9ufHGMAFx796xs4gEHTuGrkf4nk62LhLNYweQwhbGFF3PttuGLkeR2GpGPQbNaNnyRc49tzdHHhk1ljQBFTTJyJQpAN/hssvCN2KR2CorKze6/eyzsN9++sJVDNTlKI22ciWpiYcrGDo0dhqRuvXrBz/5CVxzDcDuseNIDqmgSaNdeCEsXgwwQsP0Ja/98Y/h/Ej4q1a2TjAVNGmUZ56BG26AkSMBnosdR2SLOnWqGSByEBMnxk4juaKZQiRtn38O++wDK1bAq69C+/ZfzKUnkq/Wr4eSkqfo2PHbzJ0L228fO5E0lGYKkZz5059CIbvhhppuHJH816wZwE/49NOangVJGhU0Scu8eXDJJXDccXDUUbHTiKTrDS6+GO66K1wkWVTQpMHWroXBg0Or7LrrYqcRaZxzzoG994YRI+DDD2OnkWzSeWjSYJdcAs8/D3feqeMPUphKS0tp2dKAvYDn6NHjET755AhN15YQaqFJgzz7bBglNnhwWOpepBBVVlamVlafzRVXtOTTT4+gWbMhmrw4ITTKUer16afQpw989hm8/DJ06LDx/bVXDBYpFNXVYXWIigqYPRu++lX9LuerJh3laGaVZvaymc02s4rUtm3NbLqZzU/91MRIBerMM8Ny9rfcEopZWVnZhm+0ZkZpaWnsiCJpKymBiROheXM46SRYty52IslUNrscD3T3PrWq6AXADHfvBcxI3ZYCM3EijBsXZgU56KCwraqqKtVtEy6bzp0nUii6dw8L086cCf/v/8VOI5nK5aCQAcD3UtcnAI8D5+fw9STL5s6FM86A734XLr44dhqR3Pjxj+Gxx+CKKwB+EDuOZCBbLTQHHjaz581seGpbV3dfnLr+PtC1rgea2XAzqzCziqVLl2YpjmRq1apwrlnbtrBgwb60aKEuRkmua66B8nKACbz5Zuw00ljZKmjfcvd9gMOBEWb2ndp3ejjKWueRVncf7e7l7l7euXPnLMWRTLjD0KHhJOrJk2HRogp1MUqitWoFU6cCVPOjH8Hq1bETSWNkpaC5+6LUzyXAPUBf4AMz2wEg9XNJNl5Lcu+SS+COO0IXzMEHx04j0jTCiP2TmDMHhg1Ds/IXoIwLmpm1MbN2NdeBQ4FXgPuBU1K7nQLcl+lrSe7deSdcdBGcfDKce27sNCJN7SEuuywsXPu738XOIunKxqCQrsA9qSXPmwOT3f0hM5sF3GFmw4Aq4PgsvJbk0AsvhEK2//5w001o9gQpSuefD2+8EQpar14waFDsRNJQGRc0d3+LMI/Mpts/BNRhlYfKysqoqqracLu0tJRHH63kiCPCulF33x2OKYgUI7MwlP+tt8Kx5NJS+Na3YqeShtDUV0Vo0/PIqqpWc+ihYfLhf/9b8zSKtGwZZuMvLQ2rSrz0UuxE0hAqaEVuxQqAf7F4MTz4IOy+u2YCEQHYbjuYPj2cunLooTB//pf/NjT3Y37RbPtF7JNP4MgjAb7OXXfBfj+9e3YAABG3SURBVPuF7TUtOJFiV1oaitq3vx3mfXznneqN/jZMB5rzilpoRWrlSujfH/77X+jU6WwOP1wtMpG67LZb6Ir/+GOAR3nnndiJZHNU0IpSOw47LMxfN2UKLF36V504LbIF++wTihp05jvfQbOJ5CkVtCITZhebzqxZcPvtYXorEalb7WNm++9vbL/9ID75BL7zHXjttdjpZFM6hlZE3nwzdDNCb+66C36geVhFtqiu48mvvBKOp4Wh/N+MkkvqphZakZg1K5wwvWwZwMEqZiKNtOee8PTTsO22ADO4/fbYiaSGCloRmDw5LAHTpk0YBALPxI4kUtB23hmeeQbgOU44AS69VHM/5gMVtARbuxbOPjtM3VNeHgaB7Lpr7FQi+au0tLTB52Butx3AIQwaFBYHPeYY6NGjt85Ti0gFLQHqOtlz4UI45JCwztPPfw4zZkDXOlekE5EalZWVaY34LS3dgUmTDDiLe+9dy+LF9/Pii7Vn4ana4uMlu1TQEuDLU1ntS+/e4bjZxIlw7bXQokXslCLJ80UBvJb//KcFXbuWsd9+cOWVUF0dO13xUUErQJubmuqjj+DUUwGm0qsXzJ4NgwfHTCpSPL75TXjxxTCS+Nxz4cADAXrGjlVUVNAK0KYtsrffrmTy5DCjwT/+AfB7nn46LH0BmptRpKl07gz33AO33FIzofEc/vKXcDxbck8FrcDNmweHHRYGfpSWQkUFwEUbdTFuWgA1E4hI7pjBKaeE89Xgcc47D/r0gccfjxysCKigFaj33oPTTgvnxDz7LFx/fRhGvPfesZOJCED37lBaOhI4irlz3+bAA6F163+nCp3kggpagVmyBOAydt4ZJkyAkSNhwYLws6QkdjoRqS0MGvknn37ak0sugdWr96N377Ay/FtvbbyvlqbJnApagXj7bRgxInQrwvkMGBDmkrv22tBvLyL5a+utw7lq8BXOOw+mTg3HuI8/PvSwQF2jlTXkP10qaHmsuhr+9S84+mj4ylfW8be/fc5nn93Mjjt+nylT4Ctfqftx6ZwcKiJNaRlXXBHmVT3vPHj44bAO4be/DTCQzz6Lna+wWT4t5FheXu4VYVRDUXv9dejX72qWLz8aKAM+oH37u5k79wx22ilyOBFpNDPbaLLjlSth7Fi47rrQC7PNNmGA15AhUF5uWmg3xcyed/fy+vZTCy0PuIfuw0svhb32CsPvly8/i4MPLuOOO2DNmq4sX65iJlLoNu09ad/eOPts4+23m9Gly4kcfjiMGROmqispeQezv2D2Tcya6ZhaA6iFFkmPHr15991dgEOBw4DQLXjAAaFf/ayzdsT9vZgRRSSCZcvg3nvhrrtg+vRwDluXLrBkyWTGjTuRgw+GHj1ip2xaDW2hqaA1gfXrQwts5swwtH7mTHjllfVAM9q3h4MPhieeuJBly/4BLATCNzmdLyZS3JYvhwceCMfSp0xZyvr1NSPA5tOmzctcfvkx9O0bena22ipq1JxSQYtg/fpwfti8eeGkyprLq6/CqlVhn222CQeBH3rotzz11O/p10/zLIpI/dzDZ8kjj4STtB94YAnV1V1S966hZcs3OP74r7PHHmy49OwJzROwjHPeFDQz6w9cC5QAY9z98s3tm88Fbd06WLo0nAf2wQfhsnBhOJBbWRl+vvMOfP75F4/p0iWc+LznnjBp0nl8+OE/gTcAVwtMRDLiHj6DnnuuZnKFJ/nss55A9w37NG9ec4I3zJp1J6tWvQJUAe+z/fYlzJr1AJ07h9ZdWVnZRqcK5NNnVF4UNDMrIXyCH0LoS5sFDHT3uXXtn42C5h6Kz2efffmyZs2Xt61aBStWhKZ97Z/33vsYn33WAugIdAG2o64xNF26QFlZuPTsGX7uuiucfPI3WLjwhQ375dMvh4gk14oV4RBHv36nMmrULVRVQVUV/Oc/CzHrVudCpB06wPLlb7D//rvQoUO4ffvtoznnnOF06ADt24dL27bhnLpWrbZ8adkyTPTQrFmYCixT+VLQ9gcudvfDUrdHAbj7H+vav02bct955wqqq2n0JZO306pV+E/r0AHmz5/FQQftyzbbhKLVtesXly5d4Mc//h7vvTcL+LTO51IBE5GY6mpxvfFGJQsXhh6mAQNOY+lSCF/Yu9K6dU/23/8oVqwIRXHBgiVUV7cB2mSYZC1m1bRu3YrmzdlwKSn58vVmzb4ogjXXmzWDior8KGjHAv3d/Sep24OBfu4+stY+w4HhqZt7AoU+01kn4H+xQ2RI7yE/FPp7KPT8oPeQL3Z193b17RT9cKG7jwZGA5hZRUOqcD7Te8gPeg/xFXp+0HvIF2bWoGNRuT6xehG1j1BCt9Q2ERGRrMp1QZsF9DKznmbWEjgBuD/HrykiIkUop12O7r7OzEYC/yYM2x/n7q9u4SGjc5mnieg95Ae9h/gKPT/oPeSLBr2HvDqxWkREpLE0ObGIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCSCCpqIiCRC89gBauvUqZOXlZXFjiECwPufrtvo9vat8+rPRaRoPP/88/9z98717ZdXf6FlZWVUVFTEjiECwOUv/m+j2xfs3SlSkuKm/wcxs6qG7JdXBU1ECo8KjuQLHUMTEZFEyLigmVl3M3vMzOaa2atmdlZq+7ZmNt3M5qd+dsw8roiISN2y0UJbB5zj7nsA+wEjzGwP4AJghrv3AmakbouIiORExgXN3Re7+wup6yuBecBOwABgQmq3CcAPM30tERGRzcnqoBAzKwP2Bp4Furr74tRd7wNdN/OY4cBwgB49emQzjiRAPg84yOdsIsUoa4NCzKwtcBfwC3dfUfs+d3fA63qcu49293J3L+/cud7TDEREROqUlYJmZi0IxWySu9+d2vyBme2Qun8HYEk2XktERKQu2RjlaMBYYJ67X1XrrvuBU1LXTwHuy/S1RERENicbx9AOAAYDL5vZ7NS2C4HLgTvMbBhQBRyfhdcSERGpU8YFzd2fBmwzdx+c6fOLiIg0hGYKERGRRNBcjiJZomH8InGphSYiIomggiYiIomgLkfJWO2uNnWz5adNu0M3pf83SQK10EREJBFU0EREJBFU0EREJBF0DE2Kio73iSSXWmgiIpIIaqFJVDoZWUSyRQWtCKmISFPS75s0FRU0KVr6oP2C/i0kCVTQRCTRNBCoeKigFSDN+iASRy5bsmolZ04FTXJKf6Qi0lRU0EQaqL6WsYjEpYImIpIl6pGISwVNRNKWzdZqzCKgApQsKmgiUi91t0ohyHlBM7P+wLVACTDG3S/P9WtKPPrgEylMSWit5rSgmVkJcANwCLAQmGVm97v73Fy+biFIwi9P0qgYSz7RZ0T6ct1C6wsscPe3AMzsNmAAkPcFLfa5XvplFpFsyuZnSr5+Ppm75+7JzY4F+rv7T1K3BwP93H1kXfuXt2vnFd/4RtZe/51P1m7x/h5tW+TksQ2x6fNv+nxbuj/TbPW9drrqy7Ol18rk36Ex0smaqUyzZlO677u+/4dMHtvUvwO1H5/pc+cya7Z/N3P579zUf7f2xBPPu3t5fftFHxRiZsOB4QC9t9oqq8+d7i9nth7bEJn8AmSaLZMPq8bk2dJrZfpHl66mLDJNWTzrk0kRSPfxuS5gmRTnTJ87l1mz/buZy6KT6wLWWLkuaIuA7rVud0tt28DdRwOjAcrLy53HH89xpC/0qHU93Sb05DSPt9T3fD22eG966su2aZZ030t9z7clm75WfVliZs1UplmzKd33ncnvY7b/jze9f9Ns6RweyPT3K5dZs/27mW62TGTzuepk1qDdcl3QZgG9zKwnoZCdAJyY49dslHR/mfKlz1hE4sn0cyCXnyPF+BmV04Lm7uvMbCTwb8Kw/XHu/mouX1NEZHOK8UO+mOT8GJq7TwOm5fp1pPE2/SPX8HURKUTRB4WIZELFuG5qiUgxUkGTJpXrD9okFbgt/Vsl9X3lm2xnLaT3XohU0KReSSoSTUn/biJNq1nsACIiItmggiYiIomggiYiIomgY2hFQCeNi0gxUEETkUTRF7LipS5HERFJBLXQEqipz/USEckHKmgiUtD0BUtqqMtRREQSQS00STR9excpHmqhiYhIIqiFJgUlyS2uJL83kaaggiYieUWFXRpLXY4iIpIIKmgiIpII6nKUvKLuJhFpLLXQREQkEdRCEylAWg07N9RDUNgyaqGZ2Z/N7DUzm2Nm95jZNrXuG2VmC8zsdTM7LPOoIiIim5dpC206MMrd15nZFcAo4Hwz2wM4AfgasCPwiJnt4u7VGb6eiNRBLQuRDFto7v6wu69L3ZwJdEtdHwDc5u5r3P1tYAHQN5PXEhER2ZJsDgoZCvwrdX0n4N1a9y1MbfsSMxtuZhVmVrF06dIsxhERkWJSb5ejmT0CbF/HXb929/tS+/waWAdMSjeAu48GRgOUl5d7uo8XERGBBhQ0d//+lu43s1OBI4GD3b2mIC0CutfarVtqm4iISE5kNCjEzPoDvwK+6+6f1rrrfmCymV1FGBTSC3guk9cSKXTFMtReA1QklkxHOf4V2AqYbmYAM939p+7+qpndAcwldEWO0AhHERHJpYwKmrvvvIX7LgUuzeT5RUREGkpTX4mISCKooImISCLYFwMT4zOzlcDrsXNkqBNQ6Ef79R7yQ6G/h0LPD3oP+WJXd29X3075Njnx6+5eHjtEJsysQu8hPr2H+Ao9P+g95Aszq2jIfupyFBGRRFBBExGRRMi3gjY6doAs0HvID3oP8RV6ftB7yBcNeg95NShERESksfKthSYiItIoeVvQzOwcM3MzK7iJ4czsktQq3rPN7GEz2zF2pnRtaTXyQmBmx5nZq2a23swKaoSXmfVPrfS+wMwuiJ0nXWY2zsyWmNkrsbM0lpl1N7PHzGxu6vforNiZ0mVmrczsOTN7KfUefhc7U2OYWYmZvWhmD9S3b14WNDPrDhwKvBM7SyP92d17u3sf4AHgt7EDNcJ0YE937w28QViNvJC8AhwDPBk7SDrMrAS4ATgc2AMYmFoBvpDcAvSPHSJD64Bz3H0PYD9gRAH+P6wBDnL3vYA+QH8z2y9ypsY4C5jXkB3zsqABVxNm8S/IA3zuvqLWzTYU4PvYwmrkBcHd57l7IZ6k3xdY4O5vufvnwG2EFeALhrs/CSyLnSMT7r7Y3V9IXV9J+ECtc5HifOXBJ6mbLVKXgvosMrNuwP8BYxqyf94VNDMbACxy95diZ8mEmV1qZu8CgyjMFlpttVcjl9xq8Grv0jTMrAzYG3g2bpL0pbrrZgNLgOnuXmjv4RpC42Z9Q3aOMlPIllbBBi4kdDfmtfpW8nb3XwO/NrNRwEjgoiYN2AC5Xo081xqSXyQTZtYWuAv4xSY9LwUhtWxXn9Qx8HvMbE93L4hjm2Z2JLDE3Z83s+815DFRCtrmVsE2s68DPYGXUuurdQNeMLO+7v5+E0asV30redcyCZhGHha0Rq5GnjfS+D8oJFrtPU+YWQtCMZvk7nfHzpMJd//YzB4jHNssiIIGHAD8wMyOAFoB7c3sVnc/aXMPyKsuR3d/2d27uHuZu5cRulv2ybdiVh8z61Xr5gDgtVhZGqvWauQ/2GQ1csmtWUAvM+tpZi2BEwgrwEsTsvCNeiwwz92vip2nMcysc83oZDPbGjiEAvoscvdR7t4tVQtOAB7dUjGDPCtoCXK5mb1iZnMI3acFN+SXsBp5O8Jq5LPN7MbYgdJhZkeb2UJgf+BBM/t37EwNkRqIMxL4N2Egwh3u/mrcVOkxsynAM8CuZrbQzIbFztQIBwCDgYNSv/+zUy2FQrID8Fjqc2gW4RhavUPfC5lmChERkURQC01ERBJBBU1ERBJBBU1ERBJBBU1ERBJBBU1ERBJBBU1ERBJBBU1ERBJBBU1ERBLh/wO0WhcdRszyBgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"values, edges = minkit.data_plotting_arrays(data, bins=100)\n",
"\n",
"gx, pdf_values = minkit.pdf_plotting_arrays(g, values, edges)\n",
"\n",
"fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(7, 6), gridspec_kw={'height_ratios': (3, 1)}, sharex=True)\n",
"\n",
"centers = 0.5 * (edges[1:] + edges[:-1])\n",
"\n",
"ax0.hist(centers, edges, weights=values, histtype='step', color='k');\n",
"ax0.plot(gx, pdf_values, color='blue');\n",
"ax1.plot(x.bounds, (0, 0), color='red');\n",
"\n",
"binned_ds = minkit.BinnedDataSet.from_ndarray(edges, x, values)\n",
"\n",
"pdf_r = values.sum() * g.evaluate_binned(binned_ds).as_ndarray() # evaluate_binned is normalized to unity\n",
"\n",
"res = values - pdf_r\n",
"\n",
"pos = res > 0\n",
"\n",
"w = edges[1] - edges[0]\n",
"\n",
"ax1.bar(centers[pos], res[pos], width=w, color='skyblue');\n",
"ax1.bar(centers[~pos], res[~pos], width=w, color='skyblue');\n",
"ax1.set_xlim(*x.bounds)\n",
"ax1.set_ylim(-25, +25);"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}