{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# iMinuit for LSQ\n",
"\n",
"\n",
"In this notebook we will be using iminuit to fit with LSQ.\n",
"\n",
"iMinuit: \n",
"https://iminuit.readthedocs.io/en/latest/index.html#\n",
"\n",
"Here below a quick summary of: \n",
"https://iminuit.readthedocs.io/en/latest/tutorials.html \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"from iminuit import Minuit"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD5CAYAAAAp8/5SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWmklEQVR4nO3df2xdZ33H8fcH1xTTdrgQ8yNO2jCtRJQfbbqrAAqiLbA4rYAUVI1UUAqqFonBBqyyRJhEp/IH26wxDQRkGY1KES0w6poItXWzUVZ+LF2vm1C3DYasFBo7UkxTtwU8SMJ3f9zjcute+55rH99rP/fzkq5y7nOee+/3adxPjp/z3HMUEZiZWbqe0+oCzMxsaTnozcwS56A3M0ucg97MLHEOejOzxDnozcwSd0q9DpKeB9wNnJr1/2ZEXDurz6nAjcCfAo8B746IR7J9O4CrgZPAX0fEcL3PXLVqVaxbt66hgZiZtbORkZFfRkRPrX11gx74LfDmiPiVpE7g+5Juj4h9VX2uBh6PiD+RtA34B+Ddks4FtgGvAlYD/yHpFRFxcr4PXLduHeVyOUdpZmYGIOnnc+2rO3UTFb/KnnZmj9nfstoKfDnb/ibwFknK2r8WEb+NiJ8Bh4CNDdZvZmaLkGuOXlKHpAPAUWBvRNwzq0sv8ChARJwAngBeVN2eOZy11fqM7ZLKksqTk5ONjcLMzOaUK+gj4mREnA+sATZKevWsLqr1snnaa33GrogoRUSpp6fmNJOZmS1AQ6tuImIK+C6wZdauw8BaAEmnAC8AjlW3Z9YAEwus1czMFqBu0EvqkdSdbXcBbwV+PKvbHuCqbPty4DtRuVraHmCbpFMlvRw4B/ifooo3M7P68qy6eRnwZUkdVP5h+EZEfFvSdUA5IvYA1wNfkXSIypH8NoCIeFDSN4CHgBPAh+qtuDEzW8mG9o8zMDzGxNQ0q7u76O9bz2Ubap6abBotx8sUl0ql8PJKM1tphvaPs2NwlOnjfzie7ers4NPves2Sh72kkYgo1drnb8aamRVkYHjsGSEPMH38JAPDYy2qqMJBb2ZWkImp6Ybam8VBb2ZWkNXdXQ21N4uD3sysIP196+nq7HhGW1dnB/1961tUUUWeVTdmZpbDzAnX5bbqxkFvZlagyzb0tjzYZ/PUjZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJq3sJBElrgRuBlwK/B3ZFxL/M6tMPvKfqPV8J9ETEMUmPAE8BJ4ETc10Y38zMlkaea92cAK6JiPsknQGMSNobEQ/NdIiIAWAAQNLbgY9FxLGq97g4In5ZZOFmZpZP3ambiDgSEfdl208BB4H5rthzBXBzMeWZmdliNTRHL2kdsAG4Z479zwe2ALdUNQdwp6QRSdvnee/tksqSypOTk42UZWZm88gd9JJOpxLgH42IJ+fo9nbgB7OmbTZFxAXAJcCHJL2p1gsjYldElCKi1NPTk7csMzOrI1fQS+qkEvJfjYjBebpuY9a0TURMZH8eBW4FNi6sVDMzW4i6QS9JwPXAwYj4zDz9XgBcCHyrqu207AQukk4DNgMPLLZoMzPLL8+qm03AlcCopANZ2yeAswAiYmfW9k7gzoj4ddVrXwLcWvm3glOAmyLijiIKNzOzfOoGfUR8H1COfjcAN8xqexg4b4G1mZlZAfzNWDOzxPnm4GaWnKH94wwMjzExNc3q7i76+9Yvuxt2N5OD3sySMrR/nB2Do0wfPwnA+NQ0OwZHAdo27D11Y2ZJGRgeezrkZ0wfP8nA8FiLKmo9B72ZJWViarqh9nbgoDezpKzu7mqovR046M0sKf196+nq7HhGW1dnB/1961tUUev5ZKyZJWXmhKtX3fyBg97MknPZht62DvbZPHVjZpY4B72ZWeI8dWNmS8LfTl0+HPRmVjh/O3V58dSNmRXO305dXhz0ZlY4fzt1eXHQm1nh/O3U5SXPrQTXSrpL0kFJD0r6SI0+F0l6QtKB7PHJqn1bJI1JOiTp40UPwMyWH387dXnJczL2BHBNRNyX3f91RNLeiHhoVr/vRcTbqhskdQCfB/4MOAzcK2lPjdeaWUL87dTlJc+tBI8AR7LtpyQdBHqBPGG9ETiU3VIQSV8DtuZ8rZmtYP526vLR0By9pHXABuCeGrvfIOlHkm6X9KqsrRd4tKrP4ayt1ntvl1SWVJ6cnGykLDMzm0fuoJd0OnAL8NGIeHLW7vuAsyPiPOBzwNDMy2q8VdR6/4jYFRGliCj19PTkLcvMzOrIFfSSOqmE/FcjYnD2/oh4MiJ+lW3fBnRKWkXlCH5tVdc1wMSiqzYzs9zyrLoRcD1wMCI+M0efl2b9kLQxe9/HgHuBcyS9XNJzgW3AnqKKNzOz+vKsutkEXAmMSjqQtX0COAsgInYClwMflHQCmAa2RUQAJyR9GBgGOoDdEfFgwWMwM7N5qJLHy0upVIpyudzqMszMVgxJIxFRqrXP34w1M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwSl+dWgmsl3SXpoKQHJX2kRp/3SLo/e/xQ0nlV+x6RNCrpgCTfTcTMrMny3ErwBHBNRNwn6QxgRNLeiHioqs/PgAsj4nFJlwC7gNdV7b84In5ZXNlmlsfQ/nEGhseYmJpmdXcX/X3ruWxDb6vLsiarG/QRcQQ4km0/Jekg0As8VNXnh1Uv2QesKbhOM2vQ0P5xdgyOMn38JADjU9PsGBwFcNi3mYbm6CWtAzYA98zT7Wrg9qrnAdwpaUTS9nnee7uksqTy5ORkI2WZWQ0Dw2NPh/yM6eMnGRgea1FF1ip5pm4AkHQ6cAvw0Yh4co4+F1MJ+jdWNW+KiAlJLwb2SvpxRNw9+7URsYvKlA+lUmn53bHcbIWZmJpuqN3SleuIXlInlZD/akQMztHntcCXgK0R8dhMe0RMZH8eBW4FNi62aDOrb3V3V0Ptlq48q24EXA8cjIjPzNHnLGAQuDIiflLVflp2AhdJpwGbgQeKKNzM5tfft56uzo5ntHV1dtDft75FFVmr5Jm62QRcCYxKOpC1fQI4CyAidgKfBF4EfKHy7wInIqIEvAS4NWs7BbgpIu4odARmVtPMCVevujFFLL/p8FKpFOWyl9ybmeUlaSQ7wH4WfzPWzCxxDnozs8Q56M3MEuegNzNLnIPezCxxDnozs8Q56M3MEuegNzNLnIPezCxxDnozs8Q56M3MEuegNzNLnIPezCxxDnozs8TlvpWg2Uo2tH/c12W3tuWgt+QN7R9nx+Do0zfKHp+aZsfgKIDD3tpCnlsJrpV0l6SDkh6U9JEafSTps5IOSbpf0gVV+66S9NPscVXRAzCrZ2B47OmQnzF9/CQDw2MtqsisufIc0Z8AromI+7L7v45I2hsRD1X1uQQ4J3u8Dvgi8DpJLwSuBUpAZK/dExGPFzoKs3lMTE031F40TxtZq9U9oo+IIxFxX7b9FHAQmP1TuhW4MSr2Ad2SXgb0AXsj4lgW7nuBLYWOwKyO1d1dDbUXaWbaaHxqmuAP00ZD+8eX/LPNZjS06kbSOmADcM+sXb3Ao1XPD2dtc7XXeu/tksqSypOTk42UZTav/r71dHV2PKOtq7OD/r71S/7Znjay5SB30Es6HbgF+GhEPDl7d42XxDztz26M2BURpYgo9fT05C3LrK7LNvTy6Xe9ht7uLgT0dnfx6Xe9pinTJ62eNjKDnKtuJHVSCfmvRsRgjS6HgbVVz9cAE1n7RbPav7uQQs0W47INvS2ZF1/d3cV4jVBvxrSR2Yw8q24EXA8cjIjPzNFtD/C+bPXN64EnIuIIMAxslnSmpDOBzVmbWVto5bSR2Yw8R/SbgCuBUUkHsrZPAGcBRMRO4DbgUuAQ8BvgA9m+Y5I+Bdybve66iDhWXPlmy9vMbxFedWOtpIiaU+YtVSqVolwut7oMM7MVQ9JIRJRq7fO1bszMEuegNzNLnIPezCxxDnozs8Q56M3MEuegNzNLnIPezCxxDnozs8T5DlPWNL4uu1lrOOitKXw7P7PW8dSNNYWvy27WOg56awpfl92sdRz01hStvJ2fWbtz0FtT+LrsZq3jk7HWFL4uu1nrOOitaVp1Oz+zdlc36CXtBt4GHI2IV9fY3w+8p+r9Xgn0ZHeXegR4CjgJnJjrovhmZrZ08szR3wBsmWtnRAxExPkRcT6wA/ivWbcLvDjb75A3M2uBukEfEXcDee/zegVw86IqMjOzQhW26kbS86kc+d9S1RzAnZJGJG2v8/rtksqSypOTk0WVZWbW9opcXvl24Aezpm02RcQFwCXAhyS9aa4XR8SuiChFRKmnp6fAsszM2luRQb+NWdM2ETGR/XkUuBXYWODnmZlZDoUEvaQXABcC36pqO03SGTPbwGbggSI+z8zM8suzvPJm4CJglaTDwLVAJ0BE7My6vRO4MyJ+XfXSlwC3Spr5nJsi4o7iSjczszzqBn1EXJGjzw1UlmFWtz0MnLfQwszMrBi+1o2ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpa4ukEvabeko5Jq3gZQ0kWSnpB0IHt8smrfFkljkg5J+niRhZuZWT55juhvALbU6fO9iDg/e1wHIKkD+DxwCXAucIWkcxdTrJmZNa5u0EfE3cCxBbz3RuBQRDwcEb8DvgZsXcD7mJnZIhQ1R/8GST+SdLukV2VtvcCjVX0OZ201SdouqSypPDk5WVBZZmZWRNDfB5wdEecBnwOGsnbV6BtzvUlE7IqIUkSUenp6CijLzMyggKCPiCcj4lfZ9m1Ap6RVVI7g11Z1XQNMLPbzzMysMYsOekkvlaRse2P2no8B9wLnSHq5pOcC24A9i/08MzNrzCn1Oki6GbgIWCXpMHAt0AkQETuBy4EPSjoBTAPbIiKAE5I+DAwDHcDuiHhwSUZhZmZzUiWTl5dSqRTlcrnVZZiZrRiSRiKiVGufvxlrZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeLqXo/e0jK0f5yB4TEmpqZZ3d1Ff996Ltsw5618zSwBDvo2MrR/nB2Do0wfPwnA+NQ0OwZHARz2ZgmrO3Ujabeko5IemGP/eyTdnz1+KOm8qn2PSBqVdECS7yTSYgPDY0+H/Izp4ycZGB5rUUVm1gx55uhvALbMs/9nwIUR8VrgU8CuWfsvjojz57rziTXPxNR0Q+1mloa6QR8RdwPH5tn/w4h4PHu6D1hTUG1WsNXdXQ21m1kail51czVwe9XzAO6UNCJp+3wvlLRdUllSeXJysuCyDKC/bz1dnR3PaOvq7KC/b32LKjKzZijsZKyki6kE/RurmjdFxISkFwN7Jf04+w3hWSJiF9m0T6lUWn53LE/AzAlXr7oxay+FBL2k1wJfAi6JiMdm2iNiIvvzqKRbgY1AzaC35rhsQ6+D3azNLDroJZ0FDAJXRsRPqtpPA54TEU9l25uB6xb7eSnwWnYza6a6QS/pZuAiYJWkw8C1QCdAROwEPgm8CPiCJIAT2QqblwC3Zm2nADdFxB1LMIYVxWvZzazZFLH8psNLpVKUy2kuu9/0999hvMZyxt7uLn7w8Te3oCIzS4GkkbmWsftaN03mtexm1mwO+ibzWnYzazYHfZN5LbuZNZsvatZkXstuZs3moG8Br2U3s2by1I2ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klzkFvZpY4B72ZWeIc9GZmiXPQm5klLlfQS9ot6aikB+bYL0mflXRI0v2SLqjad5Wkn2aPq4oq3MzM8sl7RH8DsGWe/ZcA52SP7cAXASS9kMqtB19H5cbg10o6c6HFmplZ43IFfUTcDRybp8tW4Mao2Ad0S3oZ0AfsjYhjEfE4sJf5/8EwM7OCFTVH3ws8WvX8cNY2V/uzSNouqSypPDk5WVBZZmZWVNCrRlvM0/7sxohdEVGKiFJPT09BZZmZWVE3HjkMrK16vgaYyNovmtX+3YI+c1GG9o/7Lk9m1haKOqLfA7wvW33zeuCJiDgCDAObJZ2ZnYTdnLW11ND+cXYMjjI+NU0A41PT7BgcZWj/eKtLMzMrXK4jekk3UzkyXyXpMJWVNJ0AEbETuA24FDgE/Ab4QLbvmKRPAfdmb3VdRMx3UrcpBobHmD5+8hlt08dPMjA85qN6M0tOrqCPiCvq7A/gQ3Ps2w3sbry0pTMxNd1Qu5nZStaW34xd3d3VULuZ2UrWlkHf37eers6OZ7R1dXbQ37e+RRWZmS2dolbdrCgz8/BedWNm7aAtgx4qYe9gN7N20JZTN2Zm7cRBb2aWOAe9mVniHPRmZolz0JuZJU6VL7UuL5ImgZ836eNWAb9s0me1gse3snl8K1ezx3Z2RNS89O+yDPpmklSOiFKr61gqHt/K5vGtXMtpbJ66MTNLnIPezCxxDnrY1eoClpjHt7J5fCvXshlb28/Rm5mlzkf0ZmaJc9CbmSWubYJe0hZJY5IOSfp4jf2nSvp6tv8eSeuaX+XC5Rjf30h6SNL9kv5T0tmtqHOh6o2vqt/lkkLSsljWllee8Un68+zv8EFJNzW7xoXK8bN5lqS7JO3Pfj4vbUWdCyVpt6Sjkh6YY78kfTYb//2SLmh2jURE8g+gA/hf4I+B5wI/As6d1ecvgZ3Z9jbg662uu+DxXQw8P9v+YGrjy/qdAdwN7ANKra674L+/c4D9wJnZ8xe3uu4Cx7YL+GC2fS7wSKvrbnCMbwIuAB6YY/+lwO2AgNcD9zS7xnY5ot8IHIqIhyPid8DXgK2z+mwFvpxtfxN4iyQ1scbFqDu+iLgrIn6TPd0HrGlyjYuR5+8P4FPAPwL/18ziCpBnfH8BfD4iHgeIiKNNrnGh8owtgD/Ktl8ATDSxvkWLiLuBY/N02QrcGBX7gG5JL2tOdRXtEvS9wKNVzw9nbTX7RMQJ4AngRU2pbvHyjK/a1VSOMFaKuuOTtAFYGxHfbmZhBcnz9/cK4BWSfiBpn6QtTatucfKM7e+A90o6DNwG/FVzSmuaRv//LFy73GGq1pH57HWlefosV7lrl/ReoARcuKQVFWve8Ul6DvDPwPubVVDB8vz9nUJl+uYiKr+NfU/SqyNiaolrW6w8Y7sCuCEi/knSG4CvZGP7/dKX1xQtz5Z2OaI/DKyter6GZ/96+HQfSadQ+RVyvl/HlpM840PSW4G/Bd4REb9tUm1FqDe+M4BXA9+V9AiVedA9K+iEbN6fz29FxPGI+BkwRiX4l7s8Y7sa+AZARPw38DwqFwRLRa7/P5dSuwT9vcA5kl4u6blUTrbumdVnD3BVtn058J3IzqSsAHXHl01t/CuVkF8p87sz5h1fRDwREasiYl1ErKNyDuIdEVFuTbkNy/PzOUTlhDqSVlGZynm4qVUuTJ6x/QJ4C4CkV1IJ+smmVrm09gDvy1bfvB54IiKONLOAtpi6iYgTkj4MDFNZBbA7Ih6UdB1Qjog9wPVUfmU8ROVIflvrKm5MzvENAKcD/56dY/5FRLyjZUU3IOf4Vqyc4xsGNkt6CDgJ9EfEY62rOp+cY7sG+DdJH6MypfH+FXSQhaSbqUyprcrOM1wLdAJExE4q5x0uBQ4BvwE+0PQaV9B/TzMzW4B2mboxM2tbDnozs8Q56M3MEuegNzNLnIPezCxxDnozs8Q56M3MEvf/tCDLT+gt98gAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Generate random numbers on a straight line\n",
"def line(x, a, b):\n",
" return a + x * b\n",
"\n",
"data_x = np.linspace(0, 1, 10)\n",
"# precomputed random numbers from a normal distribution\n",
"offsets = np.array([-0.49783783, -0.33041722, -1.71800806, 1.60229399, 1.36682387,\n",
" -1.15424221, -0.91425267, -0.03395604, -1.27611719, -0.7004073 ])\n",
"data_y = line(data_x, 1, 2) + 0.1 * offsets # generate some data points with random offsets\n",
"plt.plot(data_x, data_y, \"o\")\n",
"plt.xlim(-0.1, 1.1);"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Definition of the loss function you want to minimize (here the LSQ)\n",
"\n",
"def least_squares(a, b):\n",
" yvar = 0.01\n",
" return sum((data_y - line(data_x, a, b)) ** 2 / yvar)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Minuit"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td/>\n",
"<th title=\"Variable name\">\n",
"Name\n",
"</th>\n",
"<th title=\"Value of parameter\">\n",
"Value\n",
"</th>\n",
"<th title=\"Hesse error\">\n",
"Hesse Error\n",
"</th>\n",
"<th title=\"Minos lower error\">\n",
"Minos Error-\n",
"</th>\n",
"<th title=\"Minos upper error\">\n",
"Minos Error+\n",
"</th>\n",
"<th title=\"Lower limit of the parameter\">\n",
"Limit-\n",
"</th>\n",
"<th title=\"Upper limit of the parameter\">\n",
"Limit+\n",
"</th>\n",
"<th title=\"Is the parameter fixed in the fit\">\n",
"Fixed\n",
"</th>\n",
"</tr>\n",
"<tr style=\"background-color:#FFFFFF;\">\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"a\n",
"</td>\n",
"<td>\n",
"5.00\n",
"</td>\n",
"<td>\n",
"0.10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"yes\n",
"</td>\n",
"</tr>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td>\n",
"1\n",
"</td>\n",
"<td>\n",
"b\n",
"</td>\n",
"<td>\n",
"5.00\n",
"</td>\n",
"<td>\n",
"0.10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"------------------------------------------------------------------------------------------\n",
"| 0 | a | 5.00 | 0.10 | | | 0 | | yes |\n",
"| 1 | b | 5.00 | 0.10 | | | 0 | 10 | |\n",
"------------------------------------------------------------------------------------------"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Minuit uses introspection to detect the parameter names of your function - it means \n",
"#\n",
"# A complete instance of Minuit contains the starting point of each parameter (<name>) and \n",
"# the step size (<err_name>) for the dimension corresponding to that parameter.\n",
"#\n",
"# To correctly compute the uncertainties, it needs to know if you are doing a max likelihood \n",
"# or a least squares fit:\n",
"# errordef = 0.5 for negative log-likelihood functions\n",
"# errordef = 1 for least-squares functions\n",
"#\n",
"# You can specify the limits on the parameters to be fit (\"limit_varname\")\n",
"# lower limit: use limit_<name> = (<value>, None) or (<value>, float(\"infinity\"))\n",
"# upper limit: use limit_<name> = (None, <value>) or (-float(\"infinity\"), <value>)\n",
"# two-sided limit: use limit_<name> = (<min_value>, <max_value>)\n",
"#\n",
"# You can fix some of the parameters (ignore that dimension in the fit) setting <fix_name> = True/False\n",
"# fix_a=True,\n",
"\n",
"m = Minuit(least_squares, \n",
" a=5, b=5,\n",
" fix_a = True,\n",
" error_a=0.1, error_b=0.1,\n",
" limit_a=(0, None), limit_b=(0, 10),\n",
" errordef=1)\n",
"\n",
"m.get_param_states()\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"# Once Minuit is constructed you can still fix/release parameters as:\n",
"m.fixed[\"a\"] = False\n",
"m.fixed[\"b\"] = True\n",
"\n",
"# Trick to run over all parameters:\n",
"for key in m.fixed:\n",
" m.fixed[key] = False\n",
"m.migrad()\n",
"\n",
"# To change the value of a fixed parameter (or reset it to a different value) you can access it as:\n",
"m.values[\"a\"] = 0.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# MIGRAD"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr>\n",
"<td colspan=\"2\" title=\"Minimum value of function\">\n",
"FCN = 10.39\n",
"</td>\n",
"<td align=\"center\" colspan=\"3\" title=\"No. of calls in last algorithm and total number of calls\">\n",
"Ncalls = 42 (106 total)\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td colspan=\"2\" title=\"Estimated distance to minimum and target threshold\">\n",
"EDM = 4.32E-06 (Goal: 1E-05)\n",
"</td>\n",
"<td align=\"center\" colspan=\"3\" title=\"Increase in FCN which corresponds to 1 standard deviation\">\n",
"up = 1.0\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td align=\"center\" title=\"Validity of the migrad call\">\n",
"Valid Min.\n",
"</td>\n",
"<td align=\"center\" title=\"Validity of parameters\">\n",
"Valid Param.\n",
"</td>\n",
"<td align=\"center\" title=\"Is EDM above goal EDM?\">\n",
"Above EDM\n",
"</td>\n",
"<td align=\"center\" colspan=\"2\" title=\"Did last migrad call reach max call limit?\">\n",
"Reached call limit\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td align=\"center\" colspan=\"2\" style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td align=\"center\" title=\"Did Hesse fail?\">\n",
"Hesse failed\n",
"</td>\n",
"<td align=\"center\" title=\"Has covariance matrix\">\n",
"Has cov.\n",
"</td>\n",
"<td align=\"center\" title=\"Is covariance matrix accurate?\">\n",
"Accurate\n",
"</td>\n",
"<td align=\"center\" title=\"Is covariance matrix positive definite?\">\n",
"Pos. def.\n",
"</td>\n",
"<td align=\"center\" title=\"Was positive definiteness enforced by Minuit?\">\n",
"Forced\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td align=\"center\" style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"</table>\n",
"<table>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td/>\n",
"<th title=\"Variable name\">\n",
"Name\n",
"</th>\n",
"<th title=\"Value of parameter\">\n",
"Value\n",
"</th>\n",
"<th title=\"Hesse error\">\n",
"Hesse Error\n",
"</th>\n",
"<th title=\"Minos lower error\">\n",
"Minos Error-\n",
"</th>\n",
"<th title=\"Minos upper error\">\n",
"Minos Error+\n",
"</th>\n",
"<th title=\"Lower limit of the parameter\">\n",
"Limit-\n",
"</th>\n",
"<th title=\"Upper limit of the parameter\">\n",
"Limit+\n",
"</th>\n",
"<th title=\"Is the parameter fixed in the fit\">\n",
"Fixed\n",
"</th>\n",
"</tr>\n",
"<tr style=\"background-color:#FFFFFF;\">\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"a\n",
"</td>\n",
"<td>\n",
"0.99\n",
"</td>\n",
"<td>\n",
"0.06\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td>\n",
"1\n",
"</td>\n",
"<td>\n",
"b\n",
"</td>\n",
"<td>\n",
"1.95\n",
"</td>\n",
"<td>\n",
"0.10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"------------------------------------------------------------------\n",
"| FCN = 10.39 | Ncalls=42 (106 total) |\n",
"| EDM = 4.32E-06 (Goal: 1E-05) | up = 1.0 |\n",
"------------------------------------------------------------------\n",
"| Valid Min. | Valid Param. | Above EDM | Reached call limit |\n",
"------------------------------------------------------------------\n",
"| True | True | False | False |\n",
"------------------------------------------------------------------\n",
"| Hesse failed | Has cov. | Accurate | Pos. def. | Forced |\n",
"------------------------------------------------------------------\n",
"| False | True | True | True | False |\n",
"------------------------------------------------------------------\n",
"------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"------------------------------------------------------------------------------------------\n",
"| 0 | a | 0.99 | 0.06 | | | 0 | | |\n",
"| 1 | b | 1.95 | 0.10 | | | 0 | 10 | |\n",
"------------------------------------------------------------------------------------------"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Migrad performs Variable-Metric Minimization. It combines a steepest-descends algorithm along with \n",
"# line search strategy. Migrad is very popular in high energy physics field because of its robustness.\n",
"m.migrad()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"FMin(fval=10.387015571126001, edm=4.319932643302712e-06, tolerance=0.1, nfcn=42, ncalls=106, up=1.0, is_valid=True, has_valid_parameters=True, has_accurate_covar=True, has_posdef_covar=True, has_made_posdef_covar=False, hesse_failed=False, has_covariance=True, is_above_max_edm=False, has_reached_call_limit=False)\n"
]
}
],
"source": [
"from pprint import pprint\n",
"\n",
"# To understand the results of the fit there are two dict-like objects:\n",
"# The first one is\n",
"pprint (m.get_fmin())\n",
"\n",
"# The most important one here is is_valid. If this is false, the fit did not converge and the result is useless !\n",
"#\n",
"# When is_valid = False it can be that the fit function is not analytical everywhere in the parameter space \n",
"# or does not have a local minimum (the minimum may be at infinity, the extremum may be a saddle point \n",
"# or maximum). Indicators for this are:\n",
"# is_above_max_edm=True\n",
"# hesse_failed=True\n",
"# has_posdef_covar=False\n",
"# has_made_posdef_covar=True\n",
"#\n",
"# Migrad reached the call limit before the convergence so that has_reached_call_limit=True. \n",
"# The used number of function calls is nfcn, and the call limit can be changed with the keyword argument ncall \n",
"# in the method Minuit.migrad\n",
"#\n",
"# Migrad detects converge by a small edm value, the estimated distance to minimum. This is the difference between\n",
"# the current minimum value of the minimized function and the next prediction based on a local quadratic \n",
"# approximation of the function (something that Migrad computes as part of its algorithm). If the fit did not \n",
"# converge, is_above_max_edm is true.\n",
"#\n",
"# To have a reliable uncertainty determination, you should make sure that:\n",
"# has_covariance = True\n",
"# has_accurate_covar = True\n",
"# has_posdef_covar = True\n",
"# has_made_posdef_covar = False \n",
"# hesse_failed = False\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[Param(number=0, name='a', value=0.9908538664000521, error=0.05876818861364508, is_const=False, is_fixed=False, has_limits=True, has_lower_limit=True, has_upper_limit=False, lower_limit=0.0, upper_limit=None),\n",
" Param(number=1, name='b', value=1.945147699139382, error=0.09907833701567537, is_const=False, is_fixed=False, has_limits=True, has_lower_limit=True, has_upper_limit=True, lower_limit=0.0, upper_limit=10.0)]\n"
]
}
],
"source": [
"# The second one is a list of dict-like objects which contain information about the fitted parameters:\n",
"pprint(m.get_param_states())\n",
"\n",
"# Important fields are:\n",
"# index: parameter index.\n",
"# name: parameter name.\n",
"# value: value of the parameter at the minimum.\n",
"# error: uncertainty estimate for the parameter value."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True\n",
"True\n"
]
}
],
"source": [
"# if you need just the basic\n",
"print (m.migrad_ok())\n",
"print (m.matrix_accurate())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Hesse"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td/>\n",
"<th title=\"Variable name\">\n",
"Name\n",
"</th>\n",
"<th title=\"Value of parameter\">\n",
"Value\n",
"</th>\n",
"<th title=\"Hesse error\">\n",
"Hesse Error\n",
"</th>\n",
"<th title=\"Minos lower error\">\n",
"Minos Error-\n",
"</th>\n",
"<th title=\"Minos upper error\">\n",
"Minos Error+\n",
"</th>\n",
"<th title=\"Lower limit of the parameter\">\n",
"Limit-\n",
"</th>\n",
"<th title=\"Upper limit of the parameter\">\n",
"Limit+\n",
"</th>\n",
"<th title=\"Is the parameter fixed in the fit\">\n",
"Fixed\n",
"</th>\n",
"</tr>\n",
"<tr style=\"background-color:#FFFFFF;\">\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"a\n",
"</td>\n",
"<td>\n",
"0.99\n",
"</td>\n",
"<td>\n",
"0.06\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td>\n",
"1\n",
"</td>\n",
"<td>\n",
"b\n",
"</td>\n",
"<td>\n",
"1.95\n",
"</td>\n",
"<td>\n",
"0.10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"------------------------------------------------------------------------------------------\n",
"| 0 | a | 0.99 | 0.06 | | | 0 | | |\n",
"| 1 | b | 1.95 | 0.10 | | | 0 | 10 | |\n",
"------------------------------------------------------------------------------------------"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# The Hesse algorithm numerically computes the matrix of second derivatives at the function minimum \n",
"# (called the Hessian matrix) and inverts it.\n",
"# Pros:\n",
"# (Comparably) fast computation.\n",
"# Provides covariance matrix for error propagation.\n",
"# Cons:\n",
"# Wrong if function does not look like a hyperparabola around the minimum\n",
"\n",
"m.hesse()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr>\n",
"<td/>\n",
"\n",
"<th>\n",
"a\n",
"</th>\n",
"<th>\n",
"b\n",
"</th>\n",
"</tr>\n",
"<tr>\n",
"<th>\n",
"a\n",
"</th>\n",
"<td>\n",
" 0.003\n",
"</td>\n",
"<td style=\"background-color:rgb(140,140,250)\">\n",
"-0.005\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<th>\n",
"b\n",
"</th>\n",
"<td style=\"background-color:rgb(140,140,250)\">\n",
"-0.005\n",
"</td>\n",
"<td>\n",
" 0.010\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"---------------------\n",
"| | a b |\n",
"---------------------\n",
"| a | 0.003 -0.005 |\n",
"| b | -0.005 0.010 |\n",
"---------------------"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Covariance - colored table\n",
"m.matrix()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"---------------------\n",
"| | a b |\n",
"---------------------\n",
"| a | 0.003 -0.005 |\n",
"| b | -0.005 0.010 |\n",
"---------------------\n",
"[[ 0.00345475 -0.00490941]\n",
" [-0.00490941 0.00981865]]\n"
]
}
],
"source": [
"# to access the matrix\n",
"cov = m.matrix()\n",
"print (cov)\n",
"\n",
"# or as a numpy array\n",
"npcov = m.np_matrix()\n",
"print (npcov)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr>\n",
"<td/>\n",
"\n",
"<th>\n",
"a\n",
"</th>\n",
"<th>\n",
"b\n",
"</th>\n",
"</tr>\n",
"<tr>\n",
"<th>\n",
"a\n",
"</th>\n",
"<td>\n",
" 1.00\n",
"</td>\n",
"<td style=\"background-color:rgb(140,140,250)\">\n",
"-0.84\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<th>\n",
"b\n",
"</th>\n",
"<td style=\"background-color:rgb(140,140,250)\">\n",
"-0.84\n",
"</td>\n",
"<td>\n",
" 1.00\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"-----------------\n",
"| | a b |\n",
"-----------------\n",
"| a | 1.0 -0.8 |\n",
"| b | -0.8 1.0 |\n",
"-----------------"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Correlation - colored table\n",
"m.matrix(correlation=True)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-----------------\n",
"| | a b |\n",
"-----------------\n",
"| a | 1.0 -0.8 |\n",
"| b | -0.8 1.0 |\n",
"-----------------\n",
"[[ 1. -0.84293683]\n",
" [-0.84293683 1. ]]\n"
]
}
],
"source": [
"# to access the matrix\n",
"corr = m.matrix(correlation=True)\n",
"print (corr)\n",
"\n",
"# or as a numpy array\n",
"npcov = m.np_matrix(correlation=True)\n",
"print (npcov)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Minos"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Minos implements the so-called profile likelihood method, where the neighborhood around the function minimum \n",
"# is scanned until the contour is found where the function increase by the value of errordef. \n",
"# Pros:\n",
"# Good for functions which are not very close to a hyper-parabola around the minimum\n",
"# Produces pretty confidence regions for scientific plots\n",
"# Cons:\n",
"# Takes really long time\n",
"# Result is difficult to error-propagate, since it cannot be described by a covariance matrix\n",
"#\n",
"# The results contain information as:\n",
"# At Limit: Whether Minos hit a parameter limit before the finishing the contour.\n",
"# Max FCN: Whether Minos reached the maximum number of allowed calls before finishing the contour.\n",
"# New Min: Whether Minos discovered a deeper local minimum in the neighborhood of the current one."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr>\n",
"<th title=\"Parameter name\">\n",
"a\n",
"</th>\n",
"<td align=\"center\" colspan=\"2\" style=\"background-color:#92CCA6;\">\n",
"Valid\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Lower and upper minos error of the parameter\">\n",
"Error\n",
"</td>\n",
"<td>\n",
"-0.06\n",
"</td>\n",
"<td>\n",
" 0.06\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Validity of lower/upper minos error\">\n",
"Valid\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Did scan hit limit of any parameter?\">\n",
"At Limit\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Did scan hit function call limit?\">\n",
"Max FCN\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"New minimum found when doing scan?\">\n",
"New Min\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"</table>\n",
"\n",
"<table>\n",
"<tr>\n",
"<th title=\"Parameter name\">\n",
"b\n",
"</th>\n",
"<td align=\"center\" colspan=\"2\" style=\"background-color:#92CCA6;\">\n",
"Valid\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Lower and upper minos error of the parameter\">\n",
"Error\n",
"</td>\n",
"<td>\n",
"-0.10\n",
"</td>\n",
"<td>\n",
" 0.10\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Validity of lower/upper minos error\">\n",
"Valid\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"True\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Did scan hit limit of any parameter?\">\n",
"At Limit\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"Did scan hit function call limit?\">\n",
"Max FCN\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"<tr>\n",
"<td title=\"New minimum found when doing scan?\">\n",
"New Min\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"<td style=\"background-color:#92CCA6;\">\n",
"False\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"-------------------------------------------------\n",
"| a | Valid |\n",
"-------------------------------------------------\n",
"| Error | -0.06 | 0.06 |\n",
"| Valid | True | True |\n",
"| At Limit | False | False |\n",
"| Max FCN | False | False |\n",
"| New Min | False | False |\n",
"-------------------------------------------------\n",
"-------------------------------------------------\n",
"| b | Valid |\n",
"-------------------------------------------------\n",
"| Error | -0.10 | 0.10 |\n",
"| Valid | True | True |\n",
"| At Limit | False | False |\n",
"| Max FCN | False | False |\n",
"| New Min | False | False |\n",
"-------------------------------------------------"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.minos()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td/>\n",
"<th title=\"Variable name\">\n",
"Name\n",
"</th>\n",
"<th title=\"Value of parameter\">\n",
"Value\n",
"</th>\n",
"<th title=\"Hesse error\">\n",
"Hesse Error\n",
"</th>\n",
"<th title=\"Minos lower error\">\n",
"Minos Error-\n",
"</th>\n",
"<th title=\"Minos upper error\">\n",
"Minos Error+\n",
"</th>\n",
"<th title=\"Lower limit of the parameter\">\n",
"Limit-\n",
"</th>\n",
"<th title=\"Upper limit of the parameter\">\n",
"Limit+\n",
"</th>\n",
"<th title=\"Is the parameter fixed in the fit\">\n",
"Fixed\n",
"</th>\n",
"</tr>\n",
"<tr style=\"background-color:#FFFFFF;\">\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"a\n",
"</td>\n",
"<td>\n",
" 0.99\n",
"</td>\n",
"<td>\n",
" 0.06\n",
"</td>\n",
"<td>\n",
"-0.06\n",
"</td>\n",
"<td>\n",
" 0.06\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"<tr style=\"background-color:#F4F4F4;\">\n",
"<td>\n",
"1\n",
"</td>\n",
"<td>\n",
"b\n",
"</td>\n",
"<td>\n",
" 1.95\n",
"</td>\n",
"<td>\n",
" 0.10\n",
"</td>\n",
"<td>\n",
"-0.10\n",
"</td>\n",
"<td>\n",
" 0.10\n",
"</td>\n",
"<td>\n",
"0\n",
"</td>\n",
"<td>\n",
"10\n",
"</td>\n",
"<td>\n",
"\n",
"</td>\n",
"</tr>\n",
"</table>\n"
],
"text/plain": [
"------------------------------------------------------------------------------------------\n",
"| | Name | Value | Hesse Err | Minos Err- | Minos Err+ | Limit- | Limit+ | Fixed |\n",
"------------------------------------------------------------------------------------------\n",
"| 0 | a | 0.99 | 0.06 | -0.06 | 0.06 | 0 | | |\n",
"| 1 | b | 1.95 | 0.10 | -0.10 | 0.10 | 0 | 10 | |\n",
"------------------------------------------------------------------------------------------"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.get_param_states()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([0.99085387, 1.9451477 ])\n",
"array([0.05876843, 0.09907874])\n",
"array([[0.05866313, 0.09929038],\n",
" [0.05888831, 0.09888435]])\n",
"array([[ 0.00345475, -0.00490941],\n",
" [-0.00490941, 0.00981865]])\n"
]
}
],
"source": [
"# Other ways to access the fits results\n",
"pprint(m.np_values())\n",
"pprint(m.np_errors())\n",
"pprint(m.np_merrors()) # The layout returned by Minuit.np_merrors() follows the convention \n",
" # [abs(delta_down), delta_up] that is used by matplotlib.pyplot.errorbar.\n",
"pprint(m.np_covariance())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Plot contours"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.contour.ContourSet at 0x116874950>"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3hU1daH3z0z6b1CKCH0HlooCoKCgAVERRRUUBFRUQTFgt6riAVEwY7YsCBNVJoo0nvvhNAhhAQSSO/T1/fHgNcP6TmTkHDe58mTMnvWXpNMzu/stddeS4kIOjo6OjrXL4aydkBHR0dHp2zRhUBHR0fnOkcXAh0dHZ3rHF0IdHR0dK5zdCHQ0dHRuc4xlbUDV0p4eLjExMRc+RPFDJIPznyQIs39AkB5gvI98+EDyhtQ7plLR0dH5wrYtm1bhohEnO+xcicEMTExbN269bLGij0JKfwGLCvBmQd4gakVeHVCeXUCY+WzI8/5fPbrf17Ez72gq/+Nsx8D2y7Etgtsu8GZDthcYzwag0czlEcDMNYCUwzKEHwFr1hHR0en5Cilki70WLkTgstBnHlIwRdQ9BMoE3jd7Lrwe3ZEGcO1n9AYBV43oAARAWeqSxisO8G2C4pmIFj+558KAVNNlygYXZ8x1gRTDZTy0t4/HR0dnYtQoYRAxA5FM5GCT0Fywac3yn84yhhZaj4opcBYBYxVUN63n/HLBo4UsCeCIxGxJ7pWEZa1iHP2P55tQjzboby7gVcXlPG8qzgdHR0dTakQQiAiYF2F5I0DxxHwbIsKeBXl0aisXQNAKY8zK4Caru//8Zg4C8BxDOzHEFsCWJYieW8AoxCPFmdEoSvKVL0sXNfR0bkOUOWtxERcXJz8c49AbAeR/PfAuhaMNVABr7juplX53KQVEbAfAstixLwY7PtdD5gaul6X1w2uPQflWbaO6ujolCuUUttEJO58j5XrFYEU/47kvgTKHxXwGvg+WO4vkEop8KgHHvVQ/s8i9uNgWYKYl0DhRKTwc8AL8WyB8mwNnm3PCIO+t6Cjo3N1lNsVgTizkPTbwFQTFfIlyhBS1q65HXHmgHUrYt0M1s1g34cru8kTPJqDZ2uUd1cwNSy3KyIdHR33UCFXBJL3HkghKujd60IEAFfaqfetKO9bARBnLli3/U8YCichhRPBVBe87wKfu1DGqDL2WkdH51rHbSeLlVLVlVIrlFL7lFIJSqlh5xnTSym1Wym1Uym1VSnV4XJsi2UjmOeC3+MoUx3tnS8nKEMQyrszhsCRGMJnoyI3oAJHgwpACiYg6TfjzOqPFP2COPPL2l0dHZ1rFLeFhpRSUUCUiGxXSgUA24C7RWTvP8b4A4UiIkqpWGCWiDS4mN24uFayeWE4YEeF/4FS3m7xv7wj9uNg/h0pnufKSsLTdZDOu6vrXIV+qE1H57qiTEJDIpIKpJ75Ol8ptQ+oCuz9x5iCfzzFj/9/tPf8ODPAkY8KmayLwEVQpmjwfwb8hoBtN2KeD+ZFiGUJYEQ8W6O8urgykUzVytpdHR2dMqRUNouVUjHAaqCJiOSd89g9wFggErhTRDac5/mDgcEAsY39W+1YdSeGsJnudrvCIeIE+x7EvBQsy1xpqgCmBuB1K8qnB8pUq2yd1NHRcQsXWxG4XQjOhH9WAe+KyOyLjOsIvCEit17MXlwzX9m85mUMgW9q6+h1iNiPgWWZSxhs2wFxHcbzeQC8u5X7VFwdHZ3/UWZZQ0opD+A3YNrFRABARFYrpWorpcJFJOPCI50oY21tHb1OUaYYMD2O8nsccZyG4tlI8Swk9wXIC0F87kX53o86cyJaR0enYuLOrCEFTAb2iciHFxhT58w4lFItAU8g85LG9fCF5ihjJMr/KVT4UlTId+DZBop+RDK648wagFjWlrWLOjo6bsKdK4L2QH8gXim188zPXgOiAUTkS6A3MEApZQOKgQfkcmJVJn1F4C6UMoBXB5RXB8SRDsW/IUUzkeyBiGcHVMDLrpLaOjo6FYbyd7K4mbds2VmIUsayduW6QcQKRVORgkkgeeBzL8p/GOrvfg46OjrXOhfbIyifrSodqWXtwXWFUp4ov4GoiKXg+xgUz0fSu+HMn4DYdruEQkdHp9xSPktMOBJBz30vdZQhCBU4EvF9CCn4EAq/Qgq/AjwQj0bgEYvyiAWPWFclWFU+7zN0dK43yqcQ2BPB66ay9uK6RZmqo4I/QhyvnGnRuRusu87sJ/x0ZlAg4tEU5dkSfB9CGULL1mkdHZ0LUg6FwIjYj+ot4a8BlLEyGCujvLsDIOIA+xHXSWbbbtfngolQ+D34DQa/R/XT4Do61yDlTwiUJziOlrUXOudBKeP/eilwHwBiP4Lkj3eFkopmQMBw8O6lh410dK4hyt9/o/JyhYZ0ygXKVBtDyCRU6FQwhCG5ryCZ9yCWVZS3jDUdnYpK+RQC5ym9rHI5Q3m2QYX9igqaAM48JPsJJPNO1xkFKS5r93R0rmvKoRD4uT6bl5StHzpXjFIGlE9PVMQiVNA4wAPJewM53cmViupIK2sXdXSuS8rfgbK4ONm8MAwMYRjCppe1O5qSYc4jzZwL8Pdm+OW0nDSgiPYLx9dUvvoWiwjYtiKFP4BlKWAE79tQ/s+56iDp6OhoRoVrVal8eiMF4xF7YrkuiGZ3OojPOc769IOszzjIofyrPyhnQFEroBJNgqrTOLg6TYKqEeMfifEa3pRVSrn6LHu2RuzHkaKpUPwLYl6E+A1G+T+pZxnp6JQC5XJFsGXTn0h6J/B7HEPAi2Xt0hVTZLcwJXE1vyRtIN9uxqgMNAuuwY0R9akdUAmFq0PP2b/Npf5CdqedQ/lp7MlNZm9uCnk2V8zdz+hFg6AqxPhFEu0XTg2/cKL9wonyCblmBUIcp5H8cWD+HYzRqMA3UF4dy9otHZ1yT5n2I9CauLg42bp1K87sJ8GWgIpYiVLlY2HjECd/nNjOpENLyLTk07lSE7pFxdImrA7+Htrc+YoIx4sy2ZNznITcZPbnniSpMJ18u/nvMR7KSDXfUGL8I2kXXpeOkQ0J8wrQZH6tEMt6JG+06xS5122ogBf0cJGOTgmokEIg5sVIzrOo4K9Q3reUtVuXZEvmET7e/yeH8lNpGhzN8AZ30DQ4ulTmFhFybIUcL8wkqTCd44UZJBVmcCg/lZPF2SgUzUJqcHOlRtwc2ZgqviGl4telELFC4WSk4AvAAsaarr7LXh1dISVVvvZEdHTKkoopBGJF0juCRxyGkM/L2q0Lklqczfi9v7MmfT9R3sE8W/82bq3c9LI2gd2NiHCk4BQrTyWw4tTev/co6gVE0TqsNi1Da9E8pAYBHj5l66cjDcyLEcsqsG4CrKB8wLPdGVHopPdd1tG5BBVSCACcee9B0RRUxAqUsVIZe/Zv7E4HAzZM5GRRFo/VvoW+NW7Ey+hR1m5dkJSiTFae2sua0/vYk5OMTRwYUNQLjKJFSE1ahdWidWhtfExl18JSpBismxDLarCsAkey6wHP9qjgCXpNIx2dC1BhhUDsSUhGV/B7FkPAc2Xs2b+ZmriGTw8sZFyLh7ilUuOydueKMDtsJOQksz07ke1ZR4nPScbqtONr9OTmSo25vUpz4sJql+mms4iA4xiYFyEFn4MxEhU8CeVRv8x80tG5VqmwQgDgzBoE9n2uVcE11Gw9tTibB9Z+TOvQ2oxv2f+aCAWVBIvDxu6c4yxJ3c3StHgK7GbCvQLoFtWM26KaUT+wSpm+RrHuQnKeASlABb2P8u5WZr7o6FyLVGghEMsqJPsJVNCHKJ8eZejZ/2fEtilsyTrCzx2GE+VzbWy+aoXFYWNd+gH+PLmD9ekHsYuDyt7BdIxsSKdKjWgREoPJUPod5MRxyiUGtt3g9xTK7zGUoWL97nV0rpYyEQKlVHVgClAZcAJfi8gn54x5CHjlzLcFwNMisutidv8lBOJEMrqBIQJD2AwtX8JVk27O486V7zG4ThcG1elS1u64lRxrEWtO72PV6b1syjiExWknwORNh8gG3Fq5KTeE1ytVURCxILmvg3ku4AFenVE+97r6MKtrd39GR8fdlNXJYjswQkS2K6UCgG1KqSUisvcfYxKBTiKSrZS6HfgaaHslkyhlAN+HkPyxiG0vyqORdq/gKkkuygCgeUhM2TpSCgR7+tKzWit6VmtFsd3KpszDrDq1lzXp+1h4cieRXoHcVS2OXtXiqOQT7HZ/lPJCBb+P2B5DiueAeT5iWQSGcMT7LpTPvSiPem73Q0enPFFqoSGl1DzgcxE5b7U4pVQIsEdEql7MzrkrAgBx5rpSSb17YAh6VzOfS0Kh3YKnwYiHoXwcdtMau9PB2vQDzEnezMaMQyjgxoj63BfdlhvC65XafoKIDSyrkeLZYFkJ2MDUBBX4GsrzvDdHOjoVkjLfI1BKxQCrgSYikneBMS8CDURk0HkeGwwMBoiOjm6VlJT0r+c7c1+H4rmo8HkoUy0NvdcpKSeLspmXsoX5J7aRacmnRUgMwxrcQaOg0s39F2cWFC9Ain4AxwnwfQQV8DxKle05CR2d0qBMhUAp5Q+sAt4VkdkXGHML8AXQQUQyL2bvfCsCAHGkIhl3gyEUFfYLyuCvgfc6WmJ3OpibsoVvDy8jy1pIt6hYhtTtXuonmcVZiBSMh6JpYKyBChqD8mxdqj7o6JQ2ZSYEyrU7twBYJCIfXmBMLDAHuF1EDl7K5oWEAEAsG5Hsx8CrCyr4s3KfsllRKbRb+ClxNdMS1+IQJw2DqtI4qBqNgqrROKg61XxDS+VvJ5ZNSN6rZ1YHj2EIHHl5zxOn3mpTp9xRVllDCvgRyBKR4RcYEw0sBwaIyPrLsXsxIQCQwu+Q/PdQ/i+g/J+6Cs91SotT5lxmJW1gd04S+3NPYnHaAAj08KFRUDXiQmvTO7otfm7ss+BaHUwAxykMIRMvPE7kb3ESMYPjBMpU221+6ehoTVkJQQdgDRCPK30U4DUgGkBEvlRKfQv0Bs4G/e0XcvQslxQCESR3BJj/QIV8o5cwLifYnQ6OFpxmb24KCbnJJOSkcLggjVBPf56seys9q7ZyaxqqiP2yqthK8QLEuh7Mf7lWnV7t3eaTjo6WlPlmsZZcSgjAVY9GMu8HRxoq7DeUqXSqfF7cJ+FkcTZVffVaOJfLnpxkPjnwJ7uyk6jpF8Gz9W+nQ0T9Mgn5iXk52Pchtj2uRjrWDa4TzPqBNZ1ywsWEoEIGOpXyQQW7lvmS8wziLCpjj2B7ViIjd1as1prupklwdb5uM5j3WzyMXZyM2D6Fh9Z9xvyUrVgctlLxQcSJFM1wVT717oHyfx5xZqI826MMIYg4SsUPHR13UiGFAECZolHBE8B+EMkbVapzn7vKcoiTVmG1iPIJ5tfjG0vVl/KOUoqbKzXi5w7Deb1Jb1Dwzp7Z3LXqfb46tIQM83mzkbVDchDLcrDtBOUFzlSwHwef+84McEU9xZGGFP+BM/tJpPBH9/qko6MxFVYIANf+gN8QMM9DzEtLZc48WzEZlvy/v19wYjsf7fsDgD7RN7Az61ip+FHRMBmM9KzWimk3DuWL1o/TNDia746s5O7V4/n60FK3rRCUIRRDyDco7+5IZj8kbzTKsy3K4I+IDaU8ELEh+R8itu2uftrmhUjhT27xR0fHHVT4Y6/K/2nEssy1KvCMQxncW+Zg0cmdJBVm8GKjngDcHNmICXt/J9DDh7XpB+hVTT/NWhKUUsSF1SYurDbJhZl8dXgJ3x5Zzl+pO3mp4V3cEKFt+Yizm8jKfwioQKRoGmLbCdIHOLN5XfgV4ET5P+Pqh6D8Ees6+FsorNdUZVwdnXOp0CsCAKU8UUHvgTMbyXN/+YkulZuyI/sYVqcdgH15J2gXUZdQL3/6RLeje1Qzt/twvVDdL4x3mvXl87iBGDAwbNsPvLxjKitP7aXQbtFkDqVMrkw0+3HEsgIVNBpD8HhXTSNlQpw5iPkvlN+jfzfFEctScOa5RMC2H8kfizN7CGK9aD1FHZ0yo0JmDZ0PZ/4nUDixVHocf3FwEUmFGQR5+OJlNBEXWptOlcq+GF5Fxuq081PiaqYmrqHQbsGkjDQPqUH7iPrcEFGPmn6RJc42EtshMNWB4hngdTPKWAUpmonY96MC3kApA+LIRLIeQIX+BLb9iHkuyqMpqGDEvAAV/AnKEKTRq9bRuXyuu/TR8yFiRTLvBWcOKvxPlCHQDd65sDsdbMw4RJo5Bz+TF+0jGhBg8sYhzjKp0389YXPa2Z19nPUZB1mffoAjBacAqOYbyhtN79OkIqxYd4LyRHk0Qop/R+yHMAS8AIAz91VXmROf3kjhZJTnDX/3yXBm3ocK+A/Ks4XLzmWeXdDR0QJdCM4gtnjX+QKfuzEEjdXYs3+zKzuJZiE1/vXzPTnJhHj64WUwEe7tPkHSgVPFOazPOMjUxDWkFufwcqO7uLu6dnWFxLwYKZrmCj/a4pG8t1HhC6FoBmAB77tRpmqIdRtSMBEV8jVKmVwrVFs8mGqhAkbqJSt03E5Z9SO45lAeTRG/QVD4FeLVFeXd2a3zbcs6Sg2/cII9/QBYmhbPr0kbCfXyx6SMnCzO4o2m9xHtF+5WP0pKkdmKxe7AqBTKoDAqhdFowGQ0YDRc2xewSj7B3FO9DV0qN+W/u2YyJmEOh/JTGVr/dryNGjSq8eoCtn1I1gDwvh0V+CbK4I/TmYHybIkyuSqsimUFyqsDOJJwFs8D215U4Juu8hZFU8FvQMl90dG5Sq6rFQGcDRHdB44kVMhkt9akz7EW4mfyIstSwKRDSzhZnE3v6m1pG14XX5MnM46tY29uCuNaPOQ2H64Gh9PJxn3H2bz/OFsPJrM/+TQXepsoBR5GIx4mI9UjgmkSU5mmNSvTtGYU0ZEhGAzXTuE/u9PBxIOLmHZsLX4mL7pVjuXOqq1oGly95PsHzkKUwe/v752ZfVGB/3HdfFhWI+bFKK8uIIWIZTXKfyjKVB0p/BZxpGIIfL2kL09H56LoK4J/oJQnhHyHZD2EZD8BoT+iPGLdMtfZlcCmzMOEevkzvMGdBHv6/v14FZ8QzA6rW+a+GkSElbuOMHH+Oo6mZuFhMtK0ZmUG3d6WEH8fHE5BRHA4BYfTid1x9sOBxebgaGomC7fs59c1uwHw9/YkpnIoMZVCXJ8rh9I0pjIRwWVTItxkMDKswR3cFNmAeSnb+PPkTuakbCHaN5weVVvSo2rLqw7V/VMExFkIxmpg24sYqyL5E1B+j4AxAjFvQXm2comAswAw6sXrdMqc625FcBZxpCFZD4IzHxU6BeXRUAPv/s1pcy5v7JrF8w3vpH5glb9/Pi1xLfNTtjCiYU/ahNdxy9xXwraDKXw6dy3xianEVArhqZ43cFPTWvh4Xln4xOF0ciwti/hjaexLOk3SqSyOncrmdE7B32Nia0XRuXkdOjevQ7UI97evvBAFdjPL0/bwx4nt7Mg+hp/Ri2EN7qBXtTgNMowOIHmvgbEGGKtjCHgese5CCr9ABY1DGYIR2x6keB7KIxbl01OjV6Wjc370zeILIPZkJOshEAsqbBrK5J4L8lObvqFbVCx3Vm3J0rR4ZiVtINTLn2frdad2QGW3zHklfDBrJTNW7CAy2J8n72xHzxsaYzJqG/svNFtJTMti474klu84zP7k0wDEVAqhVb1qxNWrTqu61QgP8ruEJfeQVJDOe3vnsS3rKG3D6vCfJvdSWYMey+IsAuWDUgop+AxxpGEIehdxZIBlMWI/jPJ/HmUI0OBV6OhcGF0ILoLYE11igAEVOg1l+neWT0nZnpXInOTNFNot5FgLeapuV2oHVCIhN4XT5lyOF2ZQL7AKPaq21HzuS7FpXxJPfzqbezs05cU+N+PtWTrRwpOZuazYeYRN+4+z4/AJCs2uEFmD6pHc0aYB3VvXJyKodENITnEyO3kznx34CwOKe6Lb0K1yLPUDq2hS8dSVOfQ5KnA0UvgdKBPKqxPK66b/1+9AR8cd6EJwCcR2EMl62HXnFjYdZayqqf2znCzKpopvCNuzjrLo5C58TF6EeroudmvT93NHlRaapjZeimKrjfvfmoLJaGTmfx/Gy6NstozsDicHkk+z5UAyy3YcIiHpFAalaNswmq4t69GxaS1CA30vbUgjThZl8/GBP1hzej8OcVLdN4xuUbF0jYqlln+lEtmWgi8Qy3ow1Ub5P4syRmjktY7OxdGF4DIQ214kqz8YglGhM1DGSM3nAEgsOM2spA00DKpKw6CqVPcNw9voydrT+1mWtodRsfdd2ohG/LpmN2OmL+Pr5+8jrl71Upv3UiSmZfHnpn38tWU/JzLzUAqa1apC5+Z1uL1NA8ICSyd8lGMtYuWpBJak7WZb5lGcCE2CqvNOs74l6rPs+p+zoZTneVcC4sxG8t5yhYyugV4aOhUDXQguE7HuQrIHgKkuKnQqSnlraj/fVszo+F9pGFiVPjVuINDDB4Aj+Wm8lzCPu6rF0bNaK03nvBjvTl/Kkm0HWTH+6WsyLCEiHDqRwfKdh1m56wgHU9IxGQzc3Lw2vTs0pXX96FJLT82w5LMsNZ6vDy/Fw2Dk/Rb9iQ1xz0VaLOuRnGdAHCj/oeD3KK723zo6V48uBFeAmBchOUPBuycqaLymF8hcaxFDt37HlBufBVwZRcvT9pCQm0KMXwSP13HvAbdzGTj+Z5RSTB5xf6nOe7UkpmUxZ208v2/cS26hmbpVwxl+b0duaKT9vs6FOFZwmhe2T+G0OY/XGt/D7VWau0VExZGK5L0FlmVgqo8KfAflqRcs1Ll6yqpncXVgClAZV/eOr0Xkk3PGNAC+B1oC/xGR8Zey624hAJCCSUjBR66luf/Tmtp+Y9csVwaJCGanlVr+lagfEEX9wKpU9gniYF4qNfwi8DG5v2zxvaN/pGblUCY8Wb5SFy02O0u2HeSrBRs4kZnHDY1qcH/HZjSpWblUwkY51iJe2TGVHdnHqOUfSZ/oG7i9SnN8TV6azyXmxS5BcGaiAl4E34HX5OpN59qnrIQgCogSke1KqQBgG3C3iOz9x5hIoAZwN5B9zQiBCJI7AswLUMGfo7y7aWbb7LCxK/sYu7KTuCmyIU5xkmHJ51hhOqtP78NDGfEwmBjboh/+Jm1DU+fy/KR5pKTn8ssb5bO8gdVmZ9aqXXy7cBN5Ra6y01GhgTSJqURcvercdWNjt22A25x2Fp3cxazjG9ifdxJ/kzc9qraiT3Q7qvuFaTqXOPOR3FfBshi8uqKC3tPTTXWumGsiNKSUmgd8LiJLzvPYm0DBtSIEACJmV1qp/Yhr89hNB85+SdrAwfxUIrwC6RDZgEZB1XgvYS6BHj4MqdfdLXOe5ZM5a5i2bDvrPxmq+bmB0sRstbM/+RR7EtPYc8z1cTIzj6jQQJ7t1Z7ucfXdtpcgIsTnHGfW8Q0sS9uDQvFU3a48XLMDBg0LyYkIFH2P5H8Axmqo4M9QHg00s69T8SnzEhNKqRigBbCpNObTAqW8IXgSktkbyX4Kwn5DGbUrDmd3Ohgd/ysZlnyeqtv1/1Up9TJ44Gfydntuec3KodgdTpLTc6hZOdRt87gbb08TzWtXpXnt/6X9btp/nI9/W81/vl/Ij0u28lDnFnSPq4+nxisEpRSxITWIDanBsPp5TNj3O58f/IuNGQd5M7YPkd7a9B5QSoHfQPCIRXKGI5l9EL/+KJ8H/y5sp6Nztbj9NlAp5Q/8BgwXkavqNK6UGqyU2qqU2pqenq6tgxeb1xiJCpnk6m6W8ywiZs1sZ1sLsThtTGoz6G8ROF6YwbyUrWRY8ty2CflPzl44V+064tZ5yoK2DaKZ9upDvP3obTgcTkZNWcyd/53Mt39uotjqnv7GEd6BjG3+IP9pci97cpN5cN2nzDy2jlxrkWZzKM84VNhc8O4Chd8jGbfizH4GsWykvCV+6Fw7uDU0pFw5bwuARSLy4UXGvck1Fhr6J2JeiOQMB49WqJAvNWtqc/uKsYxr/hCCkGHJ43B+GtnWQjpGNuTGiPqlctr08QmzyMorYvabj1TYTUgRYfP+40xbvoO1exKJCg3k5QduplOs+4q9JRWk886e2ezKScLTYOLmyEbcVT2OuNBamoWMxJGKFE2Hop9BcsBUD+XbH3zuQikfTebQqTiU1WaxAn4EskRk+CXGvsk1LAQAUvwHkvuyq5FIyGRNDpytPLWXP05sp8hhoX5gFYI9/OgU2ZAa/qV32nT+hgTenLKY7168//+FVioq2w+lMHbGco6kZtK8dhU6xtaiQ+Oa1K4S5hYhPJiXyvyUrfyVupM8WzFVfELoF9Oe+6LbYdRKEMQMxQuQoh/BfgAMEajAt1DeXTSxr1MxKCsh6ACsAeJxpY8CvAZEA4jIl0qpysBWIPDMmAKg0cVCSGUlBABiWec66GMIRYV8hzLFlNhmgd2Mr9GTPJv5/5WoLi2KzFa6jfyaLi3rMnqAezenrxVsDgc/r9jJgk37OJjiCjVWCvGnd4dYBt7Wxi0byxaHjZWn9vJb8iZ2Zh+jYWBVXmtyz/+rSFtSRASsm5H8d1yC4H2PqyeCG9uy6pQfromsIa0oSyEAENtuJOsJUAZUyLcoj8Zl5otWjJm+jPkbEvhr7BME+19fIYXTOQWsTzjGsh2HWJdwjJua1uTdx27H30f7MwHgulgvTYtnwr4F5NqK6FvjRgbXuVXTcyMiVqRgIhR+7VodBI1xdUfTua65mBCU35zBMkJ5xKLCZgBeSNbDiGWj2+Y6UZTFmD1zsDsdbpsD4P5OzbDaHcxdv8et81yLRAb7c3f7Jnz6zN2M7NuZDQlJ9Hl7CpN+X09Keo7m8yml6BoVy6ybnqdn1VZMO7aWvms/Zn36AQ3n8MQQ8Dwq7GdQfkj2QJy5b5xphKOj8290IbgKlKkWKmwmGKsg2Y8j5kVumWdXdhJzU7Ywft/vbs0IqVM1nLh61fhl1W4cTueln1ABUUpxf6dmfDm8NzUrh/Htwk3c9cb3DPpwFlsOJGs+X6CHD681uYev2jyBt9GD4dt+5L87Z5JpyddsDuURiwqfC76PQ/HPSOZdiHWzZvZ1Kg66EFwlylgZFToNPJoiOcOQolmaz3FH1amF3VgAACAASURBVBYMqNmR2cmb+fW4+1YeAH1vaUFqVh6TF17fF4qWdavxxXP38ue7g3i2V3vSsvJ58uNfGTtjGUVm7duKtgityU/thzK4ThdWnErg/jUfsSwtXjP7SnlhCHzF9V7FgGQ9jDPvXUSKNZtDp/yj7xGUEJFiJHsoWFejAseifHtrat8hTl7ePpX1GQf5uNUjtA2vq6n9s4gIr//wFwu37Ofjp3txU9NabpmnvFFstTFp/nqmLd9OiL8v3ePq0S2uPrE1ozTPMkoqSGd0/K/syU2mX432DK1/GyaDUTP74ixCCsZD0VQwxqCCx7utX7fOtYe+WexmRCxI9tNgXYcKmoDy6aGp/UK7hSc2fsUpcw7f3TCEGn7anXD+J8VWG4+Pn0VKRi5TR/YjOvLqa+5XNHYdPclPS7axdk8iVruDqNBAnrizLXff2ETTeWxOO5/sX8is4xuIDY7mrdgHStT74HyIZYOrdpEzAxU0Vu+XfJ2gC0Ep4FoZDAbrVlTwxyhvbVMxTxZl8+iGiQR4+PD9DUP+7mWgNSczc3l47HTCAv34aeSDpda6srxQUGxh5e4jzF4Tz84jJ3moS0uG33sTRoO2UdbFqbt4d88crE47nSIb8UCNG2geEqPZKkScWUj2c2DbDH6DUf4voDSsjaRz7aELQSkhzkIkeyDY9qBCJqK8btbU/s7sYwzZPJk6AZX5uNUjhHq5p6fv+r3HePazOQzo2orh93Z0yxzlHbvDyUe/rWLGip00qB5Jv1ua0y2uvqbVTk8V5zDr+EbmpWwhz1ZMvYAoHqhxI92rNMPTUPJ5RKxI3ttQ/DN4dXb13zCUbp9ondJDF4JSRJz5SNYjYD+ICvkK5dVeU/trT+/n1Z0ziPQO5LO4gZqHDc7y9rQlzFuXwPcvPUDTmlFumaMi8MemfXz312YS07II9vfhgU7NGHhbGzxM2sX2zQ4rf53cxc9J6zlScIqafhGMiu1Do6CSF5tzVTWdiuSPAUMYyu8x8HlAF4QKiC4EpYw4s5GsAWBPQoVORnlq25B+d3YSz2/7ES+jB5/GPUadgMqa2gdXCOT+t3/C19uD6a8+pHnVzoqEiLD5QDIzV+xg1e6jNKgeyZiBtxOjcUVXEWFd+gHe2zuPTEs+A2p2ZFCdznhosTqwbkMKPgHrRlCB4PsQyncAyqhtbwWdskMXgjJAHBlI1sPgPIUK+Qbled7f/1VzJD+N57b+gNlhZUKrATQPidHUPvwvRNSzXSNG9e9Wav2BL4SIsOWvnaz+ZQN7Nx7EYbOjDAaXX0oRUimIBm3qUr91bao3qErVOpXx9HZ/p7d/snLXEUb/tBiLzU6Pto24u30TGtWopOkc+bZiPtz3B3+c3E7dgChGNb2PeoHarNrEugsp/AYsSwBP8O2N8n0cZaquiX2dskMXgjJCHGmulYHjBCrobZTPvZraTy3OZuiW7zllzmFM837cFKl985xv/tzIpN830KdjLCP7di6zCqUHtx3hyxE/Er96H/7BfjTt2BDfAB8cDieI4HA4SU/O5PD2o9htrpPYBoOicYcGdLinLa26NSO6QdVS8f90TgGfz13Lku0Hsdgc1K8WwUO3tuTONg01nX/16X2M2TOHPFsxj9e5hUdqdtIs3VTsR5HCb6F4HqBQASNdq4QKWqH2ekAXgjJEnNlIzjDXkttvEMp/BEppFz/OthYwfOuPHMxP5d1mfelcWdt0RhHh83nr+H7RFh7q0pIXencs9YvBicOpPNnsRXwCfBgwqg+3D+qC6QKhKqvZStLeFJIPnCRxdxIbF2zjWILrZHBIpSBu7NWGLg/dROP29TFonOlzLvlFZhZuOcDstfEcTEknrl51/vNgF2pU0m5fJ8daxPi981mctpuGgVV5M7YPNf1LXhn3LOJIQ3L/C9bV4NXFVbfIoKcVl0d0IShjRGxI/rtQNB28bkYFfajpZlyh3cKwrd+TkJvC+y0e5qZIbVsYigjjf1nJjBU7eeKOtjzd80ZN7V+K0feNZ/vS3UxO+Ijwqlcesz55JI1dKxPYvnQ3G3/fhrnIQmjlYG64qzXdH7uFBm3qXJa4/fntMqJqRdKic9Mrmt/pFOasi+eTOWux2uwMur0tj3SL03RDeWlaPOMS5mFz2hkdez+dKjXSzLaIE4p+RPLHuyrvBo1HebXVzL5O6aALwTWCFE13peuZaqKCv0SZojWzXWAzM2TLZI7kp/FhqwGan0AWEd6eupS56/fw1iPd6dFOuwvNpehTeRDt7mzJiMlDSmyruKCYDfO3sn7+Fjb9sR1zoYX6rWvzwbJR+Fyk8mpxQTHzv1jM9DG/0ePJbjwx7uErnjs9t4Dxv6xiybaD1IoK5fWHu9KslnZlqNPNeby0Yyp7c1N4um43Hq3VSdPVm9gSkJwXwHEM/J5G+T+LUnoSQXlBrz56jaB8H0SFfA+ODCTzPk0rl/p7ePNp3GPU8Ivgxe1T2Z51VDPb4CrK9uqDnWldvzpvT1vKriMnNbV/McwFZnwDtenV4OPvQ+cHb+K/M19g5omveW7iIOq3rnNRETj7vAde7kXDdvWoXt918b7Sm6iIIH/GDbqTT4b0oshiY9CEWfy0dJtmBQUjvAP5ss0TdIuKZdKhxby0YypH8tM0sQ2gPBqjwmaDz71Q+AWS9RBiT9HMvk7Zoa8IygCxH0eynwRHEirwdZRvP81sZ1kKeGrzN5w25/Jp3EBiQ7RbdQDkFprpP246RWYbP43sR1So+5uePNniRcKqhDLmj9fcPtfF+Ou75SyespIPV771/35+NS1F84stjJ6ymOU7D3NT05oMvqMdjWO0SQMWEX5KXMPkI8spdljpENGAAbU6appZJsULkLw3AIUKfBvlc4dmtnXcg74iuMZQpmhU2Czwao/kjcKZ944rDqsBoV7+TGz9OKFe/gzf9oPmK4MgP28+froXVpudZz6bQ1qWdmWTL0SjdvXYvTKBwjztmsBfKfnZBSycvIw+I+4CwOH4X4+IsyKQtPfyy1UH+HjxweAevNC7I9sPnaD/uBk8+sFMFm05gNNZspszpRQDanVkXqeXGVynC/E5xxm86Wue2PgVu7KTSmT77zl8eqDC5oKpFpI7HGfua3q/g3KMLgRlhDIEoIK/BN9HoWgKkvsyIjZNbEd4B/JF60GEewXwzJbvmJOsbWnpWlFhfDykFxk5BQwc/zOJaVma2j+Xbo/egqXYyooZ69w6z8X44+ulBIYFcENP1w2V0Wj8O6STtC+Fn976hc+enczAhsNYN/fyft9KKR6+tRULxwzi5ftvJqegmFe/+5PBH/+iSVOcYE9fBtXpwvxOLzOiYQ/SzDk8uelrvj60VJNmR8oUjQqdDn5PQvFvSMadiGV1ie3qlD66EJQhShlRAa+i/F8A83wk5xlXI3INqOwTzOR2T9EmrA5jE+Yyfu/vmnY6a1m3Gt+80Aebw8Hj438mPjFVM9vn0qBNHWrF1uD3Lxe5tUHPudhtdv76bjkZJ7PY8PtWHvyPq8S4w+H4OxyUl5XPz+/PxWF3MOq3F3nqw0fZvnT3Fc3j7+NF31taMHvUo4zq342Dyek88O5Ufl29S5PX62Py5IEaNzKzw3C6V2nOt0eWM2TLZNKKSy42SnlgCBiBCv0ZlC+SPQhn7kjEmVti2zqlh9uEQClVXSm1Qim1TymVoJQadp4xSin1qVLqsFJqt1Kqpbv8uVZRSqH8n0IFjgbLKiRrIOLUJtwS4OHDh60G8GBMB2Yd38DwbT+SZ9OuIUn96pF89+ID+Pt48eTHv7JpnzZhh3NRSnHXkO4c3ZV02XfbWmAptrJ18U4erTuUvIw8GrZ1ZWIZjUacZzq5zZ+4CL9AX+4Y1IWAEH+iG1Zl9+q9f4exctJzSU/JvKz5DAZFrxsbM+v1ATSrFcWYGct55rPZmoXf/ExejI7tw+jYPhzMO8lD6z5leZo27UmVZ3NXNzS/p6B4HpJxB2JeqoltnVJARNzyAUQBLc98HQAcBBqdM+YOYCGggHbApkvZbdWqlVRUnEULxJHaUBzpd4nTnq6p7fnJW+WGv/4r964aL8fyT2tqOz2nQO5/e4rcMPRT2XogWVPbZ7EUW+TpuJflrqD+knzghFvmuBAJGw7I03Evy6wP5klq4qm/f37qeLq8cfc42blij9jtdhERmfTCDzLp+e9FRGTjH9vk5a6jZWCjYfJOv4+kuNB82XM6nU75ZdVOuXHYZ3LT8xNl7Z5ELV+SHC/IkAHrPpfWC1+VMXvmSLHdopltpzVBHOk9xZFaVxzZw8TpyNTMts7VA2yVC1xX3bYiEJFUEdl+5ut8YB9Q9ZxhvYApZ/zcCAQrpa7bUpfK505UyJfgOIZkPYg4Tmhmu2e1VnzR5nEK7GYe2ziJrZlHNLMdHuTHl8N6ExUWyHNfzGXXUe1TSz29PRn164t4eJp4p+9H2G12zee4EI3a1eOLLeO46b527FyRwPLpawBIT84kvGooEdXDMBqNnDicyvH9J7ilXwdW/7qBVbPW0/2xzkxO+BiDQbF3w8HLnlMpxX0dm/Hzfx6mSlggwybOZcqSrZr1lK7uF8a37Z6kf82OzEnezCMbvuBAnjZ/N+XRCBX2G8p/OJiXIOm3u7KMylmG4vVEqewRKKVigBbApnMeqgr8M9UihX+LBUqpwUqprUqprenp6e5y85pAeXV0nTVwZiGZfRH7Yc1sNw+J4YcbhhDpHcjwbT+y8tRezWyHBPjy5fDehAf6MfSzOexNOqWZ7bNUqhHB8K+e5MjOY8x8b67m9i9F5ZhIbnvsFpp0cJ3c9g30IeVQKlVqu9I+p7w5i6YdGmKz2Di07SgN2tSlc78OAGScyKIw1xUustvsnDh8eXsq1SKC+W7E/dzcrDYfz17DvW/+yK9rdmPRQAg9DCaG1r+NT+MeI89WzKMbvmDSwcVYnSW3rZQHyn/ImcyiakjuC0jOEMRxusS2dbTH7UKglPIHfgOGi0jeuQ+f5yn/um0Qka9FJE5E4iIiItzh5jWF8myJCp0KOJHMBxGLdtkyUT4hfNVmMPUCohi5YxoLUrZpZjsiyJ8vh/cmwNebpz/5jY1u2DPocE9bbunXnilvzmL5jLWa278cIqNd70FPbw+sxVY+feZbfn5/Hkd3J9Hz6W6kJp7G6RRa3uoqRbF/8yFqNommev0qFOYVMebBj/li+PcMqPMsWxfvuuR8vt6evP9EDz54ogcBvl6Mmb6Mnv+dzIqd2twktAuvy8z2w7gtqjnfH11J//Wfsyfn8lNhL4byqIsK/RkV8ApY1iIZPRBL2fzddC6MW4VAKeWBSwSmicjs8wxJAf5Z37YaUHpHVq9hlEcDVOgMMEYg2QNx5n+KiDZZP0GevnzeeiCtw2rz1p7fmJao3T9mVGgg37zQh0oh/gz9fA4zV+zQPCTwwjdP07RjQ8YN+Iz187doavtKqFonitFzX8ZoNBAaFcxr04fjH+xH8v4T+Ph7U61eFZxOJ3vXH8Q/xA9lUPz4xs94+3nz7oLXeGLcwyyZsvKy5jIYFF1a1uWnV/rx5bDehAf5M+Kr3xk7Yzlma8nv4IM8fRkVex8ft3qEQruFQRu/5JP9f2J2WEtsWykTyu9xVPg8MEYi2Y8jBZM0OzujU3LcmTWkgMnAPhH58ALD5gMDzmQPtQNyRcR9eYjlDFee9i/g3QsKP3f9AzkuLwPlUviavJjQagCdKzXhkwN/MungYs0u2FXCAvn+pb50aFKT92et5J3pS7HZtUtd9fb14u35I6nXqhZjH/yEwzsTNbN9pQSGBvDMpwPp2r8TNZu4TnEfS0imcfv6AKybu4XTx9OJaVyd9ORMju8/wZCPHwPAw8sDHz9vbFbbZf/ulVK0aRDNDy89QP9bW/HL6l0MGDedIyczNHk9N0bUZ2aHYfSq3pppx9by0LrP2J6lze9XmWqhQmeBdw+k4CNXqMh5bpBApyxw54qgPdAf6KyU2nnm4w6l1FNKqafOjPkTOAocBr4BSl5VrIKhDL6ooHGowDFg3YZk9kKs2twFexpMvNu8L72qxfH90ZW8lzAXmwbxYQA/b08mPHkXj3VvzZy1exjy6Wzyi7Q5IwHgG+DD6LkvExDqzxt3jSPl0LVx/yAitL6tBZOe/4Hfv1zMNy//RP02dYnr3pxVP6+jXY9W+Af7YbXYsFlsKINCnFdeosLTw8TzvTvy2bP3kJVfzMPvTeeXVbtKfCoZwN/kzauN72Zi68dxiJOnNn/Dh/sWaPLecL2fx6MCXgfLaiTzXsS2v8R2dUqGXmuoHCG2fa7eBo5k1yE0v8dRquRaLiJMOrSYH46uoklQdcY070dln2ANPHbx5+Z9vDllMbWiwvj82XsID/LTzPbhnYmM7PY2ymBgzJ+vUbdlLc1sl4QVM9eRcuAk9VrXpu0dLTm0/Sjv9P2IHw58ilKKU0np/DT6Fxq2q8udg7teVb2is2TkFvLGj4vYuC+JuHrVeePhW6kWoc3fr9hu5fODf/HL8Y00DY5mbPN+RHoHaWJbrNuRnOfAmXemcVMvTezqnB+9DHUFQpwFSO5rYPkLvG5xrRYM2vzTL02L59342RgNBt6KvZ8bI+prYhdgw94kRnw1n/AgPyY915uq4dpcTACSD5xgZPd3KMgp5P2lo6gfV1sz21qxdfEu1s/dzHNfPEH26Vw2zN/KjmW7GfnTcxg16EsgIsxZt4ePfluNw+nk2V7teeDm5hg1ar6zJHU37+yZjY/Rg3eb9aNVmDaCK450JOd5sG0G34dRASNRqnTbi14v6EJQwRARKPoJyR8HhkhU8Mcoz2aa2E4qzODVndM5nJ/Go7VuZnCdLpq1P9x9NJXnJs7By8PExKH3UqdquCZ2AU4nZzCi0xsU5hbxwfI3qd0sRjPbWpCVls1LXUbT8tZYcjPyCIkMpv09bYjt2KhEq4FzOZWdz7vTl7F2TyLNakUxqn83YiqHamL7aMEpXtkxjZSiLJ6t150HYzpo4reIHcmfAEWTwaOF6/1svG6PE7kNXQgqKGLd5QoVOTNQQWNRPj01sWt22Jiw73fmpWylVWhN3m7Wl3CvAE1sHz6RwZDPZlNktvL6w13pHqfdqiM18RQjOo2iMLeIV6cNo12PVprZ1oK8rHzmfrqQoIhAbu3fET+Neiyci4jw5+b9fDBrBXaHk1H9u9G1VT1NbBfYzbwT/xvLTyXQuVITXm/aGz+Tlya2xfwXkvsqYEQFvqWXttYYXQgqMOLMRrKfBdsW8HsW5T9Us7vLBSe2My5hHn4mL8Y270eL0Jqa2D2Vnc/Ib/+gQXQlXnngFk1sniU9JZNRd4/j8I5jDP6gP/e9oI04ugMtVwLn43ROAS99vYD4xFT639qKoXd3wGTUZk9p6rE1TDywiGi/cN5p1pd6gdrcwYs9Ccl9EWy7wLsXKnCUpm1dr2d0IajgiFiR3NfBPAe870QFvYdS2tylHclPY+TO6Zwoyua1JvfQo6o2dQFtDgcImvbtPYu5yML7j37Oml830nfkPQx8t59bL7haUpBTyI5l8XS4t60mPtvsDib8uopZq3bRqm413ht0B2GB2mzWb808wuu7Z5FnLWJIvW70i2mPQZPkBTtS8AUUTgJjFCrofZTnea9fOleA3pimgqOUp+vi7z8CzH8gWf0RhzZ55bUDKjO53dO0CI3hrfhf+eLgIpwaHATyMBovKQJOp7B+7zHmrI2/Itvevl78Z8Zw7nziVma+N4dxAz6juFC71FV3MvezhbzVZwIfPDYRc5GlxPY8TEZG9u3MW492J+FYGv3GTNPsxHdcWG2mtx/KjRH1+eTAQp7Z8h2nNCltbcIQ8Jyr1wEKyXoYZ/6HmvXr0Pk3+oqggiHmRUjOS2AIQ4V8hfLQJjZsdzr4YO985qRsoXOlxrwZ2wdvo3uzOwrNVnYcPsGUJVvxMBl597HbCb5Eb+F/IiJMHzObH9/4meiGVXn9lxHUaFjNjR6XHIfdwbR3fmPq279SMzaaN34ZQdU62oRdDqak8+rkP0lMy+LBzi0YencHvDxK3nxeRJh/Yhsf7luASRkY2fhuukbFauDxmSy5/DFQ/CuYmqCCx6NM10aKcHlDDw1dZ4gtHsl+GqTQlYHh1UkbuyLMSFrHJ/sX0iCwChNa9ifc2/09iwE+m7uWID9vBnS98hDB9mXxjH3oE6zFVv778/O0vq2FGzzUls0Ld/Dew59gs9p55pOBdH/sFk1CRWarnU/mrOHnlTupUyWMdx67nXrVtKnfdbwwg1G7Z5GQm8LtVVrwUqOe+Ju8NbEt5kWu8CcWVMCr4PNAuQn3XSuUODSklPJWSr2glJqtlPpNKfW8Ukqbv7CO5iiPpqiwX8EYjWQ/iRT+qEn5CKUUD8Z0YHzLhzlWmM6jG75gb26KBh7/G7vDFX6ynekNXFBsIfMqexa37NKUiVveo3KtSP7bYyy/fXTtl0Ruc3sLvtw5ngZt6jBh0CTefuBDTcJb3p4mXnngFj575m6yC4rpP24Gs1ZduvDd5RDtF843bZ9kUO3OLDq5k4fWfcbBPG1OfCvv7qjw38GjJZL3BpIzFHGWXQ/risZlrQiUUrOAfGDqmR/1A0JEpI8bfTsv+org8hFnoSsDw7LMdfgscAzKGKaJ7YN5qYzYPoUMSz5P1rmV/rU6YtRgo/BcUrPyWBufyNqERB7t1poWdVxVyq8m46a4oJhxAz5j3dwtxHVvxojJQwivok2OvbtwOp38Mv53vnttGrWaxfD2/FcIr6rN3zA7v4hRUxazdk8iD9zcnBH3ddIkqwhgd/ZxXts5nQK7mbeb9eWmyAaa2BVxQtH3SP4HYGqACvkSZaysie2KTolDQ0qpXSLS7FI/Kw10IbgyXP84U5H898EQ4BIDb21SNnOtRYzbO4+lafE0D4nhzaZ9qOIbUiKbxRYbq3cfJel0NvuTT5OWlU+dKmHc06EpzWpVwWD4/xf/wycyruhgmoiw4MvFfPXiFLx8vXht+jBadS31t/EVs+nP7Yzp9zFevp68Nn04zW9pooldh9PJp3PW8tPSbdzYKIb3Bt2Bv482GWfp5jxGbP+JA3knGdbgdvrVaK9ZOEcsK10nkpWfSww8tPl9VGS0EIIfgC/PdBFDKdUWeERESr1InC4EV4fYDiK5I8B+AHweRAW+glKXv/F6QbsiLDy5kw/2zgcFLze6i9uiml/VP7zTKXz1xwa+XbiJR7u1pvdNsXgYDUQEu/LIHU7n3yUT1iUcY31CIodOZFJksfJm/25XJAjH95/g7fsncHxvCo+81Ze+I+/GoFE5BndxLCGZt/pM4MTBkzzyVl/6vXqPZhfW2WvjeW/GcmpUCuHjIb00KwFSbLfyZvwvrDiVwD3VWvNSo7s0O6kutgNI9pMg2aigD1De3TSxW1G5aiFQSsXjahTjAdQHjp/5vgawV0RKXYZ1Ibh6RKxI/keuo/zGmqjgCZrdSZ0symZU/Cx2ZSfRtXJTXml8N4EeVyc03/21mbV7Ennyzna0bVjjjO+CUwSjwcCa+KMs3LyfZrWrcH+nZkxbvoPs/CKG3t3hiuYpLjTz4RNfsnLmOpp0aMCIb5+mWr0qV+VzaVFcUMyHg79i5cx13Dm4K0MnPo7RqM2FddP+47z89QKUgrcfvY2bmmqTneMUJ5MOLuHHxFW0CavD2Ob9CLjK98a5iCMDyXkabLtQ/i+C3xP6JvIFKIkQ1LiYYRHRvgXVJdCFoOSIZQOS+zI4M1H+z5355yn5xcQhTqYcXc3Xh5cS5hXA27H3X/Vp5C0Hkvn6j4080+tGmtf+X/fSzLxC3pq6hG6t6tOhSU2C/LxZHX+UOWvj+WBwzyuOcYsIS6asYtLzP2A1Wxn47oPcM+yOa3p1ICJ8958ZzHxvDjfd146Xf3gWb19twjnHT+fwyjcLOJCSzmPdW/N0zxs12zf4PWUbYxPmUtU3hPdbPExN/0hN7IqYkdyRYP4TfHqjAkfrhevOg54+qvMvxJmD5I0C80LwaOVaWpu0ybHfm5vC67t+5mRxNkPr33bVsWERYX/yaXYfTaVzizpEBPkzbdl2svKL6HdLi7/LWQ/6cBY92zWm142N2Zt0igPJp/H18qR768uvY5SZms0nT3/NhvlbaX1bc176/hlCKmlXitsd/PbRAr56cQp1WtZk9JyXiaimzSay2Wpn/C8rmb02npZ1qzL28TuICNKmzMP2rKP/1959R0dVdQ0c/u3JpPdCDb333rsdFMSOKKLwIioq9gKvHXt7BRWsiAiioogUEQtVkN67tECAFNJ7puzvjwnKp5QAd5JAzrMWa0Hm5pzDJcyee8rePLn+SwrcTkY1vYbeVVtZ0q6qG81+F3LeB98OSOS7iO3c1qsuNCYQGCekqpA/E818ARBP4rqAKyxpO9uRz/Obv2Vx0jYurdycp5pdR9BZJifbfegoAX52qlWI4MM5fxAeEsj13Zrja/fhk3krSU7P5vbL27N53xE+/WkVl7auz/w1O7m0TX3u6del2P2oKnM+/IUPHp6EX4Af/e/tzQ2P9CMkwrr6CVZbOXctL98ylsDQAN5e/AJV61q3g2buyu289OWvhAT48d7911l23iApP4P/bviKjelxXFe9Aw81ugp/H19L2ta8WZ7EdT6VkIhxZhH5OCbFhHFCIoIEXotEzwJ7bTT9ftxZb6B67pWoQnwDeL31rdzb4AoWJGxhyB/j2Z+ddFZt1YuN+avQSkJaFtGhQfjafViwYTcrtsfRr3NT1uw6yNb9CQy9oj3Dr+rEhAeuZ8/hFPIKip+WQETod/fljF/zGq0ubsaXL8/gvo6j2L/VmkLu3tDxqra88/sYHAVOHrvkeRL2n909PpGrOjZm8uMD8fGxMezt6WzYc8iSdisGhDOhwzAG1e7OjIOruHPlhxzKTbWkbQm8GomaCupEUwagudPK/JmRssAEAgOxNcnTcgAAIABJREFUV0eipkHgzZDzMZo21JLayCLC7XV6Mq79ENIKc7jjj/EsSNhyTm12aVKLN6cv4rWvFzJ+1jJuu7QtAb529iWkUrNSJL3be/ar/7hqO/5+dgL9fVHVM6qZXLNJdZ799lHeXvw8ORm5jGj3BN+88QMul3V1l61Uu3lNXvv5aXIz83i4xzPEbbMucNWLjWHiowOIDg1ixNgZLNtiTf1iu82HkQ378Gab2ziUm8pty99jUeI2S9oWv1ZIzEzw64RmPotmPIq6cyxp+0LltakhEZkI9AWSTrS7SEQigYlAXSAfGKqqp32XMFND3qW5MzxrB7ZIJOJdywreJOal8+SGL9maEc8ttbpxT/3Lzno6YF9CKgeT06kSFUb92BiWbd3Poo27+U/vjlSOCmVXfDJzV26nSY1KXNG+Ia9MW4DD6WJrXAJj7uh9RlMcaYnpjB3xMcu+X0X9tnW46dGr6XZdR+wW5Oix2p6N+xnd5yUcBQ5emf+0pZXaUjNzue+979l96CjP3X45V3ZobFnbh3JTGb1hGtszD3F77Z6MaHC5RQVv3JDzIZo9FnxqeX6efetbMOLzU2lNDU0Cep/i9dHABlVtAQwGxnpxLEYxSdB1SPTXIHY09RbLHq0rBUbwYcfhXF+9I1/u/53blr/HprQDZ9VW7cpR9Gheh/pF5wYWrP+TShGhVI4KJSMnn9+37CPAz071ihF88etadsYnMfLabtx2WVve+2EZ6dl5xe4rslIEz377KKOmjCQ7LYeXBr7DQ92fJjEu+azG7k11W9binWUvEhwexNP9XiHpoDUZaAGiwoL48KEbaFm3Kk999hPjZy3H7bbmQ2RsUBQfd7qLa6u15/N9i3lxywyc7nN/+hKxISH3IJGTQDPR1BvQvJnnPuALkNcCgaouAU418dcE+K3o2h1ALRGp5K3xGMUnvk2Q6Bng17no0foJVIv/5nkyfjY7TzTtz9i2d5DvcnDnyg95e/tc8pyFZ92mqlI/NoaDyemoKu/MWEJGTj6t6saSkpnDiu1xf2UtrVM5GruP7YxPzooIF9/Snc92jmXUlJEc2HGIO5s/zIyxc8vcdFGV2pV4cc4oCvIKebrfq2SlZVvWdmigP+NHXkf/Lk35ZN5KHvt4Drn5Z/9vdzw/m50nm17DsLoXM/vQWkZvnEaBy5q00+LfCYmeCfbmaMbjuDP+i+r5kZa8pJTmGsFG4DoAEemA55DaCfcvishwEVkjImuSk8veJ7ELkdgikMiPkJD7If8Hz8Kb8+w+wf9T5woNmNbtAa6v0YGv4pZxy7JxrEnZc3bjFOGmnq3IK3Tw0IRZ5BU4GHZlR1rWqcI3izdyVcfGxMaEU+hwkpKZQ9XoMPILz+4NxsfHh4tv6c4H69+gWbdGTHhoEvd3Gs3ib5ZTeAaL0t5Ws0l1npn+CAe2xzP6ypfJOctkfSfia/fhmUGX8cgNPVm8cQ9D3/qGwymZlrQtIgyvfykPN7qKRYnbeGjt5+Q4z70mA4D4VESiJkHw3ZA3HU25ybKf5wuBV7ePikgtYM5J1gjC8EwHtQY2A42AYap6ylSIZo2g5GnBYjT9UcCJhD4KgQMRixLMrUvdy4tbZhCfm8p11TtwX8PeZ526ODMnn6AAP+w+Nib+tIo1u+IZP/I6AOKT0/n0p1XUj43hlovPvcqaqrLo6+V8OmoqiXHJxMRG8dhn99LmUmvy8Fth+Q+reeHGt2h3RUvGzHrS8hO3y7ft58lPfiTI35dPHr7xr51dVph7aB0vbplBw7CqvNP2DiL8rKvv7Pl5fgzEjkR+ivhat95RlpXJ7aOqmqmqQ1S1FZ41ggqANVsSDEuJf0/PLgzf1mjm854KaE5rDpW3iarDl11Hcmutbsw8uJqBv48966eDsOCAv07BxoQHc1FLz2Lp3iMp/Lr+T1wutyVBADyfXi+6uSuf736Xl38cTWBoIE9cPoaHez7D2l+sSet8rrr0b8/wN25j5dx1zJ+0yPr2m9Tik4dvpMDh5O6x35GQmmVZ21fFtuG11reyOyuBu1Z+SIIFlc+OEf+eSPQ0wI6mDkILzQfLUgsEIhIhf58DHwYsUVVrnjENy4lPrOfTU9gr4NyBHu2H5nyK6rnPkQf4+PFAoyv5uNNd+Nt8uXf1RD7485dzWjCsHBnKB3P+YOz3Sxn3/e+kZeUxtHcHAEv3lfv4+NC+d2veX/0qw1+/jeSDRxnV+yU+fmIK8bsOW9bP2brm/j4079GY8Q9MtHRb6TENqlXg/fuvIzMnn7ve+ZbkDOvWJHpUbMy4dneQVJDJHX+8z6Y06zLaiL0uEv0V2GLQ1KFowSLL2j4feXP76DSgFxADJALP4kleh6p+ICKdgcmAC9gG/EdV007XrpkaKn3qSvRsMS1YAL4tPSeS7fUsaTvXWcBb2+cw+9BamkfUYEyLAWed2jouMY35a3bSqm5VGtesRGig/1nVMTgTBXkF/O+uD1kw9XdUlQ5XtuaBCcOpWL34mVGtlhyfwr3tnyAoLIj3V71CcLj1J6U37jnMiHdnUCUqlAkPXG9ZSgqAfdlJPLruCxLy0hnV9Br6VmtrWdvqSkHThoFzJxL+OhLY17K2yxqTYsKwnCc9xVxPegrNQULug+BhiFiTKuDnIxt5ZetMBGGUhTVwj+dwunhnxlLuuKKdpW9cAEcPpfDL5CV8+dJ3OAqddOrbhpsev4YmnaypIX2mNi/dzsM9n+GOF27m1qeu90ofa3Yd5IHxPxAdGsT7I6+juoVrBhmFuYzeOI3VKXu4tVY37mvY27JCSOrORtPuBsdqJOwZJOhWS9ota0wgMLxGXSlo1hhP5kd7Y8/TgW8TS9o+lJvK0xu/ZkvGQa6u1o5HGvUl0G5dVsmNew5z99hv8fe188SAi+ndvqHlTwtH9iUye/x8fp2yhMyUbG567Gq6XtuR+m1ql3iG0yd7v8i+TXFM2T8eXz9rAvY/bdmfwMj3vsdms/He/dfSqLo1GUYBnG4X/9sxl+kHVtC1QkPGtBhAiK9FNZE131PopuA3JOQBCB5xwaWzNoHA8DrN/xnNfB7caUjIPRB8DyLnfvrW6Xbx0e7f+HzvYqoFRTG83qVcUrmZZcVN9iek8uzkn9m87wit68Vyy0Wt6dmyrmWpl4/JycjhjaHjWfb9KgBqNa3Ofe/9h5Y9m1raz6msnLuWp/q9ylNfP0zPGzt7rZ99CancO24GWXkF/O+eq2nXoLql7c84sJI3ts+mRlA0b7e9ndgga8qNqjrRjNGQPxOCBiOhoyxJz15WmEBglAh1Z6CZYyB/Fvi2QyLeRHysKfSyOmUPb2ybxf6cZKoGRjKodnf6xrYlwIKslS63m28WbWTKb+s4kppJ5chQbr+8HTf1bGn5p8K0pAxWzl3HlBemkxiXTNOuDbnpsf507tfO659AXS4Xg+veR7WGVXlt/tNe7SsxLYsR42ZwJCWTiY8NsPTJAGBtyl6e2DCVYHsAH3W4k0qB1kxDqbrRrFchdxL490LC30Zs1k4blhYTCIwSpXk/oJnPAXYk/EXLUlu71c3SpB18vncxWzIOEukXzNA6F3FTzc6WvIm63G6WbNrL1AXrWPfnIXq1rMvzgy8nNMia6Yfj5ecW8ONHv/L9uz+SsC+JJp0b0O6KVjTqWJ82lzTHx+6dT6LvPzCRnz5dwMyMzy2rbHYyRzNyGPzaNBTliydu+at+hFW2ZxxixOpPiPIL4cMOdxITEGZZ25r7pedDjb02EvEBYq9hWdulxQQCo8SpMw7NeBgcmyHwZiRslCU1ksGzUL0+bR8T9yxiVcpueldpxX+bXWtdTntVvlywnrEzllIhIoRhfTpwVcfG+Hkh0ZzT4WTeJ7/xw/s/cWD7IVSV6o1iuf35AXS9pr3lye1+mriAt4ZNYNKuccTWq2Jp2yey82ASQ978mnpVY/j44Rvxt/jvsyktjvvXfEblwAg+6DCMSD/rPr1rwR9o+khAPLUN/DtZ1nZpMIHAKBWqhZ7Mjzkfg70eEv4/xLf4VcNO377y2d5FfPDnLzQNr8YbrQdZ+qlw457DvP7NQrYfSCImPJjbLm3LwItaW75+cExOZi5rf97IpKe/4uDOw0RWCqfXgK506d+e5t0bW/KUMPejX3jn7o+YGjehxLa0/rb+Tx77aA692zXkxSF9sNmsnQJbm7KXB9dOokZwBca3/w/hVp5CdsZ5dhS54op2FN1sWdslzQQCo1Rpwe9FNZIzkbBREHiLpfPhCxO38tym6QTZ/Rnd9Fq6V2xkWduqyqodB/hs/mpW7TxI9+a1efU/VxHo751dN+CZy189bwPzPv2N1T9twFHgIKpyBL0GdKV9n9Y07dKAwJAzf7pSVR7u+QyH/jzCtPgPvT41dLzP5q/i3ZnLuKF7C0YNvNjy9ZAVR//kkbWTqR9WhffaDbVsNxGAurM8T7cFiyHoViR0tGXbpEuSCQRGqVNXCprxBBQuAf9LkPCXEJs1uz0A9mQl8PSmb9idlUDf2DY83KivpW8GANMXb+S1rxfSqEZF7rqqE52b1PLa08Exedl5rP5pAwu+XMrKuetwOlz42H1o0bMJXa/pQIN2danesOppy2nG7zrMeyMnsvbnjTz4wXCuGn6ZV8f9T6rKez8s47P5q7n1kjY8fH0Py4PB0qTtPL5+KvVDq/BO29uJ8rdwmkhdaNabkPsp+HVGIsYitrJd0/qfTCAwygRVN+RORrPeAPGDoFuQoNsRH2t2lDjcTj7ZvYDJ+5YQ7R/KU82uo1OMtYVIFm3cw4tTfyU1K5eKESH069SEm3q1tPxA2onkZeexdfku1v+2mT9mrebgzr9TWMTERtGoY31qNa2O7bjg5Ha52bl6N+t+3Yx/kB9Dxgyk/329S2WPvKry5vRFTFu4gWF9OjLi6uLXky6u35N2MGrDNCoGhPFuu6FnfSr9ZDyFm54Gn6pI1OeW7YorCSYQGGWKOnahORMgfx7gA4HXIsHDEHstS9rfmn6Q5zd/y/6cZK6r3oGRDfsQZD+zGgSn4nC5WLp5HzOXbWH51v2EBvnz9KDLuLiVNWk2iuvQ7iPEbYsnfudh9mzcz/YVf3Jkb+K/rqtarzLdrunADY/0I7JS6X6KVVVe+vI3Zvy+mReH9La00tkxm9LieGjt5/j7+DKu3RDqhVa2tH0tXIumDQdbFBI1BfE5P8qomEBglEnqPIDmfAp53wEO8L8CCbnbkpPJ+S4HH/75C1/uX0aVwAieaX49baLqnPug/2F/Qir//Wwe2w8k0bdjY+64oj11qkRb3k9xud3uf32tpE8wn47L7WbYW9+wPzGN6U8PtnxbKXimCkeumUS+q5C32g6mVWQtS9vXwvVo2hCwVSkKBqX3b15cJhAYZZq6ktHcyZA7FTQPCXvKsnwvG9L288Lmb4nPTaVtVB0G1OxC94qNLMtTA56cRRNmL+fLBespdLro0LA6N1/Ump4t6lxwaQqssi8hlYEvTaFrs9q8ObyvV+7Tkbw07l/9GYn56bzcaiDdK1r79KGFq9DUYWCviUR9UebXDEwgMM4L6s5EMx6FgkWe4jdhT1myOyPPWcj0Ayv49sAKEvLTqRIQwQ01O3F1bDtLtxqmZeXy/bItfLtkEwlpWfTt2Jinbr3UK+cPLgSTfl7NuO9/5+Whfejd3rqdXsdLK8zmwTWfsyvrCP9re7vla0ZasAxNuwvsDTxrBrZQS9u3kgkExnlD1YVmv+05e+DX0XOQx2bNgp/T7WJp0na+jvuDdWn7CPbx56nm13FJ5eaWtP9XPy43E39axQdz/qBV3aqMvLYbLWpXtXz//PnO6XLzn7e+Ye+RFCY/MZDala3bRXa8HGcBdxYVt5nY6W5qhVib7kLzF6Lp94JvC0/NDpv1U11WMIHAOO9o3g9oxn/BpyIS/ibiZ01lsWN2ZR7h1a0z2ZJxkJtrdmVkw96WJbI7Zv7qnTw/5WfyC51UjgzlsrYN6NupCfVjS682QVmTkJrFra9MJTw4gMlPDCQk0LpF/eMdyUtjyB/jCbYHMLHTPZY+CQJo/nw0/UHwa++p9S3WpyU5VyYQGOclLdyIZjwErkMQNBQJfRARC3f/uJ2M3TGPbw78QYuIGrxwDkVwTiYnv5DFm/Ywf81O/tgWh8vtKZc54uouBHopFfT5Zu2ueO4e+y3dm9fhzeH9vPbktCktjntWfUKrqFqMbXuH5YFf82Z7pjb9uiIR7yE2a4PNuTKBwDhvqTsbzXod8r4Cn7pIxGuIr7VFan45somXtszAjXJXvUsZULOL5W8SABk5+YyftZzpSzZSo2IEw6/qxEWt6pmAAEz9bR1vfbuYJ2++mJt6tvRaP3MOreOFzd9yQ42OPN6kv+Xte84ZjC6aJvqoTC0gm0BgnPe0YKknV7z7KATfhYSM4O+S1+cuIS+dN7bNYmnyDhqGVWVU02toEl7NsvaPt3rnQV768lcOJKUTHODH5W0bcE3XZjSv7f0kcGWVqjL0za9Jycpl5nNDvLqeMm7nPKbsW8ojjfoyoJb1h9o0/xdPkRt7DSRyIuJj7TmGs2UCgXFBUHcmmvkS5H8P9kZI+GuIr3VbAlWVhYlbeXP7bFILsrmpZhfurn+ppYfRjnG7lfW7DzFrxVZ+XfcneQUO7rqqE3de2ancLirPX72TURN/5N17r6Frs9pe68elbp5c/yWLk7bxdLPr6WdhDeRjtHCVJ1mdhCJRn1pW0/tclEogEJGJQF8gSVWbneD1cGAKUAOwA2+q6mena9cEAkPzf/Mc83dnFNVKvtOSamjHZDvyeX/XfGYcXEXFgDAeb3K15XvQj5ebX8irXy1gzsrtXNSyLoMubUvz2lW8nseorHE4XVz1309oVKMS4+69xqt9FbqdPLL2C1al7OaFFjdxRVXrp6PUsQ1NGwbqRCI/Rvy8N+VVHKUVCHoA2cDkkwSC0UC4qj4hIhWAnUBlVS08VbsmEBgA6k7zlMbM/xHszTxJ7Cx8OgDP4uLLW79nb3YSPSo2ZkDNLrSNqo3NwsNox6gqU35bx7szf8fpchMa6E+nxjXo3rwOV7RviG8JZgotTRNmL+eTeSuZPWYoVaPDvdpXvquQB9d8zsb0OF5uNZCLKllfNlSdBzwnkN1HPQvI/t0t76O4Sm1qSERqAXNOEghGAdWBe4FawC9AA1X99xn545hAYBxP8370VEPTdE9WyKDBnhKDFtWadbidTNn3O1P2LSHLmU/VwEiurtaOvrFtqBhg/RtVVm4+K3ccYNnW/Szbup+jGTnUqRLFqIGX0La+d9YsypLDKZn0e/pThl/Zibv6eq+u8jE5zgLuXz2RHZmHeaPNILpWsK5exjHqSkbT/gPO3Uj460hgX8v7KI6yGghCgVlAIyAUGKCqc0/SznBgOECNGjXaxsXFeWvIxnlI3WmQ+w2aOxXcCeBTHQkaBIHXIzZrCtXkuxwsStzKrPg1rEndiw2ha4WGPNT4KqoFeSfPjKqyZNNeXv9mEUdSM+ndriGDL29nef3fsmbEuO+IS0xn9pihJbJekuXI497Vn7I3O4l32t5Ou+i6lveh7iw07R5wrPKsbQVea3kfp1NWA8ENQFfgYaAunieClqqaeao2zROBcTKqTij4Fc2ZDI41IEEQcA0SOtLS2gfxuSnMil/LdwdW4FI3jzW5miurtvZaXqG8QgefzlvFtIXryStw0KZ+LLde3IZeLetekLmMji0aTxh5HR0b1yyRPtMLc7lr5YdkOPL4qtuDRFh84AxAtcCTtbRwDRI1GfGzfpH6VE4VCEpzNWoIMEM9dgP78DwdGMZZEbEjAb2xRX+JRM+EgD6QNx09eg1auN6yfqoFRTOiweVM7TqShmFVeX7ztzy96WuyHHmW9XG8QD9f7uvflZ9eHsbD1/fgSEoWj3w4m8GvTWPtn/Fe6bM09WpVF7uPjZU7DpRYnxF+QYxpOYAMRy7/2zHHK32I+CMR48AnFk2/F3WWnX+70gwEB4BLAESkEtAQ2FuK4zEuIOLbBFv4K0j0NyB2NPVWNGcyVj4BVw6MYHyHYdxd/zJ+S9jC1Ytf563ts4nLOWpZH8cLDQpg0KVt+eGFITw/+HKSM3K48+3pjBj3HRN/WsXKHQfIyivwSt8lyd/XTo2KEexLSC3RfhuEVWVInV7MO7yBpUnbvdKH2MKRyA9AHWj6Pag7xyv9nClv7hqaBvQCYoBE4FnAF0BVPxCRqsAkoAogwKuqOuV07ZqpIeNMqTsDzXgSCn6DgD5I2EuIzdqKYjszDzNl31J+S9iCU110jmnAjTU60aVCA6/sMgLPlNGXv61nzoptxCWl/fX1WpWjuK5rM27q2fK8zXz62Eez+fPQUWY+P6RE+3W4ndy+fDzpjhy+6vYgYb5nXhu6OLTgdzTtTs/Ghoj3ES/9jBzPHCgzyj1VN+R84sls6lMVCXveK1v5jhZkMfPgKmYcXMXRgixaRdbi9da3EuHn3YyUmTn5bDuQyJb9CazcfoC1f8YTGx3Gff27cVnbBufdIbXxs5Yz8adVLB97X4kHsx0ZhxiyYgJ9qrbimeY3eK0fzZmMZr0IwcOxhT7qtX6OKatrBIZRYkRsSMhwJOoLwBdN+w/u9EdQV4ql/cT4hzKs3iX80PMxRje9lm0Z8QxdMYH92UmW9vNPYcEBdGpck2F9OvLxwzcyfuR1BAf6M2rij1z9zEQ+nbeK5Ixsr47BSgUOJ3YfW6lUV2sUHsvg2j2Yc2gdCxK2eK+joNsg8GbI+QjNm+m9forBPBEY5Y5qIZr9AeR8CBKEhD4OgTd4ZQfO5vQDPLbuCwrdLh5vcjW9KjUlwKdkksy53G4WrN/Nt0s3sXrnQXxsQucmtbi4VT16NK9DVFjZyo55vPve/Z6jmTl89d9BpdJ/odvJ3Ss/Zm92IhM730OdEO/UJVZ1oGlDoXAdEjUJ8WvvlX7ATA0Zxgmpczea8Yxnq6lvByT8BcRufV3jI3lpPLLuC3ZnJRDk40fPSk24rHILOsbUw9dWMtMeB5LS+H7ZFn5es4sjqZmIQMs6VbmkdX2uaNfQK3WDz1ZWbj5XP/MZXZvW4sUhfUptHIn5Gdy+/D1C7IFM6jyCEF/v1BhQdxqacjO4k5DISV5LRWECgWGchKrbs8U06w1PveSQeyB4uKWZTcFTHW1d6j5+TtjEwoQtZDnzCfMN5JJKzRha72IqeeGU8omoKrvik1m4cQ+LNu5hV3wyNhE6Nq7BlR0a06VpLSJDvLNAWlyvTFvAd0s3MfmJgTSp6Z1P4sW1LnUf967+lC4xDXijzSCvLfyrKwFNvdWTPytqMuLbxPI+TCAwjNNQVzKa9ZInd5FPHSR8jNce0x1uJyuO7uaXI5tYmLgVX5sPDzW6ir6xbUr8gNjeIyn8uGoH81bt4Eiq5yxnvarRtKlfjXYNqtGhYQ3Cgkuu2tayLfsYOX4mN/dqxWM3XVRi/Z7K13HLeWv7HIbXu4Rh9S7xWj/qjPcEA81HoqdanrHUBALDKCYtWIxmPAfuQxBwNRJ0K/i28tobdHxuCmM2f8f6tP10jmnA6GbXltjTwfHcbmXL/gTW7DrIml3xbNx7mLwCBz42oV2D6lzUqh4XtaxLhQhrt90eE5+cztjvl/Lb+t3UqBjBlCdv8VrZyjOlqjy/+VvmHd7AO21vp3OFBt7ry7nfEwwAiZqK2GtZ1rYJBIZxBtSdi+a8B7lfguaCvR4SeAMEXmNpqopj3Opm+oEVvL9rPm5VelRsTN/YNnSIrueVSmnF4XC52BaXyJJNe1m4YTf7Ez3nFKpEhVGzUgQ1KkZSs2IkNStH0rh6RSJDz3zh2eF0cSQ1k5nLtjB1wXp8bMKQK9oz6NK2Za5qW77Lwa3LxhHmG8hnnUd4tS917kZTbgUJKAoG1iQbNIHAMM6CurMhfx6aNx0cGwBfCLjUExT8uliW4fSYQ7mpfBW3jJ8ObyTDkUuMfyi9q7aib2wbr+1aKa69R1JYsmkvfx4+SlxiGgcS08jO/ztjfGx0GE1qVaZZrcrExoQjgIggAAL5hU4Op2QQn5xB/NEMDh3NICE1C3fR+0+/Tk24r39Xrz1xWOGbuD94c/tsPut0D00jqnu1L3VsQ1MHgy3cEwwsqHJmAoFhnCN1/OkJCHkzPSmv7Q2R8DcRX+vTFjvcTn5P3sncQ+tYlrwTl7oZWLMr9zfsXWpPCP+kqqRm5bL3SCrb4hLZGpfA1v2Jf60znExUaBCxMeFUqxBOtZhwYmPCaVqzEnWrxpTQyM9etjOfvgtfpVelpjzX4kav96eFG9G0O8BWCYn+9pxPw5tAYBgWUS2E/Plo1ivgzkRCn4CgQV5bQ0gtyObTPQuYfmAFHaPr8VKrgV5Le2CFlMwcjmbkoHiCheftRfGz26kaHUZQgLW7sUraG9tmMfPgauZc9ASRft5/etGClWja7RB4I7bwMefUlgkEhmExdaWgmU9CwWJPvpiwVxAf79QlAJgVv4ZXt/5AlcAIHmp0Fe2i65bYwTTjb/uykxjw+zuMqH85d9TtVSJ9ujNfg9xPkcjPEP+uZ92OCQSG4QWqCrlfoFmvgy0MCRnp2Wlk886J3Y1pcTy54UtSCrLwt9lpE1WHLhUa0CWmIdWDvReEjP/v3tWfciDnKDN7PoZPCSSLU81Hj/YHLURiZp/1FJEJBIbhRerYiWaMAucWkDBPuoqgWxB7Dcv7KnA5WJ+2n+XJO1mevIsDuZ6U19WDormuekcG1OxcZtYRLlS/Jmxm9IZpjGs3hE4x9UukTy1ci6beAkEDsYU9d1ZtmEBgGF6mquBYg+ZOgfyfAbdnyijotqIdRt755Bifm8Ly5F0sTNzC2tR91AupzBNN+9MysmQqe5VHhW4nVy585a81m5LiznwZcichkZMR/05n/P0mEBhGCVJXApr7FeR9De4UsNdHIsZaflL0//WpyuKk7bwH1xYIAAATI0lEQVS1fTaJ+Rn0i23LiAaXE+0f6rU+y7O3ts9mxoFV/HjRKMK9UNbyRFTz0KP9ADcSPeeMpyBNGmrDKEHiUxlb6INIhcVI+BvgTkNTBqAFy7zXpwi9KjXh624PMqh2d348vJ7+i99gzObv+DPriNf6La/6xbbDoS7mH9lYYn2KBCLhr4DrEJr9pqVtm0BgGF4i4ocE9keip4NPVTRtGJo7zat9Btn9GdmwD191e5B+sW35JWETty57lxGrPmFp0nbc6vZq/+VFg7AqNAqLZVZ8yc5OiF97CBoMuVPQ/J8sa9cEAsPwMvGJRaKmgX93NPNZ3JljLC+I8081g2N4oml/Zvd8gvsa9OZAbgqPrPuC25a/z6qju73ad3lxWZXm7Mo6QmpByRb8kdDHwLcVmvEk6vjTkja9FghEZKKIJInICUv8iMhjIrKh6NcWEXGJiPWJXAyjDBBbCBIxAYJu92w5Te6GO3Uomvsd6j71adxzEe4XxOA6PZjZ41Geb3Ej2c587lszkQfWTGJPVoLX+i0P6oV40j7E5SSXaL8ifkjEuyDBaPoIS35+vFm8vgeQDUxW1WanubYf8JCqXny6ds1isXG+U8cuNH8u5M8F1wHAF/y7IwFXgf/FiM17RWIK3U6mx/3BxD0LyXEW0LlCA5pH1KBZeHUah8cSWoZPLZc1R/LS6L/4DUY3vZZrqnuvstjJaOEaTz4i/25IxAen3Zl2qsVir5VHUtUlIlKrmJcPBLw7eWoYZYT4NkB8G6AhD4JzC5o3F/J/RAsWAAFo8CAkZCQi1tcB8LPZubV2d/rGtuXzvYtZmrydZck7/3q9ZnAMTcKr0zS8GpdXaUGEX9mpXFbWVAoIx9/my/4SfiI4RvzaQehoNOsFNPtdJPSBs2/Lm9tHiwLBnFM9EYhIEBAP1FPV1JNcMxwYDlCjRo22cXFx1g/WMEqRqhsc69DcbyB/JvjURsJfRfxae73vLEce2zMOsTUjnq0ZB9mSfpDUwmyC7f7cUacXA2p2MeksTmLQsnepEBDG/9reXir9qyqaOQryZiARE5CAkxfOKbVzBMUMBAOAQararzhtmqkh40KnBcvQjP+C+wgEDUFCH/TK08FJ+1dlT3YiE3b9zNLkHVQKCOee+pfTu2pLr5VqPF+N3jCNHZmHmNHj0VIbg2oBmjIQXPs9WUpPUne7rJ8juBkzLWQYfxH/rkjMHAgcALkT0aP90by5qBaUTP8i1AutzFttBzOh/TAi/UJ4bvN0Bv4+lq/2LyPTkVci4zgf5DgLCLaXXJA+ERF/JPI9wI5mvnxWbZRqIBCRcKAn8ENpjsMwyhqxhWALfwGJ/AxwohkPoUldcWc868lTX0IZAdpG12FS53sY02IAQXZ/3t4xl6sWvsKzm6azMS2uxMZRVu3LTqJ2cMXSHgbiUxUJvg0Kl6DOM58699pisYhMA3oBMSISDzwL+AKo6gdFl10L/KyqOd4ah2Gcz8S/K8T8DIUr0LzvIe97NG8a+NSBwGshsL8l1atOxSY2rqjakiuqtmRX5mG+P7ianw5vYN7h9dQMjuGiSk3pWbEJjcNjy9XUUa6zgIT8dGqFVCjtoXgEDoDsCWjuVCRs9Bl9q8k1ZBjnkb/LZ84Ax1rABoHXI2HPIFJyxd5znQX8krCZ+Yc3sD5tPy51E+MfSveKjelVsQlto+vgZ/Pa58wyYV3qPu5e9TGvtbqFiyqfcod8iXGnPwQFS5AKS/+Vi8gknTOMC5A696O5UyH3c/BtiUS8j/iU/DRFRmEuy4/uZHHiNv44+id5rkJC7QFcUrk5faq2omVkzQvuScGlbu5c8SGH8lL5tvvDZeb8hSdd9UAk7AUk6Ob/95oJBIZxAdP8+WjGEyAhSOR4xLdFqY2lwOVgVcoefk3YxMLEreS7HFQJiOCKqq3oU7UVtYIreK2sZ0n69sAKXt82i+db3Eifqt7f4ltcqoqmXAO4kOjZ/+9em0BgGBc4dexA00eAK8kzPxzQF7GFleqYcp0FLE7axrzDG1h1dDdulEoB4bSJqkObqNq0japNbGDUeRcYjhZkcdPS/9E4PJb32g0tc+PX3G/RzNFI1BTEr8NfXzeBwDDKAXWnoukPQOFKwAd8WyP+PcC/J9gbleob1tGCLBYlbmVNyl7Wp+0jrdCzP6RiQDhtImvTuUJ9ulVoVGamWE5mS/pBXt82iz1ZCXzZ7QFqBseU9pD+RTUfTeoBfp2wRY776+smEBhGOaHqAscGtGAJFCwG5zbPC7aK4N8D8e8Jfl3Puu6tNWNU9uUksS51X9GvvaQW5uAjNtpF1eWiSk3oUakJMWWoqE5CXjrv75rP/CMbifIL4bEm/bikcvPSHtZJeaqZTUUqLkVsnlyeJhAYRjmlriQoWIoWLoGC30GzQIKRkAcgaBAipb+zx61utmUcYmHiVhYlbuVgbgqC0CKiBs0iqlMjOIYaQTHUDI4h2j+0RJ9scp0FTN63hKn7lgJwS61uDK7Tk2B7ye3QOhvq2IWm9EVCRyPBdwAmEBiGAag6PE8L2R9C4RLPdFHY8yWSz6i4VJW92YksTNzK0qQd7M1OpMDt/Ov1IB+/vwJD/bAqNA2vRuPwapa/MSfkpbM4aRuf713M0YIsLq/Sgvsa9KZyYISl/XiTO+V60EIkehYiYgKBYRh/U1Uo+BnNfBHciRB4ExL6KGIre29ybnWTlJ/JgZyjxOUkcyD3KAdyUtifk8yRvDQABKF2SAWaFmVNbRJejRDfv9M+CKd+glA874GJeRksP7qL5ck72ZOdCECz8Oo81PgqmkfU8NLf0Hs0dxqa+SwS/R3i29wEAsMw/k3d2Wj2e55zCBLmmUIIGvDXnHJZl16Yy/aMeLZkHGRbRjxb0g+S4cg9pzbt4kOryJp0rdCQLhUantfbXdWdhSZ1gcDrsYU/ZwKBYRgnp46daNbrULgU8PekrQi6HfGtX9pDOyOqyqG8VHZkHqbQ5ZlOOvZpX///hVD05n7sLV4QQn0DaRNVu8zP/58Jd/qjULAQqbgMmy2w5AvTGIZxfhDfhkjUp6jjTzR3MuTNRPO+Qf26IsG3g1+P01a/KgtEhGpB0VQLii7toZQZEngDmj8L8n8+5XVl/1/XMIwSIb71sYWPQSouQUIeBuduNG04mtIXzV9Y7jONnpf8OoCtgmc78SmYQGAYxv8jtkgk5G6kwgIk/G1QF5p+F5p2O+rYWtrDM86AiA3s9cG175TXmUBgGMYJifgigX2RmDlI6DPg2ImmXIs7/VHUdbi0h2cUl70OOPee+pISGophGOcpEV8IHgSB/dGcjyBnEpr/E+rfHfFtCb4twLc5Yis7J4GNv4lPHU5X8sUEAsMwikVsoUjoI2jQwKJDaX+gBb/99br61PGkw/Zr4QkO9oaI+JXiiA3A80RwuktKYBiGYVxAxKcqEv48AOpOB8cWcGxCHZugcCma/33RlX5o4LVIyEjEp4xU8SqPTCAwDMObxBYB/t3AvxtC0all9xFPYChYBnnfofmzIXgYBA39V9UsowTYKoEEn/qSEhqKYRjlgIh4nhgCenu2osbMA78eaPY49OjlaO43ngypRokREfA5dYoMrwUCEZkoIkkisuUU1/QSkQ0islVEFntrLIZhlA6x18QWOQ6J+gp8YtHMp9CUqz0BwX3qBUzDQu60U77szSeCSUDvk70oIhHAeOBqVW0K3OjFsRiGUYrErw0S9RUS8S4gnoCQ3A13xrOoY1tpD++Cpu5scCec8hqvrRGo6hIRqXWKS24BZqjqgaLrk7w1FsMwSp+IQMAV4H85ONajuV9B3gw0bxrq2xwJvBkCrkRsp57PNs6Q69RnCKB01wgaAJEiskhE1orI4JNdKCLDRWSNiKxJTk4uwSEahmE1EUH82mCLeB2p+DsS+hRoHpr5X89TQvoDaN5s1J1Z2kO9MDj3nPaS0tw1ZAfaApcAgcAfIrJCVXf980JV/Qj4CDzZR0t0lIZheI3YwiF4MATdBo51aN73ULAAzZ8H2FG/Doj/pRBwKeJTubSHe15S525O91ZfmoEgHjiqniNvOSKyBGgJ/CsQGIZxYRMR8GuL+LVF1Q2OjWjBr5D/K5r1AmS9gNqbIQGXQUA/xF6ttId8XlB3KuROB792wPaTXleaU0M/AN1FxC4iQUBHTjVSwzDKBREb4tcaW+hj2CrMR2LmISGPgtjR7P+hRy/GnXILmvsV6s4o7eGWaZr1Bmi2Z/rtFLz2RCAi04BeQIyIxAPPAr4AqvqBqm4XkZ+ATYAb+ERVT7rV1DCM8knsdSGkLhIyHHUdgrxZaN4PaOYzkDkG9b8ICbwa/HuZlBbH0cLVkPcdBN+J+DY45bWmQplhGOcdVQXnNjTvB8ifA+6jIGHg1wHx6wh+7cHe6LwoqOMNqoXo0f6g+UjMXMQWdMpSlSbFhGEY5x0RAd+miG9TNPRxTwK8/HlQuNKztgAg4ahfO8SvPfh1LAoMPqU78BKgqmj2BHDtQSI+LFZaDxMIDMM4r4nYwb874t8dAHUdgcJVaOFKKFz9d4ZUCUX9uiABl3umkS7AtNnq2IZmvQKFKz1nMgIuKtb3mUBgGMYFRXyqQGB/JLA/AOpK8ASEwhVQsAgtmA/4on6dioLCpYjP+V3nWF1JaPb/IG8GSAQS9hwE3lTs7zdrBIZhlBt/bU3N/xkKfgbXQUDAt61na6pva7DX9pxvOA+o5kPOZ2jOh6AOCBqMhNyD2ML+da1ZIzAMw6Cohq9fa8SvNaqPg3NnUVD4xTOlUkQlEuy1PUHBp5bn9z61wV6zTOxMUlXIn4tmvQnuw+B/GRL6OGKveVbtmUBgGEa55FlwboT4NoLQkagzHpy7wLUPde4D534oWIq6vzvuu2yoT3VPaU7fluDbEnwbI+Lv9fGqKrj2QMFiNG8uOLeAvQkS/hri3/Gc2jaBwDAMAzynlYtOLMtxX1d3Nrj2g/NYgNjlWXPIn1N0hS9qbwh+Lf8ODj41Ldm6qpoHBSvQwsVQsARc8Z4X7A2QsFcg8BpLdkKZQGAYhnEKYgsBWzPwbfb/A4QroahE50ZwbIa8mWju1KJX/VA5m7fXf6zZaiHgAgkEvy5I8HDw74H4VD27v8xJnHeLxSKSDMSV4hBigKOl2H9ZZu7NyZl7c3Lm3pyclfempqqesHj0eRcISpuIrDnZynt5Z+7NyZl7c3Lm3pxcSd2b8nn+2jAMw/iLCQSGYRjlnAkEZ+6j0h5AGWbuzcmZe3Ny5t6cXIncG7NGYBiGUc6ZJwLDMIxyzgQCwzCMcs4EghMQkd4islNEdovIkyd4vYaILBSR9SKySUSuLI1xlpZi3J+aIvJb0b1ZJCLlosCsiEwUkSQROWGlPfEYV3TfNolIm5IeY2kpxr1pJCJ/iEiBiDxa0uMrbcW4P7cW/cxsEpHlItLSyv5NIPgH8ZzXfh/oAzQBBopIk39c9hTwjaq2Bm4GxpfsKEtPMe/Pm8BkVW0BvAC8QvkwCeh9itf7APWLfg0HJpTAmMqKSZz63qQCI/H87JRHkzj1/dkH9Cz6PzUGixeRTSD4tw7AblXdq6qFwFdA/39co8CxPK/hwOESHF9pK879aQIUVQNh4QlevyCp6hI8b2gn0x9PgFRVXQFEiEiVkhld6TrdvVHVJFVdDThKblRlRzHuz3JVTSv64wrA0qdsEwj+LRY4eNyf44u+drzngEEiEg/8CNxfMkMrE4pzfzYC1xf9/logVETO78of1ijOvTOM0/kPMM/KBk0g+Dc5wdf+ucd2IDBJVasBVwJfSPmpkl2c+/Mo0FNE1gM9gUOA09sDOw8U594ZxkmJyEV4AsETVrZrso/+WzxQ/bg/V+PfUz//oWg+T1X/EJEAPMmhkkpkhKXrtPdHVQ8D1wGISAhwvapmlNgIy67i/GwZxgmJSAvgE6CPqqZY2XZ5+RR7JlYD9UWktnhKEd0MzPrHNQeASwBEpDEQACSX6ChLz2nvj4jEHPeENAqYWMJjLKtmAYOLdg91AjJU9UhpD8oo+0SkBjADuE1Vd1ndvnki+AdVdYrIfcB8wAeYqKpbReQFYI2qzgIeAT4WkYfwPNrfoeXkiHYx708v4BURUWAJcG+pDbgEicg0PH/3mKL1o2cBXwBV/QDPetKVwG4gFxhSOiMteae7NyJSGViDZxOGW0QeBJqoamYpDblEFeNn5xkgGhgvIgBOK7OSmhQThmEY5ZyZGjIMwyjnTCAwDMMo50wgMAzDKOdMIDAMwyjnTCAwDMMo50wgMAzDKOdMIDAMwyjnTCAwDAuIyEwRWSsiW0VkeGmPxzDOhDlQZhgWEJEoVU0VkUA8aTh6Wp0PxjC8xaSYMAxrjBSRa4t+Xx1P8RkTCIzzggkEhnGORKQXcCnQWVVzRWQRnkSEhnFeMGsEhnHuwoG0oiDQCOhU2gMyjDNhAoFhnLufALuIbMJTT3ZFKY/HMM6IWSw2DMMo58wTgWEYRjlnAoFhGEY5ZwKBYRhGOWcCgWEYRjlnAoFhGEY5ZwKBYRhGOWcCgWEYRjn3f6AperxOwwDIAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"m.draw_mncontour('a','b', nsigma=4, numpoints=100) # nsigma=4 says: draw four contours from sigma=1 to 4"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gcV7n48e+rbjVLsoqLJMuyVSw7duzINXasJCSYlJsGAaeHBNMuod0L4YYEXy5cCPzgQijhGlJJhwTSSMPJOonj7rjbstwky5bVLWnVy/n9sWuuItS1s7Pl/TzPPFrNzM68c3b33dkzZ84RYwxKKaWCR4jdASillPIuTfxKKRVkNPErpVSQ0cSvlFJBRhO/UkoFGU38SikVZDTxK6VUkNHErzxKXO4XkVr39BMRkUHWvUdEykSkUUSeEZH4XsuniMiLIlInIuUi8oU+z18rIsUi0iMit3kg9htEpFREmkXkryKSNMi654rIdhFpcf89t9eyNSLSKSLOXlN2r+UXicgO9zEfFZHVfbb9FRE55l6+TUSWjfXY+ok/SUT+4j7WUhG5YZB1B31NRcS4t3P2WP/Qa1mkiPxORCrdr+PLIjKl1/IsEfmbiNSLyGkR+bWIhHn6eFUfxhiddPLYBHweKAbSgSnAfuALA6x7K3AQyABigReBx3otfwf4BRAOzAXqgAt7Lf8ycDGwDbhtjHHPApqAC9yxPAU8M8C6EUAp8HUgErjL/X+Ee/ka4IkBnhsONLjLSYAFgBOY616+CGgGznMv/yJQDYQO4xhuAx4d5vE+DTzrPtZl7phmjeY1BQwwY4DnfgvYBaQBUcAfgRd6Lf8b8Kh72URgD3CX3e/jQJ/0jD+IiMjdInJERJpEZL+IXGPBbm4FfmaMKTfGnAR+hish9edK4CFjzAljjBO4H/i0iESLSCxQBPzQGNNpjNkF/Bn47NknG2N+Y4xZB7R5IO4bgZeNMe+6Y7kXuFZE4vpZtwgIA35hjGk3xjyAK0lfNIz9JAHxwB+Ny1bgAFDgXp4F7DPGbDeuzPg4kAykjv7QPkpEYoDrgHuNMU5jzPvAS8DNAzxlJK9pX9OAN4wxlcaYNuAZXF+yvZc/Z4xpM8acBl7vs1xZQBN/cDkCLAfGA/8JPCEik/pb0V3tcWaQKXOAfczCdYZ31i4G/iCLe+r9fySQ02t+3+WzBzm+sfhI3MaYI0AHkDvAurvdifms3Xz0OK90V23sE5Ev9tpuJa6z7dtFJFRElgBTgffdq7wGhIrIIhEJxfVFtxM4PeYj/D+5QLcx5lCveYO9TsN5Td91V9W8ICJZveY/BJwvIpNFJBrXF+xrvZb/EviM+8t+CvAJXMlfWUgTfxAxxvzJGHPKGNNjjHkWKAEWDrDuU8aYhEGmsgF2E4ur2uCsBiB2gHr+14A73fW844Fvu+dHG2OagA3AvSISJSLzcZ2lRo/8yIelb9xnY+/vjH+odZ8DZgIpwOeA+0RkVa91nwbuA9qB94B7jDEn3MuagOdxfRG0A98DVvf5khmrkRxrf+v3fU1X4Pqlkg+cAl7pVU9/CCgDTgKNuMrl+722tR7Xl0gjUI6r2u6vIz4iNSKa+IOIiNwiIjvPnrXjOntO9vBunLiqMs6KB5wDJK6HcSVBB7APV50+uBIAuM4OpwEngAeBJ3stGzYRyex9oXWYcZ+NvWmk6xpj9ru/YLuNMR/gOqv9pDuWfFx167fgulYwC/iWiFzu3s6duM7yZ7mX34QrkU4e4Nh+2+v1/C3Q+5fabg8ca3/rf+Q1dVePdRhjzgBfxfWazXSv+yCu+vsJQAzwAu4zfhEJAd5wz4vB9V5MxFXlpyykiT9IiMhU4PfAvwITjDEJwF4+WpXSe/0b+7RK6TsNVNWzD9eF2LPmuuf9E/cvj+8ZY7KMMenu9U66J4wxpcaYK4wxKcaYRbiSx5aRHrsxpswYE3t2Gk7c7lY4kbjOWPtbd06fXzFzGOA4cV38PLvubKDYGPOG+/iLgVdxVXHgjuFlY8wh9/LXgQpg6QDH9qWzv8KALwG9f6nNGSCeQ0CYiOT0mjfg68QIXtN+jncurgvOdcaYduBXwEIRScZ1vSMD+LX7Wkkt8Ahw2SDbVh6giT94xOD6QFYDiMjtDFJfbox5sney7GcaqKrnceAb4mqKORn4Jq5WG//E3aRwuru5YAHwc+D7xpge9/KZIhInIhEichNwqXuds8+PEJEoXEkm3F0lNNr39JO46uWXuy9+fh9X65P+zoIdQDdwl7u54r+657/tjusqEUl0H9dCXK1+XnSv8yGQI64mnSIi04Er+L869K3A5SKS7V5+Ca46+b2jPK5/YoxpxnWW/X0RiRGR84GrcLW46c+Ar6mIzBJX09ZQ9wX5n+H64j7Q63huEZHxIhKO68vplDGmxhhTAxwDvigiYSKSgOtC8i6UtaxsMqSTb03AD3E1iazBlUDXA3d6eB8C/MS9nzr3Y+m13Aksdz/OxdVMsAVXc8hv9NnW13B9UTXjqvMu7LPcgevLrPdUNIbYb8BVH92MK1En9Vr2GvAfvf6fB2wHWoEdwLxey54Gat3HepA+zROB63El8iZcVVf3AyG9yu/77jiacCXQm4cZ/20MvzlnEq669Gb3vm7otWw5rqqcIV9TXC2Zit3bqXJvM6fXcyfg+lKtAs64X8eFvZaf634d693vyz8BqXZ/VgJ9OvviKaWUChJa1aOUUkFGE79SSgUZTfxKKRVkNPErpVSQ8Yte8JKTk01WVpbdYYxZc3MzMTExdocxqFNNpyzbdmRPJO0h7ZZt39uST9QCUJMxYdjPCbQyGKnJcZP94nNgNW+Vwfbt22uMMSl95/tF4s/KymLbtm12hzFmDoeDoqIiu8MY1BrHGsu2nefMozi22LLte9ttX3sUgEd/cduwnxNoZTBSa4rW+MXnwGreKgMRKe1vvlb1KKVUkNHEr5RSQUYTv1JKBRlN/EopFWQ08SulVJDRxK+UUkFGE79SSgWZgE7875VU81vHYbvDUEqpEWvr7OaBdSWcbmjz+LYtS/wikiEi74jIAfeA01/ts/zfRMS4R+KxxPslNfz8zUNUNXm+4JRSykpv7q/k528d4nDVQKOFjp6VZ/xdwDeNMTOBxcCX3aMsISIZwCW4BoCwzPULMujqMbyw46SVu1FKKY97dmsZ6YnjWDp9+F2CDJdlid8YU2GM2eF+fHYkoSnuxf8DfAvXiEmWmZ4Sy4KsRJ7begIdcEYp5S/KalvYcLiW6wszCAnpd1jsMfFKXz0ikoVrqLrNIvIvwEljzK6PjlX9T89ZDawGSEtLw+FwjGrfc2I7eeh4B2v/8jZ5SaGj2oanOJ3OUR+Ht+Q58yzbdmR3pKXb97Zx3dHAyMos0MpgpBwOh198Dqw2VBk8X9KBAFM6TuBweL7GwvLE7x6A+Xlc46d2AffgGjR7UMaYtcBagMLCQjPaDo0WdnTxbMk6DnUl8/miuaPahqf4Q+dU2knb8C0JbQEY0TEFWhmM1KqiVX7xObDaYGXQ1d3D3R+8w4q8BK77xEJL9m9pqx4RCceV9J80xrwATAemAbtE5DiQDuwQkYlWxRAdEcaVcyfz6p5TNLZ1WrUbpZTyiHdLqjnd2MZnFmRYtg8rW/UI8BBwwBjzcwBjzB5jTKoxJssYkwWUA/ONMaetigPgMwsyaOvs4aWd1vU1r5RSnvDs1hMkx0ZwUX6aZfuw8oz/fOBm4CIR2emeLrNwfwOakz6e/IlxPLv1hB27V0qpYaluamfdgSqunZ9ORJh16dmyOn5jzPvAoJej3Wf9lhMRPr0gg/98eT97TzYwe8p4b+xWKaVG5M/by+nqMVxfaF01DwT4nbu9XTNvCpFhITyz1dJbB5RSalR6egzPbi1jYVYSM1JjLd1X0CT+hOgILjtnEi9+eIqWji67w1FKqY/YdLSW47UtrFpk7dk+BFHiB1i1MJOm9i5e2V1hdyhKKfURT289wfhx4Xxi9iTL9xVUiX9BViLTU2J4ZotW9yilfEets5039p7mmnlTiAq3/kbToEr8IsKqhZnsKDtD8ekmu8NRSikAXthxko7uHlYtzPTK/oIq8QOuZlKhITytZ/1KKR9gjOHprWWcNzWRvIlxXtln0CX+pJgIVs6eyAs7ymnr7LY7HKVUkNtyrI6j1c2W3qnbV9AlfnBd5G1s6+JVvcirlLLZU1vKiIsK44o5k722z6BM/Iuzk8hOieHJzaV2h6KUCmK1znZe23Oa6+anMy7Ce70HB2XiFxFuXDSVHWVn2H+q0e5wlFJB6s/by+no7uGGRd65qHtWUCZ+gOvmu+7kfWqLnvUrpbyvp8fw1BbXnbq5ad65qHtW0Cb+hOgIrpgzmb/sOImzXe/kVUp514YjNZTWtnDjYu+e7UMQJ36AGxZl0tzRrd01K6W87qnNZSRGh7NytmXDkQwoqBP//MwE8ifG8eTmUh2TVynlNfVtPby5v5JPFWYQGeb9IWGDOvGLCDcunsq+U43sPHHG7nCUUkHi3fIuunuM1+7U7SuoEz+4umuOiQjliU16J69Synpd3T04TnSxPCeZackxtsQQ9Ik/NjKMa+en8/LuU9Q3d9gdjlIqwP39QBX17YabFk+1LYagT/wANy2eSkdXD3/arkMzKqWs9cSmUpKihIvzU22LQRM/kDcxjoVZSTyxqYyeHr3Iq5SyxpFqJ+8frqEoI4ywUPvSryZ+t5uWTKWsroV3S6rtDkUpFaCe3FRGeKhwQbplw50Pi2WJX0QyROQdETkgIvtE5Kvu+T8VkYMisltE/iIiCVbFMBIrZ00kOTaCJzbpnbxKKc9r7ejmz9tPsHL2JBIi7T3ntnLvXcA3jTEzgcXAl0WkAHgLmG2MmQMcAr5jYQzDFhEWwmcWZLLuYBUn6lrsDkcpFWBe2nWSxrYubrbxou5ZliV+Y0yFMWaH+3ETcACYYox50xhzto+ETUC6VTGM1A2LMhHgyc3atFMp5TnGGB77oJS8tDgWZCXaHQ5eqWgSkSxgHrC5z6LPAs8O8JzVwGqAtLQ0HA6HdQH2Mi81lD9+cIT5ERVEhIpHt+10Or12HKOV58yzbNuR3ZGWbt/bxnVHAyMrs0Arg5FyOBx+8TnwtJL6bvZXtHHbrAjWr19vexlYnvhFJBZ4HviaMaax1/x7cFUHPdnf84wxa4G1AIWFhaaoqMjqUAGIzKhl1e83cWb8DK4v9OyIOA6HA28dx2itcayxbNt5zjyKY4st2763LQl1VQmO5JgCrQxGalXRKr/4HHjan5/aQXxUNd/69IVER4TZXgaWXmEQkXBcSf9JY8wLvebfClwB3Gh8rJOcxdlJ5KXF8dgHx7X/HqXUmFU2tvH63tNcX5hBdIS9rXnOsrJVjwAPAQeMMT/vNX8l8G3gX4wxPncVVUS4Zamr/54dZfV2h6OU8nNPbi6j2xhuXmL/Rd2zrDzjPx+4GbhIRHa6p8uAXwNxwFvueb+zMIZRufrcKcRFhfHYB9q0Uyk1eh1dPTy1uYwL81KZOsGefnn6Y9nvDmPM+0B/V0f/ZtU+PSUmMoxPnZfB4xuP893LZ5IaH2V3SEopP/Ta3gpqnO3c4kNn+6B37g7oliVT6TaGJ7Rpp1JqlB7ZcJxpyTFckJNidygfoYl/AFnJMVyUl8pTm0tp7+q2OxyllJ/5sKyenSfOcOuSqYSEeLZp+Fhp4h/EbednUePs4OVdFXaHopTyM49+cJy4yDA+6eFm4Z6giX8Qy2Ykk5MayyMbjmnTTqXUsFU2tvHq7go+VZhBbKRvNOHsTRP/IESE287PYt+pRraVatNOpdTwPLGplG5juHWpb13UPUsT/xCumTeF8ePCeWTDMbtDUUr5gbbObp7aXMbF+b7VhLM3TfxDiI4I4zMLMnh972lOnmm1OxyllI97adcpaps7uP38aXaHMiBN/MNw85KpiAiPf3Dc7lCUUj7MGMMjG46TmxbL0ukT7A5nQJr4hyE9MZqVsyfy1JYymtu7hn6CUioobTxSy4GKRu5YNg1XrzW+SRP/MN2xbBpNbV38aZsOyK6U6t8f3j/GhJgIrjp3it2hDEoT/zDNz0xkfmYCD284TrcOyK6U6uNItZO3D1Zx0+KpRIWH2h3OoDTxj8Ady7Ipq2vh7wcq7Q5FKeVjHtlwjIiwEG7ygaEVh6KJfwQ+PiuNKQnjeOg9bdqplPo/9c0d/Hl7OVefO5mUuEi7wxmSJv4RCAsN4fbzs9hyvI7d5WfsDkcp5SOe2lJGW2cPdyzLtjuUYdHEP0LXL3Ddgv0HPetXSgHtXd089sFxluckkzcxzu5whkUT/wjFR4WzamEGr+6poLze5wYQU0p52Us7T1HV1M7nlvvH2T5o4h+V28+fhuDqa1spFbyMMfz+vaPkT4xjeU6y3eEMmyb+UZicMI4r5kzimS1lNLR22h2OUsomjkPVHKp0svqCbJ++YasvTfyj9LkLsmnu6ObpLTpCl1LB6vfvHmVifBRXzJlsdygjYlniF5EMEXlHRA6IyD4R+ap7fpKIvCUiJe6/iVbFYKVZk8ezbEYyj2w4RkdXj93hKKW8bO/JBj44Usvt52cREeZf59BWRtsFfNMYMxNYDHxZRAqAu4F1xpgcYJ37f7/0uQuyqWxs56Vdp+wORSnlZb9/7yixkWGsWpRpdygjZlniN8ZUGGN2uB83AQeAKcBVwGPu1R4DrrYqBqtdkJNM/sQ41r57hB7txkGpoHGiroVXdlewamEG8VHhdoczYl4ZE0xEsoB5wGYgzRhTAa4vBxFJHeA5q4HVAGlpaTgcDm+EOmIXpHaxdnc7D/x5HeemDl6cTqfTZ4/jrDxnnmXbjuyOtHT73jauOxoYWZkFWhmMlMPh8IvPwVD+uL8djKEg9DQOR9WIn293GVie+EUkFnge+JoxpnG4V76NMWuBtQCFhYWmqKjIshjH4vzuHl4tc7ChLoqvXb900HUdDge+ehxnrXGssWzbec48imOLLdu+ty0Jdd3HMZJjCrQyGKlVRav84nMwmFpnOxvWvc2189O5ZuXcUW3D7jKw9IqEiITjSvpPGmNecM+uFJFJ7uWTgJF/XfqQ8NAQ7lw+ja3H69leWmd3OEopiz22sZS2zh4+v8J/btjqy8pWPQI8BBwwxvy816KXgFvdj28FXrQqBm/59IIMEqPDedBx1O5QlFIWauno4vGNx7mkII0Zqf7RPUN/rDzjPx+4GbhIRHa6p8uAHwOXiEgJcIn7f78WHRHGLUuy+PuBSkoqm+wORyllkWe2nOBMSydfWDHd7lDGxMpWPe8bY8QYM8cYc657+psxptYYc7ExJsf9NyDqR25dmsW48FB+t17P+pUKRJ3dPTz0/jEWZiVx3lS/vP3oH/zrrgMflhQTwWcWZvDizpPaeZtSAeivH57k5JlWvnihf5/tgyZ+j/rc8mxEXLdxK6UCR3eP4cH1RyiYFE9Rbord4YyZJn4PmpwwjmvnpfPM1hNUN7XbHY5SykPe2Heao9XNfOnC6X7VGdtANPF72BeKptPZ3cPDG3SgFqUCgTGG37xzmGnJMXxi9iS7w/EITfweNi05hsvOmcQfN5Zql81KBYD1h6rZd6qRL66YTmiI/5/tgyZ+S3ypaAbO9i4e/+C43aEopcbot+8cYdL4KK6eN8XuUDxGE78FCibHc1F+Kg9vOEZze5fd4SilRmnz0Vq2HK9j9QXZftf18mAC50h8zFcumkF9SydPbCq1OxSl1Cj96u3DJMdGsmqh/3W9PBhN/BaZl5nI8pxkfv/eUVo7uu0ORyk1QttL63n/cA2rL5hGVHio3eF4lCZ+C33lohxqnB06PKNSfuhXb5eQGB3OjYum2h2Kx2nit9DCaUksmpbE79Yfoa1Tz/qV8he7y8/gKK7mzuXZxER6ZdgSr9LEb7G7Ls6hqqmdP207YXcoSqlh+tXbh4mPCuOWJYF3tg+a+C23dPoE5mcm8FvHETp1eEalfN6+Uw28tb+S28+fRpwfDqs4HJr4LSYifO1juVQ0tPFeuTbtVMrXPbCuhLioMD67bJrdoVhGE78XLM9JZn5mAq8c7aS9S+v6lfJV+0418Ma+Su5YNo3x4wLzbB808XvF2bP+ujbDc9vK7Q5HKTWAX/7ddbZ/+/mBe7YPmvi9ZnlOMjMSQvjtO4f1rF8pH7T3ZANv7g/8s33QxO81IsLVMyKoaGjjua3awkcpX/PLdSXEB8HZPmji96pZE0IonJrIb97Rdv1K+ZI95a6WPHcsyw74s33QxO9VIsI3Ls3ldGMbT23Wu3mV8hU/e6uYhOhwbl+WZXcoXmFZ4heRh0WkSkT29pp3rohsEpGdIrJNRBZatX9ftXR6MkuyJ/Bbx2FaOrR5p1J2215ah6O4ms9fMJ34AG2335eVZ/yPAiv7zPsJ8J/GmHOB+9z/B51vXppLjbODxz7QnjuVstvP3jxEcmwEty4NzLt0+2NZ4jfGvAvU9Z0NxLsfjwdOWbV/X1aYlURRXgr/++4Rmtp0lC6l7PLB4Ro+OFLLl4pmEB0ReH3yDMTbR/o14A0R+X+4vnSWDrSiiKwGVgOkpaXhcDi8EqCVnE7nP46jKKkbR3En9z3xDlfNiLA3sF7ynHmWbTuyO9LS7XvbuO5oYGRlFmhlMFIOh+MjnwM7GWP44eY2EiOF9PbjOBze+wVudxl4O/F/Efi6MeZ5EbkeeAj4WH8rGmPWAmsBCgsLTVFRkdeCtIrD4aD3cWxs3MZbh2u5b9VSEmN8I/mvcayxbNt5zjyKY4st2763LQltARjRMQVaGYzUqqJV//Q5sMs7B6s4fGYrP7h6Npcu9m41j91l4O1WPbcCL7gf/wkIuou7vX3z0jyaO7p4cP0Ru0NRKqj09Bjuf/0gUydE8+kFGXaH43WDJn4RyRxsGsX+TgEr3I8vAkpGsY2AkZsWx7Xz0nn0g+NUNLTaHY5SQePl3ac4eLqJb1ySS3ho8LVqH6qq51VcF2Sl1zwDpACpwIDjkYnI00ARkCwi5cD3gM8BvxSRMKANdx1+MPvax3J4addJHlhXwo+unWN3OEoFvI6uHn725iHyJ8Zx5ZzJdodji0ETvzHmnN7/i0gW8G1c9fL/PcRzVw2w6Lzhhxf4MpKiuXHRVP64qZTPLc8mOyXW7pCUCmjPbjtBWV0LD99WSEiIDP2EADSs3zgikiMijwKvAduBAmPMr6wMLJj860UziAwL4WdvHrI7FKUCWmtHNw+sK2FBViIX5qXaHY5thqrjn+2usnke+Dsw2xjzB2OMNj73oOTYSO5cNo1X91Sw68QZu8NRKmA9vOEY1U3tfGtlPiLBebYPQ5/x7wKWAO/haoHzPyLywNnJ8uiCyOcuyGZCTAQ/eu0AxugQjUp5Wq2znQcdR/jYzDQWZCXZHY6thrq4eweui7nKYnFR4dx1cQ7fe2kfjuJqLswP3p+hSlnhV2+7+sf69srgvYHurKES/zNAnDGmuvdMEUkFGi2LKkitWpjJIxuO8ePXDnJBbgqhQXrhSSlPK61t5snNpXx6QQY5aXF2h2O7oap6HgCW9zP/EuB/PB9OcIsIC+HfP55PcWUTz+/QIRqV8pSfvlFMaIhrCFQ1dOJfZox5oe9MY8yTwAXWhBTcLjtnInMzEvj5m4do7dDBWpQaq10nzvDK7gruXJZNWnyU3eH4hKES/2B1DcF3u5sXiAj3XDaT041t/OG9o3aHo5RfM8bww1cPMCEmgs+vyLY7HJ8xVPKu6m+wFBFZAFT3s77ygIXTkvj4rDQeXH+EqqY2u8NRym+9sa+SLcfr+PolucQFySArwzFU4v934DkRWSMiV7qn/wSecy9TFrn7EzPp7O7h53pTl1Kj0tHVw49fO0BOaiyfCcKO2AYzaOI3xmzB1X5fgNvckwCLjDGbrQ4umE1LjuHmxVk8t+0EB09rAyqlRuqPm0o5XtvCf1w+k7Ag7IhtMEP2zmmMqTLGfM8Yc517us8YU+WtAIPZXRfPIC4qnB++qjd1KTUSZ1o6eGBdCctzkinKTbE7HJ8z1NfgX88+EJHnLY5F9ZEQHcFdF+fwXkkN7xTrd61Sw/XLdSU0tXVyz+Uzg7prhoGMpFWPXhK3wc2Lp5KdEsMPXjlAR1eP3eEo5fNKKpt4fGMpqxZmkj8xfugnBKGhEr8Z4LHykoiwEO69vICjNc08vvG43eEo5dOMMXz/lf3ERITyjUv0Zq2BDJX454pIo4g0AXPcjxtFpElE9Iqjl1yYn0pRXgq//HsJNc52u8NRyme9fbCK90pq+NrHcpkQG2l3OD5rqFY9ocaYeGNMnDEmzP347P/6G8qLvnt5Aa2d3dpnv1ID6Ojq4b9e2c/0lBhuXuLdwdP9jbZx8hMzUmO5dWkWz2wtY+/JBrvDUcrnPLLhGMdrW7j3ioKgHEd3JLR0/MhdF+cwISaC+17cS0+PXnJR6qzTDW08sK6Ei/NTKQrikbWGy7LELyIPi0iViOztM/8rIlIsIvtE5CdW7T8QjR8XzrdX5rOj7AwvfHjS7nCU8hn//bcDdPYYvnflLLtD8QtWnvE/CqzsPUNELgSuAuYYY2YB/8/C/Qek6+anMz8zgR+/doCGVh0BU6lNR2t5adcpvrBiOpkTou0Oxy9YlviNMe8CdX1mfxH4sTGm3b2O3pU0QiEhwvevmk1tcwe/+Lte6FXBrbO7h++9uI8pCeP44orpdofjN7xdx58LLBeRzSKy3t3Lpxqh2VPGc+OiTB7fWMqBCm1Vq4LX4xtLKa5s4r4rCxgXEWp3OH5jqKEXrdhfIrAYWICr589s009HNCKyGlgNkJaWhsPh8GaclnA6nR47jsXRhhdDDV95bAP/sSiKEA/dlp7ntG480sjuSEu3723jul3VCiM5pkArg5FyOBwe+xzUt/Xw0/daOSc5lIiqAziqD449QC/xZC4YDW8n/nLgBXei3yIiPUAy/fTtb4xZC6wFKCwsNEVFRd6M0xIOhwNPHkdr0gn+/c+7qYqZzmcWZnpkm2scazyynf7kOfMoji22bPvetiS0BWBExxRoZTBSq4pWeexz8OUnd2CknV/fvpypE2LGHpwXeToXjJS3q3r+ClwEICK5QLzeRR4AABXfSURBVARQ4+UYAsYnz0tn4bQkfvz6QeqaO+wORymvWX+omlf3VPCvF87wu6TvC6xszvk0sBHIE5FyEbkDeBjIdjfxfAa4tb9qHjU8IsIPrp6Ns62LH/3tgN3hKOUVbZ3d3PfiXrJTYlitwymOimVVPcaYVQMsusmqfQaj3LQ47lyeze/WH+GT56WzKHuC3SEpZanfvnOY0toWnrpzEZFhekF3NPTO3QBw18UzSE8cx3f+sof2rm67w1HKMocqm3hw/RGumTeFpTOS7Q7Hb2niDwDREWH88JpzOFrdzG/eOWJ3OEpZoqfH8J0X9hAbGcZ3L59pdzh+TRN/gFiRm8LV507mQcdhSiqb7A5HKY97cnMp20vr+e7lBdrl8hhp4g8g915RQExkGHe/sEc7cVMB5XRDG/e/XsyyGclcO3+K3eH4PU38AWRCbCTfvbyA7aX1PLm51O5wlPIIYwz3vriXrp4efnjNbB1D1wM08QeY6+ZPYXlOMj9+7SDl9S12h6PUmL28u4K39lfy9Y/lapt9D9HEH2BEhB9dew4A33lhD3qbhPJntc521ry0j7kZCdy5XNvse4om/gCUnhjN3Z/I572SGv60rdzucJQate+9tA9nWxc//eQcQkO0isdTNPEHqBsXTWXRtCT+69X9nG5oszscpUbs9b2neWV3BXddPIPctDi7wwkomvgDVEiIcP91c+js7uHuF3ZrlY/yK3XNHXz3r3spmBTP57WffY/TxB/AspJjuHtlPo7iap7desLucJQaFmMM3/3rHhpaO/jZ9XN14HQLaIkGuFuWZLEkewL/9cp+TtRpKx/l+17adYq/7TnN1y/JZeakeLvDCUia+ANcSIjwk0/OQUT49z/v0hu7lE+rbGzjvhf3MS8zgdXaiscymviDQEZSNPdeMZNNR+t45IPjdoejVL+MMXz7+d20d3Xzs0/NJUyreCyjJRskri/M4OL8VO5//SDFp7UvH+V7nthUiqO4mrtX5pOdEmt3OAFNE3+QEBHu/+Qc4qPC+OozH2r3zcqnHK5q4gevHmBFbgq3Ls2yO5yAp4k/iCTHRnL/dXM4eLqJn74evOO+Kt/S0dXDV5/ZSUxkGD91X49S1tLEH2QunpnGTYsz+cP7x3i/RIc7Vvb7+VuH2HeqkR9few6p8VF2hxMUNPEHoXsuK2B6SgzfeG4ntc52u8NRQez9khr+990jrFqYwaWzJtodTtDQxB+ExkWE8sCqeZxp7eSbf9ImnsoeNc52vv7cTqanxHLvFQV2hxNULEv8IvKwiFSJyN5+lv2biBgR0UEzbTJr8ni+e/lMHMXVPLzhmN3hqCDT02P45nO7aGjt5Nc3zCM6IszukIKKlWf8jwIr+84UkQzgEqDMwn2rYbh58VQuLUjj/tcPsrv8jN3hqCDyh/ePsv5QNfdeUUD+RL0719ssS/zGmHeBun4W/Q/wLUDrF2wm4rqrNyU2ki8/tYOG1k67Q1JB4HB9Nz95vZiVsyZy06JMu8MJSl79fSUi/wKcNMbsGqrJloisBlYDpKWl4XA4rA/QYk6n0yeP47P58KMtrdz+4Do+mZOHVa3pIrsjyXPmWbNxG4zrjgYY0TEFWhmM1MtvvsOvP2wlMTKEKyc2sn79ertDsoXducBriV9EooF7gEuHs74xZi2wFqCwsNAUFRVZF5yXOBwOfPE4igCSj/KDVw8QkniG2dMqLdlPnjOP4tjAuX9gSair07uRHFOglcFIGAPHjk3H2dnCX1afzznp4+0OyTZ25wJvtuqZDkwDdonIcSAd2CEi2obLB9yxbBorZ01k+6EMKuv1dnnlebuPTuLdQ9XcMDMiqJO+L/Ba4jfG7DHGpBpjsowxWUA5MN8Yc9pbMaiBiQg/+dQcYse149g5nZZ2bWWhPOdkTTw7S6Zw1bmTuTBD31t2s7I559PARiBPRMpF5A6r9qU8Iz4qnAvnHaazKxTHzhl09+it82rsmloiWL9rOgmxrfzo2nO0SwYfYGWrnlXGmEnGmHBjTLox5qE+y7OMMdpngI9Jimtl6ezjVNXHsfVght3hKD/X1R3C2x/mgIEL5x/W9vo+Qu/cVf8ke1Ids7JOc7AsjcMnJ9gdjvJTxsAHe7OobxrHBXOPEh+t3YP4Ck38ql/n5Z5gUlIjH+zLovpMjN3hKD+07/hEjlZMYF7OSdJTGuwOR/WiiV/1KyQEVpx7mJjIDt7+cAbNbeF2h6T8yImq8WwrTidrYh1zsivsDkf1oYlfDSgqopuL5pfQ1RXK2zty6OzSt4saWn3TONbvms6E+BaWnXPMshsC1ejpJ1kNKjGujQvmHqG2MZr390zDaEcbahBtHWGs2zGD8DDXSUNYaI/dIal+aOJXQ8pIbaAw7wSllUlsP5RudzjKR3V1C+t25NDaHsFF8w4TE6V9P/kqbVulhmVWViVNLVHsPTaJuOh28jKq7Q5J+RBj4L3d2VSfiaHo3COkJDTbHZIahCZ+NSwisGhmKc7WCDbtn0pMVDvpKY12h6V8xPZD6ZRWJlGYV0bWxHq7w1FD0KoeNWwhIVB07hESY1tw7JxBTUO03SEpH3CgNJW9xyaRn1nJrCxrOvhTnqWJX41IeFgPHzuvhMjwLv6+PZfG5ki7Q1I2OlaRyOYDmWSk1rMwv0xb8PgJTfxqxKKjOrl0QTHGwJvbcmnRNv5B6VRNPO/tziY10cmKuUcI0WziN/SlUqMyPqadj51XQltHOG9tz6W9M9TukJQX1TRE8/aHM4iPaePi+SWEhWo7X3+iiV+NWkpCMxfOO0yDM4q/b8vVG7yCRH3TON7alkdkeBeXFB4iMrzb7pDUCOknVY3JlORGVsw9Sk1jDOt25NDVrZW8gayxOZI3t+USGtLDxxcWa1t9P6WJX43Z1In1LDvnGKfr4rQf/wDmbI3gja159PQIly4o1t42/ZgmfuUR0yfXsqSglPLqBBwfTtfkH2CcrRG8viWfzq5QLl1QTEJsm90hqTHQxK88Ji+zmsUFxzlRnajJP4A0tbiSfkenK+lPiG+1OyQ1Rpr4lUfla/IPKE0tEbyx9f+SfvL4FrtDUh6giV95XO/kv267dufsrxqcUby2ZaYm/QCkn0hlifzMas6ffZSK2nje2pZLh7bz9yu1jeN4bUs+PT3CyoUHNekHGMsSv4g8LCJVIrK317yfishBEdktIn8RkQSr9q/sl5Ney4pzj1DTEMMbW/Nwdmq1jz+oqo/l9S35hIQYPrHwIElapx9wrDzjfxRY2WfeW8BsY8wc4BDwHQv3r3xA1sR6LppfwhlnFL/bFU1ji/bt48vKqhJ4Y2seURFdXLboAOO19U5AsizxG2PeBer6zHvTGNPl/ncToKN6BIH0lEY+vrCY1i7hb5tmaq+ePqq4LIV3dswgMa6FyxYdIHZch90hKYvY2R//Z4FnB1ooIquB1QBpaWk4HA4vhWUdp9Pp88eR58yzZrthMGd2F/97IIw3NxewKr+V/KSuoZ/ow8Z1u77ARlJmkd2RlpXxaPUYeKs0ko3lkeQldnJDfg8RndPBgptyHQ6HX3wOrGZ3GdiS+EXkHqALeHKgdYwxa4G1AIWFhaaoqMg7wVnI4XDg68exxrHGsm3nkcclS4+wbnsuj++PpjD/BAVTK/22K98loa4LnsWxxcN+Tp4zb0TrW62rO4T3dk+jtDKenPRqFhUc55iFFcCrilb5xefAanaXgdcTv4jcClwBXGyMDt0dbKIju1i58CDv75nG1oOZNDijWFRQRmiIvhW8raUtnHU7cqhtjKYwr4xZWf77JaxGxquJX0RWAt8GVhhjtH1YkAoP66Ho3CN8WNLG7qOTOdM8jqK5R4jWDr+8prI+FsfO6XR2hXLx/BIyUhvsDkl5kZXNOZ8GNgJ5IlIuIncAvwbigLdEZKeI/M6q/SvfJgLzc09ywdwj1DVG8/LGAqrqY+0OK+AZ4xoq8fUteYSF9nD54gOa9IOQZWf8xphV/cx+yKr9Kf+UPamOxNhW3v5wBq9tyaMwr9yv6/19WWdXCBv3ZXG0YgIZKfUsm3NM+9IPUna26lEKgMS4Vq5Ysp/3d7vq/U/XxnH+OceIitCk5Cm1jdGs3zmdppZI5s0oZ870Cv1yDWLaZYPyCZHh3Vw0/zAL80s5WTOelzbMpqI2zu6w/J4xsP94Gq9unElXdwgfX1jM3Bma9IOdnvErnyECBVlVpCU6Wb9rOm9szadg6mnm55brmK6j4GyNYMOeaVTUxZOecoZl5xwjKsK/751QnqGJX/mcCeNbuHLpPrYfSmd/6URO1oxn2TnHSElotjs0v2AMHD6ZzJYDmRhg6axj5KTX6Fm++gdN/MonhYf1sLigjMzUM7y/dxqvbprJzMwq5uWWExHWY3d4PquhOZKN+7I4XRdPWmIjy845Rly0dr2gPkoTv/Jpk5MbuXrZHnYcSudAWSqllYksmllGZlq9nsH20tUt7Ds2kV1HJxMa0sOSguPkZlRrGal+aeJXPi/CffafPbmWjXuzeGfnDCYlNbJwZhmJccHdZbAxUFaZyNbiDJytkWRNrGNhfpneDKcGpYlf+Y3UhGauXLqP4hOp7Dw8hZc2zCIno5pzp58KykRX0xDDtuJ0TtfFkxDbwqULDjJ5QpPdYSk/oIlf+ZWQEJg5tYrsSbXsPDyFgydSOHIymZlTK5k9rSIo2v7XN0Xx4eF0yioTiQzvZHHBcXLTqwnRxtlqmDTxK78UGdHNooIyCrIq2Xl4MnuPTeRgWSp5GVXMyqoMyF8ANQ3R7Dk6idLKRMLDupk3o5yCrErC9WK3GiFN/MqvxUW3s3zOMWZPO83uo5PYf3wiB0rTmDGlhplTq/z+GoAxcLJmPPuPp3GqdjzhYV3Mya6gIOt0UPy6UdbQxK8CQmJcKyvmHmV+zkn2HpvI4ZPJHCpPJS2xkfzMajJS6/3qJrDW9jCOnJpAcVkqTa1RjIvs4LzcE+RlVBOh/euoMdLErwJKXHQ7S2aVMj+nnJKTKRwsS2H9rumEh3UxbWId06fUkprg9Mlmjl3dQnl1AkdOTqC8ZjzGhJCW2MT83JNkptXrmAXKYzTxq4AUGdHN7GmnKcg6zenaeA6fmsDRigkcKk9lXGQHGalnyEw9w8SkJsJC7asjb20P42TNeMqqEjlVE09XdyjRkR3Myqpk+uQaEuN0sHPleZr4VUALEddNYJOTG+nsKuVEVQJlVYkcPTWBQydSCZEeUhOdTJrQSEqCk+T4FsuqUoyBM+1CqTORyvpYKmrjqXe6xu2NjuxgxpQaMtPqmZjURIgP/iJRgUMTvwoa4WE9ZE+uI3tyHV3dwum6OCpq46mojefDkvR/rBcf00pibCvjY9oYH9tGbFQ70VGdjIvsGPI6gTHQ2RVKS3s4LW0RNLVE0tAcRUNzFHVN0bS2RwBxhIT0uKpxJpUzObmRCfHNPln9pAKTJn4VlMJCDekpjaSnNALQ3hFKTWMMNQ2uqb5pHGVViRjz0WwcGtJNeFgP4WHdXN4QA8Bf3ptNd4/Q2RVKZ1coPSakz766iY9pY/KERgqiYuhJPUZiXItfXWxWgUUTv1K4rglMSW5kSnLjP+Z19whNLZE0t0XQ0hZOS3sEHZ2hdHa7Eny4+9pAQmwroaE9RIR1ExbaQ1REF9GRHURHdRI7rp2YqI5/nM3nOfMojtVeRpW9NPErNYDQEENCbBsJsf1fYE2Idd0jcOG8I94MS6kxs3Kw9YdFpEpE9vaalyQib4lIiftvolX7V0op1T8re/d4FFjZZ97dwDpjTA6wzv2/UkopL7Is8Rtj3gXq+sy+CnjM/fgx4Gqr9q+UUqp/3q7jTzPGVAAYYypEJHWgFUVkNbAaIC0tDYfD4Z0ILeR0On3+OPKceZZtO7I70tLte9u4blcb/JEcU6CVwUg5HA6/+BxYze4y8NmLu8aYtcBagMLCQlNUVGRvQB7gcDjw9eNY41hj2bZdLVqKLdu+ty0JbQEY0TEFWhmM1KqiVX7xObCa3WXg7R68K0VkEoD7b5WX96+UUkHP24n/JeBW9+NbgRe9vH+llAp6VjbnfBrYCOSJSLmI3AH8GLhEREqAS9z/K6WU8iIxxvdvGxeRaqDU7jg8IBmosTsIGwX78YOWAWgZgPfKYKoxJqXvTL9I/IFCRLYZYwrtjsMuwX78oGUAWgZgfxno8MxKKRVkNPErpVSQ0cTvXWvtDsBmwX78oGUAWgZgcxloHb9SSgUZPeNXSqkgo4lfKaWCjCZ+DxCRlSJSLCKHReSfupoWkakisk5EdouIQ0TSey271T0+QYmI3Nr3uf5ijGXQLSI73dNL3o3cM/obf6LPchGRB9zls1tE5vdaFijvgbGUgd+/B2BYZZAvIhtFpF1E/q3PskE/Qx5ljNFpDBMQChwBsoEIYBdQ0GedPwG3uh9fBPzR/TgJOOr+m+h+nGj3MXmzDNz/O+0+Bg+UwQXAfGDvAMsvA14DBFgMbA6k98BYyiBQ3gPDLINUYAHwQ+Dfes0f8jPkyUnP+MduIXDYGHPUGNMBPINr3IHeCnANPAPwTq/lHwfeMsbUGWPqgbf458Fr/MFYyiAgmP7Hn+jtKuBx47IJSHB3VBgo74GxlEHAGKoMjDFVxpitQGefRcP5DHmMJv6xmwKc6PV/uXteb7uA69yPrwHiRGTCMJ/rD8ZSBgBRIrJNRDaJSKAOzjNQGQXKe2A4BjvWYHgPDMar7wNN/GMn/czr20b234AVIvIhsAI4CXQN87n+YCxlAJBpXLev3wD8QkSmWxapfQYqo0B5DwzHYMcaDO+BwXj1faCJf+zKgYxe/6cDp3qvYIw5ZYy51hgzD7jHPa9hOM/1E2MpA4wxp9x/jwIOYJ4XYva2gcooUN4DwzHgsQbJe2AwXn0faOIfu61AjohME5EI4DO4xh34BxFJFpGzZf0d4GH34zeAS0UkUUQSgUvd8/zNqMvAfeyRZ9cBzgf2ey1y73kJuMXdsmUx0GBcw5AGyntgOPotgyB6DwxmyM+QR9l9FTwQJlytFQ7huip/j3ve94F/cT/+JFDiXucPQGSv534WOOyebrf7WLxdBsBSYA+uawB7gDvsPpZRHv/TQAWui3blwB3AF4AvuJcL8Bt3+ewBCgPwPTCqMgiU98Awy2Cie34jcMb9ON697J8+Q1ZN2mWDUkoFGa3qUUqpIKOJXymlgowmfqWUCjKa+JVSKsho4ldKqSCjiV8ppYKMJn6llAoymviVGiUR+auIbBeRfSKy2u54lBouvYFLqVESkSRjTJ2IjMN1y/0KY0yt3XEpNZQwuwNQyo/dJSLXuB9nADmAJn7l8zTxKzUKIlIEfAxYYoxpEREHEGVrUEoNk9bxKzU644F6d9LPxzWUoFJ+QRO/UqPzOhAmIruB/wI22RyPUsOmF3eVUirI6Bm/UkoFGU38SikVZDTxK6VUkNHEr5RSQUYTv1JKBRlN/EopFWQ08SulVJD5/5o0F3yqNepHAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"m.draw_profile('a');"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3xU9Z34/9c79yshEAgQAiEBwv1mABGVaL3gvdq1iu2qrX6p3fW73e52e7OtrG73Z9vtt669rGW90LqK1a1aW1G8DnhDBCHcwyVACAEChEsmIff374+Z2DFMbiRnzszk/Xw85pGZcz7nzPtzZua8c87nc85HVBVjjDGmvRi3AzDGGBOeLEEYY4wJyhKEMcaYoCxBGGOMCcoShDHGmKAsQRhjjAnKEoQxxpigLEEYV4jPT0TkuP/xUxGRTsreJyLlInJaRJ4VkQEB83NE5E8iUi0iFSJyT7vll4pIqYi0isidfRD7bSKyX0RqReQlERnUSdkZIrJeROr8f2cEzFsiIk0i4g145AfMv1REPvHXuUxEFrdb9/8Vkb3++etE5MLe1i1I/INE5EV/XfeLyG2dlO30MxUR9a+nra6PBcxLFJFHReSI/3P8s4jkBMzPE5EVInJCRA6LyK9EJK6v62vaUVV72CPkD+BrQCkwEsgBtgH3dFD2DmAHkAukAX8Cfhcw/x3gYSAemA5UA5cEzP974HPAOuDOXsY9GagBLvbH8gzwbAdlE4D9wDeBROAf/K8T/POXAP/TwbLxwCn/dhJgNuAFpvvnzwVqgfP8878OHAViu1GHO4Fl3azvcuAP/rpe6I9p8rl8poACYztY9ttACZANJAFPAS8EzF8BLPPPGwZsBv7B7e9xtD/sCMKcRUS+KyJ7RKRGRLaJyI0OvM0dwM9VtUJVDwI/x7fjCuY64HFVPaCqXuAnwC0ikiIiaUAx8GNVbVLVEuB/ga+2Layqv1bVt4D6Poj7S8CfVXW1P5YfAjeJSHqQssVAHPCwqjao6iP4duaXduN9BgEDgKfU52NgOzDJPz8P2Kqq69W3B/09kAUMPfeqfZaIpAJfAH6oql5VfQ94GfjbDhbpyWfa3hhgpaoeUdV64Fl8yThw/nOqWq+qh4HX2s03DrAEYYLZA1wEZAD/CvyPiAwPVtB/uuVkJ49RHbzHZHz/MbYpoeMfvPgfga8TgXEB09vPn9JJ/XrjM3Gr6h6gERjfQdlN/h14m018tp7X+U+pbBWRrwes9wi+/96/IiKxIjIPGA285y/yKhArInNFJBZfQtwIHO51Df9qPNCiqjsDpnX2OXXnM13tP0X0gojkBUx/HJgvIiNEJAVfIn41YP5/Arf6/ynIAa7ClySMgyxBmLOo6vOqWqmqrar6B2AXMKeDss+o6sBOHuUdvE0avtMVbU4BaR20Q7wK3O0/D50BfMc/PUVVa4D3gR+KSJKIzML3X29Kz2veLe3jbos92BFEV2WfAyYCQ4D/A/xIRBYFlF0O/AhoAN4F7lPVA/55NcAf8SWMBuB+YHG7ZNRbPalrsPLtP9MF+I58JgCVwF8C2hF2AuXAQeA0vu3yQMC6VuFLNqeBCnynC1/qcY1Mj1iCMGcRkdtFZGPbUQC+/8az+vhtvPhOobQZAHg72ME9gW9n6QG24mtzAN+OAnz/bY4BDgD/BTwdMK/bRGRUYINxN+Nui72mp2VVdZs/Ebeo6gf4/kv+G38sE/Cd+78dX1vGZODbInKNfz134ztqmOyf/2V8O9wRHdTtNwGf52+AwCO/TX1Q12DlP/OZ+k/LNarqSeAb+D6zif6y/4WvfWEwkAq8gP8IQkRigJX+aan4vouZ+E41GgdZgjCfISKjgf8G7gUGq+pAYAufPYUTWP5L7XrhtH90dIppK74G5TbT/dPO4j+SuV9V81R1pL/cQf8DVd2vqteq6hBVnYtvJ7O2p3VX1XJVTWt7dCduf6+jRHz/AQcrO63dUdE0OqgnvkbctrJTgFJVXemvfynwCr5TK/hj+LOq7vTPfw04BFzQQd3+ru2oDvg7IPDIb1oH8ewE4kRkXMC0Dj8nevCZBqnvdHwN59Wq2gD8EpgjIln42mNygV/523KOA08CV3eybtMHLEGY9lLx/XCPAojIV+jkfL6qPh24Uw3y6OgU0++BfxJfF9URwD/j66VyFn9XywJ/N8pJwP8DHlDVVv/8iSKSLiIJIvJl4Ap/mbblE0QkCd/OKN5/Kupcv/tP42s3uMjfiPsAvt42wf6r9gAtwD/4u3He65/+tj+uG0Qk01+vOfh6Of3JX2YDME58XV1FRAqAa/nrOf6PgWtEJN8//3J8bQZbzrFeZ1HVWnz/tT8gIqkiMh+4AV8Po2A6/ExFZLL4uvzG+jsW/Bxfgt8eUJ/bRSRDROLxJbFKVT2mqseAvcDXRSRORAbiaxAvwTjLyS5S9ojMB/BjfF1Fj+Hb0a4C7u7j9xDgp/73qfY/l4D5XuAi//Px+LpP1uHrJvpP7db1j/gSWi2+c/JF7eZ78CW9wEdxL2K/Dd/58lp8O/RBAfNeBb4f8HomsB44A3wCzAyYtxw47q/rDtp12wS+iG+HX4PvlNlPgJiA7feAP44afDvav+1m/HfS/W6ug/Cd66/1v9dtAfMuwncKqcvPFF/PrVL/eqr86xwXsOxgfMm3Cjjp/xznBMyf4f8cT/i/l88DQ93+rUT7o+3DM8YYYz7DTjEZY4wJyhKEMcaYoCxBGGOMCcoShDHGmKCi6m6IWVlZmpeX53YYvVJbW0tqaqrbYXSpsqbSsXUntibSENPg2PpDLevAcQCO5Q7u9jLRtg16akT6iIj5LTglVPVfv379MVUdEmxeVCWIvLw81q1b53YYveLxeCguLnY7jC4t8SxxbN2F3kJK00odW3+o3fmPywBY9vCd3V4m2rZBTy0pXhIxvwWnhKr+IrK/o3mOnWISkSdEpEpEzrpwR0S+Jb57wwe9fYOItPhv9bBRRF52KkZjjDEdc/IIYhnwK3xXV35KRHKBy/FddNORM6o6o5P5xhhjHObYEYSqrsZ3NWV7v8A3OIhdoWeMMWEspG0QInI9cFBVS4Lf1flTSSKyDmgGHlLVDm/rK75hGBcDZGdn4/F4+jDi0PN6vRFRh0JvoWPrTmxJdHT9oZbc4rvzeE/qFG3boKc8Hk/E/BacEg71D1mC8A8Cch++G6l1ZZSqVvrvlPm2iGxW38AsZ1HVpcBSgKKiIo30Rq1IaZizRurumxdbB9CjOkXbNuipRcWLIua34JRwqH8or4MowHf/9xIR2Ydv3NpPRGRY+4KqWun/W4bvBl0zQxemMcYYCGGCUNXNqjpUfff0z8N3d8pZ6htf9lP+2x8n+p9nAfPxDX5ujDEmhJzs5roc+BAoFJEKEbmrk7JFIvKY/+VEYJ2IlOAbOewhVXUsQdQ3tfDfq8v4YM8xp97CGGMikmNtEKq6qIv5eQHP1+EbQhH1Db041am42ouNEf773TKm5mRwQUFfj6ppjDHOeuK9vXy8r5qHb51BYlxsn66739+LKT42hi8W5fJOaRWVJ8+4HY4xxnSbqvLUmv0crWno8+QAliAAuGV2Lgo8t+6A26EYY0y3rSmrZu+xWhbN6Wjo996xBAHkDkrhonFD+MPHB2hptev3jDGR4Zm15QxIiuOaacMdWb8lCL/b5uRy6FQ9q3ZWuR2KMcZ06bi3gZVbDnPTrJEkxff96SWwBPGpz03MZkh6Is98ZKeZjDHh74VPDtLY0urY6SWwBPEpX2P1SN7ecYTDp+rdDscYYzqkqixfW855ozMpHJbu2PtYgghw6+xRtKo1VhtjwtuasmrKHGycbmMJIoCvsTrLGquNMWFtub9x+lqHGqfbWIJo57Y5ozh48gyrdx11OxRjjDlLdW0jrzncON3GEkQ7l03KJistkeUfdTaekTHGuOOFTypobGnl1jm5jr+XJYh24mNjuLloJG/tqOLIaWusNsaED1XlmbXlzBo1kAnDBjj+fpYggrh1di4trcpzH1tjtTEmfHy0t5qyo7XcNnd0SN7PEkQQowencuHYLJ61xmpjTBhZvrac9KQ4rpnqbON0G0sQHVjkb6x+1xqrjTFh4ERtI69uPsxNM3NITnC2cbqNJYgOXD4pm8GpCSxfa43Vxhj3/dHfOL1orrPXPgSyBNGBhLgY/qZoJG9ur6LKGquNMS4KdeN0G0cThIg8ISJVIrIlyLxviYj6hxUNtuwdIrLL/7jDyTg7cuvsUbS0Ks+vr3Dj7Y0xBoC1/sZpp6+cbs/pI4hlwML2E0UkF7gcCHr+RkQGAfcDc4E5wP0ikulcmMGNyUrlgoLBLF9bTqs1VhtjXNLWOH3ttBEhfV9HE4Sqrgaqg8z6BfBtoKO97pXAG6paraongDcIkmhCYdGcUVScOMN7u23MamNM6J2obWTFltA2TrdxbEzqjojI9cBBVS0RkY6K5QCBFyFU+KcFW99iYDFAdnY2Ho+n74IFklqV9Hh45JX1tFYm9em6g/F6vX1eBycUegsdW3diS6Kj6w+15JYUoGfbLNq2QU95PJ6I+S04pa3+K/c10djcytiYqpBvj5AmCBFJAe4DruiqaJBpQY82VHUpsBSgqKhIi4uLexNiUIvqt/PEe3uZdN75DE13Nkl4PB6cqENfW+JZ4ti6C72FlKaVOrb+UJsXWwfQozpF2zboqUXFiyLmt+AUj8fDggULeHD9KmaOSuFvr5sf8hhC3YupABgDlIjIPmAk8ImIDGtXrgIIvNHISKAyJBEGcevsXJpblf+1xmpjTAh9vO8Ee1xonG4T0gShqptVdaiq5qlqHr5EMEtVD7cruhK4QkQy/Y3TV/inuSJ/SBrn5w/i2bUHrLHaGBMyz3y0n/SkOK4LceN0G6e7uS4HPgQKRaRCRO7qpGyRiDwGoKrVwIPAx/7HA/5prlk0ZxTl1XV8sOe4m2EYY/oJb6OyYsthbnShcbqNo20Qqrqoi/l5Ac/XAXcHvH4CeMKx4HroysnDyEyJZ/naci4cF/TSDWOM6TPvVzbT2NzKrbPdOb0EdiV1tyXFx/KFWSNZufWwXVltjHGUqrLqQBMzcgcyaUTorpxuzxJED3z5/NG0qPK0DSZkjHHQmrJqKmuV20J436VgLEH0QF5WKpcUDuXpj8ppbG51OxxjTJR68v29pMXD9dPdaZxuYwmih+68II9j3gZWbD7kdijGmCh0oLqON7cfoTg33vExp7tiCaKHLhybRf6QVJ78YJ/boRhjotBTa/YjIlySG/IbXZzFEkQPxcQId16QR8mBk2woP+F2OMaYKFLX2Myza8tZOGUYg5Pd3z27H0EEumnWSNIT41hmRxHGmD704oaDnK5v5isX5LkdCmAJ4pykJcZxc1Eur2w6xBHr8mqM6QOqyrL39zElZwDnjQ756AZBWYI4R7fPsy6vxpi+88Ge4+yq8nLnBWPo5E7XIWUJ4hzlZaVyaeFQnvloPw3NLW6HY4yJcE++v5fBqQlcO22426F8yhJEL9w5P49j3kZe2WRdXo0x527/8Vre2lHFbXNHud61NZAliF64cGwWY4em8eT7+1C1u7waY87N7z/cT6wIX5o72u1QPsMSRC+ICHfMG83mg6f4pPyk2+EYYyJQbUMzz318gKumDmdYhvOjVvaEJYhesi6vxpjeeOGTCmoamvnK/Dy3QzmLJYheSk2M44uzc3l1s3V5Ncb0TGursuyDfUwfmcHM3IFuh3MWxxKEiDwhIlUisiVg2oMisklENorI6yIS9E5UItLiL7NRRF52Ksa+8mmX1zX73Q7FGBNB3tt9jD1Ha7lzfl7YdG0N5OQRxDJgYbtpP1PVaao6A/gL8KMOlj2jqjP8j+sdjLFPjB6cyucm+O7yal1ejTHd9eT7e8lKS+TqqeHTtTWQYwlCVVcD1e2mnQ54mQpETdefOy8Yw/HaRv5SYl1ejTFd23uslndKj/KluaNIjAufrq2BQn67QBH5MXA7cAq4pINiSSKyDmgGHlLVlzpZ32JgMUB2djYej6dvA+4mVWVEqvDLlZsZdHrXOR8uer1e1+rQE4XeQsfWndiS6Oj6Qy25JQXo2TaLtm3QUx6PJ2J+C+fq6e0NxAqMaanA46k8a3441D/kCUJV7wPuE5HvAfcC9wcpNkpVK0UkH3hbRDar6p4O1rcUWApQVFSkxcXFDkXetb9L3s8PXtrCgPzpnDd60Dmtw+Px4GYdumuJZ4lj6y70FlKaVurY+kNtXmwdQI/qFG3boKcWFS+KmN/Cuaipb+Led97muukj+PzCmUHLhEP93ezF9AzwhWAzVLXS/7cM8ADBt2CYuWlWDulJcTz5/j63QzHGhLE/rq/A29DMnfPHuB1Kp0KaIERkXMDL64EdQcpkikii/3kWMB/YFpoIeyclIY5binJ5dcthDp0643Y4xpgw1Nqq/O7D/cwcNZAZYdi1NZCT3VyXAx8ChSJSISJ3AQ+JyBYR2QRcAXzDX7ZIRB7zLzoRWCciJcA7+NogIiJBANw+L49WVZ5eY3d5NcacbdWuo+w9VsudYTLmQ2cca4NQ1UVBJj/eQdl1wN3+5x8AU52Ky2mjBqfwuQnZPLO2nHsvHRtWN94yxrhv2fv7GJqeyFVTwrNrayC7ktoBX5mfR3VtIy+XnN0zwRjTf+2u8rJq51G+fP5oEuLCf/cb/hFGoAsKBlOYnc5j75bR2ho1l3oYY3pp2Qd7SYiNYdGcUW6H0i2WIBwgInxtQT47j3h5p7TK7XCMMWHgmLeB59dVcOPMHIakJ7odTrdYgnDIddNHkDMwmUdXBb18wxjTzyx7fx+NLa0sXpDvdijdZgnCIfGxMdx14Rg+3neC9furu17AGBO1vA3N/P7DfVwxKZuCIWluh9NtliAcdOucXAamxPNfnjK3QzHGuOjZteWcrm/mngUFbofSI5YgHJSSEMft8/J4c/sRdh2pcTscY4wLGptbeezdvZyfP4iZozLdDqdHLEE47M4L8kiKj+G3q+0owpj+6KWNBzl8uj7ijh7AEoTjBqUmcEtRLn/aeNBuv2FMP9Paqvx21R4mDh/AgvFD3A6nxyxBhMDdF+XTqvD4u3vdDsUYE0JvbD/CnqO13LMgPyxHjOuKJYgQyB2UwnXThrN8bTmn6prcDscYEwKqyqOr9pA7KJlrwnTEuK5YggiRry0ooLaxhafW7HM7FGNMCKzdW82G8pMsviifuNjI3NVGZtQRaOLwARQXDuHJ9/dR32TjVhsT7R5dtYfBqQncXJTrdijnzBJECN2zoIDjtY08v77C7VCMMQ7afug075Qe9fdijNw7OluCCKG5YwYxI3cg/726jOaWVrfDMcY45Ler9pCaEMvt8/LcDqVXLEGEkIhwz4ICyqvreHXLYbfDMcY44EB1HX/edIhFc0aRkRLvdji94miCEJEnRKRKRLYETHtQRDaJyEYReV1ERnSw7B0issv/uMPJOEPpiknZ5A9J5dFVe1C1W4EbE20ee7eMGIG7Lgrv8aa7w+kjiGXAwnbTfqaq01R1BvAX4EftFxKRQcD9wFxgDnC/iETWNeodiIkRvnZxPlsrT/Pe7mNuh2OM6UPHvQ38Yd0BPj8jh+EZyW6H02uOJghVXQ1Ut5t2OuBlKhDs3+grgTdUtVpVTwBvcHaiiVifn5lD9oBEuxW4MVHmdx/up76pla9F0C29O+PYmNSdEZEfA7cDp4BLghTJAQ4EvK7wTwu2rsXAYoDs7Gw8Hk+fxuqUBcNbea70OMv+9BZ5GX/t5eD1eiOiDoXeQsfWndiS6Oj6Qy25JQXo2TaLtm3QUx6PJ2J+C23qm5XHV9cxc2gsFdvWU7Gtd+sLh/q7kiBU9T7gPhH5HnAvvtNJgYJdkx70hL2qLgWWAhQVFWlxcXEfRuqc885v4tWH3uZjbyZ33jDr0+kej4dIqMMSzxLH1l3oLaQ0rdSx9YfavNg6gB7VKdq2QU8tKl4UMb+FNo+/t5fapm388G/mMqsP7toaDvV3uxfTM8AXgkyvAAKvLhkJVIYkohBJT4rny+eP5tUth9h3rNbtcIwxvdDU0srj75YxZ8ygPkkO4SLkCUJExgW8vB7YEaTYSuAKEcn0N05f4Z8WVb4yP4+42BiWvmu3Ajcmkr28sZLKU/V8PQJv6d0Zp7u5Lgc+BApFpEJE7gIeEpEtIrIJ347/G/6yRSLyGICqVgMPAh/7Hw/4p0WVoelJfGHWSP53fQWHT9W7HY4x5hy0tvpuyjdhWDrFhZF3S+/OON2LaZGqDlfVeFUdqaqPq+oXVHWKv6vrdap60F92nareHbDsE6o61v940sk43fR3xQW0tir/5dntdijGmHOwcuthdlV5+XpxQUTe0rszbrdB9Hu5g1K4uSiX5WsPUHnSBhQyJpK0tCq/eHMnBUNSuXZa0Gt+I5oliDBw76VjUZRfv2NHEcZEklc2H2LnES//eNl4YmOi6+gBLEGEhZyBydwyO5fn1h3gaJ3dxM+YSNDSqjz85k7GZ6dF7IBAXbEEESb+/pKxCMKfy2zEOWMiwcslByk7Wss3LxtPTBQePYAliLAxPCOZ2+aO4r2DzZQfr3M7HGNMJ5pbWvnPN3cxafgArpw8zO1wHGMJIox8vbiAWIFfvr3L7VCMMZ14ccNB9h2v45uXR+/RA1iCCCvZA5K4NDeOFzYcZK9dXW1MWGpqaeWRt3cxNSeDyyYOdTscR1mCCDNX5ycQHyv88i07ijAmHP1xfQUHqs/wT5ePj7rrHtqzBBFmMhKF2+fl8dLGg+yu8rodjjEmQENzC798ezczcgdG3VXTwViCCENfuzifpPhYHrGjCGPCynPrKjh4sn8cPYAliLA0OC2ROy7I48+bKtl5pMbtcIwxQH1TC79+ezez8zK5aFyW2+GEhCWIMLX4onxS4mP5TzuKMCYsPLu2nMOn6/lmPzl6AEsQYSszNYGvXjiGVzYdYsfh010vYIxxTH1TC7/27OH8/EFcUNA/jh7AEkRYu/vCfNIT43j4DTuKMMZN/7NmP0drGvjmZePdDiWkLEGEsYyUeL564Rhe23qYLQdPuR2OMf1SXWMzj67aw4Vjs5ibP9jtcELKEkSY++qFYxiQFMfDb9pRhDFu+P2H+znmbeSbl4/runCUcSxBiMgTIlIlIlsCpv1MRHaIyCYReVFEBnaw7D4R2SwiG0VknVMxRoKM5Hj+z0X5vLn9CJsqTrodjjH9irehmd+u2sOC8UM4b/Qgt8MJuU4ThIiM6uzRxbqXAQvbTXsDmKKq04CdwPc6Wf4SVZ2hqkVdVSLa3Tk/j4Ep8XYUYUyI/e6DfZyoa+KfLu9fbQ9t4rqY/wqgQGCfLgWGAEOB2I4WVNXVIpLXbtrrAS/XAH/Tg1j7rfSkeBZfnM9PXytlQ/kJZo7KdDskY6Le6fomlq4u47KJQ5meG/RkR9TrNEGo6tTA1/4d/neAy4B/7+V7fxX4Q0dvDbwuIgr8VlWXdrQSEVkMLAbIzs7G4/H0Mix3eb3eoHXIb1bSE+C7y9fw3TlJrvfDLvQWOrbuxJZER9cfasktKUDPtlm0bYOe8ng8Hf4WQuXFXY2cOtPERZmnXYnD7fpD10cQAIjIOOA+YC7wc+AfVPWcR7YRkfuAZuDpDorMV9VKERkKvCEiO1R1dbCC/uSxFKCoqEiLi4vPNayw4PF46KgOR1P38cM/baVp6ESucPke9Es8Sxxbd6G3kNK0UsfWH2rzYn3je/SkTtG2DXpqUfGiTn8LTjtyup6vv+XhmqnDueP6Wa7E4Gb923TVBjFFRJYDfwTexNd+8Fgvk8MdwLXAl1RVg5VR1Ur/3yrgRWDOub5fNFk0ZxQFQ1J56NUdNLXY0KTGOOXnr5fS0qp8Z+EEt0NxVVe9mEqAecC7+HbSvxCRR9oePX0zEVmI7xTV9aoadNg0EUkVkfS258AVwJZgZfubuNgYvn/1RMqO1fLMR+Vuh2NMVNpWeZrn11dwxwWjGTU4xe1wXNXVKaa78LUH9Jj/yKMYyBKRCuB+fL2WEvGdNgJYo6r3iMgI4DFVvRrIBl70z48DnlHV184lhmh06YShXFAwmIff3MnnZ+aQkRzvdkjGRA1V5ccrtpGRHM+9l/S/6x7a6ypBPAukq+rRwIn+toFObxCkqouCTH68g7KVwNX+52XA9C7i6rdEhO9fPZHrfvUev3lnN9+7eqLbIRkTNd4preL93ce5/7pJZKTYP19dnWJ6BLgoyPTLgV/0fTimO6bkZHDTzJE8+f4+DlQHPVNnjOmh5pZW/n3FDvKzUvny+aPdDicsdJUgLlTVF9pPVNWngYudCcl0x79cWUhMDPx0Zf/t6WJMX1r+8QF2V3n57lUTiI+1uxBB1wmis872tgVdNCwjicUX5fPnkko2lJ9wOxxjIlpNfRMPv7GTuWMGcfmkbLfDCRtd7eSrROSsLqYiMhs4GqS8CaGvLShgSHoi//bKdjroMWyM6YbfePZwvLaRH1wzyfWLUMNJVwniX4DnRGSJiFznf/wr8Jx/nnFRamIc/3z5eNbvP8GrWw67HY4xEelAdR2Pv7eXm2bmMHVkhtvhhJVOE4SqrsV3/YMAd/ofAsxV1Y+cDs507eaiXAqz03no1R00NLe4HY4xEednK0sR4FtX9t9bm3Sky7u5qmqVqt6vql/wP37kv8LZhIHYGOH710ykvLqOpz7c73Y4xkSUjQdO8nJJJYsvzmfEwGS3wwk7XZ1ieqntiYj80eFYzDlaMH4IF48fwiNv7eJkXaPb4RgTEVSVf/vLNrLSEvnaggK3wwlLPenFlO9kIKZ37rt6It6GZh55a7fboRgTEV7dcph1+0/wz1eMJy2xW/ct7Xe6ShDawXMTZgqHpXPL7FyeWrOPfcdq3Q7HmLDW0NzCQ6/uoDA7nS8W5bodTtjqKkFMF5HTIlIDTPM/Py0iNSLS6a02TOh98/LxxMfG8NCrO9wOxZiw9tSH+ymvruP710wkNsa6tXakq15Msao6QFXTVTXO/7zt9YBQBWm6Z2h6EvcsKOC1rYdZu7fa7XCMCUsnaht55K1dXDx+CAvGD3E7nLBmV0NHmf9zUT7DBiTx41e20dpqZwWNae+Rt3fhbcL0RykAABl3SURBVGjmPrvRZZcsQUSZ5IRYvr2wkJKKU/xh3QG3wzEmrOw6UsNTH+7nltmjKByW7nY4Yc8SRBS6cWYOc8cM4qFXd3DM2+B2OMaEhdZW5b4Xt5CWFMe3rhjvdjgRwbEEISJPiEiViGwJmPYzEdkhIptE5EURGdjBsgtFpFREdovId52KMVqJCD++cSp1jc38+JXtbodjTFj43/UVrN1XzfeumsDgtES3w4kITh5BLAMWtpv2Br5xracBO/GNMPcZIhIL/Bq4CpgELBKRSQ7GGZXGDk3jngUFvLjhIB/sPuZ2OMa46ri3gX9/dTuz8zK5+Tzr1tpdjiUIVV0NVLeb9rqqNvtfrgFGBll0DrBbVctUtRHfqHY3OBVnNPv7S8YyenAKP3hpi92nyfRrP16xndqGZv79xqnEWLfWbnPz8sGvAn8IMj0HCGxdrQDmdrQSEVkMLAbIzs7G4/H0YYih5/V6+7QON49p4T/WNfCdZW/x+bEJfbbeQq9zNzZLbEl0dP2hltziG/i+J3WKtm3QUx6Pp89+C9uPt/DCJ/Vcmx/Pwe3rORghZ137el9wLlxJECJyH9AMPB1sdpBpHfbXVNWlwFKAoqIiLS4u7osQXePxeOjLOhQDpU0bWLHlMN/4/GzGZKX2yXqXeJb0yXqCKfQWUpoWPSPlzYv1DQvbkzpF2zboqUXFi/rkt9DQ3MIDD7/LqEEp/MdXLiYpPrZvAgyBvt4XnIuQ92ISkTuAa4EvafBRbiqAwJOEI4HKUMQWrX54zUQS42L4wUubbWAh06886imj7FgtD35+SkQlh3AR0gQhIguB7wDXq2pdB8U+BsaJyBgRSQBuBV4OVYzRaOiAJL69sJD3dx/n5RLLtaZ/2Husll97dnPttOF2xfQ5crKb63LgQ6BQRCpE5C7gV0A68IaIbBSRR/1lR4jICgB/I/a9wEpgO/Ccqm51Ks7+4ra5o5k+MoMH/7KNU3VNbodjjKNUlR+8tJnE2Bh+dK11gjxXjrVBqOqiIJMf76BsJXB1wOsVwAqHQuuXYmN810Zc/6v3+MnKHfz7jVPdDskYx7y08SDv7z7Og5+fwtABSW6HE7HsSup+ZEpOBl+ZP4ZnPipn/f4TbodjjCNO1jXyb3/ZzozcgXxpzii3w4loliD6mW9ePp7hGUnc9+Jmmlpa3Q7HmD73k9d2cPJMk13z0AcsQfQzaYlx3H/dZHYcruHJ9/e6HY4xfWrdvmqWrz3AV+fnMWmEjUjQW5Yg+qErJ2fzuQlD+cUbuzh48ozb4RjTJ5paWrnvxS2MyEjiHy+zm/H1BUsQ/ZCI8K83TAbg/j9ZBzETHR57dy+lR2r41xumkGpjTPcJSxD91MjMFL5x2Tje3H6E17YcdjscY3rlQHUd//nWTq6cnM3lk7LdDidqWILox+66cAwThw/gvhc327gRJmK1tirf+eMmYkVYcv1kt8OJKpYg+rH42BgevmUGNQ3NfPePm+w2HCYiPfH+Xj7Yc5wfXjuJ4RnJbocTVSxB9HOFw9L59pWFvLm9ij98bEOUmshSeriGn64s5bKJ2dwy28Z56GuWIAxfnT+GCwoG88BftrHvWK3b4RjTLQ3NLXzj2Q0MSIrjoS9MRcSueehrliAMMTHCf9w8nbgY4Z+e20izXUBnIsD/e30nOw7X8NO/mUaWDSHqCEsQBoARA5N58PNT+KT8JI+u2uN2OMZ0ak3ZcZa+W8Ztc0dx6QTrteQUSxDmUzfMyOG66SN4+M1dbKo46XY4xgR1ur6Jf36uhLzBqfzgmoluhxPVLEGYz/i3G6aQlZbIN/+wkTONNo61CT/3/2krh0/X84tbZpCSYBfEOckShPmMjJR4fv7F6ew5WstDr0bI4L2m3/hzSSUvbjjI/710LDNyB7odTtSzBGHOMn9sFl+dP4bffbifVTuPuh2OMQAcOnWG+17czIzcgdx7yVi3w+kXnBxR7gkRqRKRLQHTbhaRrSLSKiJFnSy7T0Q2+0edW+dUjKZj315YyLihafzL8yWcqG10OxzTz7W2Kt96voTmVuUXt8wgLtb+tw0FJ7fyMmBhu2lbgJuA1d1Y/hJVnaGqHSYS45yk+Fh+ccsMTtQ1ct9Lm+0qa+OqJz/Yx/u7fVdLj8lKdTucfsOxBKGqq4HqdtO2q2qpU+9p+taUnAy+efl4Vmw+zIsbDrodjumnSg/X8JPXdnDZxGxutaulQypcuwAo8LqIKPBbVV3aUUERWQwsBsjOzsbj8YQmQod4vd6wqsMEVcZnxnDfCyW0HtlJVrLvf4pCb6Fj75nYkujo+kMtuSUF6Nk2i7Zt0FMejwev18sbb7/DAx/WkxTTyvXDali1apXboYVMOOwLwjVBzFfVShEZCrwhIjv8RyRn8SePpQBFRUVaXFwcwjD7nsfjIdzqMG56HVf957s8X57M03fPJS42hiWeJY69X6G3kNK06DnQnBdbB9CjOkXbNuipRcWL8Hg8fHgmmwM1ZTx+RxGfm9i/LogLh31BWLb0qGql/28V8CIwx92I+rfcQSk8cMNkPtpbzc9e7787LRNaJUebWbq6jEVzRvW75BAuwu4IQkRSgRhVrfE/vwJ4wOWw+r2bZo1kQ/lJfruqjKk5GW6HY6Lc3mO1PFrSwMRhA/jRtZPcDqffcrKb63LgQ6BQRCpE5C4RuVFEKoB5wCsistJfdoSIrPAvmg28JyIlwFrgFVV9zak4Tff98NpJFI3O5F+e38SJGrvvvnFGU3MMi3+/jliB3/7teSQnxLodUr/l2BGEqi7qYNaLQcpWAlf7n5cB052Ky5y7hLgYfvPlWVz3y/d4+5OxXDtvG4kJdjsO03dU4b3NYzhw1Mu3ipLIHZTidkj9Wli2QZjwNTQ9id986Txq6xNYvSmfVrs8wvShTWXD2X9kEN+/eiKTBtuRg9ssQZgeO290JnMnlXPw2EA27spxOxwTJSqOZrBhVw75w49z14Vj3A7HEIaN1CYyFOYe5fipFDaVjWDwgDpGDzvhdkgmgp2uTWRVST6D0uu4YMo+Gx0uTNgRhDlncyeVk5Xh5d3NYzjpTXI7HBOhmppjeHvDOGJEuWTmbuJibUTDcGEJwpyz2Ji//qDf/mQcjU12ztj0TFuj9ClvEgtm7CE9xW4MGU4sQZheSU1q4pKZu6k542u0tnv6mZ5oa5QuKjzAiME1bodj2rEEYXotO9PL3InlVBwdyMbdI9wOx0SIwEbpSXlH3A7HBGEJwvSJwtyjjM05SsmeHMqrbKQv07mzG6XdjsgEYwnC9AkROH/SfgYPqOXdknxrtDYd+mujNNYoHeYsQZg+ExerXDpzN7Gxrbyxbjy19fFuh2TCTEursKqkwNcoPd0apcOdJQjTp1KTG7n8vJ00Nsfy+seF1DfapTbGp1XhvU1jqDg6kPMn72dE1mm3QzJdsARh+tzgjDo+N2sX3jOJvLFuPE3N9jXr71RhzbbR7D08mPPGH6Aw96jbIZlusF+uccSwQV6KZ+ymuiaZtz4ZR3OLtUL2Z5/sHMnOA0OZMuYQU/MPux2O6SZLEMYxuUNPceHUvRyuHsCqkgJarS2yX9pcNozNe4czPreK88ZXuB2O6QFLEMZRBSOqmTtxPweqMnl/yxi7kK6f2Xkgi/U7c8kbdpzzJ+237qwRxloQjeMmjq6isSmWDbtHkhjfzOwJB2xH0Q/sPZTJB1vzyMk6yUXT9hJjn3nEcXJEuSdEpEpEtgRMu1lEtopIq4gUdbLsQhEpFZHdIvJdp2I0oTOt4BCTRh9m2/5hbCob7nY4xmEVRwfw7qZ8sjO9XDJzD7ExdugYiZw8xbQMWNhu2hbgJmB1RwuJSCzwa+AqYBKwSERsUNoIJwKzJxygYMQxNuwayfb9Q90OyTjkyIk03tkwloHpZ/jcrF12IVwEc3LI0dUiktdu2nagq3u9zwF2+4ceRUSeBW4AtjkSqAkZEZg/ZS+NzbF8tH00ifHN5I+odjss04eOn07mzfXjSE3yXQ+TEG9D0kaycGyDyAEOBLyuAOZ2VFhEFgOLAbKzs/F4PI4G5zSv1xsRdSj0Fp7zsmPHwrKGZt7bnM+Y5uFMGNT8mfmJLYm9Wn+4SW7xjavckzpF4jY4eiaG/92UQkos3DOpmYFNBdB0buvyeDwR81twSjjUPxwTRLDDiw5PYKrqUmApQFFRkRYXFzsUVmh4PB4ioQ5LPEt6tfy82TGcWjuBp3Yks2Ba2WdGpCv0FlKaVtrLCMPHvNg6gB7VKdK2wSlvEq9vHU+LtHD5nO0cSW2gN/dnXVS8KGJ+C04Jh/qHYzfXCiA34PVIoNKlWIxD4uNaubyolMED6nhnYwE7yoe4HZI5R1UnU1nx0QRaWmO4vKiUjNQGt0MyfSQcE8THwDgRGSMiCcCtwMsux2QckJTQwpWzS8kdcpI12/L4ZGeOXScRYcqrBrJybSEJ8S1cff52Bg8443ZIpg852c11OfAhUCgiFSJyl4jcKCIVwDzgFRFZ6S87QkRWAKhqM3AvsBLYDjynqluditO4Ky62lUtm7mbcyKNsKhvBB1vyaLEkERFKDwzhnU/Gkpl+hqvnbmdAih05RBsnezEt6mDWi0HKVgJXB7xeAaxwKDQTZmJi4ILJ+0hJbKRkTw5P1TVRdF4M8XHWPTIcqcLG3SMo2ZPDyCEnWTB9j31WUSocTzGZfkgEZo6rZN6kfew8EcdKu1V4WGpthQ+25lGyJ4dxOUe5dOZuSw5RzBKECSuFo47ypYlnOFGTwoo1E6ipS3A7JOPXNhLcroohTC84yAVT9hFjV0hHNUsQJuxMHtzMFbNLqW+KZ8VHEzl+OtntkPq9+kbfUd3BoxnMm7SPmeMq7X5a/YAlCBOWsjO9XD13OzGivPbRRA4dT3c7pH7rdF0ir6yZyImaFC6ZuZvCUTbYT39hCcKErYFp9Vx9/nZSkxt4Y914dh7Ism6wIVZ1MpUVaybS2BTHlbNLGZV90u2QTAhZgjBhLTWpiavm7iB7UA0fbB2DZ2MBDU2xbocV9VoVSvYM59WPJhIX28pVc7czNNPrdlgmxKybiAl7ifEtXFG0ky17h/HJrhyOvT+Zi6aVMWyQ7bCc4D2TwLub8jlyIp28YceZN3k/iXbTvX7JEoSJCCIwNf8wwwfXsKokn5VrJzC14BAzCiqtJ00f2ntoEB9sHY2qcOHUMgpGHLfG6H7MEoSJKFkZtVx/wVY+2j6KTXtGcOjYAC6eXka6XcXbK03NMazZNpo9lVkMyfBy0fQyuzLaWIIwkSc+rpULp+4jJ+s0H2wdzcvvT+b8yfsosLElzknVyVTeLcnHeyaR6QUHmV5wyI7KDGAJwkSwMcOrGTLQy+qSfN7dVMDBYxmcP2k/CXZlb7e0KmzeM5yNe3JISWxk4dwdZFtDtAlgCcJEtLTkRhbO2cGmshGU7B5B1Yk0Lp5WxtDMWrdDC2veMwms3pRP1Yl0xgw/zvmTrCHanM0ShIl4MTEwY2wlwwefZnVJPis+mkj+iOPMHHuQ9JRGt8MLK03NMWzdN4wte4chwEXTfA3RxgRjCcJEjexMLzfM38rmsuFs25/NvkODmDCqimkFh0hKaO56BVGspVXYeWAIJXtGUN8Yz6jsE8wuLLcEajplCcJElYT4Fs4rrGDC6CNs3J3D9v3Z7DqYxdQxh5k4+ki/u/Ooqq/r6oZdOdScSSI78zSXztrF0IF2Cs50zbEEISJPANcCVao6xT9tEPAHIA/YB3xRVU8EWbYF2Ox/Wa6q1zsVp4lOqUlNzJ+yj8l5h1m/cySf7BrJ9vKhzBhbybico8RE+T0EVOHgsQF8snMk1TWpZKbXcdl5O8nJOmXXNZhuc/IIYhnwK+D3AdO+C7ylqg+JyHf9r78TZNkzqjrDwdhMPzEwrZ7PzdrNkRNprC8dyYdb89i6L5vzxlcwaujJqNxZHj2ZyvqdIzlcPYC05HounraHMcOro7KuxllOjii3WkTy2k2+ASj2P/8d4CF4gjCmT2Vnerlq7g4OHB3I+tKRvLNhHEMyvEwec5jcoSeJjfB+/6pw7FQqW/YOY/+RQSQlNDF34n7G5x6N+LoZ94S6DSJbVQ8BqOohERnaQbkkEVkHNAMPqepLIYvQRC0RGDX0JCOzTrKnMosNu3PwbBxLQlwzecOqyR9xnOxMb0T9p32qNpGyysGUHRpMTV0ScbEtzBh7kMl5h/tde4vpe+HaSD1KVStFJB94W0Q2q+qeYAVFZDGwGCA7OxuPxxPCMPue1+uNiDoUegsdW3diS6Kj6weYOBCuOq+ePSeb2FgVz9bKIeysGMrAxFZmDGli5tAmhqb0zQ42uSUF6Nk262wb1DQKJUfj2Xg0noPeWAQlP6OFK8aeYUpWE0lx6VAf2eNneDyeiPktOCUc6h/qBHFERIb7jx6GA1XBCqlqpf9vmYh4gJlA0AShqkuBpQBFRUVaXFzsRNwh4/F4iIQ6LPEscWzdhd5CStNKHVv/Z6TD9FyY1BxDedVAyioHs6oiA09FIoMH1JI/4jhjhh8nJfHcu8nOi60D6FGd2m+DxuYYyo9kUlY5mEPH01GEQQNqmV14nDHDq0lJagJg/zlHGV4WFS+KmN+CU8Kh/qFOEC8DdwAP+f/+qX0BEckE6lS1QUSygPnAT0Mapel34uNaKRhRTcGIauoa4th7aDBllYP5eMco1u3IZWhmDVkZdWSm+x4D0+odO7evCtX1QnntQE54kzl+KpWDxzJoaY0hLbmeqfmHyB9xnIFp9Y68vzFtnOzmuhxfg3SWiFQA9+NLDM+JyF1AOXCzv2wRcI+q3g1MBH4rIq34BjR6SFW3ORWnMe2lJDYzOe8Ik/OOcNKbRFnlYA4ey2BH+VBaWn39Y0VayUitZ1D6GTLT6xiUXkdm+hmSE5t61IZR3xjLiZoUTnqTOVGTTLX/eXNLLOA7TZSeXM+4kUfJ9997KpLaSExkc7IX06IOZn0uSNl1wN3+5x8AU52Ky5ieGJhWz6zxB5k1/iCtrXC6LokTNSlU1yRzoiaFwyfSKTs0+NPyifFNJATc00iAa06lAvDCu1M+nQbQ2BzLmYaEzyybmX6GcTnHmJgwgIasvQxMO2ONzcY14dpIbUzYiYnxJYyBafWMGf7X6Q2NsZzwplB9Ojngv39oOwEVH+vbwQ9Kr+Ov6QFiY1rJTD9DZtrZRx++Ngi72tm4yxKEMb2UmNDCsEE1DBtUE3T+wLQzABTPKAtlWMb0WpTfcMAYY8y5sgRhjDEmKEsQxhhjgrIEYYwxJihLEMYYY4KyBGGMMSYoSxDGGGOCsgRhjDEmKFGNnsFEROQokX9DyyzgmNtBuMy2gW0DsG0QqvqPVtUhwWZEVYKIBiKyTlWL3I7DTbYNbBuAbYNwqL+dYjLGGBOUJQhjjDFBWYIIP0vdDiAM2DawbQC2DVyvv7VBGGOMCcqOIIwxxgRlCcIYY0xQliBCSEQWikipiOwWke8GmT9aRN4SkU0i4hGRkQHz7hCRXf7HHaGNvO/0chu0iMhG/+Pl0EbeN0TkCRGpEpEtHcwXEXnEv302icisgHnR8h3ozTboD9+BCSLyoYg0iMi32s3r9PfT51TVHiF4ALHAHiAfSABKgEntyjwP3OF/finwlP/5IKDM/zfT/zzT7TqFchv4X3vdrkMfbIOLgVnAlg7mXw28im9s0vOBj6LpO9CbbdCPvgNDgdnAj4FvBUzv8vfT1w87ggidOcBuVS1T1UbgWeCGdmUmAW/5n78TMP9K4A1VrVbVE8AbwMIQxNzXerMNooKqrgaqOylyA/B79VkDDBSR4UTPd6A32yAqdFV/Va1S1Y+BpnazuvP76VOWIEInBzgQ8LrCPy1QCfAF//MbgXQRGdzNZSNBb7YBQJKIrBORNSLyeWdDdU1H2yhavgPd0Vld+8N3oCMh/w5YgggdCTKtfR/jbwELRGQDsAA4CDR3c9lI0JttADBKfbceuA14WEQKHIvUPR1to2j5DnRHZ3XtD9+BjoT8O2AJInQqgNyA1yOBysACqlqpqjep6kzgPv+0U91ZNkL0ZhugqpX+v2WAB5gZgphDraNtFC3fge7osK795DvQkZB/ByxBhM7HwDgRGSMiCcCtwGd6YYhIloi0fSbfA57wP18JXCEimSKSCVzhnxZpznkb+Oue2FYGmA9sC1nkofMycLu/J8/5wClVPUT0fAe6I+g26EffgY50+fvpa3FOrtz8lao2i8i9+H7UscATqrpVRB4A1qnqy0Ax8P+JiAKrgb/3L1stIg/i+4IAPKCqnTXyhaXebANgIvBbEWnF94/NQ6oacTsHEVmOr45ZIlIB3A/EA6jqo8AKfL14dgN1wFf886LiOwDnvg3oJ98BERkGrAMGAK0i8o/4eiudDvb7cTRWf/cpY4wx5jPsFJMxxpigLEEYY4wJyhKEMcaYoCxBGGOMCcoShDHGmKAsQRhjjAnKEoQxxpigLEEY4zAReUlE1ovIVhFZ7HY8xnSXXShnjMNEZJD/SuhkfFdCL1DV427HZUxX7FYbxjjvH0TkRv/zXGAcYAnChD1LEMY4SESKgcuAeapaJyIeIMnVoIzpJmuDMMZZGcAJf3KYgG8ITWMigiUIY5z1GhAnIpuAB4E1LsdjTLdZI7Uxxpig7AjCGGNMUJYgjDHGBGUJwhhjTFCWIIwxxgRlCcIYY0xQliCMMcYEZQnCGGNMUP8/ef+0FUkRhgcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"m.draw_mnprofile('a');"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1f3/8ddJJvtCAtmA7EAIYVEg7JuKKCruS1kEakW0am2tXWzV+v2p/Vbb2m/VugFiXcC1WCkqosgOAQKybyEhG5CdbISsc35/JFqaEgiZ5c6d+TwfjzwIMxPue67D28u9556jtNYIIYQwHy+jAwghhOgaKXAhhDApKXAhhDApKXAhhDApKXAhhDApizM3FhERoRMTE525SSGEML0dO3aUaa0j2z/u1AJPTEwkMzPTmZsUQgjTU0rlnetxOYUihBAmJQUuhBAmdcECV0otVkqVKKX2neO5XyiltFIqwjHxhBBCdKQzR+B/B6a2f1ApFQdMAfLtnEkIIUQnXLDAtdbrgYpzPPV/wK8AmUxFCCEM0KVz4EqpG4DjWuvdds4jhBCiky56GKFSKhB4DLiqk6+fD8wHiI+Pv9jNCSGE6EBXjsD7AEnAbqVULhAL7FRKxZzrxVrrBVrrdK11emTkf41D75QNWaW8svZol35WCCGMVN/Uwoursyiqqrf7n33RBa613qu1jtJaJ2qtE4FCYJjWusju6dpszCrjL6uOUFJj/x0ghBCOtOpAMX/56ghHS2rt/md3Zhjhe8AWoL9SqlApdbfdU1zAHSPiaLZqlu087uxNCyGETT7Ynk9seABj+/Sw+5/dmVEoM7TWPbXWPlrrWK31G+2eT9Ral9k92Vn6RAYzIjGcD7cXICsICSHMIr+8jk1Hy7kjPQ4vL2X3P980d2LekR5HTtlptueeMjqKEEJ0ykc7CvBScNvwWIf8+aYp8OuG9CTYz8IH2wuMjiKEEBfU3GLlo8xCJqZE0isswCHbME2BB/pauP6SXny29wTV9U1GxxFCiPNan1VKUXU900fEOWwbpilwgOkj4qhvsrJ81wmjowghxHl9sL2AiGBfrkiNdtg2TFXgQ2K7kRoTIqdRhBAurbSmgdUHS7hlWCy+FsfVrKkKXCnFD0bEsfd4FfuOVxkdRwghzunjHYU0WzV3pDvu9AmYrMABbh7aGz+LF+9vl0kQhRCux2rVfLA9n5GJ3ekbFezQbZmuwMMCfbl2cE8+/fYEdY3NRscRQoj/kJFTTm55HTNGOfboG0xY4AAzRsZT09DMij0njY4ihBD/4b3tBXQL8OGaQT0dvi1TFviIxHD6RAbx/jY5jSKEcB3ltQ18ua+Im4f2xt/H2+HbM2WBK6WYMTKenfmVHC6qMTqOEEIAsGzncRpbrMwY6Zyps01Z4EDr8BxvL96To3AhhAvQWvPe9nyGJ4TTPybEKds0bYF3D/Jl6qAYlu0spL6pxeg4QggPt+1YBTmlpx1652V7pi1waL2YWV3fzGdyMVMIYbCl2/IJ8bcwbUgvp23T1AU+Ork7yZFBLNmaZ3QUIYQHK69t4Iu9Rdw6LJYAX8dfvPyOqQtcKcWsUQnszK/kwIlqo+MIITzUxzsKaWyxMnOUc9f9NXWBA9w6rPXOzKXb5ChcCOF8Vqtm6bbWOy9Top1z8fI7pi/wsEBfpg3pxSc7j1PbIHdmCiGca1N2GXnldcwa7dyjb3CDAgeYOSqe040tMs2sEMLplm7NJzzQh6mDYpy+bbco8GHxYaTGhLBka56smSmEcJri6npWHSjm9vQ4/CzOu3j5nc6sSr9YKVWilNp31mN/UkodUkrtUUp9opQKc2zMC2Zk1ugE9p+oZldBpZFRhBAe5IPtBbRYtdPuvGyvM0fgfwemtnvsK2CQ1noIcAT4jZ1zXbSbh/YmyNebdzPkzkwhhOM1t1hZujWfCf0iSIoIMiTDBQtca70eqGj32Cqt9XdXDDMAxyy5fBGC/SzcMiyWf+05wanTjUbHEUK4ua8PllBUXc+doxMMy2CPc+A/Ar7o6Eml1HylVKZSKrO0tNQOm+vYnaMTaGy28tEOWXJNCOFY72bk0aubP5NTowzLYFOBK6UeA5qBJR29Rmu9QGudrrVOj4yMtGVzF9Q/JoSRid15NyMfq1UuZgohHCO7tJaNR8uYOSoei7dxY0G6vGWl1FxgGjBLu9DQjzvHJJBfUcf6LMce7QshPNeSjHx8vBV3OHHiqnPpUoErpaYCvwZu0FrX2TeSbaYOjCEi2Jd3M+TOTCGE/Z1pbOHjHQVMHdSTqBB/Q7N0Zhjhe8AWoL9SqlApdTfwNyAE+EoptUsp9ZqDc3aar8WL6SPiWX2ohIIKl/p/ixDCDSzffZzq+mZmG3jx8judGYUyQ2vdU2vto7WO1Vq/obXuq7WO01pf2vZ1nzPCdtbMUfEoYMlWGVIohLAfrTVvbc6jf3QIIxLDjY7jHndittcrLICr0mJ4f3u+LPYghLCbHXmnOHCymrljE1FKGR3HPQscYO7YRCrrmli+W+ZHEULYx9835xLqb+Gmoc5btOF83LbARyd3p390CG9tzpX5UYQQNiuurmflviLuSI8j0NdidBzAjQtcKcWcsa3zo+zMP2V0HCGEyS3Zmk+L1sweY/zFy++4bYED3HRpb0L8Lby1WYYUCiG6rrG5dd6Ty/tHkdDDmHlPzsWtCzzIz8Ltw+P4fO9JSqrrjY4jhDCpL/adpKy2gTkudPQNbl7gAHPGJNCiNe/KkEIhRBe9uSmXpIggJvZz7HQgF8vtCzwxIogr+kexdGseDc0ypFAIcXG+zT/FroJK5o5JwMvL+KGDZ3P7Agf44bhEymob+dfuk0ZHEUKYzN835xLiZ+G2dGPnPTkXjyjw8X0j6BcVzJubjsmQQiFEpxVX1/PZnpPcnh5HsJ9rDB08m0cUuFKKH45LZP+JajLzZEihEKJz3s3Io0Vr5o51rYuX3/GIAofWJde6Bfjw5qZjRkcRQphAfVMLS7fmMznVtYYOns1jCjzQ18L0EXGs3FfE8cozRscRQri45btPUH66kbvGJRkdpUMeU+AAs8ckoJTi7c25RkcRQrgwrTVvbsolJTqYsX16GB2nQx5V4LHhgUwdFMPSbfmcbmi+8A8IITzSluxyDp6s5u7xSS4x62BHPKrAAe4en0RNfTMfZcrCx0KIc1u08Rg9gny58dLeRkc5L48r8GHx4QyLD2PxplxaZOFjIUQ72aW1fHOohDtHJ+Dv4210nPPyuAIHuHt8MvkVdXx9sNjoKEIIF/PmpmP4Wry40wWWTLsQjyzwqwdG0zssgDc2yJBCIcS/nTrdyMc7Crnp0l5EhvgZHeeCPLLALd5e3DUukW25FewprDQ6jhDCRSzdlk99k5W7xycbHaVTOrMq/WKlVIlSat9Zj3VXSn2llMpq+9X41T0v0h0jWm+NXSRH4UIIoKG5hbc25zKhXwT9Y0KMjtMpnTkC/zswtd1jjwKrtdb9gNVtvzeVUH8fZoyM47O9Jyk8VWd0HCGEwZbvOkFJTQP3TDDH0Td0osC11uuBinYP3wi81fb9W8BNds7lFHeNS0LROtevEMJzaa1ZuCGH1JgQJvSLMDpOp3X1HHi01vokQNuvUR29UCk1XymVqZTKLC0t7eLmHKNXWADThvTk/W35VJ1pMjqOEMIga4+UcqS4lvkTk136xp32HH4RU2u9QGudrrVOj4x0rdUsAO6ZmMzpxhbe2yYr9gjhqRauzyEm1J9pQ3oZHeWidLXAi5VSPQHafi2xXyTnGtirG+P7RvDmpmM0NluNjiOEcLJ9x6vYnF3OXeMS8bWYa2BeV9MuB+a2fT8X+NQ+cYxxz8RkiqsbWL77hNFRhBBOtnBDDsF+FmaMijc6ykXrzDDC94AtQH+lVKFS6m7gWWCKUioLmNL2e9Oa2C+C1JgQFqzPxiq31wvhMQoq6lix5yQzRsYR6u9jdJyLdsE1grTWMzp4arKdsxhGKcW9k5J5+IPdrDlcwuQB0UZHEkI4waINOXgpTHPjTnvmOuHjQNOG9KJ3WACvrcs2OooQwgnKaxv4ILOAmy7tTUw3f6PjdIkUeBsfby/mTUhie+4pduS1H/YuhHA3b23Jo77Jyr2TzHn0DVLg/+EHI+IID/Th1bU5RkcRQjhQXWMzb2/JZUpaNH2jzHHb/LlIgZ8l0NfCnDGJfH2wmKziGqPjCCEc5P1tBVTWNXHfpD5GR7GJFHg7c8cmEuDjzWvr5ChcCHfU1GLljY3HGJnYneEJppuH7z9IgbfTPciX6SPj+HTXcZnkSgg39M9vj3O88gw/vtzcR98gBX5O90xIRqnW22uFEO6jxap5dV02aT1DuSzF9ab2uFhS4OfQKyyAW4bG8v72AkprGoyOI4Swky/3F5FTepr7L+9jqkmrOiIF3oH7LutDU4uVxZtkwQch3IHWmpfXHCUpIohrBvU0Oo5dSIF3ICkiiGsH9+SdLXky1awQbmDdkVL2n6jmx5P64O1l/qNvkAI/r/sv60ttQzNvb841OooQwkavrMmmZzd/bhra2+godiMFfh5pvUK5IjWKxZuOcbqh2eg4Qogu2ppTzrbcCuZPTDbdlLHn4z7vxEF+ckVfTtU18W5GntFRhBBd9NI3R4kI9mPGSPNNGXs+UuAXMDQ+nAn9Ili4IYczjS1GxxFCXKQdeafYeLSM+ROT8PfxNjqOXUmBd8JPruhHWW2jLLsmhAm99E0W4YE+zBqVYHQUu5MC74SRSd0ZldSd19ZlU98kR+FCmMWewkrWHi5l3oRkgvwuuPyB6UiBd9JDk/tRUtPAR5kFRkcRQnTSS98cJdTfwpwx7nf0DVLgnTa2Tw+GxYfxytpsGprlKFwIV7f/RBVfHSjmrnFJhJhwubTOkALvJKUUP7syhZNV9XyYWWh0HCHEBby4OosQfws/Gp9kdBSHsanAlVIPK6X2K6X2KaXeU0qZc12iTprQL6L1KHzNUTkKF8KF7T9RxZf7i7l7fBLdAtzz6BtsKHClVG/gISBdaz0I8Aam2yuYK5KjcCHM4YWvW4++7xrnvkffYPspFAsQoJSyAIHACdsjubYJ/SIYnhAuR+FCuKh9x6tYdcD9j77BhgLXWh8H/gzkAyeBKq31qvavU0rNV0plKqUyS0tLu57URbQehfdrPQrfLiNShHA1L6zOItQDjr7BtlMo4cCNQBLQCwhSSt3Z/nVa6wVa63StdXpkpPknUAcY3zeC9IRwXl4j48KFcCV7C1tHntw9Ptntj77BtlMoVwLHtNalWusmYBkw1j6xXJtSip9flUJRdT1Lt8rdmUK4iue/OkxYoA93jU80OopT2FLg+cBopVSgal3aYjJw0D6xXN/YPhGMSe7BK2uPUtcoMxUKYbQdeRWsPVzKvRP7EOqm477bs+Uc+FbgY2AnsLftz1pgp1ym8MhVKZTVNvLWZpmpUAijPb/qCBHBvswd6553XZ6LTaNQtNZPaq1TtdaDtNaztdYetYBkemJ3Lusfyevrs6mpl1V7hDDK5qNlbM4u5/7L+hLo635znnRE7sS00SNT+lNZ18TijblGRxHCI2mtef6rI8SE+jNzlHvN930hUuA2GhzbjasHRrNwQw6nTjcaHUcIj7P2cCk78k7x4BV93W6+7wuRAreDR67qz+nGZl5dl210FCE8itWqeW7lIRJ6BPKDEXFGx3E6KXA7SIkO4Zahsfx9cy4nq84YHUcIj/GvPSc4VFTDz6ek4OPteXXmee/YQX52ZT+01ry4OsvoKEJ4hMZmK8+vOkJqTAjXD+lldBxDSIHbSVz3QGaNSuDDzEJySmuNjiOE2/sgs4D8ijp+NbU/Xl7K6DiGkAK3owev6IufxYvnVx0xOooQbu1MYwsvrs5iRGI4l/ePMjqOYaTA7Sgi2I9545P4bO9JdhdUGh1HCLe1eNMxSmsa+NXUVFpvBPdMUuB2ds/EZHoE+fKHLw6itTY6jhBup7y2gVfXZnPlgGhGJHY3Oo6hpMDtLMTfh4cm9yMjp3VeBiGEfb30Tev8Q7+e2t/oKIaTAneAGSPjSewRyLNfHKLFKkfhQthLXvlplmzN4wcj4ugXHWJ0HMNJgTuAr8WLX16dyuHiGv6xU5ZeE8Je/vTlYby9Wpc2FFLgDnPt4BguiQvjL6uOcKZRFn0Qwla7CypZseck88YnEx3q1uund5oUuIMopXjs2gEUVdezaEOO0XGEMDWtNb//7CA9gny5d1Ky0XFchhS4A41M6s7VA6N5dV02JTX1RscRwrS+3F/MttwKHp6SQoiHLNbQGVLgDvboNQNoarHyF7m5R4guaWy28uwXB+kXFcx0D5yw6nykwB0sKSKI2aMT+TCzgENF1UbHEcJ03snII7e8jt9eNwCLB05YdT6yN5zgocl9CfH34fefyc09QlyMyrpGXlydxYR+EVyWEml0HJcjBe4EYYG+PDS5HxuyylhzuMToOEKYxgurs6ipb+Kx6wZ49C3zHZECd5LZoxNIjgzimRUHaWy2Gh1HCJeXVVzD21vymDEyntSYUKPjuCSbClwpFaaU+lgpdUgpdVApNcZewdyNr8WLJ65LI6fsNG9vyTU6jhAuTWvNUysOEOTrzc+nyE07HbH1CPwFYKXWOhW4BDhoeyT3dXlqFJf1j+SFr7Moq20wOo4QLuubQyVsyCrjZ1em0CPYz+g4LqvLBa6UCgUmAm8AaK0btdYyh+oFPH5dGmeaWmTOcCE60Nhs5ekVB+gTGcTsMQlGx3FpthyBJwOlwJtKqW+VUouUUkHtX6SUmq+UylRKZZaWyux8faOCmTs2kfe357PveJXRcYRwOW9uOkZueR1PTEvzyHUuL4Yte8cCDANe1VoPBU4Dj7Z/kdZ6gdY6XWudHhkpw4AAHprcjx5Bvvzu031YZbZCIb5XVFXPi6uzmJwaxWUevNJOZ9lS4IVAodZ6a9vvP6a10MUFdAvw4ddTU9mZX8myb48bHUcIl/G/nx+kyap58vqBRkcxhS4XuNa6CChQSn03q/pk4IBdUnmAW4fFMiw+jGe/OEjVmSaj4whhuIyccpbvPsF9k/oQ3yPQ6DimYOsJpp8AS5RSe4BLgf+1PZJn8PJSPHXjIMpPN/LXr+WCpvBsTS1Wnvx0P73DAvjxpD5GxzENmwpca72r7fz2EK31TVrrU/YK5gkG9e7GrFHxvL0lj4MnZZ4U4bne3pLH4eIafnd9GgG+3kbHMQ25xGuwX1zVn24BPjz+T7mgKTxTUVU9//fVESalRHJVWrTRcUxFCtxgYYG+/OaaVHbkneLDzAKj4wjhdE+vOEBTi5Wnbhwo851cJClwF3Db8FhGJnXn2ZWHqDjdaHQcIZxm3ZFSPtt7kgcv70tCj/+6jURcgBS4C1BK8cxNg6itb+YPn8tsBMIz1De18LtP95EcGcR8WSatS6TAXURKdAjzJiTz0Y5CtuaUGx1HCId7Zc1R8srreObGQfhZ5MJlV0iBu5CHJvclNjyA33yyl4ZmWcleuK8jxTW8ui6bm4f2ZmzfCKPjmJYUuAsJ9LXw+5sHk1N6mpfXZBsdRwiHsFo1v1m2l2A/C49fN8DoOKYmBe5iJqVEctOlvXh17VGyimuMjiOE3S3ZmseOvFM8fl2aTBVrIylwF/TEtDSC/Cw8umyvjA0XbqWoqp7nVh5mfN8IbhnW2+g4picF7oJ6BPvx+HVp7Mg7xZKteUbHEcIutNY88ek+mq1Wfn/zIBnzbQdS4C7q1mG9mdAvgme/OEThqTqj4whhs3/tOclXB4p5+MoUGfNtJ1LgLkopxR9uGQzAb5btRWs5lSLMq7y2gf9Zvp9L4sKYN0HGfNuLFLgLiw0P5NFrUtmQVcZHmYVGxxGiy55cvp/a+mb+dNsQvL3k1Im9SIG7uFmjEhiV1J2nPztAUVW90XGEuGgr9xWxYs9JHprcl5ToEKPjuBUpcBfn5aV47tYhNLVYeXTZHjmVIkyl4nQjj/9zH2k9Q7lX5vm2OylwE0iMCOLRqamsPVzKB9tlxkJhDlprHv/nXqrONPL8HZfIAsUOIHvUJOaMSWRMcg+eXnGAggoZlSJc3/LdJ/h8bxEPT0lhQM9Qo+O4JSlwk/DyUvzxtiEopfjlx7vlBh/h0oqr6/ndp/sZGh/GfBl14jBS4CYS1z2QJ6YNICOngjc35xodR4hz0lrz63/soaG5hedvvwSLnDpxGNmzJnNHehyTU6N4buUhDhfJXCnC9bybkcfaw6U8OjWV5Mhgo+O4NZsLXCnlrZT6Vim1wh6BxPkppXjutiGE+lv46fvfyrSzwqUcLanhmc8OMiklkrljE42O4/bscQT+U0CWkXGiiGA/nrt1CIeKavjTysNGxxECgMZmKz99fxdBfhb+1Ha9RjiWTQWulIoFrgMW2SeO6KzJA6K5c3Q8izYeY2NWmdFxhOAvXx1h/4lqnr1lMFGh/kbH8Qi2HoH/FfgVYO3oBUqp+UqpTKVUZmlpqY2bE2d77No0+kQG8fMPd1Fe22B0HOHBNmaV8fr6bGaMjOOqgTFGx/EYXS5wpdQ0oERrveN8r9NaL9Bap2ut0yMjI7u6OXEOAb7evDhjKJVnmnjkIxlaKIxRVtvAwx/uok9kME9MSzM6jkex5Qh8HHCDUioXeB+4Qin1rl1SiU4b2Ksbj183gLWHS1m86ZjRcYSHsVo1j3y4m6ozTfxt5lACfS1GR/IoXS5wrfVvtNaxWutEYDrwjdb6TrslE502e3QCV6VF89zKQ+wprDQ6jvAgizbmsO5IKU9MSyM1Ru62dDYZB+4GlGq9SzMy2I8Hlu6k6kyT0ZGEB9iRd4o/rjzM1IEx3Dkq3ug4HskuBa61Xqu1nmaPP0t0TVigLy/NHMbJynp++dFumbVQOFTF6UYeXLqTnmH+PCdDBg0jR+BuZHhCOI9ek8qqA8Us2iDnw4VjWK2an32wi/LaRl6ZOZxuAT5GR/JYUuBu5u7xSUwdGMOzKw+xPbfC6DjCDb285ijrj5Tyu+vTGBzbzeg4Hk0K3M0opfjj7UOICw/ggSU7KamRVXyE/aw/Usr/fX2EGy/txSw57204KXA3FOrvw6t3DqemvpkHluyksbnD+6yE6LSCijp+8t63pESH8IdbBst5bxcgBe6mBvQM5bnbhrA99xS//+yA0XGEyZ1pbGH+OzvQWvP67OEy3ttFyH8FN3bDJb3YU1DJoo3HGBwbxm3DY42OJExIa82jy/ZwqKiaxT8cQUKPIKMjiTZyBO7mHr0mlbF9evDbT/bybf4po+MIE1qwPodPd53gkSkpXN4/yug44ixS4G7O4u3FyzOHERPqz/x3dnCy6ozRkYSJrD5YzLMrD3HdkJ48cHlfo+OIdqTAPUB4kC+L5qZzprGFe97OpK6x2ehIwgQOF9Xw0HvfMrBXKH++7RK5aOmCpMA9REp0CC/OuJT9J6r5hcxcKC6gvLaBeW9vJ9DPwsI56QT4ehsdSZyDFLgHuSI1mt9eM4DP9xbxxy9lJR9xbvVNrf9SK6luYMHs4fTsFmB0JNEBGYXiYeZNSCK3/DSvrcsmvnsgM+VmDHEWq1Xz8w938W1BJa/MHMbQ+HCjI4nzkAL3MEop/t8NAzleeYYnPt1HrzB/LpORBaLNc18e4vO9RTx27QCuGdzT6DjiAuQUigeyeHvxt5nD6B8dwgNLdrK3sMroSMIFvL0ll9fX5TB7dALzJiQZHUd0ghS4hwr2s/DmXSMIC/Tlh29u41jZaaMjCQOt2HOCJ5fv58oB0Tx5fZqMODEJKXAPFh3qzzt3j0QDcxZvpaRaJr7yRBuzynj4g12kJ4Tzt5lDsXhLLZiF/JfycMmRwSz+4QjKaxuZ++Z2Wc3Hw+wprOTedzJJjghm0ZwR+PvIcEEzkQIXXBoXxmt3DudoSQ0/fHMbpxvkRh9PcLiohjmLtxEW6MtbPxpJt0BZmMFspMAFABNTInlpxlD2FFYx761M6ptajI4kHOhY2WnufGMrfhYvlt4ziphu/kZHEl3Q5QJXSsUppdYopQ4qpfYrpX5qz2DC+aYO6smfbx9CxrFy7pd5xN1W4ak6Zi3MoMWqWTJvlMwuaGK2HIE3A49orQcAo4EHlFJp9okljHLz0FieuWkQ3xwq4f4lO6TE3UzhqTqmL8igtqGZd+4eSd+oEKMjCRt0ucC11ie11jvbvq8BDgK97RVMGGfWqASevnEgXx+UEncnBRWt5V19pol3541iYC9Zz9Ls7HIOXCmVCAwFtp7juflKqUylVGZpaak9NiecYPaYRJ6+adD3Jd7QLOfEzaygoo4ZC1vLe8m80QyJDTM6krADmwtcKRUM/AP4mda6uv3zWusFWut0rXV6ZGSkrZsTTjR7dALPtJX4vLdkGlqzyi6t5Y7Xt3xf3rKSvPuwqcCVUj60lvcSrfUy+0QSruTO0Qn86bYhbDpaxpw3tlFdL+PEzWT/iSrueG0LTS1W3p8/RsrbzdgyCkUBbwAHtdZ/sV8k4WpuT4/jbzOHsbuwkpkLMyivbTA6kuiEHXkVTF+Qga/Fiw/uHUNar1CjIwk7s+UIfBwwG7hCKbWr7etaO+USLubawT1ZOCedoyW13PrqZvLKZe4UV/bVgWJmLtxKjyBfPrpvDH0ig42OJBzAllEoG7XWSms9RGt9advX5/YMJ1zLZf2jWHrPaKrONHHrq5vZU1hpdCRxDku25nHvO5mkxoTw8Y/HEhseaHQk4SByJ6a4KMPiw/n4x2Px9/Fm+oIMvjlUbHQk0cZq1fzpy0M89sk+JqVE8t780UQE+xkdSziQFLi4aH0ig1l2/1iSI4OY91YmizbkoLWssWmkM40tPPjeTl5ek830EXEsnJNOoK+s1+LupMBFl0SF+PPhvWO4Ki2GZz47yG8/2Ss3/BikuLqeO17fwhf7WlfS+cMtg2VKWA8h/5VFlwX6Wnhl1jAeuLwP720rYNaiDJlT3MkycyuY9tJGckprWTg7nXsmJstiDB5EClzYxMtL8curU3lxxlD2Ha9m2ksb2ZFXYXQst6e15q3NuUxfkEGQrzf/uH8sV6ZFGx1LOJkUuLCLGy7pxScPjCXA15sfvJ7BGxuPyXlxB6ltaObhD/musToAAAjwSURBVHbx5PL9XNY/kk8fHE9qjIzx9kRS4MJuUmNCWf7geC7rH8nTKw5wz9uZnDrdaHQst7LveBXXv7SR5btP8PMpKSyYnU63AFmIwVNJgQu76hbgw8I56Tx5fRrrj5RxzQsb2JJdbnQs07NaNYs3HuOWVzZzprGFpfeM5qHJ/fDykvPdnkwKXNidUoq7xiWx7P7WUyozFmbw1L8OyCo/XVR4qo4739jKUysOMKFfBJ//dAKjk3sYHUu4ABkoKhxmUO9ufPbQeJ794hCLNx1j3ZES/nz7JQyNDzc6milorfkos5CnVhxAa82ztwzmByPiZJSJ+J5y5oWm9PR0nZmZ6bTtCdexMauMX368m6LqeuaOSeSRq1II8Zdztx3JKa3lsU/2sSWnnFFJ3fnz7ZcQ111uifdUSqkdWuv0/3pcClw4S019E3/+8jBvZ+QRHeLP/9yQxtUDY+SI8iz1TS0sWJ/D39Ycxc/ixaPXpDJjRLyc6/ZwUuDCZezMP8Vvl+3lUFENY/v04HfXp3n8MDitNV/uL+KZzw5SeOoM1w3pyZPT0ogKldXihRS4cDHNLVaWbsvnL18dofpME9NHxvPTyf2I9sDC2l1QybNfHGJLTjn9o0P43fVpjOsbYXQs4UKkwIVLqqxr5K9fZ/FuRh7eXoofjk3kvkl9CA/yNTqawx0pruH5VYf5cn8x3YN8eXhKCjNGxMk8JuK/SIELl5ZfXsdfvz7CJ7uOE+jjzazRCcwbn+SWpxD2FlbxytqjrNxfRLCvhXsmJvOj8UkE+8mgMHFuUuDCFA4X1fDymqOs2HMCi5cXtw7vzdyxiaY/R261atYdKWXxpmNsyCojxN/CnDEJzBuf7BH/2hC2kQIXppJfXsfr67P5eEchDc1WRiZ1Z/boBKakRePv4210vE4rq23gk53HeScjj/yKOiJD/PjRuCRmjY4nVIZRik6SAhemVFnXyIeZBbyTkUdBxRlC/C1MG9KLW4b1Znh8uEsOr6tvauGbQyUs21nI2sOlNFs1IxO7M3tMAlcPjMHXIue4xcWRAhem1mLVbM4uY9nO46zcV8SZphYiQ/yYkhbNlLRoRif1IMDXuCPzstoG1h0uZdWBItYfKeNMUwvRoX7cNLQ3tw6LJSU6xLBswvwcUuBKqanAC4A3sEhr/ez5Xi8FLuyhtqGZ1QeLWbW/mDWHS6hrbMHX24thCWGM6xPBsIRwBsd2c9gpCq01J6vq2VNYybZjp9icXcahohoAokP9uCothqmDYhid3ANvF/wXgjAfuxe4UsobOAJMAQqB7cAMrfWBjn5GClzYW31TCxk55WzOLmdjVhkHTlZ//1xyZBD9o0PoExlMn6ggeocFEhPqT1So3wXPo2utqWlopqS6nqKqBvIqTpNdcprs0loOnKymtKYBAF+LFyMSwxnbJ4LxfSMY3LubS57WEebWUYHbMm5pJHBUa53TtoH3gRuBDgtcCHvz9/Hmsv5RXNY/Cmg9Z76nsIo9hZXsLqzicFENqw4U02LV7X7Oi2A/C8F+lv8Yd93Q3EJtfTOnG1pobPnPNT4DfLxJjgxiQt8ILokLY0hsNwb0DDXVRVXhXmwp8N5AwVm/LwRGtX+RUmo+MB8gPj7ehs0JcWFhgb5MTIlkYkrk9481NLdQUFHHicp6iqtbv6rrm6ltaKa2vvn7ctdo/CzeBPtZCPKz0D3Ih+hQf6JD/YkND6BXtwA5uhYuxZYCP9cn+b/Ox2itFwALoPUUig3bE6JL/Cze9I0KoW+UXEgU7sWW8UyFQNxZv48FTtgWRwghRGfZUuDbgX5KqSSllC8wHVhun1hCCCEupMunULTWzUqpB4EvaR1GuFhrvd9uyYQQQpyXTbPnaK0/Bz63UxYhhBAXQe7pFUIIk5ICF0IIk5ICF0IIk5ICF0IIk3LqbIRKqVIgz2kbdJwIoMzoEAby9PcPsg9A9gE4bx8kaK0j2z/o1AJ3F0qpzHNNLOMpPP39g+wDkH0Axu8DOYUihBAmJQUuhBAmJQXeNQuMDmAwT3//IPsAZB+AwftAzoELIYRJyRG4EEKYlBS4EEKYlBT4WZRSU5VSh5VSR5VSj57j+QSl1Gql1B6l1FqlVOxZz81VSmW1fc11bnL7sXEftCildrV9mXJqYaXUYqVUiVJqXwfPK6XUi237Z49SathZz7nLZ8CWfWD6zwB0ah+kKqW2KKUalFK/aPfcef8O2ZXWWr5arwN4A9lAMuAL7AbS2r3mI2Bu2/dXAO+0fd8dyGn7Nbzt+3Cj35Mz90Hb72uNfg922AcTgWHAvg6evxb4gtYVqUYDW93pM2DLPnCXz0An90EUMAL4PfCLsx6/4N8he37JEfi/fb9Is9a6EfhukeazpQGr275fc9bzVwNfaa0rtNangK+AqU7IbG+27AO3oLVeD1Sc5yU3Am/rVhlAmFKqJ+7zGbBlH7iNC+0DrXWJ1no70NTuqc78HbIbKfB/O9cizb3bvWY3cGvb9zcDIUqpHp38WTOwZR8A+CulMpVSGUqpmxwb1TAd7SN3+Qx0xvneqyd8Bs7HqZ8DKfB/68wizb8AJimlvgUmAceB5k7+rBnYsg8A4nXrbcUzgb8qpfo4LKlxOtpH7vIZ6IzzvVdP+Aycj1M/B1Lg/3bBRZq11ie01rdorYcCj7U9VtWZnzUJW/YBWusTbb/mAGuBoU7I7Gwd7SN3+Qx0Rofv1UM+A+fj1M+BFPi/XXCRZqVUhFLqu332G2Bx2/dfAlcppcKVUuHAVW2PmU2X90Hbe/f77jXAOOCA05I7z3JgTttIjNFAldb6JO7zGeiMc+4DD/oMnI9zF3s3+mqvK33RenX9CK1XkR9re+wp4Ia2728DstpeswjwO+tnfwQcbfu6y+j34ux9AIwF9tJ6jnwvcLfR76WL7/894CStF6cKgbuB+4D72p5XwMtt+2cvkO6Gn4Eu7QN3+Qx0ch/EtD1eDVS2fR/a9tx//R1y1JfcSi+EECYlp1CEEMKkpMCFEMKkpMCFEMKkpMCFEMKkpMCFEMKkpMCFEMKkpMCFEMKk/j8CEUy8paukbQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"px, py = m.profile('a', subtract_min=True)\n",
"plt.plot(px, py);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}