diff --git a/notebooks/pumpingtest_hypothetical.ipynb b/notebooks/pumpingtest_hypothetical.ipynb index e03c83e..449ad1f 100644 --- a/notebooks/pumpingtest_hypothetical.ipynb +++ b/notebooks/pumpingtest_hypothetical.ipynb @@ -1,5 +1,12 @@ { "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calibrating TTim model to hypothetical pumping test" + ] + }, { "cell_type": "code", "execution_count": 1, @@ -9,6 +16,9 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", + "import sys\n", + "\n", + "sys.path.insert(1, \"..\")\n", "from ttim import *" ] }, @@ -21,7 +31,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 72, "metadata": {}, "outputs": [], "source": [ @@ -40,7 +50,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 73, "metadata": {}, "outputs": [ { @@ -69,14 +79,15 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 77, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "......................................................................\n" + "0.007522051457496564\n", + "..................................\n" ] }, { @@ -112,36 +123,36 @@ " \n", " \n", " kaq0\n", - " 60.000021\n", - " 2.495565e-05\n", - " 0.000042\n", + " 60\n", + " 7.357572e-06\n", + " 1.22626e-05\n", " -inf\n", " inf\n", - " 10\n", - " [60.00002072174987]\n", + " 59\n", + " [59.99999146945327]\n", " \n", " \n", " Saq0\n", - " 0.000100\n", - " 2.087595e-10\n", - " 0.000209\n", + " 9.99999e-05\n", + " 6.106179e-11\n", + " 6.10619e-05\n", " -inf\n", " inf\n", " 0.0001\n", - " [9.999968805513572e-05]\n", + " [9.999986102899141e-05]\n", " \n", " \n", "\n", "" ], "text/plain": [ - " optimal std perc_std pmin pmax initial \\\n", - "kaq0 60.000021 2.495565e-05 0.000042 -inf inf 10 \n", - "Saq0 0.000100 2.087595e-10 0.000209 -inf inf 0.0001 \n", + " optimal std perc_std pmin pmax initial \\\n", + "kaq0 60 7.357572e-06 1.22626e-05 -inf inf 59 \n", + "Saq0 9.99999e-05 6.106179e-11 6.10619e-05 -inf inf 0.0001 \n", "\n", " parray \n", - "kaq0 [60.00002072174987] \n", - "Saq0 [9.999968805513572e-05] " + "kaq0 [59.99999146945327] \n", + "Saq0 [9.999986102899141e-05] " ] }, "metadata": {}, @@ -150,23 +161,24 @@ ], "source": [ "cal = Calibrate(ml)\n", - "cal.set_parameter(name='kaq0', initial=10)\n", + "cal.set_parameter(name='kaq0', initial=59)\n", "cal.set_parameter(name='Saq0', initial=1e-4)\n", "cal.series(name='obs1', x=robs, y=0, layer=0, t=tobs, h=hobs)\n", + "print(cal.rmse())\n", "cal.fit(report=False)\n", "display(cal.parameters)" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 75, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "rmse: 5.347440894236035e-07\n" + "rmse: 5.813313744921566e-07\n" ] } ], @@ -182,7 +194,7 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, "execution_count": 6, @@ -191,7 +203,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcU9X9//HXJwyTaKtUlEXFAS1g3QuOS1gzwOBCKyiLUqxIEcQF9cdDkYqo1cFxrVVxAVzRKgqooFURZiZKIeCGWJciiKCoP/dqEQaEnO8fd6iIA7MkmZvl/Xw88kjuzZncD+cRPvfmnHPPMeccIiKSWwJ+ByAiIg1PyV9EJAcp+YuI5CAlfxGRHKTkLyKSg5T8RURykJK/iEgOUvIXEclBSv4iIjlIyV9EJAfl+R3Ajuy1116uTZs2fochIpJRXnvttS+dc81qKpe2yb9Nmza8+uqrfochIpJRzGxNbcqp2UdEJAcp+YuI5CAlfxGRHKTkLyKSg5KS/M3seDNbbmYrzWxcNe8HzeyxqveXmFmbZBxXRETqJ+Hkb2aNgDuAE4CDgcFmdvB2xYYD3zjn2gK3ANcnetydicVilJaWEovFUnkYEZGMlYyhnkcDK51zqwDMbDrQF3hnmzJ9gauqXs8EJpmZuRSsIbn4pZe4tbiY7zZvpjwvj5vvvJPDjz0WdtnFe+y6q/fcuDGYJfvwIiIZIRnJf1/go2221wLH7KiMc26zmX0L7Al8uW0hMxsJjAQoKCioVzBLnn+e6Zs2eRubNsFZZ1VfMBD48USw/YkhgdfL3n+fuatW0bW4mHA4XK9/g4hIqiUj+Vd3+bz9FX1tyuCcmwJMASgsLKzXr4JjTjiBY//6V/J++IHd8vL468SJHNS6NWzYAOvXe8/bvq5u34YN8NVX1b9fw4+VI4BDgH9fdRVf9O5Ns969oUMH+O1vYY896vNPEhFJumQk/7XAfttstwI+2UGZtWaWBzQBvk7CsX/m2K5duaWigmg0SiQS4aBkXn075/2a2MFJ47H772f2tGkc5BwdnaN1LAZz5/74923aeCeCbR/77KPmJxFpcMlI/q8A7cxsf+Bj4DTgD9uVmQMMBWLAAKA8Fe39W4XD4dQ0uZhBMOg9fvWrn71dkJ/PU48/zuObNpGfn0/Z888TPuAAeOMNWLr0x8eTT/74R82a/fyE0Lat1ywlIpIilowcbGYnAn8DGgH3OecmmtnVwKvOuTlmFgIeAjrgXfGftrWDeEcKCwtdJs7tE4vF/verY4cnoP/+F5Yt++kJ4e234YcfvPd/8Qs44oifnhAOOYTY66/X/NkiktPM7DXnXGGN5VJ4AZ6QTE3+9bZpk3cC2PaEsGwZrFsHQDwvj7e2bKECeCI/n+vKywl36uRvzCKSdpT8s0E8DitXwtKlLLrrLr5/8UW6AUHgm6ZN2WPUKBg8GA491O9IRSRN1Db5q2E5nQUC0L49nHoqVlpK3112Ye9AgJGNG2Nt28J118Fhh3mPa6+FVTttSRMR+R9d+WeQn/UnfP45zJgBjzwCixZ5hY45BgYP5tUDDmDeW2+pf0Akx6jZJ9esWQPTp8Ojj8KyZWwByoE78/MZW15OuHNnvyMUkQagZp9c07o1XHopvPEGUy66iFIzfgM8uWkTbU85Be67DzZu9DtKEUkTSv5Z6LBBg7g2FKJdIMCfGjdm1yZNYPhw7yaz0lL45hu/QxQRnyn5Z6FwOExZWRlXlpQw4sUX+cXy5TBvHhx+OFx2Gey3H1x4IXzwgd+hiohP1Oafa958E26+2eskjsdhwAC45BJiP/ygG8hEsoA6fGXn1q6F226DyZPhu+94KRDgJmB+fj5l5eU6AYhkKHX4ys61agU33AAffcT8E0+kTTzOnHicFyoreW/aNL+jE5EUU/LPdbvvzi8uv5xDQyHOMaMdMPTuu2HgQO/uYhHJSkr+QjgcZm55OQUTJ7J6/ny48kp49lk46CC44AL48suaP0REMora/KV6n34KV10F99wDv/wljBsHF13krVomImlLbf6SmL339jqD33oLunf3hoi2bw8PPABbtvgdnYgkSMlfdu6gg2DOHIhGvRPCsGHQsSPv3HILpaWlxGIxvyMUkXpQ8pfa6d4dFi+G6dOp/PJLDh4zhkMvu4wzi4p0AhDJQEr+UnuBAJx6KreNGsWlZvQAXt+4kfXXXPPjKmQikhGU/KXOuvbqxe2hEIcFAkQDAXo+9xwUFnq/DEQkIyj5S51tnTtoREkJTRcsgCeegK++gk6d4Nxz4T//8TtEEamBhnpKcvz3vzBhAtx+OzRvDn/7GwwaBGZ+RyaSUzTUUxrWbrt5Cf/ll2HffeG00+CEE3h95kyNChJJQ3l+ByBZ5sgjYckSuPNONo8bx0Fz5zLDjN7BIC9owjiRtKErf0m+Ro1g9GjuGj2a54FS5yirrORfjz/ud2QiUkXJX1KmsG9fhoRCDDbjAOCsO++EG2/UHcIiaUDJX1ImHA5TVl7O4RMnsvqZZwj06QNjx0KXLrB8ud/hieQ0jfaRhuMcTJ8O558P69ezesQIprdoQfcePdQXIJIktR3tow5faThmMHgwFBXx9aBBtLn9djoDw4JB7q+o0AlApAGp2UcaXsuWTD7+eIaacTjw8saNfHnLLd4vAxFpEEr+4otIUREzQiE6BAIsCwT4/YwZcOqp8PXXfocmkhPU7CO+2DpFRDQaJa9rV1i40LtDeOFC3rn0UmZ//z2RSERNQSIpklCHr5k1BR4D2gCrgUHOuW+2K/Nb4C5gd2ALMNE591hNn60O3xz0+uus79+fXVev5lYzrgwGeU43honUSUNN7zAOKHPOtQPKqra3tx44wzl3CHA88Dcz+1WCx5Vs1LEjdwwbxiQzLnSOispK3pwxw++oRLJSosm/L/Bg1esHgX7bF3DOveecW1H1+hPgc6BZgseVLNWluJixoRD9AgH2A8666y64/351BoskWaLJv4Vz7lOAqufmOytsZkcD+cD7CR5XstTWvoBjSkpYM3s2jcJh+NOfYMgQ+PZbv8MTyRo1tvmb2XygZTVvjQcedM79apuy3zjn9tjB5+wNRIGhzrlqV/0ws5HASICCgoIj16xZU5t/g2SzLVvg+uvhiiuobNGCx046ifZnnKF+AJEdqG2bf6IdvsuBiHPu063J3Tl3YDXldsdL/KXOuVo14qrDV7b1rylT2H3UKPZxjivy8jgpGiXcubPfYYmknYbq8J0DDK16PRSYXU0g+cCTwLTaJn6R7T3z1Vd0NGM2ULp5M3uddZa3epiI1Euiyf86oNjMVgDFVduYWaGZ3VNVZhDQDTjTzN6oevw2weNKjolEImwIBjktEGBM48b8+v33oUMHWLTI79BEMpImdpOMEYvFiEaj3s1f+fneMpFr1rDm7LN5ZJ99iGiCOJGGafNPJSV/qdG33/LlKaewV3k5zwIjQiFm6qYwyXFaw1eyX5MmTO3Zk/PN6AksrKzk3w8/7HdUIhlByV8yWqSoiPtCIboHAgTMGDp1KkyerJvCRGqg5C8ZbetNYX1LSvjs2WcJ9OgBo0bBmWfC+vV+hyeSttTmL9llyxYoKYG//AUOPRRmzYJ27fyOSqTBqM1fclOjRnDllfDcc/Dxx1BYCE884XdUImlHyV+y03HHwdKl8JvfQP/+cPHFLF6wgNLSUmKxmN/RifhOzT6S3TZuhDFj4M47+WcgwGnA18EgZWVlGhIqWUnNPiIAwSDccQezBw2iQzzOK/E4hRs3Eo1G/Y5MxFdK/pITml90Ed2DQdYBZfE4p371lYaDSk5T8pecEA6Hub2igqcnTGBdp04ccPPN3joBGzb4HZqIL9TmL7knHoerr/aGgx55pDcctHVrv6MSSQq1+YvsSCAAV10Fc+bAihXecNDycr+jEmlQSv6Su37/e3jlFWjeHIqLWT16NKXXXquhoJIT8vwOQMRX7dvD4sV81bcvbSZN4gAzfh8M8rRmB5Uspyt/kd12Y0qvXlxmxkDnmF9ZydInn/Q7KpGUUvIXwZsd9G+hEL8PBGgDjJgyRf0AktWU/EX4cXbQLiUlfPDYYzTeZx/o3RtuvVX3A0hW0lBPkep89x388Y/eiKChQ1l85plUxGLeEpLqC5A0pmUcRRK1zf0Ar5rR34wvNC+QpDmN8xdJVNX9ADNPP50DnWNxPM5vN2xg2rRpfkcmkjAlf5Ea7HvuuXTLy2M9UAHE77lH9wJIxlPyF6lBOBzm2LPO4mjgJWDy5s3kXXwxbN7sd2gi9abkL1ILZ5xxBht22YXfBQJMatSIoxYtghNO4JW5c7VAjGQk3eErUgtbh4JGo1GOjETg3XeJn302TcvKeMSMa9QRLBlGyV+klsLh8I/JPRzm4ZdfpvfkySx0jtOrFohR8pdMoWYfkXpqN3QoXYNBVgJPxeMM/uQT3RAmGUNX/iL1FA6HmVZRQfkLL9D6xRdpM2kS/Pe/MHmyt3ykSBpT8hdJwP+aguJxuOYab52A996DJ56Ali39Dk9kh9TsI5IMgQBceSXMmAFvvAFHHw1Ll/odlcgOJZT8zaypmc0zsxVVz3vspOzuZvaxmU1K5JgiaW3AAFi40Gv779LFWyJSJA0leuU/DihzzrUDyqq2d+Qa4MUEjyeS/jp08FYIO/xwGDCAD0eO1AphknYSbfPvC0SqXj8IRIFLty9kZkcCLYDngRonHBLJeC1bQkUFn59yCgVTp/JrM34XDPKMVgiTNJHolX8L59ynAFXPzbcvYGYB4GbgkgSPJZJZQiHu7dKFS80Y4BzPV1by6pw5fkclAtTiyt/M5gPVDVsYX8tjnAs865z7yMxqOtZIYCRAQUFBLT9eJH1FioroGQqxfONGHorHOfzee2HgQOjY0e/QJMclNJ+/mS0HIs65T81sbyDqnDtwuzJ/B7oCceCXQD5wp3NuZ/0Dms9fskYsFiMajXJiq1Yccfnl8MUX8NBD0L+/36FJFmqQxVzM7EbgK+fcdWY2DmjqnBu7k/JnAoXOufNr+mwlf8lKn30GJ58MsRiUlMBll0ENv4hF6qKhFnO5Dig2sxVAcdU2ZlZoZvck+Nki2adFC29h+CFD4PLLvaUiKyv9jkpykJZxFPGDc1BaCuPHw7HHwlNPeScGkQRpGUeRdGbmNfnMmgVvvglHHQXLlvkdleQQJX8RP51yCixY4M0N1Lkz/77hBi0OIw1CE7uJ+K1jR3j5Zdb17En7Sy/lP2b0DAYp0w1hkkK68hdJB/vsw12nncZMM653jkkbN/LS/Pl+RyVZTMlfJE106d2bYcEgE834k3OMeuop+Pprv8OSLKXkL5ImwuEw88vLCUycyIorrqDJW295I4Hee8/v0CQLaainSLpauBD69YMtW7zFYSIRvyOSDKChniKZrnNnWLIE9t4biovhvvv8jkiyiJK/SDo74ABYtAh69IDhw2HsWG9YqEiClPxF0l2TJvCPf8A558CNN3oTwn3/vd9RSYZT8hfJBHl5cMcdcNttMGcOdO0KH3/sd1SSwZT8RTKFGYweDU8/DStWeIvEv/aa31FJhlLyF8k0J57o9QPk5UG3bvDkk35HJBlIyV8kEx12GLz8svd8yilUHH88sUWL/I5KMoiSv0imatGCJdddx4xGjSiaO5f3unVj8Usv+R2VZAglf5EMVh6LMdg5rgaGbtlCy2HD4Jtv/A5LMoCSv0gGi0Qi5AeDXN2oEcMbN6bgo48gHIaVK/0OTdKcpnQWyWDhcJiysjKi0SiRSITA5s3eGsFbVwfr0sXvECVNaW4fkWyzciX06QOrV7Ni3DhmhkJEIhGtDZAjaju3j678RbJN27YQi/FtcTHtrr6azVocRqqhNn+RbNS0KXf368f9Zkxwjns3bmTBvHl+RyVpRMlfJEt169WL84JBLjNjsHOMmjULvvjC77AkTSj5i2SpcDhMWXk5u02cyPKSEnZ/7z2vI/jf//Y7NEkD6vAVyRWLF0PfvrBpE8ya5U0TLVlHi7mIyE8de6y3OMw++8Bxx2lxmByn0T4iuaRNG29SuIEDYfhwPo5GmXbggUR69NBIoByjK3+RXFO1OMxnffuy70MP0XbCBPr06EEsFvM7MmlASv4iuahxY+47+mjGmtHfOZ6rrOTlp5/2OyppQEr+IjkqUlTEpFCIgYEAhwFnP/AAvPOO32FJA1Gbv0iO2nZeoJXNm3P45ZdDp068c/XVzP7+e00JkeUSGuppZk2Bx4A2wGpgkHPuZ/PJmlkBcA+wH+CAE51zq3f22RrqKdLAPvyQ9UVFNF61ivPMeDgUoqysTCeADNNQQz3HAWXOuXZAWdV2daYBNzrnDgKOBj5P8LgikmwFBdx1+umUA1Oc48rKSqLl5X5HJSmSaPLvCzxY9fpBoN/2BczsYCDPOTcPwDm3zjm3PsHjikgKdDr+eAaEQtxtxqXOMaK8HDZs8DssSYFEk38L59ynAFXPzasp0x74j5k9YWZLzexGM2uU4HFFJAXC4TAvlJfzTUkJq0ePZq+KCigqgs/1Yz3b1Njha2bzgZbVvDW+DsfoCnQAPsTrIzgTuLeaY40ERgIUFBTU8uNFJJnC4fCP7fxFRTBkCBxzDPzjH3Dwwf4GJ0lT45W/c66Xc+7Qah6zgc/MbG+AqufqLg/WAkudc6ucc5uBp4COOzjWFOdcoXOusFmzZvX/V4lIcpx8Mrz4IlRWQqdOvH377ZSWluqGsCyQaLPPHGBo1euhwOxqyrwC7GFmW7N5D0CDiUUyxVFHweLFrN9zT9pfcAHvjx9Pz549dQLIcIkm/+uAYjNbARRXbWNmhWZ2D4BzbgtwMVBmZv8CDJia4HFFpCG1bs3dp59OBXCPc1yhkUAZT1M6i0itxGIxjuvRg5s2bmSkc3zZsyd7PfMMhEJ+hybb0JTOIpJU4XCYueXlfFVSwprzzmOvsjLo2VOrg2UoTe8gIrX2k5FAkQj88Y/eOgHPPgsHHuhrbFI3uvIXkfoZMACiUVi3DsJhb1SQZAwlfxGpv2OO8ZaHbNkSiovhoYf8jkhqSclfRBKz//6wcCF07QpnnAFXXQVpOpBEfqTkLyKJ22MPeO45GDYM/vIXvjjhBG645hrdC5DGlPxFJDny8+Hee/nw7LNpNncu4SuuYKCWh0xbSv4ikjxm/L11a/5gxlFAWWUly2bN8jsqqYaSv4gkVSQS4alQiN6BAHsCw++9FxYt8jss2Y6Sv4gk1dblIU8oKeGjxx+ncbNm0KMHPPaY36HJNnSTl4gk3U9uBuvRA/r1g9NOY000yiP77UekqEjLQ/pMyV9EUmvPPWHePL486SRa3303zc04Lhhkbnm5TgA+UrOPiKReKMTUSISJZgx3jicqK4k9/7zfUeU0JX8RaRCRoiImhkIMN6M7MOqhh2DNGr/DyllK/iLSILZ2BLedOJH3bruNXb/+2psUTlO3+0Lz+YuIP955B0480ZsS+pFHoG9fvyPKCprPX0TS28EHw5IlcMgh3lrBt97qd0Q5RclfRPzTooU3LXS/fnDRRXDBBbBli99R5QQlfxHx1667wowZMGYM3H67dyJYt87vqLKekr+I+K9RI7j5ZrjjDm9VsO7d4ZNP/I4qqyn5i0j6OPdcmDMHli/3RgL9619+R5S1lPxFJL306QMLFnht/507884tt1BaWqqpoZNM0zuISPrp0AGWLOH7oiLajxnDajN6hkKUlZVpSogk0ZW/iKSnVq24e8gQ5gOTneOqykqi5eV+R5U1lPxFJG11Ou44BoZC3G3GWOcYUVEBGzb4HVZWUPIXkbQVDod5obycb0pKWHPeeexVVgY9e3p3BUtC1OYvImntJ2sDRCLwxz9COOwNCW3f3tfYMpmu/EUkcwwYAOXl8N133glgwQK/I8pYSv4iklnCYVi8GJo1g169vEnhpM7U7CMimeeAA7xF4U8+GYYM4cNolL+3aaPlIesgoeRvZk2Bx4A2wGpgkHPum2rK3QD0wfulMQ+40KXrXNIikhmaNoUXXuCLvn0pmDqVlloesk4SbfYZB5Q559oBZVXbP2FmnYDOwOHAocBRQPcEjysiAsEg93TrxjVmDHOOJysriT33nN9RZYREk39f4MGq1w8C/aop44AQkA8EgcbAZwkeV0QE8JaHLA2FGGZGN2DUww9rechaSDT5t3DOfQpQ9dx8+wLOuRhQAXxa9ZjrnHs3weOKiAA/Lg/ZXstD1kmNbf5mNh9oWc1b42tzADNrCxwEtKraNc/MujnnXqqm7EhgJEBBQUFtPl5E5Kf3AvTs6S0P2b07PPoonHSSv8GlqRqv/J1zvZxzh1bzmA18ZmZ7A1Q9f17NR5wMLHbOrXPOrQOeA47dwbGmOOcKnXOFzZo1q/+/SkRy18EHe0NBDznEWxjmttv8jigtJdrsMwcYWvV6KDC7mjIfAt3NLM/MGuN19qrZR0RSp2VLb3nIvn3hwgu9JSK1PORPJJr8rwOKzWwFUFy1jZkVmtk9VWVmAu8D/wKWAcucc08neFwRkZ3bdVeYOdNL/LfeCv37s6S8XGsDVLF0HW5fWFjoXlWHjYgkw6RJuAsv5HXnOMmMb4LBrF0bwMxec84V1lRO0zuISPY7/3xmDhnCb5xjYTxO240biUajfkflKyV/EckJrc45h97BIEHgpXick375S79D8pWSv4jkhHA4zE0VFTwxdiz5++/PIWPGwAMP+B2Wb5T8RSRnhMNhzrv+enZdutRbG2DYMJgwAdK07zOVNKuniOSeJk28xWDOOQdKSuCDD1g8YgQVixYRiUSysiN4e0r+IpKbGjeGqVO96aHHj+eHRx/lr8A1WTwSaFtq9hGR3GUGl13G7FNP5eh4nAXxOPvmyEggJX8RyXnNL7yQPvn5NAMWxuP8bq+9/A4p5ZT8RSTnhcNhrolGmTlmDLvvuy+HXXABzJrld1gppeQvIoJ3Ajj75psJLV0KHTrAwIFw001ZOxJIHb4iIttq1gzKymDoULjkEv7/okU80LEj3Xv2zKpOYF35i4hsb5ddYPp0Ph4yhJZPPsnhEyZwUo8eWTUhnJK/iEh1AgGmHXII55jRG3ihspJXZ1c3a31mUvIXEdmBSCTCg6EQ/QIB2gIj778f3nzT77CSQslfRGQHtq4P3LmkhFUPPkiwcWPo0gXmzvU7tIRpPn8Rkdpauxb69IG334a77oIRI/yO6Gc0n7+ISLK1agULFkBxMYwcCZddBvG431HVi5K/iEhd7L47PP00nH02lJbCH/4AlZV+R1VnGucvIlJXeXles8+vfw1jx8Latbwyfjzz33gjY2YFVfIXEakPM7jkEmjThvjpp7NHnz48YJYxs4Kq2UdEJBEDB/LwsGH8qmp94CMzZFZQJX8RkQS1GzqUSDDI18C8eJz+W7b4HVKNlPxFRBIUDoeZWlHBs5dfzqYjjqD9hAlwww1pPSmcxvmLiCRTZaW3NvD06d5w0Dvu8DqIG0htx/mrw1dEJJlCIfj7373lIa+9Fj78EB5/HHbbze/IfkLNPiIiyRYIwMSJ3hrB8+bx/ZFHMmncuLSaFVTJX0QkVc46i3duuoktK1Zw8vXX8/+KitLmBKDkLyKSQrM3bKB7IEAcmLdxIx9Onux3SICSv4hISkUiEZYHg3QKBFhlxqCHH4Y0OAGow1dEJIW2TgsdjUapPOoo7JZbYNQoWLXKmxso4M81eELJ38wGAlcBBwFHO+eqHZtpZscDtwKNgHucc9clclwRkUwSDod/nO4hEoHRo737AD74gCXnnkt5LNbgcwIleuX/FnAKsMPfMGbWCLgDKAbWAq+Y2Rzn3DsJHltEJPPk5cGdd3qTwl1yCW7WLG6FBp8TKKHfG865d51zy2sodjSw0jm3yjm3CZgO9E3kuCIiGc0MLr6YJwYP5oh4nH/G47Ru4DmBGqKxaV/go22211btExHJaXuPHs0J+fk0Af4Zj/P7PfZosGPX2OxjZvOBltW8Nd45V5ul7K2afdXOKWFmI4GRAAUFBbX4aBGRzBUOhymNRpk5axZnzpjBoRdeyHuff86sxo1T3gdQY/J3zvVK8Bhrgf222W4FfLKDY00BpoA3t0+CxxURSXv/6wz+85/5rmdP2l95JevM6BkMUlZenrITQEM0+7wCtDOz/c0sHzgNmNMAxxURyRx77snd/fvziBmFzvHDpk0p7QNIKPmb2clmthYIA/8ws7lV+/cxs2cBnHObgfOBucC7wOPOubcTC1tEJPt07dWLs4JBhgQCNA4GiUQiKTuWpnQWEUkjsViMaDRa7zZ/TeksIpKBfnJDWAppbh8RkRyk5C8ikoOU/EVEcpCSv4hIDlLyFxHJQUr+IiI5KG3H+ZvZF8Ca7XY3Ab6tpnh1+7ffV9P2XsCX9Qq2ZjuKO1l/t7Nydamz6vZna73VVCaR71p1+7bdTmWd7SieZP1NKustW79rNZVL9v/Rds65JjVG5JzLmAcwpbb7t99Xi+1XGzruZP3dzsrVpc5yqd5qKpPId62mektlndW33pLxXUu03rL1u1ZTuVT/H93RI9OafZ6uw/7t99W0nUr1PVZt/25n5epSZ9Xtz9Z6q6lMIt+16vale70l47u2s/f1f7Tu76X0/2jaNvs0NDN71dXilmj5KdVb3anO6kf1llyZduWfSlP8DiBDqd7qTnVWP6q3JNKVv4hIDtKVv4hIDlLyFxHJQUr+IiI5SMm/lszsF2b2mpn9zu9YMoGZHWRmd5vZTDM7x+94MoWZ9TOzqWY228x6+x1PpjCzA8zsXjOb6XcsmSLrk7+Z3Wdmn5vZW9vtP97MlpvZSjMbV4uPuhR4PDVRppdk1Jlz7l3n3ChgEJATw/OSVG9POedGAGcCp6Yw3LSRpHpb5ZwbntpIs0vWj/Yxs27AOmCac+7Qqn2NgPeAYmAt3iLzg4FGQOl2H/En4HC8W8tDwJfOuWcaJnp/JKPOnHOfm9lJwDhgknPukYaK3y/Jqreqv7sZ+Ltz7vUGCt83Sa63mc65AQ0VeybL+mUcnXMvmVmb7XYfDax0zq0CMLPpQF/nXCnws2YdMysCfgEcDGwws2edc/GUBu6jZNRZ1efMAeaY2T+ArE/+SfquGXAd8FwuJH4JDV6XAAABLElEQVRI3vdN6ibrk/8O7At8tM32WuCYHRV2zo0HMLMz8a78szbx70Sd6szMIsApQBB4NqWRpbc61RswGugFNDGzts65u1MZXBqr6/dtT2Ai0MHM/lx1kpCdyNXkb9Xsq7H9yzn3QPJDyRh1qjPnXBSIpiqYDFLXersNuC114WSMutbbV8Co1IWTfbK+w3cH1gL7bbPdCvjEp1gyheqsflRv9aN6S7FcTf6vAO3MbH8zywdOA+b4HFO6U53Vj+qtflRvKZb1yd/MHgViwIFmttbMhjvnNgPnA3OBd4HHnXNv+xlnOlGd1Y/qrX5Ub/7I+qGeIiLyc1l/5S8iIj+n5C8ikoOU/EVEcpCSv4hIDlLyFxHJQUr+IiI5SMlfRCQHKfmLiOQgJX8RkRz0f+tyPwWh+4rQAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmUU0X6//H3E5pO1FFGlEXFBh3AcR+0XcKaBtqNGVFZlEFFRBAX1B9HkVERRxsbt3FUVBZXdBQFVHBF6O6oAwE3xHVQRFDUn/uoCA1C6vvHbUbEhl6S9M3yeZ2Tk9yb6tyHOuG5N1V1q8w5h4iI5JaA3wGIiEjDU/IXEclBSv4iIjlIyV9EJAcp+YuI5CAlfxGRHKTkLyKSg5T8RURykJK/iEgOUvIXEclBeX4HsDW77rqra9Omjd9hiIhklNdee+1r51yzmsqlbfJv06YNr776qt9hiIhkFDNbWZtyavYREclBSv4iIjlIyV9EJAcp+YuI5KCkJH8zO8bMlprZMjMbXc37QTN7pOr9RWbWJhnHFRGR+kk4+ZtZI+B24FhgP2CAme23RbEhwHfOubbAzcB1iR53W2KxGKWlpcRisVQeRkQkYyVjqOfhwDLn3HIAM5sG9Abe3axMb+CqqtczgAlmZi4Fa0gufPFFbiku5ocNGyjPy+OmO+7goCOPhO228x7bb+89N24MZsk+vIhIRkhG8t8D+GSz7VXAEVsr45zbYGbfA7sAX29eyMyGAcMACgoK6hXMoueeY9r69d7G+vVw1lnVFwwEfjkRbHliSOD1kg8/ZM7y5XQpLiYcDtfr3yAikmrJSP7VXT5veUVfmzI45yYDkwEKCwvr9avgiGOP5ch//IO8n39mx7w8/jFuHPu2bg1r18KaNd7z5q+r27d2LXzzTfXv1/Bj5WBgf+A/V13FV0cdRbOjjoIOHeBPf4Kdd67PP0lEJOmSkfxXAXtutt0K+GwrZVaZWR7QBPg2Ccf+jSO7dOHmigqi0SiRSIR9k3n17Zz3a2IrJ41H7r2XWVOnsq9zHOIcrWMxmDPnl79v08Y7EWz+2H13NT+JSINLRvJ/BWhnZnsBnwKnAH/dosxsYBAQA/oC5alo798kHA6npsnFDIJB7/H73//m7YL8fJ549FEeXb+e/Px8yp57jvDee8Mbb8Dixb88Hn/8lz9q1uy3J4S2bb1mKRGRFLFk5GAzOw74J9AIuMc5N87MrgZedc7NNrMQ8ADQAe+K/5RNHcRbU1hY6DJxbp9YLPa/Xx1bPQH9+CMsWfLrE8I778DPP3vv77ADHHzwr08I++9P7PXXa/5sEclpZvaac66wxnIpvABPSKYm/3pbv947AWx+QliyBFavBiCel8fbGzdSATyWn8/48nLCHTv6G7OIpB0l/2wQj8OyZbB4MQvuvJOfXniBrkAQ+K5pU3YePhwGDIADDvA7UhFJE7VN/mpYTmeBALRvDyefjJWW0nu77dgtEGBY48ZY27YwfjwceKD3uPZaWL7NljQRkf/RlX8G+U1/wpdfwvTp8NBDsGCBV+iII2DAAF7de2/mvv22+gdEcoyafXLNypUwbRo8/DAsWcJGoBy4Iz+fUeXlhDt18jtCEWkAavbJNa1bw6WXwhtvMPmiiyg144/A4+vX0/akk+Cee2DdOr+jFJE0oeSfhQ7s359rQyHaBQKc2bgx2zdpAkOGeDeZlZbCd9/5HaKI+EzJPwuFw2HKysoYW1LC0BdeYIelS2HuXDjoILjsMthzT7jwQvjoI79DFRGfqM0/17z5Jtx0k9dJHI9D375wySXEfv5ZN5CJZAF1+Mq2rVoFt94KkybBDz/wYiDAjcC8/HzKyst1AhDJUOrwlW1r1Qquvx4++YR5xx1Hm3ic2fE4z1dW8v7UqX5HJyIppuSf63baiR2uuIIDQiHOMaMdMGjiROjXz7u7WESykpK/EA6HmVNeTsG4cayYNw/GjoVnnoF994ULLoCvv675Q0Qko6jNX6r3+edw1VVw113wu9/B6NFw0UXeqmUikrbU5i+J2W03rzP47behWzdviGj79nDffbBxo9/RiUiClPxl2/bdF2bPhmjUOyEMHgyHHMK7N99MaWkpsVjM7whFpB6U/KV2unWDhQth2jQqv/6a/UaO5IDLLuOMoiKdAEQykJK/1F4gACefzK3Dh3OpGd2B19etY8011/yyCpmIZAQlf6mzLj17clsoxIGBANFAgB7PPguFhd4vAxHJCEr+Umeb5g4aWlJC05degsceg2++gY4d4dxz4b//9TtEEamBhnpKcvz4I4wZA7fdBs2bwz//Cf37g5nfkYnkFA31lIa1445ewn/5ZdhjDzjlFDj2WF6fMUOjgkTSUJ7fAUiWOfRQWLQI7riDDaNHs++cOUw346hgkOc1YZxI2tCVvyRfo0YwYgR3jhjBc0Cpc5RVVvLWo4/6HZmIVFHyl5Qp7N2bgaEQA8zYGzjrjjvghht0h7BIGlDyl5QJh8OUlZdz0LhxrHjqKQK9esGoUdC5Myxd6nd4IjlNo32k4TgH06bB+efDmjWsGDqUaS1a0K17d/UFiCRJbUf7qMNXGo4ZDBgARUV8278/bW67jU7A4GCQeysqdAIQaUBq9pGG17Ilk445hkFmHAS8vG4dX998s/fLQEQahJK/+CJSVMT0UIgOgQBLAgH+Mn06nHwyfPut36GJ5AQ1+4gvNk0REY1GyevSBebP9+4Qnj+fdy+9lFk//UQkElFTkEiKJNTha2ZNgUeANsAKoL9z7rstyvwJuBPYCdgIjHPOPVLTZ6vDNwe9/jpr+vRh+xUruMWMscEgz+rGMJE6aajpHUYDZc65dkBZ1faW1gCnO+f2B44B/mlmv0/wuJKNDjmE2wcPZoIZFzpHRWUlb06f7ndUIlkp0eTfG7i/6vX9wAlbFnDOve+c+6Dq9WfAl0CzBI8rWapzcTGjQiFOCATYEzjrzjvh3nvVGSySZIkm/xbOuc8Bqp6bb6uwmR0O5AMfJnhcyVKb+gKOKClh5axZNAqH4cwzYeBA+P57v8MTyRo1tvmb2TygZTVvXQ7c75z7/WZlv3PO7byVz9kNiAKDnHPVrvphZsOAYQAFBQWHrly5sjb/BslmGzfCddfBlVdS2aIFjxx/PO1PP139ACJbUds2/0Q7fJcCEefc55uSu3Nun2rK7YSX+Eudc7VqxFWHr2zurcmT2Wn4cHZ3jivz8jg+GiXcqZPfYYmknYbq8J0NDKp6PQiYVU0g+cDjwNTaJn6RLT31zTccYsYsoHTDBnY96yxv9TARqZdEk/94oNjMPgCKq7Yxs0Izu6uqTH+gK3CGmb1R9fhTgseVHBOJRFgbDHJKIMDIxo35w4cfQocOsGCB36GJZCRN7CYZIxaLEY1GvZu/8vO9ZSJXrmTl2Wfz0O67E9EEcSIN0+afSkr+UqPvv+frk05i1/JyngGGhkLM0E1hkuO0hq9kvyZNmNKjB+eb0QOYX1nJfx580O+oRDKCkr9ktEhREfeEQnQLBAiYMWjKFJg0STeFidRAyV8y2qabwnqXlPDFM88Q6N4dhg+HM86ANWv8Dk8kbanNX7LLxo1QUgJ//zsccADMnAnt2vkdlUiDUZu/5KZGjWDsWHj2Wfj0UygshMce8zsqkbSj5C/Z6eijYfFi+OMfoU8fuPhiFr70EqWlpcRiMb+jE/Gdmn0ku61bByNHwh138O9AgFOAb4NBysrKNCRUspKafUQAgkG4/XZm9e9Ph3icV+JxCtetIxqN+h2ZiK+U/CUnNL/oIroFg6wGyuJxTv7mGw0HlZym5C85IRwOc1tFBU+OGcPqjh3Z+6abvHUC1q71OzQRX6jNX3JPPA5XX+0NBz30UG84aOvWfkclkhRq8xfZmkAArroKZs+GDz7whoOWl/sdlUiDUvKX3PWXv8Arr0Dz5lBczIoRIyi99loNBZWckOd3ACK+at8eFi7km969aTNhAnub8ZdgkCc1O6hkOV35i+y4I5N79uQyM/o5x7zKShY//rjfUYmklJK/CN7soP8MhfhLIEAbYOjkyeoHkKym5C/CL7ODdi4p4aNHHqHx7rvDUUfBLbfofgDJShrqKVKdH36A007zRgQNGsTCM86gIhbzlpBUX4CkMS3jKJKoze4HeNWMPmZ8pXmBJM1pnL9IoqruB5hx6qns4xwL43H+tHYtU6dO9TsykYQp+YvUYI9zz6VrXh5rgAogftdduhdAMp6Sv0gNwuEwR551FocDLwKTNmwg7+KLYcMGv0MTqTclf5FaOP3001m73Xb8ORBgQqNGHLZgARx7LK/MmaMFYiQj6Q5fkVrYNBQ0Go1yaCQC771H/OyzaVpWxkNmXKOOYMkwSv4itRQOh39J7uEwD778MkdNmsR85zi1aoEYJX/JFGr2EamndoMG0SUYZBnwRDzOgM8+0w1hkjF05S9ST+FwmKkVFZQ//zytX3iBNhMmwI8/wqRJ3vKRImlMyV8kAf9rCorH4ZprvHUC3n8fHnsMWrb0OzyRrVKzj0gyBAIwdixMnw5vvAGHHw6LF/sdlchWJZT8zaypmc01sw+qnnfeRtmdzOxTM5uQyDFF0lrfvjB/vtf237mzt0SkSBpK9Mp/NFDmnGsHlFVtb801wAsJHk8k/XXo4K0QdtBB0LcvHw8bphXCJO0k2ubfG4hUvb4fiAKXblnIzA4FWgDPATVOOCSS8Vq2hIoKvjzpJAqmTOEPZvw5GOQprRAmaSLRK/8WzrnPAaqem29ZwMwCwE3AJQkeSySzhELc3bkzl5rR1zmeq6zk1dmz/Y5KBKjFlb+ZzQOqG7ZweS2PcS7wjHPuEzOr6VjDgGEABQUFtfx4kfQVKSqiRyjE0nXreCAe56C774Z+/eCQQ/wOTXJcQvP5m9lSIOKc+9zMdgOizrl9tijzL6ALEAd+B+QDdzjnttU/oPn8JWvEYjGi0SjHtWrFwVdcAV99BQ88AH36+B2aZKEGWczFzG4AvnHOjTez0UBT59yobZQ/Ayh0zp1f02cr+UtW+uILOPFEiMWgpAQuuwxq+EUsUhcNtZjLeKDYzD4Aiqu2MbNCM7srwc8WyT4tWngLw596KlxxhbdUZGWl31FJDtIyjiJ+cA7Gj/eu/I88Ep54wjsxiCRIyziKpDMz+NvfvJvA3nwTDjsMlizxOyrJIUr+In466ST497+9XwKdOvGf66/X4jDSIDSxm4jfOnSAl1/mx549aX/ppfzXjB7BIGW6IUxSSFf+Iulgt92YePLJTDfjOueYsG4dL86b53dUksWU/EXSROfiYs4MBikx40znGP7EE/Dtt36HJVlKyV8kTYTDYeaVl9No3Dg+GDuWJm+/7Y0Eev99v0OTLKShniLpav5874awDRu8xWEiEb8jkgygoZ4ima5TJ1i0CHbbDYqL4Z57/I5IsoiSv0g622svWLAAuneHIUNg1ChvyUiRBCn5i6S7Jk3g6afh3HPhhhu8CeF++snvqCTDKfmLZIK8PJgwAW69FWbPhi5d4NNP/Y5KMpiSv0imMIMRI+DJJ2HZMm+R+Nde8zsqyVBK/iKZ5rjjvJFAeXnQtSs8/rjfEUkGUvIXyUQHHggvv+w9n3QSFcccQ2zBAr+jkgyi5C+SqVq0YNH48Uxv1IiiOXN4v2tXFr74ot9RSYZQ8hfJYOWxGAOc42pg0MaNtBw8GL77zu+wJAMo+YtksEgkQn4wyNWNGjGkcWMKPvkEwmGvQ1hkGzSls0gGC4fDlJWVEY1GiUQiBDZs8KaE2LQ6WOfOfocoaUpz+4hkm2XLoFcvWLGCD0aPZkYoRCQS0doAOaK2c/voyl8k27RtC7EY3xcX0+7qq9mgxWGkGmrzF8lGTZsy8YQTuNeMMc5x97p1vDR3rt9RSRpR8hfJUl179uS8YJDLzBjgHMNnzoSvvvI7LEkTSv4iWSocDlNWXs6O48axtKSEnd5/3+sI/s9//A5N0oA6fEVyxcKF0Ls3rF8PM2d600RL1tFiLiLya0ce6S0Os/vucPTRWhwmx2m0j0guadPGWxymXz8YMoRPo1Gm7rMPke7dNRIox+jKXyTXVC0O80Xv3uzxwAO0HTOGXt27E4vF/I5MGpCSv0guatyYew4/nFFm9HGOZysrefnJJ/2OShqQkr9IjooUFTEhFKJfIMCBwNn33Qfvvut3WNJA1OYvkqM2nxdoWfPmHHTFFdCxI+9efTWzfvpJU0JkuYSGeppZU+ARoA2wAujvnPvNfLJmVgDcBewJOOA459yKbX22hnqKNLCPP2ZNURGNly/nPDMeDIUoKyvTCSDDNNRQz9FAmXOuHVBWtV2dqcANzrl9gcOBLxM8rogkW0EBd556KuXAZOcYW1lJtLzc76gkRRJN/r2B+6te3w+csGUBM9sPyHPOzQVwzq12zq1J8LgikgIdjzmGvqEQE8241DmGlpfD2rV+hyUpkGjyb+Gc+xyg6rl5NWXaA/81s8fMbLGZ3WBmjRI8roikQDgc5vnycr4rKWHFiBHsWlEBRUXwpX6sZ5saO3zNbB7Qspq3Lq/DMboAHYCP8foIzgDuruZYw4BhAAUFBbX8eBFJpnA4/Es7f1ERDBwIRxwBTz8N++3nb3CSNDVe+TvnejrnDqjmMQv4wsx2A6h6ru7yYBWw2Dm33Dm3AXgCOGQrx5rsnCt0zhU2a9as/v8qEUmOE0+EF16Aykro2JF3bruN0tJS3RCWBRJt9pkNDKp6PQiYVU2ZV4CdzWxTNu8OaDCxSKY47DBYuJA1u+xC+wsu4MPLL6dHjx46AWS4RJP/eKDYzD4Aiqu2MbNCM7sLwDm3EbgYKDOztwADpiR4XBFpSK1bM/HUU6kA7nKOKzUSKONpSmcRqZVYLMbR3btz47p1DHOOr3v0YNennoJQyO/QZDOa0llEkiocDjOnvJxvSkpYed557FpWBj16aHWwDKXpHUSk1n41EigSgdNO89YJeOYZ2GcfX2OTutGVv4jUT9++EI3C6tUQDnujgiRjKPmLSP0dcYS3PGTLllBcDA884HdEUktK/iKSmL32gvnzoUsXOP10uOoqSNOBJPILJX8RSdzOO8Ozz8LgwfD3v/PVscdy/TXX6F6ANKbkLyLJkZ8Pd9/Nx2efTbM5cwhfeSX9tDxk2lLyF5HkMeNfrVvzVzMOA8oqK1kyc6bfUUk1lPxFJKkikQhPhEIcFQiwCzDk7rthwQK/w5ItKPmLSFJtWh7y2JISPnn0URo3awbdu8Mjj/gdmmxGN3mJSNL96maw7t3hhBPglFNYGY3y0J57Eikq0vKQPlPyF5HU2mUXmDuXr48/ntYTJ9LcjKODQeaUl+sE4CM1+4hI6oVCTIlEGGfGEOd4rLKS2HPP+R1VTlPyF5EGESkqYlwoxBAzugHDH3gAVq70O6ycpeQvIg1iU0dw23HjeP/WW9n+22+9SeE0dbsvNJ+/iPjj3XfhuOO8KaEfegh69/Y7oqyg+fxFJL3ttx8sWgT77++tFXzLLX5HlFOU/EXEPy1aeNNCn3ACXHQRXHABbNzod1Q5QclfRPy1/fYwfTqMHAm33eadCFav9juqrKfkLyL+a9QIbroJbr/dWxWsWzf47DO/o8pqSv4ikj7OPRdmz4alS72RQG+95XdEWUvJX0TSS69e8NJLXtt/p068e/PNlJaWamroJNP0DiKSfjp0gEWL+KmoiPYjR7LCjB6hEGVlZZoSIkl05S8i6alVKyYOHMg8YJJzXFVZSbS83O+osoaSv4ikrY5HH02/UIiJZoxyjqEVFbB2rd9hZQUlfxFJW+FwmOfLy/mupISV553HrmVl0KOHd1ewJERt/iKS1n61NkAkAqedBuGwNyS0fXtfY8tkuvIXkczRty+Ul8MPP3gngJde8juijKXkLyKZJRyGhQuhWTPo2dObFE7qTM0+IpJ59t7bWxT+xBNh4EA+jkb5V5s2Wh6yDhJK/mbWFHgEaAOsAPo7576rptz1QC+8XxpzgQtdus4lLSKZoWlTeP55vurdm4IpU2ip5SHrJNFmn9FAmXOuHVBWtf0rZtYR6AQcBBwAHAZ0S/C4IiIQDHJX165cY8Zg53i8spLYs8/6HVVGSDT59wbur3p9P3BCNWUcEALygSDQGPgiweOKiADe8pCloRCDzegKDH/wQS0PWQuJJv8WzrnPAaqem29ZwDkXAyqAz6sec5xz7yV4XBER4JflIdtrecg6qbHN38zmAS2reevy2hzAzNoC+wKtqnbNNbOuzrkXqyk7DBgGUFBQUJuPFxH59b0APXp4y0N26wYPPwzHH+9vcGmqxit/51xP59wB1TxmAV+Y2W4AVc9fVvMRJwILnXOrnXOrgWeBI7dyrMnOuULnXGGzZs3q/68Skdy1337eUND99/cWhrn1Vr8jSkuJNvvMBgZVvR4EzKqmzMdANzPLM7PGeJ29avYRkdRp2dJbHrJ3b7jwQm+JSC0P+SuJJv/xQLGZfQAUV21jZoVmdldVmRnAh8BbwBJgiXPuyQSPKyKybdtvDzNmeIn/llugTx8WlZdrbYAqlq7D7QsLC92r6rARkWSYMAF34YW87hzHm/FdMJi1awOY2WvOucKayml6BxHJfuefz4yBA/mjc8yPx2m7bh3RaNTvqHyl5C8iOaHVOedwVDBIEHgxHuf43/3O75B8peQvIjkhHA5zY0UFj40aRf5ee7H/yJFw331+h+UbJX8RyRnhcJjzrruO7Rcv9tYGGDwYxoyBNO37TCXN6ikiuadJE28xmHPOgZIS+OgjFg4dSsWCBUQikazsCN6Skr+I5KbGjWHKFG966Msv5+eHH+YfwDVZPBJoc2r2EZHcZQaXXcask0/m8Hicl+Jx9siRkUBK/iKS85pfeCG98vNpBsyPx/nzrrv6HVLKKfmLSM4Lh8NcE40yY+RIdtpjDw684AKYOdPvsFJKyV9EBO8EcPZNNxFavBg6dIB+/eDGG7N2JJA6fEVENtesGZSVwaBBcMkl/P8FC7jvkEPo1qNHVnUC68pfRGRL220H06bx6cCBtHz8cQ4aM4bju3fPqgnhlPxFRKoTCDB1//05x4yjgOcrK3l1VnWz1mcmJX8Rka2IRCLcHwpxQiBAW2DYvffCm2/6HVZSKPmLiGzFpvWBO5WUsPz++wk2bgydO8OcOX6HljDN5y8iUlurVkGvXvDOO3DnnTB0qN8R/Ybm8xcRSbZWreCll6C4GIYNg8sug3jc76jqRclfRKQudtoJnnwSzj4bSkvhr3+Fykq/o6ozjfMXEamrvDyv2ecPf4BRo2DVKl65/HLmvfFGxswKquQvIlIfZnDJJdCmDfFTT2XnXr24zyxjZgVVs4+ISCL69ePBwYP5fdX6wIdmyKygSv4iIglqN2gQkWCQb4G58Th9Nm70O6QaKfmLiCQoHA4zpaKCZ664gvUHH0z7MWPg+uvTelI4jfMXEUmmykpvbeBp07zhoLff7nUQN5DajvNXh6+ISDKFQvCvf3nLQ157LXz8MTz6KOy4o9+R/YqafUREki0QgHHjvDWC587lp0MPZcLo0Wk1K6iSv4hIqpx1Fu/eeCMbP/iAE6+7jv9XVJQ2JwAlfxGRFJq1di3dAgHiwNx16/h40iS/QwKU/EVEUioSibA0GKRjIMByM/o/+CCkwQlAHb4iIim0aVroaDRK5WGHYTffDMOHw/Ll3txAAX+uwRNK/mbWD7gK2Bc43DlX7dhMMzsGuAVoBNzlnBufyHFFRDJJOBz+ZbqHSARGjPDuA/joIxadey7lsViDzwmU6JX/28BJwFZ/w5hZI+B2oBhYBbxiZrOdc+8meGwRkcyTlwd33OFNCnfJJbiZM7kFGnxOoIR+bzjn3nPOLa2h2OHAMufccufcemAa0DuR44qIZDQzuPhiHhswgIPjcf4dj9O6gecEaojGpj2ATzbbXlW1T0Qkp+02YgTH5ufTBPh3PM5fdt65wY5dY7OPmc0DWlbz1uXOudosZW/V7Kt2TgkzGwYMAygoKKjFR4uIZK5wOExpNMqMmTM5Y/p0DrjwQt7/8ktmNm6c8j6AGpO/c65ngsdYBey52XYr4LOtHGsyMBm8uX0SPK6ISNr7X2fw3/7GDz160H7sWFab0SMYpKy8PGUngIZo9nkFaGdme5lZPnAKMLsBjisikjl22YWJffrwkBmFzvHz+vUp7QNIKPmb2YlmtgoIA0+b2Zyq/bub2TMAzrkNwPnAHOA94FHn3DuJhS0ikn269OzJWcEgAwMBGgeDRCKRlB1LUzqLiKSRWCxGNBqtd5u/pnQWEclAv7ohLIU0t4+ISA5S8hcRyUFK/iIiOUjJX0QkByn5i4jkICV/EZEclLbj/M3sK2DlFrubAN9XU7y6/Vvuq2l7V+DregVbs63Fnay/21a5utRZdfuztd5qKpPId626fZtvp7LOthZPsv4mlfWWrd+1msol+/9oO+dckxojcs5lzAOYXNv9W+6rxfarDR13sv5uW+XqUme5VG81lUnku1ZTvaWyzupbb8n4riVab9n6XaupXKr/j27tkWnNPk/WYf+W+2raTqX6Hqu2f7etcnWps+r2Z2u91VQmke9adfvSvd6S8V3b1vv6P1r391L6fzRtm30ampm96mpxS7T8muqt7lRn9aN6S65Mu/JPpcl+B5ChVG91pzqrH9VbEunKX0QkB+nKX0QkByn5i4jkICV/EZEcpORfS2a2g5m9ZmZ/9juWTGBm+5rZRDObYWbn+B1PpjCzE8xsipnNMrOj/I4nU5jZ3mZ2t5nN8DuWTJH1yd/M7jGzL83s7S32H2NmS81smZmNrsVHXQo8mpoo00sy6sw5955zbjjQH8iJ4XlJqrcnnHNDgTOAk1MYbtpIUr0td84NSW2k2SXrR/uYWVdgNTDVOXdA1b5GwPtAMbAKb5H5AUAjoHSLjzgTOAjv1vIQ8LVz7qmGid4fyagz59yXZnY8MBqY4Jx7qKHi90uy6q3q724C/uWce72BwvdNkutthnOub0PFnsmyfhlH59yLZtZmi92HA8ucc8sBzGwa0Ns5Vwr8plnHzIqAHYD9gLVm9oxzLp7SwH2UjDqr+pzZwGwzexrI+uSfpO+aAeOBZ3Mh8UPyvm99Z2tKAAABKUlEQVRSN1mf/LdiD+CTzbZXAUdsrbBz7nIAMzsD78o/axP/NtSpzswsApwEBIFnUhpZeqtTvQEjgJ5AEzNr65ybmMrg0lhdv2+7AOOADmb2t6qThGxDriZ/q2Zfje1fzrn7kh9KxqhTnTnnokA0VcFkkLrW263ArakLJ2PUtd6+AYanLpzsk/UdvluxCthzs+1WwGc+xZIpVGf1o3qrH9VbiuVq8n8FaGdme5lZPnAKMNvnmNKd6qx+VG/1o3pLsaxP/mb2MBAD9jGzVWY2xDm3ATgfmAO8BzzqnHvHzzjTieqsflRv9aN680fWD/UUEZHfyvorfxER+S0lfxGRHKTkLyKSg5T8RURykJK/iEgOUvIXEclBSv4iIjlIyV9EJAcp+YuI5KD/A3ecPwMr21qgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -216,8 +228,8 @@ "output_type": "stream", "text": [ "correlation matrix\n", - "[[ 1. -0.86793949]\n", - " [-0.86793949 1. ]]\n" + "[[ 1. -0.86628016]\n", + " [-0.86628016 1. ]]\n" ] } ], @@ -226,6 +238,190 @@ "print(cal.cormat)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calibrate parameters in multiple layers\n", + "Example showing how parameters can be optimized when multiple layers share the same parameter value." + ] + }, + { + "cell_type": "code", + "execution_count": 118, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "self.neq 1\n", + "solution complete\n" + ] + } + ], + "source": [ + "ml = ModelMaq(kaq=[10., 10.], z=(-10, -16, -18, -25), c=[10.], Saq=[0.1, 1e-4], tmin=1e-5, tmax=1)\n", + "w = Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, 788)], layers=1)\n", + "ml.solve()\n", + "hobs0 = ml.head(robs, 0, tobs, layers=[0])[0]\n", + "hobs1 = ml.head(robs, 0, tobs, layers=[1])[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 119, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "............................................................................\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
optimalstdperc_stdpminpmaxinitialparray
kaq0-110.0001611.545898e-050.0001550.0000030.010.0000[10.00016139633568, 10.00016139633568]
Saq00.1000002.346147e-070.0002350.000010.20.0100[0.09999972370676695]
Saq10.0001009.612878e-110.0000960.000010.20.0001[0.00010000052418362365]
c110.0000555.080174e-060.0000510.10000200.01.0000[10.000054570919623]
\n", + "
" + ], + "text/plain": [ + " optimal std perc_std pmin pmax initial \\\n", + "kaq0-1 10.000161 1.545898e-05 0.000155 0.00000 30.0 10.0000 \n", + "Saq0 0.100000 2.346147e-07 0.000235 0.00001 0.2 0.0100 \n", + "Saq1 0.000100 9.612878e-11 0.000096 0.00001 0.2 0.0001 \n", + "c1 10.000055 5.080174e-06 0.000051 0.10000 200.0 1.0000 \n", + "\n", + " parray \n", + "kaq0-1 [10.00016139633568, 10.00016139633568] \n", + "Saq0 [0.09999972370676695] \n", + "Saq1 [0.00010000052418362365] \n", + "c1 [10.000054570919623] " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "cal = Calibrate(ml)\n", + "cal.set_parameter(name='kaq0-1', initial=10., pmin=0., pmax=30.) # layers 0 and 1 have the same k-value\n", + "cal.set_parameter(name='Saq0', initial=1e-2, pmin=1e-5, pmax=0.2)\n", + "cal.set_parameter(name='Saq1', initial=1e-4, pmin=1e-5, pmax=0.2)\n", + "cal.set_parameter(name='c1', initial=1., pmin=0.1, pmax=200.)\n", + "cal.series(name='obs0', x=robs, y=0, layer=0, t=tobs, h=hobs0)\n", + "cal.series(name='obs1', x=robs, y=0, layer=1, t=tobs, h=hobs1)\n", + "cal.fit(report=False)\n", + "display(cal.parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": 121, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 121, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FPX5wPHPs0kgBAgknIEEwhGRK3KEAApI5bRawFuLFVCKteLP1uqrtFWBVi2ttlprW8UTK1brwWWrIggClisRRE5BzgBCCAjEECDZ5/fHLDGEzblJ9nrer9e+dmb2m/k+OyzPd+Y7M98RVcUYY0x4cfk7AGOMMbXPkr8xxoQhS/7GGBOGLPkbY0wYsuRvjDFhyJK/McaEIUv+xhgThiz5G2NMGLLkb4wxYciSvzHGhKFIfwdQmqZNm2pycrK/wzDGmKCSmZl5RFWblVcuYJN/cnIyGRkZ/g7DGGOCiojsqUg56/YxxpgwZMnfGGPCkCV/Y4wJQ5b8jTEmDFVL8heRkSKyTUR2iMgUL5/XFZE3PZ+vFpHk6qjXGGNM1fic/EUkAvgbcCXQBbhFRLqUKHYHcExVOwJPAn/wtd6ybF27iJWzfs3WtYtqsppar7sy687cc4y/LdlB5p5jFVp3TZUvq1xpn1XX8srGXNltYEwwq45LPdOBHaq6E0BE3gBGA5uLlRkNTPNMvw08IyKiNfAMyS/+91/y33+YOigndyxk0ZZraJrQrrqr8erIwV002j7nvLqbtGoPyHeFxJnWomVStPzCZUV/RPb+r2i87S1cqhzZvoQ5X+2iWfseuCPqUhgRjTsiGrerDojw1eFc3v3oYzrqPuZJEtcOu4IOzRuUGvdXh3P500fbKChUIiOEXwzrVC3lyypX2mcXLB/eiY7NGrAjO5c/LSx/+f3DO9GxeYNzm9mz9YTth3N5/MOtReV+OfJiTzlBgO2Hc5nx/paiz3/9/c50atGw6N9h+6Fclm3PRoBBFzXj4paxiMCXh06y6cAJurduROeEWFzi1Cfi/FO7RNhy8ARf7D9Oj6TGdG/dCJcILnHKuFyCy1POJcIX+4/z2Z5jpLeLp3fbOCJcQoRLiHQJUvxLFZO55xirdubQr30TereNq/BnpX1e3t+Y0CC+5l8RuR4YqaoTPfM/Avqq6uRiZTZ6ymR55r/ylDlSYl2TgEkAbdq06b1nT4UuVz3PxzPv54oDz1f16wS9fI0inzrkU4cTGsMuTWCHtuIrdyt2aGt2agK5xPg7TFMFIhDpchqJc42CKuSeLigq0zI2mpi6EUS5XJwpdLM751tUwSWQmtiYpg3qEBXhIirCxYn8syz/MptChQgXjO7RmqgIF+9kZlHoViIihLsu78DFLWOJjnIRHRXB7pxv2fb1Sfokx5PeLp6YOhHE1IkkwnV+w1RaA2INS80TkUxVTSuvXHXs+XvbHSnZolSkDKo6E5gJkJaWVqVWqXn3IWze9wERuCnExZGek2nVoVtVVlVpB77aSNN1zxSr+25ate/63ZcvamjVy7wzLaWUObhrM002zCSSQhTheOebaZmQhBTm4yrI97yfRgrzOb33M+od3UpLOUYbOcww+QwX7qI4z9ZvyenGHTndqCOn4zryZWECL63YSUs9zBbpyG1jrqJzQmyp33PLwRP8Zu4XFBQokZHCo2O6ey1fVrnSPiu5/JEx3encMpYtX5/gwfKWRwi/G9Odi1s2LIrh3NbcevAED8/bWLRnP31UNy5q2bCoxNaDJ5m+YFPR5w9e3YVOLRqiwJx1+3lz7b7zvtsNaYmg8HZmForTfzqmV2uu6p6AW0FVcSu8v/Eg89cfQHH+E3y/ewJDuzTH7Qa3KqrOe6EqS7YeZvGWw0VlB13UjL7t4yksdD4vdJ//KnArn+/7hnX7vnF+O0Bc/SjaN2tAQaGbHYdzi35iboVDJ/I5W+j2vJSc3NMUej4vdDvfs/i+YEGh8tePd3j9Dby68vwdszqRLqchiIpARDhw/FRRozMwpRltm8SQm1/A/M8PUOh2tvFvrupCenI8cfWjiIupQ3RUhNe6rMGoGdWR/LOApGLzicCBUspkiUgk0Ag4Wg11X6Bb/yvZGhnF0c0fE9flCgb1GVoT1XjV8ZIBbE3qWCN1t08bztbkbmR71p1W1rr3rcH9yg+g8CxEROG6bT7ENIUjX8KRbUQd2U5U9jYa7HgHzpykFTA40tMESSSuRh2g9bBSV9+tdSPaN2tQ7n/IssqV9llpy7snNqJDJZaX1COpMSktGpZarnfbeC5OiPX6eVSEizmfZXHGkynrRLq4uU8bABZsOMDZAjdRkS7G9m17wXqbNazLh5u+Lipz+4B2pcZ4cctYVuw4UlT2/4aklJvsMvccY+wLq4r+5pEx3c/rvin+2TM/7HXBXnjxz2dP7EehW/nRi6s4W6BERQp/vrEnHZs3IP9sIW+s3csba/YVNU5Du7Sgb7t48s4Uel4F5J0p5Ius4+c1Ouv3fcPnWd9w/NTZouVnC5Vp8zed912io1zExdShcUwdGteLIq5+FAWFyuKth3C7KWq0B13UlOYNo6kT6Trvu1gDUTnV0e0TCXwJDAH2A2uBH6rqpmJl7ga6q+pPRORm4FpVvbGs9aalpakN7+CDfWtg93JIHghJ6d7LqMLJg7DkUVg3m6L9ZFck9J4A/e+G+No5XxLoMvcc453PshDg2l6Jleofr0xiqkoSq60+f2+NRWnr9FYuY/dRz3Ln6O033+9C84Z1OZZ3lmN5Z/gm7wzH8s4WvR/LO8PXx/PJO1Po9Xs3qV+H5rHRREe52LDvG9zqNBC/G92N4V1bEhcTVep5klBW0W4fn5O/p7LvA08BEcBLqvqoiPwWyFDV+SISDfwT6Imzx3/zuRPEpbHkX4v2rYFZo6DwjJP42w+Grz4GLYQuY+Cy/4NWPf0dpQkAFW2cqqvPP3PPMX74/MqiBuOB4RcTWy+SQydOc+hEPodOnGbj/uN8fSL/gr9tWDeSNk1iaNskhjbx9WkTH8PpgkIOnchnaOcWpCXHV20jBLhaTf41wZJ/LSt5pHDiAKx+FjJehtMnoN0guOxe6DAEwnBvyvhPRY5exr6wijNn3URGCvcN7URUpIu9Od+y52gee3Py2Hcsj7OF5+e6zgmxpLWNo1PLhlzcsqHTZbX/eNB3HVnyN9Uj/zhkzoJVf3e6iFp0g0v/Dxonwd6VZXcrGVNLymsgCt3KHz/YysxlO4vOWSQ0jubkqQJOFrtaCpyT1OMvbccPLkmgS6tYNu4/EVTnEyz5m+pVcAa+eAv+9zRkb8X57yMQWRfGzbcGwAQ8b+cierVpzIHj+fxp4TbmfLb/gksQI11CoeeqrKgIYdaEdC7t2NQv8VeUJX9TM9xumHcXfP6GZ4HAkIdg4C/8GpYxFVHRk9lP39wTtyovfbqbNbu+uzAxwiX0ax/PwJRmDOjYlC4JsbhcgdUNWpvX+Ztw4nJB2h2wcS4U5gMKeTVy1a4x1a532zivXTe928Yxe2K/CxqGZg2ji84nREQII7u15Muvc5nx/lbAueLoso5NGZjSlEb1oth+ODd4uodsz99Uyb41zhVBO5c6ff+XT4HBU+xksAk53o4WDp3IZ8X2I6zYcYTl27M5knumqHykS5j5o95c0bmFX+K1bh9TOwoLYMG9sP416DMRrnzcOTowJky43cr0BZt4deWeonMGIjAopRmje7RieNeWbPv6ZK2dNLZuH1M7IiJh9DMQE++cDM47Ctc8B5F1/B2ZMbXC5RJG9WjNmxn7OFvgdA9d3b0Vq3cd5b5/f06dyC8oKHSjCnUihdd/3D8guoUs+RvficDw30H9pvDRw5D/Ddz0GtSp7+/IjKkV3s4ZuN3KZ3uP8ch/trDeM/7S6QLl70t38PTNPalf17/p17p9TPVa9xrMvwda9YKxbzlHBMaEsXN3KZ8pUBBnVJWG0ZHckt6G2/q35dCJ09XaJWR9/sZ/tv4H3poAcckwdKpzX4DdDGbCWPGTxiLw0opdvL/xa87lX1WoG1X6eEmVYcnf+NfuFTD7Bjh7CrsZzJgLHfjmFD97c/159xFMuCyZq1Nb+XQkUNHkb5dlmJqRPAB6jsUZKdQNBaedsYOMMQC0alyPX468mLqRUvTMj1c+3c0Nz/6PJz7cxtgXVtXoI0Ut+Zua0/1GiKjrTItA2wH+jceYANO7bRyv/7g/94/oxCsT+tCnXbzzMCDgzFk3q3bm1FjdlvxNzUlKh/HvOSOBaiFkb/F3RMYEnN5t47j7ex0Z3Kl50ZGAC6gT5aJf+yY1Vq/1+Zua53bDa9c6dwJPWgrNO/s7ImMClq9PJbM+fxM4XC7nxq+6DZ2rgM6e8ndExgSsc0cCNX0jmCV/UzsatnAagOwt8MGv/B2NMWHPkr+pPR2HOE8Dy3wZNs3xdzTGhDVL/qZ2XfEQtE6D+ffCsT3+jsaYsGXJ39SuiCi4/kVA4fWb4JM/OsNDG2NqlU/JX0TiReQjEdnuefd6hkJEPhCRb0TkPV/qMyEiLhkuvcfp/1/yGMwaZQ2AMbXM1z3/KcBiVU0BFnvmvXkc+JGPdZlQ4orAeQ6w2t2/xviBr8l/NDDLMz0LGOOtkKouBk76WJcJJckDi93965k3xtQaX5N/C1U9COB5b+57SCYsJKXD+AWeu3/dcPqEvyMyJqyUm/xFZJGIbPTyGl3dwYjIJBHJEJGM7Ozs6l69CTRJ6XDLvyC+Pbw/BQrOlP83xphqUW7yV9WhqtrNy2secEhEEgA874d9CUZVZ6pqmqqmNWvWzJdVmWARWRdG/B5ytsOamf6Oxpiw4Wu3z3xgnGd6HDDPx/WZcHTRCOg4DD75A+T6tP9gjKkgX5P/DGCYiGwHhnnmEZE0EXnhXCERWQ68BQwRkSwRGeFjvSaUiMDI3ztj/iye7u9ojAkLPj1BWFVzgCFelmcAE4vN26UcpmxNU6DfXfC/p6H37ZDY298RGRPS7A5fEzgGPQANWsC8n8KyJ+zGL2NqkCV/EziiY6HXbc4D3z9+xO78NaYGWfI3gSUy2jNhd/4aU5Ms+ZvA0m4QuOo4064Iu/PXmBpiyd8ElnPP/W2UBNGNIOESf0dkTEiy5G8CT5u+MOqvkHcEPnvV39EYE5Is+ZvA1H4wJPWDFU86ff/GmGplyd8EJhEY/Es4sR/W/dPf0RgTciz5m8DV/nuQmA7Lbe/fmOpmyd8ELhEYPAVOZMG61/wdjTEhxZK/CWwdroDEPrD8zzbkszHVyJK/CWzF9/7X296/MdXFkr8JfB2GQOs02/s3phpZ8jeBTwQG/wqO74O3x9t4P8ZUA0v+JjjUbQgIbP0PvPIDawCM8ZElfxMc9qwAxJkutAHfjPGVJX8THJIHQkSd7+bbDvBfLMaEAEv+JjgkpcP4BXDRlYD6Oxpjgp4lfxM8ktLh+hchujGsfMbf0RgT1Cz5m+BSpz6kTYCt78HRXf6OxpigZcnfBJ/0SSAuWP2cvyMxJmj5lPxFJF5EPhKR7Z73OC9leojIShHZJCIbROQmX+o0hthW0O06Z7TPU9/4OxpjgpKve/5TgMWqmgIs9syXlAfcpqpdgZHAUyLS2Md6Tbjr91M4k2sPezGminxN/qOBWZ7pWcCYkgVU9UtV3e6ZPgAcBpr5WK8Jd616OJd/rn4OCgv8HY0xQcfX5N9CVQ8CeN6bl1VYRNKBOsBXpXw+SUQyRCQjOzvbx9BMyOt/tzPg25Z5/o7EmKBTbvIXkUUistHLa3RlKhKRBOCfwARVdXsro6ozVTVNVdOaNbODA1OOlBEQ3wH+9wyoXftvTGVElldAVYeW9pmIHBKRBFU96Enuh0spFwv8B3hQVVdVOVpjinO5oP9P4T+/gPd+Bj3GOvcCGGPK5Wu3z3xgnGd6HHDB8beI1AHmAK+q6ls+1mfM+ZqkOO+Zr8CsUTbgmzEV5GvynwEME5HtwDDPPCKSJiIveMrcCAwCxovIes+rh4/1GuPYn0HRgG8FNuCbMRVVbrdPWVQ1BxjiZXkGMNEz/Rpgj2AyNePcgG+Fp50bv5IH+jsiY4KC3eFrgltSOox/D5p0gHqNoFVPf0dkTFCw5G+CX1I6DHsE8nLgyw/8HY0xQcGSvwkNKcMhtjWsfdHfkRgTFCz5m9AQEQm9x8POJZDj9R5CY0wxlvxN6Oj5I5AIyHzZ35EYE/As+ZvQEZsAF18F62bD2Xx/R2NMQLPkb0JL2u1w6ihstvF+jCmLJX8TWtpd7oz3k/GSvyMxJqBZ8jehxeVy9v73rYJDm/wdjTEBy5K/CT09fggRdW3v35gyWPI3oScmHrpdC5+/Cadz/R2NMQHJkr8JTWm3w5mT8IUNJGuMN5b8TWhK7AMtukPGi/agF2O8sORvQpMIpE2Ar7+A/Zn+jsaYgGPJ34Su1BuhTgMb78cYLyz5m9BVt6HTAGx6F/KO+jsaYwKKJX8T2tJuh4J8+Pxf/o7EmIBiyd+EtpbdITHduebfTvwaU8SSvwl9fe6AnB2wa5m/IzEmYFjyN6GvyxioF2d3/BpTjE/JX0TiReQjEdnueY/zUqatiGSKyHoR2SQiP/GlTmMqLSoaeoyFLQtg0TTYt8bfERnjd77u+U8BFqtqCrDYM1/SQeBSVe0B9AWmiEgrH+s1pnJa9wYthBVPwaxR1gCYsOdr8h8NzPJMzwLGlCygqmdU9bRntm411GlM5R3b5ZlQKDgNu5f7NRxj/M3XRNxCVQ8CeN6beyskIkkisgHYB/xBVQ/4WK8xlZM8EFxRzrQr0pk3JoyVm/xFZJGIbPTyGl3RSlR1n6qmAh2BcSLSopS6JolIhohkZGdnV/xbGFOepHS4bR5ExUBSH2femDBWbvJX1aGq2s3Lax5wSEQSADzvh8tZ1wFgE+B1t0tVZ6pqmqqmNWvWrPLfxpiyJF8GfSbC3lVw8mt/R2OMX/na7TMfGOeZHgdc8OBUEUkUkXqe6TjgMmCbj/UaUzW9xjknfte95u9IjPErX5P/DGCYiGwHhnnmEZE0EXnBU6YzsFpEPgc+AZ5Q1S98rNeYqmna0env/2wWuN3+jsYYv4n05Y9VNQcY4mV5BjDRM/0RkOpLPcZUq97j4Z07YOcS6HjBz9eYsGCXXZrw0/kHUC/e2fs3JkxZ8jfhJ7Ku85D3rf+B3DKvUTAmZFnyN+Gp93hwF8D62f6OxBi/sORvwlPTFGg7ADLtxK8JT5b8TfjqPd4Z9mG3DfVswo8lfxO+Ov/AGeo58xV/R2JMrbPkb8JXVDRc8kPY8h7k2nAiJrxY8jfhrfc4cJ+Fz1/3dyTG1CpL/ia8NesEbS51un7sGb8mjFjyN6b3eDi608b4N2HFkr8xXUZBdGM78WvCiiV/Y6LqwSW3wOb5sPi39ohHExYs+RsD0LqXc+J3+Z/tGb8mLFjyNwbg+D7PhD3j14QHS/7GgD3j14QdS/7GgOcZv3PtGb8mbFjyN+ac5AHQ907YuxKOZ/k7GmNqlCV/Y4rrPcG52csu+zQhzpK/McXFtYWLRjhDPRec8Xc0xtQYS/7GlNRnInx7GLYu8HckxtQYS/7GlNRhCMQlw9oX/R2JMTXGp+QvIvEi8pGIbPe8x5VRNlZE9ovIM77UaUyNc7kg7Q7Y8ykc2uTvaIypEb7u+U8BFqtqCrDYM1+a3wGf+FifMbWj560QUdf2/k3I8jX5jwZmeaZnAWO8FRKR3kALYKGP9RlTO2Liodt1sOFNyD/h72iMqXa+Jv8WqnoQwPPevGQBEXEBfwIeKG9lIjJJRDJEJCM7256sZPysz0Q4k+s0AMaEmHKTv4gsEpGNXl6jK1jHT4H/quq+8gqq6kxVTVPVtGbNmlVw9cbUkNa9IKGH0/VjD3oxISayvAKqOrS0z0TkkIgkqOpBEUkADnsp1h8YKCI/BRoAdUQkV1XLOj9gjP+JOHv/8yfDnv9B8mX+jsiYauNrt898YJxnehwwr2QBVR2rqm1UNRm4H3jVEr8JGt2ug+hGsPYFf0diTLXyNfnPAIaJyHZgmGceEUkTEfvfYoJfnRjocStsngeLpts4/yZkiAZoX2ZaWppmZGT4OwxjYOO78PYEQCAyGsbNt1E/TcASkUxVTSuvnN3ha0x5ju3yTNiDXkzosORvTHmSB4KrjjMtQL0mfg3HmOpgyd+Y8iSlw5V/cKbVDe//0vr+TdCz5G9MReQfw9ntBwrPWNePCXqW/I2piOSBzlg/4LQB9oxfE+Qs+RtTEUnpMH4BtBvkdP1EN/Z3RMb4xJK/MRWVlA7Xv+xc7rnyr/6OxhifWPI3pjLqN4UeP4TP34CTh/wdjTFVZsnfmMrqPxkKz8Kamf6OxJgqs+RvTGU16QAXX+WM93PmW39HY0yVWPI3piouuxfyv4F1r/k7EmOqxJK/MVWRlA5J/WDl36CwwN/RGFNplvyNqapL74Fv9sCcSXbHrwk6lvyNqar6TQGBje/AKz+wBsAEFUv+xlTVnk+/my600T5NcLHkb0xVFR/yAaBNf//FYkwllfsM30By9uxZsrKyyM/P93coJgRER0eTmJhIVFRU1VZwbsiHtS/AhjfhyJfQ9tLqDdKYGhJUyT8rK4uGDRuSnJyMiPg7HBPEVJWcnByysrJo165d1VeUlA6JfeDoLlj6B0i9CaLqVV+gxtSQoOr2yc/Pp0mTJpb4jc9EhCZNmlTPUaQIDJ0KJw/Yg95N0Aiq5A9Y4jfVplp/S8kDoMMQWP5nyD9Rfes1pob4lPxFJF5EPhKR7Z73uFLKFYrIes9rvi91BqLdu3fTrVu3Kv3tK6+8wuTJk6s5orLNmjWLlJQUUlJSmDVrVq3WHdKGPAynjsLKZ/wdiTHl8nXPfwqwWFVTgMWeeW9OqWoPz2uUj3WaSigoOP/u06NHjzJ9+nRWr17NmjVrmD59OseOHfNTdCGmVQ/oMsa56zc329/RGFMmX5P/aODcruMsYIyP6wt4f/7zn+nWrRvdunXjqaeeKlpeUFDAuHHjSE1N5frrrycvLw+AKVOm0KVLF1JTU7n//vvLXPeCBQvo27cvPXv2ZOjQoRw6dAi3201KSgrZ2U4ycbvddOzYkSNHjpCdnc11111Hnz596NOnD59+6lx3Pm3aNCZNmsTw4cO57bbbzqvjww8/ZNiwYcTHxxMXF8ewYcP44IMPqnMThbcrHoSzp2D5n/wdiTFl8jX5t1DVgwCe9+allIsWkQwRWSUitdpAZO45xt+W7CBzj+97t5mZmbz88susXr2aVatW8fzzz7Nu3ToAtm3bxqRJk9iwYQOxsbH8/e9/5+jRo8yZM4dNmzaxYcMGHnzwwTLXP2DAAFatWsW6deu4+eab+eMf/4jL5eLWW29l9uzZACxatIhLLrmEpk2bcu+99/Lzn/+ctWvX8s477zBx4sTzYp03bx6vv/76eXXs37+fpKSkovnExET279/v87YxHk1TnPH+M16Eb/b6OxpjSlVu8heRRSKy0ctrdCXqaaOqacAPgadEpEMpdU3yNBIZ5/Z0fZG55xhjX1jFnxZuY+wLq3xuAFasWME111xD/fr1adCgAddeey3Llzt3dSYlJXHZZZcBcOutt7JixQpiY2OJjo5m4sSJvPvuu8TExJS5/qysLEaMGEH37t15/PHH2bRpEwC33347r776KgAvvfQSEyZMAJyGYPLkyfTo0YNRo0Zx4sQJTp48CcCoUaOoV+/CSw5V9YJldhK9mg2eAohz6acxAarc5K+qQ1W1m5fXPOCQiCQAeN4Pl7KOA573ncBSoGcp5WaqapqqpjVr1qyKX+k7q3bmcKbAjVvhbIGbVTtzfFqft8R5TskEKiJERkayZs0arrvuOubOncvIkSPLXP8999zD5MmT+eKLL3juueeKLkNMSkqiRYsWfPzxx6xevZorr7wScLqAVq5cyfr161m/fj379++nYcOGANSvX99rHYmJiezbt69oPisri1atWpX/5U3FNUqEPhPh89che5u/ozHGK1+7feYD4zzT44B5JQuISJyI1PVMNwUuAzb7WG+F9GvfhDqRLiIEoiJd9GvfxKf1DRo0iLlz55KXl8e3337LnDlzGDhwIAB79+5l5cqVAPzrX/9iwIAB5Obmcvz4cb7//e/z1FNPsX79+jLXf/z4cVq3bg1wwVU4EydO5NZbb+XGG28kIiICgOHDh/PMM99dWVLe+gFGjBjBwoULOXbsGMeOHWPhwoWMGDGi4hvBVMzA+yAqBj5+xN+RGOOVr8l/BjBMRLYDwzzziEiaiJy726UzkCEinwNLgBmqWivJv3fbOGZP7Md9wzsxe2I/erf1eiVqhfXq1Yvx48eTnp5O3759mThxIj17OgcxnTt3ZtasWaSmpnL06FHuuusuTp48ydVXX01qaiqXX345Tz75ZJnrnzZtGjfccAMDBw6kadOm5302atQocnNzi7p8AJ5++mkyMjJITU2lS5cuPPvss+V+h/j4eB566KGik8QPP/ww8fHxVdgapkz1mzqPe9wyH/Z/5u9ojLmAlNWV4U9paWmakZFx3rItW7bQuXNnP0XkXxkZGfz85z8vOsdgqkeN/qbyT8BfLoGES+C2uTVThzEliEim5xxrmYLuDt9wNGPGDK677jp+//vf+zsUUxnRsTDwF7BzCcybbOP9m4BiyT8ITJkyhT179jBgwAB/h2IqK6EHILDun/bAFxNQLPkbU5OyVgOeK8HsgS8mgFjyN6YmJQ+EyLo4DYBC3Vh/R2QMYMnfmJqVlA7j5sPlv3Su/1/+Z8g76u+ojLHkb0yNS0qH7/0KbnoNvs2GBf8HAXqVnQkflvyrQbAN6Txy5EgaN27M1VdfXav1hr1WPZ1hn7csgMxX/B2NCXOW/ENcySGdAR544AH++c9/+iEaQ//J0P578MGvbOgH41eW/Csp2Id0BhgyZEjRGECmlrlccM2zUCcG3r4DzlbDYySNqYLQT/771jhjq1fD9dWhMKSzCQANW8Lov8OhL2DxdH9HY8JUaCf/fWtg1ij4+FHn3ccGIBSGdDYBotNISL8TVv0dtn/k72hMGArt5L97ORSeAS103n28wSYUhnQ2AWTYb6F5F5h7F+R6HQ3dmBoT2sk/eSD4z2FaAAATWklEQVRE1AGJcN6TB/q0ulAY0tkEkKhouP4lOH3SaQDcbn9HZMJIaCf/czfYXPEb5z0p3afVhcKQzgADBw7khhtuYPHixSQmJvLhhx9WckuYatO8M4x4FHYsgtX/8Hc0JozYkM5BwoZ0rhkB8ZtShTd+6PT9/3ixMwS0MVVkQzqHEBvSOcSJwKhnnAfAvDEWlv7eRv80Nc6SfxCwIZ3DQP0mMOBncHwfLJ1hwz+bGmfJ35hAceZbvhv+OR82vOnXcExos+RvTKBIHgiR0Tj/LQUyZzljAAXoeTkT3CL9HYAxxuPc1Wm7l0PLVOcGsAX3wr61cNUTEGU37ZnqY8nfmECSlP7dJckdrnBO/i57HL7eADe+CvHt/BufCRk+dfuISLyIfCQi2z3vcaWUayMiC0Vki4hsFpFkX+oNFcnJyRw5cqTCZRo0aFCp9U+bNo0nnnjiguXjx4/n7bffrtS6KqK04an9MWz1rFmzSElJISUl5YIb5oKGKwKueBBueRO+2QMzL4cvF/o7KhMifO3znwIsVtUUYLFn3ptXgcdVtTOQDti97KbalBy2+ujRo0yfPp3Vq1ezZs0apk+fzrFjx/wUXTXoNBImfQKN28DrN8CSx8Bd6O+oTJDzNfmPBs7tVs0CxpQsICJdgEhV/QhAVXNVNc/Hev1i9+7dXHzxxUycOJFu3boxduxYFi1axGWXXUZKSgpr1jiX5h09epQxY8aQmppKv3792LBhAwA5OTkMHz6cnj17cuedd543VtBrr71Geno6PXr04M4776SwsOz/3I8//jh9+vQhNTWVqVOnFi1/9NFH6dSpE0OHDmXbttLHi1+2bBmXXnop7du3P+8ooLT1jhkzht69e9O1a1dmzpxZtPzll1/moosu4vLLLy8aUrostTFs9YcffsiwYcOIj48nLi6OYcOG8cEHH5QbW0CLbwd3fAQ9xsInf4DZN9jjII1PfE3+LVT1IIDnvbmXMhcB34jIuyKyTkQeF5EIbysTkUkikiEiGecSQaDZsWMH9957Lxs2bGDr1q28/vrrrFixgieeeILHHnsMgKlTp9KzZ082bNjAY489VpScpk+fzoABA1i3bh2jRo1i7969gHOX6Ztvvsmnn37K+vXriYiIKBrC2ZuFCxeyfft21qxZw/r168nMzGTZsmVkZmbyxhtvsG7dOt59913Wrl1b6joOHjzIihUreO+995gyZUqZ6wVnNNHMzEwyMjJ4+umnycnJ4eDBg0ydOpVPP/2Ujz76iM2bN5e7/Wpj2Or9+/eTlJRUNJ+YmMj+/fvLjS3gRdWD0X+Dq59yTgo/dzns/8zfUZkgVe4JXxFZBLT08tFvKlHHQKAnsBd4ExgPvFiyoKrOBGaCM7xDWSudvmATmw+cqGAIFdOlVSxTf9C1zDLt2rWje/fuAHTt2pUhQ4YgInTv3p3du3cDztDP77zzDgBXXHEFOTk5HD9+nGXLlvHuu+8CcNVVVxEX55wiWbx4MZmZmfTp0weAU6dO0by5t3bUsXDhQhYuXFg0rlBubi7bt2/n5MmTXHPNNUVDR48aNarUdYwZMwaXy0WXLl04dOhQmesdNGgQTz/9NHPmzAFg3759bN++na+//prBgwfTrFkzAG666Sa+/PLLMrdfVlYWN910EwcPHuTMmTO0a+ecwLz99tsZPXo0P/vZzy4Ytrp4o1KRYau9DVlSctTVoCUCaRMgIRX+PQ5eGgH97oa6DaDdIJ/HrzLho9zkr6pDS/tMRA6JSIKqHhSRBLz35WcB61R1p+dv5gL98JL8g0HdunWLpl0uV9G8y+Uq6nsuK/l4S0Kqyrhx4yo8fIOq8qtf/Yo777zzvOVPPfVUhZNc8e9xLt7S1rt06VIWLVrEypUriYmJYfDgwUXDTVc2qd5zzz3cd999jBo1iqVLlzJt2jTgwmGrzx0FnBu22luSL23Y6sTERJYuXVo0n5WVxeDBgysVZ8Br3ds5D/D6DfDpk4BARF0Yv8AaAFMhvl7qOR8YB8zwvM/zUmYtECcizVQ1G7gCyPBSrlLK20P3p0GDBjF79mweeughli5dStOmTYmNjS1a/uCDD/L+++8XnYQcMmQIo0eP5uc//znNmzfn6NGjnDx5krZt23pd/4gRI3jooYcYO3YsDRo0YP/+/URFRTFo0CDGjx/PlClTKCgoYMGCBRck8rKUtt7jx48TFxdHTEwMW7duZdWqVQD07duXe++9l5ycHGJjY3nrrbe45JKyByWryLDVP/rRjy4YtvqBBx4AnGGre/ToUe73+PWvf120fRcuXBia4yLVbwKdrvR0/ahzV/Dcu5wrhDp9HyLrlrsKE758Tf4zgH+LyB04XTo3AIhIGvATVZ2oqoUicj+wWJzdxEzgeR/rDWjTpk1jwoQJpKamEhMTU5Tkpk6dyi233EKvXr24/PLLadOmDQBdunThkUceYfjw4bjdbqKiovjb3/5WavIfPnw4W7ZsoX///oBzCehrr71Gr169uOmmm+jRowdt27YtetZARZW23pEjR/Lss8+SmppKp06d6NevHwAJCQlMmzaN/v37k5CQQK9evco9UX1u2OrWrVvTr18/du3aVfTZqFGjmDBhwgXDVt99992kpqZSUFDAoEGDyh26Oj4+noceeqioG+3hhx8mPj6+UtsiaLS7HCL/BAWnQVyQfxzeGg/14qD7jZDQA3IPOncP2xGBKcaGdDYBwx/DVofEb2rfGucEcPJApzto5xJYNxu2LoDCs04ZiYBLJ0PXa6FFV4iIuvBvrXEICRUd0tnu8DUBYcaMGfzjH/8o8yonU4ridwUDdBzqvD5+BJY9AajzKNNP/+K8IupCy+7QqDVs/a9zz0BEFNw2F9peWnaDYI1FyLA9fxPWQvo3tW8NzBrlPL/aFQXXzgQtcM4RHFgHWWu+OzI4JyYe8o4B6hwtdP4BNOsE0Y3h1FFY8RdwFziNxehnIKkv1KkPUTHOpajnLgAorZGwxqPG2Z6/MeGu+EBxxZNtt+uc9z0r4dXRTgPgioBLboKvN35385gWwrb3YfPcC9ddeBre/XGJheI0BBFRcOobQJ1lrXpCwwQoyIedS0HdTn19fuwcgUTHQt2GUDfWeUV73iPrOo2JNRg1wpK/MaGsZJdQcW37w/j3zk+sxY8WIuo4jUfr3s6J5F3L4N1JnsYiEi5/wEnqZ/OcZxGczYMzebB3JZxa56lE4dvDzgnpkwecBgWco4fynlnsioI6MZB/wlmPuJzB7pp2cq50iin2OrEfsr90urva9K2urRfSLPkbE85KNg6lHS3ExEPXMRDbqvy98JINyPUvX9iwuKLg5tnQpAOcPgmnTzhJvmj6uDO9e/l3dzGrG7LWOkcsZ7/1XveyPzoxxrVzGqaGLZ35hi2d+dxsOLIN2g8O+6MIS/7GmPOVdbRQ1mfFy3hrQEpbXpaSDcnYt52/O3vK6Z7Ky4E1zzlXN53rZqoX5zwAZ38mnDzodDeVtOQx5y7plt0hLtlpLM586zxGM2V4WDQMlvz9KDk5mYyMDJo2bVqhMg0aNCA3N7fC6582bRoNGjTg/vvvP2/5+PHjufrqq7n++uurHLs3r7zyChkZGTzzzDMVWl6TRo4cyapVqxgwYADvvfderdVrPEprJCrSeJQs763BiKrnXK3UqDX0GgdfvPNdA3H1U9+VU4X8b+DEQVj5V1j/L5xGAqfx2L4Icr8+v85lj0PCJdA6zTnZrZ6uq4tGhlSjYMnfBL2CggIiI8//KT/wwAPk5eXx3HPP+SkqU23KazDKOqIQz5FAvTjoPQE2zinWHfWSU/ZMHiz+Lax+DnAD4hxRbHzb6X46Z/mfnUthO3zPOYmd0NM59xCk7Bm+lWBDOgfHkM7gDJnRsGHDcuMxISIpHQb+omKNxBW/cd7Pla0TA92u9VxdFOE8R/n6l+GXe2DAfZyXJrO3OvdPvHYdPN4enuwOb94Ky/8Eq/4BS37vdFUFA1UNyFfv3r21pM2bN1+wrDbt2rVLIyIidMOGDVpYWKi9evXSCRMmqNvt1rlz5+ro0aNVVXXy5Mk6bdo0VVVdvHixXnLJJaqqes899+j06dNVVfW9995TQLOzs3Xz5s169dVX65kzZ1RV9a677tJZs2apqmrbtm01OztbVVXr16+vqqoffvih/vjHP1a3262FhYV61VVX6SeffKIZGRnarVs3/fbbb/X48ePaoUMHffzxxy/4HuPGjdPrr79eCwsLddOmTdqhQ4cy16uqmpOTo6qqeXl52rVrVz1y5IgeOHBAk5KS9PDhw3r69Gm99NJL9e67776gvpdffrlo+dGjR9Xtdquq6vPPP6/33XefqqpOmzZNn3zyyaI4rr32WlVVveWWW3T58uWqqrpnzx69+OKLVVV16tSp2qtXL83Lyyv132vJkiV61VVXlfq5qv9/UyaA7F2tuuwJ5734st+1UJ0W57zvXa166hvVnctUV/xF9a0Jqn/poTo1ttirseqcu1R3fqJ65lTp664hQIZWIMcGb7fP+1Pg6y+qd50tu8OVM8osYkM6B/6QzsZUibfupdK6lNoNdF7nfPyo525qt/Na/zqsn+3cTd28MxzaCG63090UICOvBm/y9xMb0jnwh3Q2plpV5CR1yjD431+/O59wy7+cext2L4cv3nbuawBn5NUPfuV0PSUPdG6I89NNbMGb/MvZQ/cnG9K5dLUxpLMxta60I4ROI6HLaHjlB85d0SLOUcA/r4G6jSAxDXYtdxqHyLrnn4uoYcGb/AOYDelc9rap6SGdAQYOHMjWrVvJzc0lMTGRF198kREjRlRqexhTKWVd3jp+wXcNQ8vu8NUS2PoebHoX3GeccgX5sO6fzs1sez6t8SMBG9jNBAwb0tmEnT0rnZvY3Gcpuv8AT1dqZHSVjgQqOrCbXeppAsKMGTO47rrrQvOJW8aUpm1/mPAfGPIQ3DYful6D0wioc/5gd83tCFm3jwkIU6ZMYcqUKf4Ow5jaV7y7KKoebPvAOVkcUcfp+qkhlvyNMSZQVGX8oyoKuuSvqpW+vNAYbwL1fJcJc5Ud/6iKgqrPPzo6mpycHPtPa3ymquTk5BAdHe3vUIzxC5/2/EUkHngTSAZ2Azeq6rESZb4HPFls0cXAzarq5fFAZUtMTCQrK6toDBhjfBEdHU1iYqK/wzDGL3zt9pkCLFbVGSIyxTP/y+IFVHUJ0AOKGosdwMKqVBYVFVU0HIAxxpiq87XbZzRw7jbNWcCYcspfD7yvqnk+1muMMcYHvib/Fqp6EMDzXvpoZI6bgX/5WKcxxhgfldvtIyKLgJZePvpNZSoSkQSgO/BhGWUmAZOAoqEPjDHGVD+fhncQkW3AYFU96EnuS1W1Uyll7wW6quqkCq47G9hTYnEj4LiX4t6Wl1xW3nxT4EhFYquC0uKurr8rq1xltpm35aG63cor48tvzduy4vM1uc1Ki6e6/qYmt1uo/tbKK1fd/0dTVLVRuRFVZND/0l7A48AUz/QU4I9llF0FfM/H+mZWdHnJZRWYr9ADEKoz7ur6u7LKVWabhdN2K6+ML7+18rZbTW6zqm636vit+brdQvW3Vl65mv4/WtrL1z7/GcAwEdkODPPMIyJpIvLCuUIikgwkAZ/4WN+CSiwvuay8+ZpU1boq+ndllavMNvO2PFS3W3llfPmteVsW6NutOn5rZX1u/0cr/1mN/h8N2FE9a5uIZGgFRsIz57PtVnm2zarGtlv1Cqo7fGvYzPKLGC9su1WebbOqse1WjWzP3xhjwpDt+RtjTBiy5G+MMWHIkr8xxoQhS/4VJCL1RSRTRK72dyzBQEQ6i8izIvK2iNzl73iChYiMEZHnRWSeiAz3dzzBQkTai8iLIvK2v2MJFiGf/EXkJRE5LCIbSywfKSLbRGSHZ0TS8vwS+HfNRBlYqmObqeoWVf0JcCMQFpfnVdN2m6uqPwbGAzfVYLgBo5q2205VvaNmIw0tIX+1j4gMAnKBV1W1m2dZBPAlzo1pWcBa4BYgAij5BPHbgVScW8ujgSOq+l7tRO8f1bHNVPWwiIzCufP7GVV9vbbi95fq2m6ev/sTMFtVP6ul8P2mmrfb26p6fW3FHsyC7jGOlaWqyzx3GBeXDuxQ1Z0AIvIGMFpVfw9c0K3jeSBNfaALcEpE/quq7hoN3I+qY5t51jMfmC8i/wFCPvlX029NcO6Ufz8cEj9U3+/NVE7IJ/9StAb2FZvPAvqWVlhVfwMgIuNx9vxDNvGXoVLbTEQGA9cCdYH/1mhkga1S2w24BxgKNBKRjqr6bE0GF8Aq+3trAjwK9BSRX3kaCVOGcE3+3p4AX27/l6q+Uv2hBI1KbTNVXQosralggkhlt9vTwNM1F07QqOx2ywF+UnPhhJ6QP+FbiiycgebOSQQO+CmWYGHbrGpsu1WNbbcaFq7Jfy2QIiLtRKQOzhPG5vs5pkBn26xqbLtVjW23GhbyyV9E/gWsBDqJSJaI3KGqBcBknKeKbQH+raqb/BlnILFtVjW23arGtpt/hPylnsYYYy4U8nv+xhhjLmTJ3xhjwpAlf2OMCUOW/I0xJgxZ8jfGmDBkyd8YY8KQJX9jjAlDlvyNMSYMWfI3xpgw9P9EdKyoawHuqQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.semilogx(tobs, hobs0, '.C0', label=\"obs layer 0\")\n", + "plt.semilogx(tobs, hobs1, '.C1', label=\"obs layer 1\")\n", + "\n", + "hm = ml.head(robs, 0, tobs)\n", + "plt.semilogx(tobs, hm[0], 'C0', label=\"modelled head layer 0\")\n", + "plt.semilogx(tobs, hm[1], 'C1', label=\"modelled head layer 1\")\n", + "\n", + "plt.legend(loc=\"best\")" + ] + }, { "cell_type": "code", "execution_count": null,