Stochastic simulation and probability-entropy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import os\n",
    "os.environ[\"THEANO_FLAGS\"] = \"mode=FAST_RUN,device=cuda\"\n",
    "import theano"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "if '/' in sys.path[0]:\n",
    "    sys.path.append(r'/home/elisa/PycharmProjects/gempy')\n",
    "elif '\\\\' in sys.path[0]:\n",
    "    sys.path.append(r'C:/users/elisa/documents/git/gempy')\n",
    "\n",
    "import gempy as gp\n",
    "#import vtk\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for plotting\n",
    "SMALL_SIZE = 8\n",
    "MEDIUM_SIZE = 10\n",
    "BIGGER_SIZE = 12\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "plt.rc('axes', linewidth=1)\n",
    "import matplotlib\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "matplotlib.rcParams['ps.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Regular model setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Active grids: ['regular']\n",
      "Active grids: ['regular' 'sections']\n",
      "Active grids: ['regular' 'topography' 'sections']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Grid Object. Values: \n",
       "array([[ 166.66666667,  166.66666667,  166.66666667],\n",
       "       [ 166.66666667,  166.66666667,  500.        ],\n",
       "       [ 166.66666667,  166.66666667,  833.33333333],\n",
       "       ...,\n",
       "       [1000.        ,  200.        ,  974.6835443 ],\n",
       "       [1000.        ,  200.        ,  987.34177215],\n",
       "       [1000.        ,  200.        , 1000.        ]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAEPCAYAAADyGVe/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXmwbVdd5z/v3jcmUCFBhkAIESM/KKBsG8oGtVMSK0oXplVoUShsg4AhEImQqEEkQGTSJgzpZjQBtMTSBqFbCAG6m1abQRxQG9BaEEOHDECEJJC8IW/sP/ZZeev93u+39trn7nPv2e/9PlWv9jl7WHudfc9b57t+09p06NAhgiAIgnFY2egOBEEQHEvEoBoEQTAiMagGQRCMSAyqQRAEIxKDahAEwYjEoBoEQTAiMagGoyIi9xaR+2x0P4Jgo4hB9ThDRM4TkVtFZNvs/Rkisl9EHthz3Ued/e9R174YeMR4PW5HRH5ERC7diHsHQWbzRncg2BD+Gvgp4I+BZwF/BiAiPw+cDWwF3gMk4HXAN4DvmZ1zGXBf4BTglWWjIrIC/CjwABG5CbgC+Bdgd0rpV0Xky8BVwOnAq4BHA+cBe4H/BXwQeAdwG/ADwK8AvzG7//uBnwN2zvryzFn7X5/tux64GThXRL571r/zgatSSk8RkWcC+1JKf7zmpxcEFUKpHp+8H/hpEVkFHgDcNNv/8ymlZ9ENdBcDFwCvSym9CLhFRB4JPAnYDdwJ/EjZaErpIPAJ4PeBX55deyGwTUQeBdyQUvptuoH1l+gGzV+Y3fNZwNOB96WUXgj8w6zZewEXAZ8Efg/IivnRs+1/SSn9FvC02fv/k1I6H/gOcCJws4gI8FTgT+Z+YkHQSAyqxye76dTf+cA1xvFDwKbZNrOf7vtyXUrpUuDtwBeca5mde3D2+uCsvfx920GnTlfUPbYV7w/MtntTSrvoBtELgNuBL83ay/cpt3fMtvuAVeBNwH8G/iqltN/obxCMSkz/j1+uBj4AZBUH8Aci8ruz168Dvgy8RkS+CTwwpfRFEbldRN4OnAo832j3OuDXgBfNrv0asDOl9AUROVNE3gDcB/h14HPAu0TkTjr1ei3wFhE5C/ghuoE7853ZdU8BHgX8xWz/ZSKyb/Z5jiKl9M8ishl455CHEwTzsikKqgTrhYh8NKX0pMrxhwGXAHcB24EXpZQOVM5/D3BpSunrlXPeClyfUnr93B0PggHEoBoEQTAiYVMNgiAYkYXZVEXkTDpP7veLyK8CZwAn0dnattKFw9wGfDGl9BZ9TkrpXxbVtyAIjh2WbaxZiFKdBYM/B9gpItuBs1JKL6BzJjyXzut8ZUrp+cCTReTexjlBEARVlnGsWYhSnTkOLp1l4ZwC3Do7dBPwILpfjxtn+24HTjbO6UVEwiAcBOtASmlT/1lH8rCHPezQli1b5r3lDSmlM/pOWq+xZgjrEVJ1K10GDsBpwC10Cvk0ug97ymyfPqeJyy+/nO3btwOwdevWe/Zv3rz5iO2mTUd+J/T7FkqnXr4+b1dWVo643+rq6hHb8rV3jnWNvo8+J9/X2qevraE/h35fvtbneH2Fw88sbw8c6Jz5Bw8ePGK/fm21a52XX+f29u/vQlHvvvtuAHbv3n3Pubt27Tpim4/l9zt37jxiC3DXXXcdsdXX7t27F4B9+/bdc82ePXuO6IPuSz5e9je3k99n8nel/G7ngUof27ZtG8A9/x/KfXmbr83vTzjhhCO2ADt27ADgxBNPBOA1r3kN87Blyxa+9KUvzXXtwx/+8IfOcdlCx5pWFj6oppT2i8j/noW2nAw8jy74+w0ich7wAeecJrZs2XLUYJX3gz34zIs1qOZ28330fa1BVQ9OtT56A6IerODwwJLR7bUMermNfK01gPVFjFjHc7t6a507JCJFD9Z5UNJbODxw6UFOD7bloKoH2rzNA2NuM2/L9vM+3W450OfBWP/Y6B/pfLx8nQfG8jum8X789XeubEMP2mthjP93rSx6rGlloYNqjklMKb1ZHfo28Ax1rj4nCIKJM8+McB6WaayZfEbVli1bzF9WrRyt6Sy0TT9raJU8RKl6fSrp+1JqdVrrYw2tQmuqNO/TKsjqq1akWTnW+q3v2fJ38RRqqSC1qsxqMytHva0d01P88j5aAedz9fvyOu+Z5ql9+Wy1KUir3BI9K8rb/D3N7Ze2z7yvNDnMy3oq1WVh8oNqEATLy3op1WUiBtUgCBZGDKoTZPPmzfdMU8rpv/aw6yiAvJ3XUeJ543XUQTn17vOet9xvniiGlun0kCn3kOiCPCXNU9TaVFWbHPTfqNbv3G52/OgtHG0S8Lz/VsRAnrL3efbLa/RWO7fKPunPrJ1E5TP2nKPaCVi+ztdrU5k2A5T71hAOZfb7eGHyg2oQBMtL2FQniBdSpX+JtXJscarUfmU9B4AXi2qd26LEMp6ythwYHrX2vT7UQqq8/zBWmFfe6hAiq09eXKp1rlbAuf2sIEulml97Dqt8TelI0gpVX6vDpsrXemvFtOrnk5+p/i6Wzzp/p71naKlaHaKlnVGWUm1xcPYRg2oQBMGIxPR/gqyurpqq0FOqNXtgX3hSTQF4irUWUpWphTH12XetL62n9PJ+63N6wfiW0uhT1uV+T0nWbKv6c9TUs2dL1YH+4IdB6XPLa3S7OmNLZ1iV+/L9dJhX+Zm1TdWyi8KRf4fcX21LtUL0PFu/9721rlkLMagGQRCMSAyqE2R1dbVqU/VqAMzjgS/xAvhrStW7p1YpFp5HfMi52lNuHavhpbBqBVx+Dq1Ua6qtNf3Vat9TkqXqzN59rSQ9W2v52ksqsPL5tT1Xf9YyddZT7HqGYCnJnOPfYlPVCrVWb2LM1O6wqQZBEIxIDKoTRceMgm9LrcWK9n0BagVVdHu1ilNWe3202Bf7ogl0zGjt3FofPWVdU5Ja6WkFq19bfasVkdGq0PLKa/XqefJL77zuv7a/6miDcp/uk6XOPaWqZwXls6kVowE7IqTv+1oy5hJLMf0PgiAYkRhUJ8imTZuqtT/7Sp7VrmnJ5PEUqqUE+sr4ZVqKpNTa8bz8NZtnX5tlu14cr6WEtbc/q7W8v6baPCVsRUdoNWjZVD1FqhWlFdvqKdVsU7Xuoz+zZVP1SiLm70++fznj0c+5JXpEzyL0s7aU8BjEoBoEQTAiMagGQRCMyBhZWVPjmBhUrUD+1ornQ6b/1rRqyJIlnpOpxczgTbmHpJ5ajqQhjiqvfqoO1SqnktoxpcOXLEeVt+SK9Qz09FU7oSwHkjYVaBOF1Sc9/dehVtY1nunDKnjiTeUtR5XXf6sv3nS/9jerpRAPJZRqEATBiMSgOkEOHTpUXUTOU1eWUvVUZ+0X2yuE4fW11n4tZKilLy3t6fdWQkDZx5Y0WO1kKR0xXkV+qzK/VlpeCqf1WVpK/3lKVas1SxXqECqtuGv3qbXfp8atv2G+XqtlrYyte3qOsVCq4zH5QTUIguUlBtUJcvDgwSYlM4SWoHPPTqkVhtWnIXZY7361c/qurYUxtTwvT2lbSkyrJ0/pwdHl9XSfan30kgysItWecqwF5+sQrVphGG/FWOv59xVQ0X2z+qnfW8/fmynUZnl94XYtxKAaBEEwIov0/ovI44FLgRuBvwTOmR16LPAm4Drg1bPt51NKVyysMwWTH1QPHDhg/tr2qQTrF1SnqXpKyWq3ZWVRrxBJyxIp3ucYolxrSrUv9bSl31qNlq89u5+1Cqm2Fdb63xfUbtlsPZVWW46kz3Zb9rXPpj0kUqRmh/USKmrRC57duhaRsBYWrFSfDrw8pfQPIvJHwC8C9wdemlK6WkReBnwNWAE+vciOlEx+UA2CYHlZ8KD6RuAyEbkNuBdwH+CVwGWz49cAVwJ3Ax8Bzl5kZzKTH1QPHjxo2rX67KJDPO2WgvFUiFa7Q4pIe/2o9WVIpEDNk+wp4JpNWH9WrYbK155CtcrsZdtqX7yqtc9TcVZfrFTZodRs3d7sZcgSOPo+1j1bVKfn7a89Jz1TmIcFD6qnA69KKV0vItcA+4GDKaUbZsefAHw5pbRHRPa5rYzM5AfVIAiWl3lL/zU6mG8CXi8idwB/CDwO+Nvi+PXAO0RkF/C2uToyBzGoBkGwMOYdVFvsuSml64GnVI5fC1w7VwfWwOQH1f3791eDznM9Ve2gsuqTesdq0/++Kbf1peoLqRricBuS8mh9Zu+aGl5t11oguT5mpZF69U5bnHGe2aI2re1Li7U+s36fv1+WI0nX1M33bXFE1taQ6ktoqeF9f2p/s7WwyEF1WZn8oBoEwfIScaoT5MCBA+YvqzbAe+qwpiRrDhJP1XipqFb7LemwnmPNUqpagelzautbDVEGXuGZlkB4HfRfru2kVzvVjqSWdOGWRAft0BmS+OCtTmrNeLznbzmAvDTnlrXOaisEZzxFX/vOjaEWYzmVIAiCEQmlOkEOHjxo2s+yyvHsTrWQIU0tlMc7t+XLpH/Fa+rBC263PnNfmqTVbostWPezpfSfp1Atm6ouTmKVB/T6ot9bn9lT8n3B+nC0Qt22bdsR51phUrogjLUahLa3eokC1krBuoxlfl8rgdmSIt2SEt1KDKpBEAQjEkWqJ8jBgweP8vTD4V9+z/Nq/bH7bJzWr64XVaCPl+31ef0tRex5a6015C3F2Hf/lgSEjPecrCLJnkK1SubpcoC6Pav/fTOGlkIn+nOVKj1/p7Zu3XpEnzzbZ3m9tq1a35H82bwUZt0POLxScO5T3ub9Vl/6CrBbs4ExlGrYVIMgCEYkpv8TpPT+1wpJeL/ULcyzZIpWHlZ7LQWo+7zataVL+pboKPviFUtpsb3V1LNXNLqloLJXrHqekonla0+hWp72rP68IiO1ONLcTlbe+X2t9KJHqVSzMs123bzVirV8na9vWTF4zIIqoVSDIAhGJJTqBNm/f391wTPvV7dWAjDTUghjHhXbZwesXePFG5b7PLvfWot/eH1rsfN6sayWzdmjxRas97dEbGj7eu3voPti2ei9WNaW76nXx1Lx7dixA4Dt27cf8b5FqeZtzX8wZkZVDKoLQkROAy4Hbgc20dU4PAM4CXgRsBW4ArgN+GJK6S3r0a8gCBbL8ej9Xy+DxyPoahmeDnwTOCul9ALgauC5wPnAlSml5wNPFpEtbktBEEyGlZWVuf5NmfWa/t9IN6h+Bfg4XckuZtsH0SnVG2f7bqdTsN9safjAgQP3OD1KY74OjPaKpNT+gF4aYNm+xjpXo6d6tTCsTF+Qe42Wab9XfKX2fDwTS0vd2draSPreVipoxnPG1Z6hnrK3mDy8z5Gn12USQ+6vFxpWC4Pz7ld+t/O0/8QTTzzife6LNf33CrTUvq8TqPy/lKzXT8KFwEkppUPAt4GHzvafBtwCfHX2GuAU4I516lcQBAtk06ZNc/2bMuulVN8NXC4iXwU+C+wVkbcCJwPPA3YAbxCR84APpJSaS47v37/f/LX1VKZWqOU18xSq0Az5QrQUVNHneqFP5Wtdji5vrXAyzwFjKT6vvzVHlXaeecrS6sOQNN6aA89DK9aWWYa+Jn8eKzg/hzq1VNn3lGrGUqr3ute9gMOOKh1iVV6nEwN0qms58xpzUJv6VH4e1mVQTSl9Dji3csq3gWesR1+CIFg/jkdH1eRDqqC+3nn5Cw9H25LKP7oOjq+tN1Wzt1rvS/I957GTaqVtpVRqZaqvrd1zSCqrp1TL+3irjtYKkfTdr5bGO0QBe3+jWuic/t7kz1zaMbNS1KFt1tpPrSnF5fdYK9K81bbV8jpPmVpFWMZk6lP5eTgmBtUgCJaTGFQnjpWyWXpl4WjlauF5xNcaNN+qxKx9nmfcUtr6vVbcLavOWuj+5+etVWhZJEUXVNHqrZaa60UXWAkPLYrb+xvVPOHeud4yMeVr7xyryEtfCm35N8x/c227tYL/8z7PtmrZVMecsodNNQiCYERCqU4cqyScVgm1X+EhXv6+a1s8yB4tSsxKj9Q203kKcteKaeglP7StUKtSOLxEii71Z9kXtb1Vt299jj6F2rIMif5OWEpVF0PxtlZfaunDfbMV6/uk43d1CcDS+6+P9UUDlO2PoTIXqVRF5PHApXQx7p8GzgFWgUN0S1LvBl4C7AI+nFL64MI6U3D8afMgCNaN1dXVuf418nTg5SmlXwZ+EngMsBO4C/gicAlwcUrp2cAFC/h4JpNXquUfwVpyoi/WdEiRlJqd1LO9WfS1bxW30MWMM+W5Wknq91bUgWc/rilWfc6Qxfy0QrU84fpYSym6ltKF2gPunVs+H61QtZJsKe6sn5c1Y/BmL5bNWStJbUvN7619LQVW1jJj0yx4+v9G4DIRuQ24N/DGlNIfishPAC8E7pdSumV2bn/g8kiEUg2CYGEsOKPqdOBVKaVLgBOA75rt/xZd6vvNInJq7sqoH6zC5JVqEATLy4KV6k3A60XkDuAq4N+IyJuB+wC/BtwXuEJE9gJXLrIjJZMfVLdu3XqU4R6OdkIMCRPR05+WtdeHmAH6plfl9NBbQaDWhj5Xt1H2SU9Fa+tc6ZAmvaaUtUJqnv57Ff+HhCJZSQ3eM6xV5G/925WftfV9rU+Z2vTfM0VYZh79vbQKqmjHlDYDWN/TMeqoZhY5qKaUrgeeUux6rzrlG2xApubkB9UgCJaXCKmaIGX4iOWU6EtJHBLQ33KOF1BuXev1pZYOW8NzqnhpseW5WhnVgtp1CFVWppZS1SrWW0ervKe3ram2lpmCV/auJQwr4yVwWPcZMqDooi4tSRhaqWrnU/laK1V9P0udjqFYI/g/CIJgREKpTpBt27aZv+p95dyGlIYbi9a+WCFPuo0aXnu1kCp9bq1IjbalekVTynO8Qs1l+17okVZMls3TCoeCNqWq2225RodlWdf0FTIvr/fCvGozHa+Mn1VkR9td9f2sIjWxmup8TH5QDYJgeQmlOkF27NhhFqpoDWBuKWJSO2cteCrLUmKe7atmf9WKyWq/LznCWqNeK1LPa1++7rOTWvu8gPuWFOBaeqenCvV55bnaXqlt9i2JJ9bSKH3K11KseZ+nUGufuZaSm9HRHWshBtUgCIIRien/BNm6dWu1UEVfoWNLKfWVk+s7pvuiqaVFeudqaiq6b3G3Fpth7T+DZ3e1+qTv2WKn6ysiM6QAdS2mtSVioC/OUx8vr9eK1Xq2nsqslSMcsoifbl9HF9RmYWMMiKFUgyAIRiQG1SAIghGJQXWClEZ/axVP732Lg0qnRZZT17xPB89nakHcfU60FjPDkFVDa2FGui/5eWZnVG1lAau9si2r/ZY0T4+aI0/fb54A/Jp5RE/z9XpQlqOqZT2ovml/y2qntfathIYS63s0RiiV19fjgckPqkEQLC8xqE4Qb80prwiHl5ZZvm6pf9kXyK9rmpbHhqSeemFGNUfYkDTPIUpxzJTcWkiVh1X5f57VbPuK09QceVoV6sr64DuoautB6c+h97esYDDk76BnY1bq7xiE9z8IgmBEQqlOkFIhlGj1oQOZa/Y5rxxeiV43yVMLlqqqFbEoj5ftt9hSvXM8W3H5uiXhwVNTejZQqir97MYIJ1srnqqthTx5wfnWNV5IlU4VLc/RFfj1NS1hcC1hZPo7oovweO3MSwyqQRAEIxKD6gRZXV017apacWnbTq0MXkZ7/VsCpWsKwwtmr/WjJUJA7/PWmbeSJLxC0BaeLVh/5rKNPi98ub8vLdVqw7OH1lRbn6e9Fh3RUpjHU5tWaT4dPeDZUlvssLXn481WrMI5Yw6EMagGQRCMSAyqE2R1ddW0hWpb0RDvsP71tjyYWr3q+2ibVdmOPrdW0MNTGpaq04UwdDFpq9yeVwyl1n5LBIJmiOLWx2p/O88TXiuo0lfGz1K3fc/fek6eJ79Uql75Pm9/rb+14igZb6Xd2nduLYT3PwiCYERCqU6QzZs3myXV+gqqtKgsT+XC0crO82pb3v++0nO1ItJe3GrZJ61U9bZUqn1qtiVSwCusYn2mIYVaNDUlOSQmtEXh6T7N85kzXjRA2Sd9jlaqLUVYagpTq+fad9uaZc1LDKpBEAQjEoNqEATBiCxyUBWRxwOXAjcCnwZ+iG5MewBwMfAQ4NXAdcDnU0pXLKwzBZMfVFdXV03HSV8d1ZYQIn2tFerUF4ZVW7uoJd1TT5dr005dmV+/t8Jn9HS/xVE1T1EU73nV6tnqa2v1Qr3K/GWgfd+6ZVaftPPPW03BMtnobS2Q3yu+UktTHVL71nP6WY4q657zsmBH1dOBl6eU/kFE3g/8fkrpT0XkqcA5wP2BrwErdIPuujD5QTUIguVlwdP/NwKXichtwHbgb0TkTOBngWcD3wtcCdwNfAQ4u69BEfmxvnNSSh+vHZ/8oFqGVNXoU67WsZoDqa+UYC0Qe0jKo1Y7tRRaT6nqbZmy661yWns+HlZywJDwHM9x1KLkPaU6JGh+HsfVPJ/H6r93Tkvl/5YCPUPC4MZQqJkFK9XTgVellK4XkWuAHwb+LXBeSmmXiDwB+HJKaY+I7Gts8/eAawHvS/DjwINqDUx+UA2CYHlZsFK9CXi9iNwB/BXwFuBjwDtF5H3A9cA7RGQX8LbGNq9KKb3MOygiv9XXwLoMqiJyBvAy4FbgTuA+wDbgJOB8QICXALuAD6eUPtjadmkza6FW8sz7Fa+pk3nWU+oriGEpPV302lIeXjhUvkaHTZXXa4Vq2Zz1Z+vrm0XNxtenIK2CJLmgjmeLtOyvmSGprUOKXvcVp5knqL5WJGUepV1Tt2Mq1UUOqiml64GnFLteaZx27cA2XwYgIt8NnEtnVsjHfqc24GbWK93hYrpflVOBbwAnpZQuAv4ceCpwCXBxSunZwAXr1KcgCBbMpk2b5vq3BPx34BQ6e2z+18R6Tf/PBH4D+ALwceATs/03Ad8H3C+ldMts36Cf8T6bqlYL2ls8pGBzS8pmi3fbU6o1+1lfuiT4NlTtwc6KteVz1LzzLfZFT0XVCil726xQa6XzWuzTnjKt2WxblkbReIq1JaJCU7PRt6Rge7OKWjrvcZ6memNK6RXzXLheg+rXge+klPaJCMCDZ/tPA24BbhaRU1NKX8M3EAdBMDGWRHXOw4dE5HXAP+YdKaXfb7lwvQbV3wFeKyK3AlcBjxKRN9PZVJ8L/D1whYjspQuBaGZlZeWomMHytVZXNVXo/eJbqtBLV9TqsLRJahXbV+Cj3Kfvm7Ha99JJawVV+uyAZV9aFlHMaIWkZwhWaUSvFF/No++ldVrPUp/j2W6tfS1KdYj670ujtlSj97xrA1jfLKBmR14LEx5Ufw74J+CRs/fNsn1dBtWU0j8BT6uc8o/AM9ajL0EQrB8THlTvTinN5d+JkKogCBbGhAfVG0TkJcDnmKnUvqD/zOQH1dJbOCSo2rrGS4+sJQjokCRdXd+annthLbXK8979W9I89fFan1qqLnl9siof6VVH9X6r4nyfWWRIlafyXK8ClBdqZfVJt2t9z/oqWlkrL/Q5rGrfwdq1reFX1mceY0AcMzxrndkCPHz2D7qB9fgYVIMgWF6mplRF5FdSSm9KKT2rdrzWxuQH1b64Nk/9tPxC19Sa56Dyappa12RaaoBqLHXitavvb4VhDQmP6kvJraWG5veW06yv4Iz1N+sLpbL6oteFqjlkdF9aHFX6OdXWg6o5sYYyT0pxTZWOMSBObVAFXiQi93eObaIr4nJsD6pBECwvExxUL+s5/vK+Bo6JQbVFdbb8cfM5np2rFrKlFeqQlNCMZePT68F797c+R+6Td+1Q+kJ5LPWmbZ1amVmlFz27Zc1m2JdQUfahzy7aYlP1iuLo60usmUJfQkXtu+cx5Nyh5w1laoNqSun31trGMTGoBkGwnExtUB2DY2pQHcsu5AVvWwqgz5ZapoR69tdMPrdMw8z7tDKykgtagtk1+RqtaluUi066sOyY5bphJbVIhb6SdlaxcE8tWxEhur0h9sW++1mfQ1NLU/Vs27UklZY0bW//kGLk8zDhNNW5OaYG1SAIloupKlUReQidU6qsUnV5y7WTH1QPHTpUTSMdw5taK17iRQFYNtW+0nzaplvua1k6wysQYpXZ09d4+6042L5CJJYd00u3bYmt7EvVta5tKdhSs6Wuhb7Y0JJWpdqiJIco1jGKvLQw1UEVeB/wP+nWvxrE5AfVIAiWlwlP/+9MKf3mPBdOflA9cOBAU/yfVpYtHn3Pw1++9myoe/fuPeJ9eY5nq60p1ZZF6zx7n7W4W6ZPrVmefL0ufM0GamU2lZSfdR5bp9WOR+2zQT1rSV/jFZcpr6k9d32u992rZc7pPtWUcN99a/1fCxNWql8QkZ8D/o7Daapfarlw8oNqEATLy4QH1X81+5c5RMPCgRCDahAEC2Sq0/+U0hNF5L7A9wDXp5S+2Xrt5AfVAwcOmCFKejqup+u1lE3P+VQ6nXK7erqvp/35fbnPm+pZIVA65Kll+q+LlQwJM6pN16002j76KuQPcSjp4yVeQZiW8CWvrZK+NNIWM5KVjKH/vtpk4yVwlAwZuIYkFxzP038R+RngVXQ1VR8tIq9IKf1By7XuX0NEHjNS/4IgOE7JtTmG/lsCXgw8NqX0U8D3Axe1XlhTqu8Tkbf3VWTZaEqlaoUveWmjVsETL+RJq1I4Wpnm7d13333E1rqmbwXWmqOnJZDf21opr/l1TjjQqseqgt8XhmX1qYV5/jP1OZ9qeI4kKzRPJzpoZ51Fn/OpPNa6LWkpcNMaSmWlzlopxENZ5AApIo8HLqULe/o03fJMZ9CtKPIiYCtwBXAb8MWU0lsGNH8wpXQXQErpThHZ03ph7Rv/2K7f8jEReeCAzgRBEAALV6pPB16eUvpl4CeBH0spvQC4mm6ZpvOBK1NKzweeLCJD1rP/ZxG5QkR+UkSuAP659UJXqaaUdgIXiMhZwKdE5LPFsaVZ+mT//v2DbJ6eYrX21VJO+5RqzaaqC6rUlIZX2s5SjX3l77RdtmzfKxwyJDyqplS99bksPPulVQRbn6PvZ53bFxZlKVXdvrfOVYkXslULw9L219paat4zrYVJeQrVChccQ6kuuEj1G4HLROQ24P7ADbP9NwEPolNWwR5VAAAetklEQVSqOXj/djoF2+pw+kW6QfkcuuWeLm3tVHVuJiKPAF4L/BnwjuJfEARBLwtWqqcDr0opXQLsBh4w259Xaf7q7DXAKcAdfQ2KyONmL88Gvgz8KXAd8MTWTrlKVUR+HXgecGFK6ZrWBteb/fv3m6rTU6aegi3P8bblL3dWpH0K1VK3nlK1FEaf+rSKMHve82w3tdSDVm21vmh12/KfYIgHXyuuFg94633LdobYF7XKzc9ySHHpWrqwRy1KpaX9Pm+/LkdZvl52myqdIn29iNwB/CHwXSLyVuBkurFrB/AGETkP+EBKab/b0mF+FPgbOtNCySFGWE7lccDjUkrfamkoCIJAs8hBNaV0PfCUyinfZuAqzSml3569/FRK6aq8X0Re2NpGzab6M0M6s1GUSrVmv/Tel9d46tZSI15kQEucat5qe1atuIgXBVCW1vPWvNfK1EoN9Yq8WH3x7K8t8aM1G6RXcrElzdO7Zp7/1DWlqhV9S+k8r/BMrZ/aJlyqRs+jbz2nvpjZmhIeQ6lOLfhfRJ4O/HvgiSKSM6hWgMcAV7a0Mfng/yAIlpcliTkdwkeBrwH35bD/6CBjeP+DIAjWytSUakrpdjrH/J+JyKl0S1VvAh5K5/zqZfKD6oEDB0zDunYy6am2Ve80T9X37NlzxDXWmlJ5nw72z9daCQN6n1c5y8JzOm3duvWec3I72QzghWqV6kFP9VqqeGU851lLDVbdRnmunmJr00Ftqt3yn1hPibXjqmxDT+H19Dl/D0ozTJ/D0AqPynipoeVn1s+5lj7spV7XzBf6nLUwQaUKgIhcDTwBOBE4gU6pPr7l2mn9jARBMClWVlbm+rcEPBJ4FPCx2evmjKrJK9Uy+L8WwOwVVCmv8ZxNtYItWaHu3r37iPc6KaC8xutbLSRGO6ys8J+sWj3FaqXBeumQua+l+inXzrKohVi1OJC0QtXntjjRWoqBeIrY6r+ntGphWPpaLwmg794eXuiZ1QevOFBtlhQFVYCuSPUhETkxpfRNEdnaf0nH5AfVIAiWlwkPqn8rIpcAt4jIHwHNqWGTH1T3799frfzvhdxYQc/azprto1aaqmdT1ckAtaBqr1hGrcp7Vp+5jVI99tlQLVXi2VCtvmn7oS4qYiltL63WSrX0Au31Z2+pVl8r36fDyPTx8r7aTtmSxODZKy07tZcIou3VVhicfoa1kCpPoVrqf0hCQx8THlQvo1v0bw/w74C/ar1w8oNqEATLy4Jz/xfJ54EPAVellD405MLJD6oHDx48yk5avvbWkqqpT21b1SrUOsfz+tfWtepbP0u/hsNqTadJlmhvcIuKa7ER9l1jKUxPodbSYLWCHPIfsxa9ULPJ6n5n9Ll9RWWse9dK6eliN95swypc7hUwt9AzM28mZ12zFiasVL+PLgngDSKyHXh3Sum9LRcuhZstCIJjkwUXVFkYKaW9KaX3A6+jK8TSvLLq5JWqV6S6z/tvxZF6x7TdFHyFqoum1GJnPe+/VXouf9G0TbVUNvqcvNWKda22spYCzbpP3n8UKxJBt+up3ZK+Yizldd5qpJaS9NJrh6TBDlkSxZtNWKmt+u/b8nxa0pHH+p54fZkCInIZ8DTgc3Q1Wf+i9drJD6pBECwvSxJzOg+3Az+UUvr20AsnP6i2LvzXt1CftU/bUnMsqnXMKwFYKuGsQvQCgDXVoD272vNekhVL7ou2/1nec0+NWEpM22q9BQVbbKpW/70iyy1K1aNUvV7RlXKGo+/flx1l2Vi9uNRagRJPoVp2av131Gra6osX1VH7zo1hU53woPpB4C0icj/g/cD/TSl9tucaIGyqQRAskKnaVOmKqbyLbvWAvwDe3HphDKpBECyMCaepbk8pfQI4lFJKLGuaqoi8ly726yGMtOphX5qqV09VT9fLfV5N1NJUoB1UOk3VqqfqOQlqaar6C6ZDtazpuXZU6fCrIdM6a2UBvdUrsbakeQ5RI0PMC3q/5fTrq9NatuWtA+WtHVaj5ojU72uOKt03b/0y/bpsN9JUe7lbRH4cWJ2t2jrKaqqjIiIvBu6avT1rxFUPgyBYUiY8/f8l4FnAdwGXABe0XrguSlVEzqWL9foM3UB+6+zQmlc97FtNVatMr/AJHFaf+pjeX7af29PnWO23pqe2pHnqsm9wtDLN276g97Ld2qoBuv1cwKVW+b8W7K/Rn7/POVTu85SqlRLqMUTd1oqkZDxFbIV59T0fK7HFC/NqcZrprZXifTwq1aJwyq3Af5ynjfWa/j+TbrCU2fusWPOqhyuz1zfSuOphEATLzwTTVBNwiK4wdf6ly68f1tLAugyqKaWfBZitargHeMAIqx4CnWK0bKqeQs2KUr8Hv+C0ZVP12tNtWCFVfQrV+nXXytFal2vbtm1HHOsr7FGiVWBWo2XxZS/ovKZUPaxn4RUr8Wy55TlDQre8QiRWyJlW+V54VE2p1v6+fSFU1jVeEWkrpdX7m9QK5xzPSjWl9N3lexE5Bbg9pdScCbGujqqU0nucQ4NXPQyCYPlZ5KAqIqcBl9PNgk+hU5QAjwXeBFwHvHq2/XxK6YoBbZ8FvJWu5N/7ROSGlNLVLddOPvi/VKpWIH9Wkt7W8v579ljLpqqVqd5atiqvAIalSrTNrWbLq5WYs+4Dvh1We7fBD4SfJ3Uzt9VS+s9TxrW+ZCz7on6mNW+3l0LsrapaHmt5/vocnarbElWgP7Nls/X+JjWlOsZqqgtWqo8Azgb+GvhMSukNInIG8NKU0tUi8jK6RfxWgE8PbPtVwFnAnwCvAT5F51jvZSkCwoIgODZZcJzqjXSD6tOAHxeRE4BX0g2IANcAz6Hz5L96YNcPppRuo4tT3QPc2Xrh5JXq3r17Tftin42zpj61YrXKBOZz+gq1lNd4JfK0wqspMb1/yLmW91yrP60KLZuqp1Ct9r0+WerQU0a1lNC+YiJWwRmvJKIVs6kVtY4PtuKFPfVZK16iaVlmRQ8+tdRf/b6m0iekVC8E3jVb9uROuvCngymlG2bHnwB8OaW0R0T2ua3YXCcirwXuKyKXAjf0XZCZ/KAaBMHysuBB9d3A5SLyVeCzwKOBvy2OXw+8Q0R2AW8b2Pbz6FTuJ4GddPH0TcSgGgTBwljkoJpS+hxwbuX4tcC1c7a9H3j7PNdOflDdtWuXOT33wqK0g6q2RpWe2pfTIa8+Za3yv54GequdWqmnXgqqtd68nsJr51N5jVd7teao0tNwL3XTOkc7icrno6fa3npNlvmiZfqvYya9lM2a8yx/F7x1orx2vD7NsxqBR81B6ZlsamuRTWD6v5RMflANgmB5iUF1guzevbuqVL3QKu2MKl97wfNWqFOtmEV5rX5dopVeTv8s9+ljOdB/+/bt95yb9+Vz8lZfW7av9+lrLCXsOY7057CO6edTnuspRh1aZSnhlv+8norVCtlynnmqTTusav231lJrDYOz8MK9aiuv9rVVvi7/b8xLDKpBEAQjsiRl/NaVyQ+qu3fvPkphwtEhVFq51irzZ7RCKtvXtjzPFlna8bxyblpJZhtouU+r0KxQ836AHTt2HLHPU67lNfocLwmg7K9W6Z5d1noOtf9ktXJ01n2tPmVqdkVNS5qwR/5OWGmk3n2svvWlw9ZK82nbvxVG1le4xbKp9hWgCWwmP6gGQbC8xPR/guzatcu0AXlrUnl20xKv/F1pi8z31MpUp32WaLWgFV5WjTWlmtWoVqXl66xivXNrNtV879raS55ttRYxkI+1FI3xFKulSrXybbFJev/Rrb7pYHnv2hZPuU46KNHfQ/3ZrciTfI5OfmkphKLbt+zIYyjVGFSDIAhGJAbVCVIqVav0n44M0OmltZVFtRIrj3uFSPL9LKWqVZWnUEv1qW2n3ns4Wpnmrefht/a1xLQO8f57sa0WXgxlxlKLfemvNfuixuqbF82hlXd5bV+KaS1d1SvTWIsj1d9xS9Xqe+rZQG0porUQg2oQBMGIhPd/gtx1112m3anvV7xWvFi/t5RqVnS5XW/pkppXWC9LYinJrES1QtU21vL1CSeccMQ1eb8VMdDn9a8VX/FiUGs21VoREG1z1s+95tH3zmlZQsZ733d9eU35t9XPpSz009e+/hyW7b8ve8/6f6BjWPPWUrcRp7o2jr+fkSAIggUyeaUaBMHycjwq1ckPqjt37rznda3KvlfkolYHU2M5SPSUWDuoyulznk5pR08tOF9P4b2wKfCn/Xm/NgNY965V2fccVF6olXWNLmpSCxVqwQuW95xdJV5gfNlHz1mmsRxVuS/6O2EVVPH6osOmymt0kL5lKvCO1ZxbVtr3vMSgGgRBMCLhqJogu3fvriocL0Slth5UxkuThKMdMKVzybumDDWyrtVOKDjaUaXVZqlUvXNrwf/agVRb88lTprXye167mXlUoYX+u1qOwr50VEtVed+JWtqnPqbVuuWc0+17KajlPq80ZS20sCUMy1pJY15CqQZBEIxIDKoTZN++fdWiwC0KNTNkPSV9TVZ8XmqldY22qVoB/bqAimdjLV/rZIJamNSQlM0+ZVorM9dSrKSmkktqCRs1+j7rkOLOLd8JL9SsNjvqU6FwdCiVPqe0hWpF6ilWq1h72FTnY/KDahAEy0sMqhPkwIEDZkC/t6Jli5ry3tcKTnuJAtrbXfZBK8lacL4O0rcSBTxFaq2M2ocVyO8V5m5xRgyxaQ8pktKX3ml9pnnw1LMVZeA9J6vYuYcX1VC+HqJq+xSqVYxojOVUjkcmP6gGQbC8hPd/gpRK1fo1r3nw9XFPKVnXeopIx2paS1voYiwtcapesepSfWrl26dc9euyj1aR7b6lUWq20CFppF67VpGUITZVr/0xrrWKket+15bl8bCem6dQLa99q0KtRReshUVO/0XkNOBy4HZgE3AfYBU4RLck9W7gJcAu4MMppQ8urDMFx9/PSBAE68amTZvm+tfII4CzgdOBm4DHADuBu4AvApcAF6eUng1cMP6ns5m8UoV6+TJPIdUybbR6aym04WUMWe14irJWzk97/60ygVr5emX8WkrzaQ9/eUxTs1N7XvIhC/ZpWjKuajZb757eMi61dmvef2+JbctO3Wevtvqsv/fWYpZeDKtnay3bHZLZ5rFgR9WNdIPqV4CPAi9PKX1ERH4CeCFwv5TSLbNzh09p5iSUahAEC2PBSvVC4KSU0iE6hXrmbP+3gK3AzSJyau7KqB+swjGhVIMgOC55N3C5iHwV+Axwpoi8mc62+mvAfYErRGQvcOV6deqYGFSteqdWKIr13sKreN6SsqnTPq3iItpRpYP2rel/y8qoXthVbbVTb/pvOe1aQ6lqle1r5/RRu6YlDKsvmaBl2u+tDmG15zlFLeeoV6PWMsNok4q3hf6wq5rJY4yp+yK9/ymlzwHnVk75BvCMhXXAIab/QRAEIzJ5pbqystL0i9qy3k6fo8JK7/ScQC2q0FOdllL1Vgeorbyqw7FqIVUtStVLt6ylefYVIJlHUVrn9oXOWe30bedlnnb6Un7LtvoSHaxzdem/liSJMZRqZFQFQRCMSAyqE2Tz5s1NYS3a3jSPirAUXt8aT1Zwvlc2UIdNldf0FUspX3tq1gro10rVC62qnTOEISmafWnDQ/EUaU299fWhps7HUr591PrYpzprz7QWHjhv344HJj+oBkGwvMSguiBE5AeBi+gyHW4A7gVsA04CzgeEOdPJSqVa/gG1+tC2KiuNVRcR1p5XSxVqtam35a+9trNqZaltoHC0QtVRAJZN1VPPlifZU6qW0tb912XwMi32zBp9tlpr6ZKW9j2F6rVVnqv70lJsp8V+OQ/6+XhbsGccVltW5EzYVOdjvbz/JwPPmaWL/TBdwO5FwJ8DT2WD0smCIAjGZl2UakrpGhHZJCIvBT4F5J/1m4DvYw3pZJs3bzZ/WXU8qv7F9gqulO14MYTlPq0GW5Sqjh/VdtKW2FMrTrXPrmt9Di+KwSoXaEU09OEpulrx7oyn7CwlpmcrQ1JO9bUtcbY1ldsXP2qd6/XT2q//NtmzrwulW+23zNjGJJTqghCRewNXAX8JvAt48OzQacAtbFA6WRAEi2XBaapLyXo5qt4MfC/wLOAA3SD6Zjqb6nOBv2cD0smCIFgsUx8g52G9pv+/2HPKPzJnOtmWLVvMKWWeEnlpfy3Bz7VAeC8syqsQVe7rq+ZvrTvlOcKslVG9OqpWeukQR9VaQqnmYYhTSzsZaxX5tWmoZdrfF2hvrTrRlxpqnasrT1mmLe9vZfVRX6/v4zkbvX1DiUE1CIJgRGJQnSCrq6vmr23el1WCdq60OKhaAqO100mrzvK+Xh1VvTaV5ajSzqxa5X/P6WSlqXpKvuVcL/SppTL/kCSA2jVaNdfqnHoqVtc9LZVsXxhWS2X+FkeV57CqFbbR30Hdp7Jd3U6tT2OGgB2Pg2oUVAmCIBiRY0Kp1sJOsnobUt6t5ddVl/HTtlVL6WklqVWnLvNXHsv7vAQCqy96W1t3qq+EYXmsb42qFtVTe9ZrWXeqZbbhrULQku6s+2gpzL41pFpWRtXK0vo75O9GSyiaXpPKC60qrw+lOh+TH1SDIFheYlCdIF4guqemaitbDrEV6l96z+NuBc97UQA6XRWOtqW2eOf7FGqtnN+iV1Ot4RXB0c/fKhjSsgKu16faWkyeQtXXWuozq0JPsZavre8j2Eq1/Jtb1GYKLc90nqJDHjGoBkEQjEgMqhNk8+bNpkKybIJw9C92rUxai1fb85pbSrJ1lVOrcMuQlVG9PlkK0yuObNFqt5z3P5JXOKUlprjlfWt0R80Oq5WppT51YeiWOFX9vmVFU20LtmYXOnXVU6FDUmeHcDwOquH9D4IgGJHJK9WtW7dWF9nz7H8tUQBDlv7QatBSnV78aM3Tq9VHLQ5Wq8yWotJ92WM1pTEkrtejRQ1pu19NfY6BNSPx1Ge2m1pK1VvKxLK/emq8FkeqqSnMvudsFSMPpTofkx9UgyBYXhY5qIrIacDlwO10hZi20o1pDwAuBh4CvBq4Dvh8SumKhXWmIKb/QRAsjAVXqXoEcDZwOnAH8NGU0vOAPwDOAc4CvkY3zn16/E9nM3mlWk7/S6eNNzVtWR8+0xJqo78AOjDbchp4Afe1wi162q9DrKx2awVhNC3Pqy9NVbdVu09LhX6vPSskrK+N2jk155Cesufpft7u3bvXvcYzGVhTec+ZZTm3tMmglnLqraJaq2s7ZuGcBU//b6QbVL8CfBR4vYicCfws8Gy6ynhXAncDH5mdu3BCqQZBMFUupFtF5BBwJ/Af6JZtOi+ldCfwBOBgSmkPsG+9OnVMKVXL2K7f1xww+tc8/8oPSXEdUjnfc27VVjBtCany+jYk/bYl/GqIUh2iWLSaHRImlamFzHn3sf7OrQ6q/B78gipWlX1d+s/blo6wvvZr4VGeY6z87ugZ1FpYsFJ9N3C5iHwV+CxwBfAx4J0i8j7geuAdIrILeNsiO1Iy+UE1CILjk5TS54Bzi13/yTjt2nXqzj1MflBdWVk5KkAejk7lq61nn9F2M10gxLKxevZRy+ap97Woz761oyylOkQdeOFKVnhXn1Jda3hTLXTK2z+PWu6z69aKo2jlmG2qpVLV6amePdPap1Wn3pb39pRqLbxrSLGUMVRmhFQFQRCMSAyqE2R1dfUoz7h+DX5RkfKPnn/p8zlZhWjbqoVnU20J/veKsZSvW1Y71eq7xdPufY6aKuw7dx6lbPXXu18tiaElqsDbtnj/tRrMqvTuu+++5xptb63ZPL3+1pSqp0xrNtXMEPU/hk31eGTyg2oQBMtLKNUJsrKyUvX+a9tgVnYt9jn9i11TGFpFtSxH4hWgtpZI8WJZa3GqmpaCJHp/eXwtpf88StXZtxhd7W/WklLcuoifpQr1NTXvvLazasXaUsS7Fnuqbf+ePdZq14uNtrz/wXxMflANgmB5OR4H6BhUgyBYGDGoHoN4U9aW9Dx9rjUF62vDqtjeVyO1Fvyvp20t5ovMPGE0Q9JUW/4D1abp3rTfc1y1tF8LX5onzVM7fqz9OoVVO7Vq0/++1VvLfV7IVs3M45m/an/nYBjH/KAaBMHGEUp1gqyurpp/OP2r7RVSqSk9/ateqlOtqrwiJlZBlb41pFrWhRqSRqqPD1lR1mq/z2FV0qeOW/7TzfMf07pvXyB8zVHlpXfqNNbyteewshxJNbXsfbaWNba8tcdqRXHGLKhyPBJPLQiCYESOCaXakp6Zf81rVc294HPL7uSpvlqapxfsr22pln2rZgvW9+5TdrXjawn+r+Ep1vLaIbbTvnat4312y1rwvFakWn1aSlUnBlil/7wiPnr/kNmApTozLTMea8Y0LzH9D4IgGJEYVCfIli1bBpXBy1hK0kvpsxSkZ2PT11h90kpV98UKxPa88mN9aVtSTluLl1ipp5ohK3a2FEJpKY7iRQboYHpLdeZtVp15m+2meWvt01vL+++V87PoK7Ruqc78nfMKCllFfMoU66CdyQ+qQRAsL6FUJ8jq6qoZ3+mpP/2+ZjOsFdqwUlfL9q1yhHqJFS8KoMWjv+gv6zz2Uusafcx7P6QPQ+ywtZVRvSVSrDJ+e/bsAWD37t1HbPN+S6lmNZvP0SUBrb7ka2vPtM8rb0Wc6IgWL2qlfB0FVeZj8oNqEATLSyjVCVJ6/2sxoX3FUkq0vc/K9PFiBPXCf7XSf1oR6P3l6xbvv+5bS0bVmF/6mmocshBfX59abKq1pVH6yviVqjMr0l27dpnv87l5PxxWplrFWgWttc1WZ+pZz8KLdsnfDUsJ52P5PrqIe3j/xyPiVIMgCEZk8ko1CILl5XhUqpMfVEtHVW3l0kwtxVJPn2oFKjwHlXZU1WqjekH/tWlXrSDJkAr/3jUtiQF91w6Z6g+ZylvtD6nir6f9XshTOZXP0/y83blzJ3B4aq+dUbVjLdN/vcqErgMMflEdK4zPM3nUkgvWyxl6rLI0g6qIPJhuidnbgC+mlN6ywV0KgmCNLHJgFpHTgMuB24FNwNeAM4CTgBcBW9mAMWVpBlXgfODKlNKnReQjIvLOlNK+vos2b95cVYVeGTPL8TMkPMdzEmjVWTqqvJVWver+5WsvHMjC63+LkvSSGbx9Q2m5t9fPlmr+tZVFdXiU3monFBxWpnmrHVY6fKp8nc/RiQKlUs3HLCcZ1FeDsByb1vuSltVVhzhF+1iw2n0EcDbw18DfAWellM4VkScCzwW2M8eYslaWaVB9IHDj7PXtdL8239y47gRBsOTcSDeofgX4OHDTbP9NwIPolOq6jynLNKh+FTiN7iGcAtzRctErXvGKBXYpCII1cIOIPHTeaxvOuRB4V0rpkIh8G8j3Og24hS66afCYslY2zePcWAQi8kDgDcCdwN+klH53g7sUBMESIyL/GnglnSD7f8BeQICTgecBO9iAMWVpBtUgCIJjgQj+D4IgGJEYVIMgCEYkBtUgCIIRiUE1CIJgRGJQDYIgGJEYVIMgCEZkmYL/m1nmOgEi8oPARcBddAHM9wK20WVznE8XR/cSYBfw4ZTSBzeoq0cgIu8FPgQ8hCXIn3b6eAbwMuBWutjD+7DEz3ZZc9MtRORM4H0ppe8XkV/t66c+J6X0LxvS8SVkqko11wl4PvBkEVmmFcpOBp6TUno28MPASSmli4A/B54KXAJcPDt+wcZ18zAi8mK6HwHo8qdfAFxNlz+9TM/6YroUxFOBb7D8zzbnpp9Olx65lM92lnjzHGCniGxv6Oe9jXOCGVMdVK06AUtBSuka4C4ReSnwKY7OR75fSumW2b4Nz7wQkXPp0vc+Q/d9uHV2KPd3mZ71mcB/o/tP/EyW/NlyODf9acATOZx3vlTPNqX09ZTSpXQ/rKfQ/x042TgnmDHVQTXXCYB1zOltYfYrfhXwl8C7gAfPDuV85JtF5NTZvmUoWPlM4AeAX6AbrO4/25/7u0zP+uvAd4pKQ8v+bC+kU9OHACs3fZmebeZW4L6z114/bzHOCWZMMk11mesEiMi7gO+ls6ceAG4GTqRTIc+dHftNujzl96eUPrxBXT0CETkP2AM8gCXIn3b6+Ei6XO9b6ZT1o1jiZ7usuekeIvLRlNKTROQievqpz0kpfXuj+r1sTHJQDYIgWFamOv0PgiBYSmJQDYIgGJEYVIMgCEYkBtUgCIIRmWRGVTAtROQc4PXA41NKu0XkQcDHgCellG7e2N4FwbiE9z9YF0Tkd+hSdi8CPgFcnlL6HxvbqyAYnxhUg3VhloL5SeBbwGdTSq/c4C4FwUIIm2qwLsyyoH4XOAd49wZ3JwgWRijVYF2YLVX8ceDtwE8DT0wpHdjYXgXB+IRSDRaOiGwF/itdibg30qVtvnxjexUEiyEG1WA9uAL4ZErpI7P3zweeLiI/snFdCoLFENP/IAiCEQmlGgRBMCIxqAZBEIxIDKpBEAQjEoNqEATBiMSgGgRBMCIxqAZBEIxIDKpBEAQjEoNqEATBiPx/V5Iv2VbzR48AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "geo_model = gp.create_model('FaultModel')\n",
    "gp.init_data(geo_model, extent=[0,1000,0,1000,0,1000],resolution = [3,3,3],\n",
    "                         path_i = r\"data/faultuncertainty_interfaces_surface.csv\",\n",
    "                         path_o = r\"data/faultuncertainty_orient_surface.csv\")\n",
    "\n",
    "gp.map_series_to_surfaces(geo_model, {'fault_series':'fault','form_series':('layer2','layer1','basement')})\n",
    "geo_model.surfaces.colors.change_colors({'layer1': '#ffbe00', 'layer2': '#9f0052', 'fault': '#015482',\n",
    "                                        'basement': '#728f02'})\n",
    "\n",
    "geo_model.set_is_fault(['fault_series'],change_color=False)\n",
    "\n",
    "geo_model.set_section_grid({'s1':([0,200],[1000,200],[80,80])})#p1,p2,resolution\n",
    "geo_model.set_topography(source='saved', filepath=r'data/topo.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<gempy.plot.visualization_2d.PlotData2D at 0x1cbfa22efd0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAGtCAYAAAClVis3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAF/VJREFUeJzt3XuQlOWd6PFvAyp4AYSAigRcFR9xwEQNGhRZGdCY8pbEmLNuPCfuKguCSrygUF4SMKmjEolSx0RXzW7cyuaU2awayG68TQxlanbVzRpl4nkqwXAPEsSo3Bwuff7oRh0ZBuJMd89v+H7+Gebtt2d+1dX217ffp98pFItFJEmKoFutB5AkaU8ZLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhg9KvWDU0pHAz/KOZ+QUpoOHAH0Aa4B9gXuAtYBTTnnez+8T875j5WaTZIUU0WOtFJKhwKXAxtSSj2BsTnnqcBDwERgEjAv5zwFOCeldFAr+0iS1EJFopVzXp1zngGsB/oBa8o3rQAGAYcCy8vb3gQObmUfSZJaqNjbgx+wBuhf/vdgYBWlWA6mFK5+5W0f3me3UkrFnHOHDitJqqhCe+5c8WjlnLemlH6eUvoOpSOqyUAvYG5K6VLgX3exjyRJLRQiX3vQIy1JCqddR1oueZckhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQ60LPPPsuKFStqPYYk0dzczI9//ONaj9HhjFYHaWxs5LOfGcfJJx7L5s2baz2OpL3c6aecyhe/+EVuv/32Wo/SoYxWB1i3bh3nfvYvufV/wjGDNjB2zKdqPZKkvdjVV11FfmkRUziFr8+8hYULF9Z6pA5jtNpp+/btnHrK8Yz7xBZm/DX86yxYsaSJKVdMqvVokvZCjzzyCH//f77LtZzGKA7n8wzn3Alns3bt2lqP1iGMVjt96aLPU9y8kodnQKEA/XrDE3fAw//49/zwhz+s9XiS9iKLFy/m0osv4SucyBD6AnAWR3Psln6MGnkC27dvr/GE7We02uHuu+/mqSd+ws/ugP17vr995JHwvRvg7y77Mjnn2g0oaa/R3NzMqSeM4rTtQxjNx9/bXqDA5XyK5tVvcf4559Vwwo5htD6ixsZGbp55DT/+OvzFYTvf/qUzYPJ5RcadfpILMyRV3OmnnEqfd7rxV4zc6bZ96c61nMbPf/ZU+IUZRusjeG/hxf+CCSfter/bJ8KwQRv4yzGjqjecpL3OjoUXV/Fpuu/iZb0/+3Mlnw6/MMNo/Zk+uPBi+v9oe9/u3eHRWbB8ySKmTr2iOgNK2qt8cOHFQezX5r7DGRB+YUaPWg8QzUVfvKDFwguAhl/BjQ/Ab5bApmbotS8cdwTcMRHqTywtzBh95X2MGTOWiy++uJbjS+pCfvvb33LpX7VcePEqa/gRi1jJOzSzjX3pzuEcxEWMYDgDOYujWbxlHaNGnsDilUvp1i3WsUusaWts7ty5PP3EghYLL2Y/DOfdDOM+CS8/CM1Plr6O+2Rp++yHXZghqeNt3ryZ0046hTHF9xdePM6r3E0jxzKA2xjPg3yO2xjPsQzgbhp5nFffW5jxbtCFGYVisVjrGT6ylFKxWhFobGzkzPpTeey2989jNfyqFKan58Doulbu0wQTpsP8b5SOuK7/Lvzzswfw2rK19OzZc+c7SNIeOuWET/HmS0u4kdPpTjdeZQ1308h0xnA0/Xfa/3e8wRye46uMZjgDWctGbuUZbv3fs5kxY0Y1Ry+0584eae2BDy+8WPc2/OerMOnbMPWC1oMFpe1TLoCv3gub3oU7/s6FGZLa74MLL4rAatbzMC8xniNbDRbA0fSnniP5J37Nm2yiP724KuDCDI+09sAnRhxF06uvkT4OS1+HLVuhz0E9WL9hK688BEcdvuv7Ll4JIy+H5i3Q90D4WB/Iy+GGG67njjvmVHx2SV3LU089xWfP+gx96cU2tvM277JfoQdbi9v4BhMYyIG7vO8a1nMzz1AAtlPkYHqxgWa67dODtzavr9b5rXYdabkQYw9M+MznOHn025x44omMHTuW4cOH061bN7p1KzD00LbvO+QQ2NwMmze/y/PPP89zzz3Hr3/9a0aO/ER1hpfUpfTv35/zP/856urqGD16NGPGjKF3794UCgX6s3+b9+3H/mxlG9uLRZYtW8bChQt54YUXeP3119m+fXuIRRkeabXDAT0LvPzg7o+0jr8cNmyO+zhL6vz2K/TgNsbv9kjrFp7h3eLWKk62E89p1cpxdXXcv6Dtfe5bAHUjRlRnIEl7rePq6niWJW3u8yy/D/96ZLTa4Y4587j38dIqwdY0NsF3Hofb77ynuoNJ2uvcNe/bPMNifscbrd7+O97gGV7jW/fMrfJkHctotUN9fT03zpzFhOkw/f7SW4Fbtpa+Tr+/tNz9xpmzqK+vr/Wokrq4+vp6bpp1K3N4jkd4hTWsZyvbWcN6HuEV5vAcN826Nfzrkee0OkBDQwMzbphG06JF710Ro27ECG6/857wTxBJsTQ0NHD9tGtpWrSILWxjH7pTN2IE37pnbmd5PWrXOS2jJUmqJhdiSJL2DkZLkhSG0VLFHdC7D0OOHMY3v/lNNm7cWOtxJAVmtFRxgwYdzvI1b3DznHkc0LsvRx6TuPPOO/2LzpL+bEZLFfej//vPsHk9nHYJnHUlv99vKDd+81v0OrA3w449jrlz59Lc3FzrMSUF4OpBVcWnThnNf63eROHkC9/bVly/DlYsgiX/DZve5phjjuGqKZOZPHkyPXp4WUypi3LJuzq/pUuXcsRRw2D8JAq9B+50e/GdtbD8FVj6EoV3NzB8+HCunnoFl112mQGTuhajpc6voaGB879wIRvWr4dtW6H7PnDQADj+LAoDj2yxb/HtP74fsC2buO1rt3DTTTfVaHJJHczPaalzmz17NuPPOpsNh9TBmVfCF74GZ06FgX8Bv/wBxd/8vMX+hd4DKNTVw4TJFItFhg4dWqPJJXU2HmmpohoaGhh/1tkw9isU+g/Z6fbiG8tg4ffhtC/vfMT14mMc06uZ/JtXqjWupMrzSEud19XXXAdHndxqsIDS9qNOhleearG9+MZyWLGIJ/99N3/7RdJexWipopqamuDIUW3vdOQoeHvNe98Wi9vhxUeZ+LeX+tagpBaMlipr2xbYv0/b++zfp7TfDr97nv0L27jvu9+t7GySwjFaqqzu+8DGt9reZ+Nbpf2A4uZ3oOlpfvhP/0i3bj49JbXkq4Iqqq6uDl57oe2dXnsBdnx266V/54RPfpLzzz+/8sNJCsdoqaLmffsuWPx8aZVgK4pvLIPFz8PIMymu+T2F13/Lz346v8pTSorCaKmi6uvrmXXrzbDw+xRffoLi+nUUt28rfX35idJy9zQGPjYEXnyU6796NQMH7nzFDEkCP6elKmloaGDatdezaNGi0qKL7vtAsQgnf4HC4BEUX11In9Uvs27NHzyXJXVt7fqclhd1U1XU19fzyku/arHtwH4D2dCtB8UNf4L/9wsee/pJgyWpTb5CqGZGHnsMrFsB/z2f08ecxhlnnFHrkSR1ckZLNTNh/Dh47QW6rVvBTx57tNbjSArAaKlmLrroImjexDe+fgt9+/at9TiSAvCclmrm+OOPZ/78+Zx77rm1HkVSEB5pqaYMlqQ/h9GSJIVhtCRJYRgtSVIYRkuSFIbRkiSFUZUl7ymlwcBs4E1K1536A3AE0Ae4BtgXuAtYBzTlnO+txlySpFiqdaR1LFAPDAHWAmNzzlOBh4CJwCRgXs55CnBOSmmfKs0lSQqkWh8uXk4pWr8HngRWlLevAAZROtJaXt72JqUjsLVVmk2SFES1jrSuBPrknIvAW8DQ8vbBwCpgWfnfAP2AP1VpLklSIFX5e1oppROBWZTitARoBhJwMDAZ6AXMBd4BXsw5P7CHP9e/pyVJsbTr72n5RyAlSdXUrmi55F2SFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYfSoxi9JKR0B3AKsAd4B+gL7AX2ASUACZgIbgQU550erMZckKZZqHWldB6wADgNeB/rknKcBvwAuBK4Hrss5XwZcUaWZJEnBVCtaRwOPAROBSygFjPLXQcCAnPOq8rZilWaSJAVTrWitBt7OOW8pf394+etgYBWwMqV0WHlboUozSZKCKRSLlT+wSSkNB2ZROqfVCNQBB1A6pzURGAbcDDQD/5JzXrCHP7eYc67IzJKkimjXgUlVolUpRkuSwmlXtFzyLkkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQpjl9FKKTWmlD7dyvafVHYkSZJa19aR1seA76eULvvQ9oMqOI8kSbvUVrSWA6cDX0kpfSeltGPfYuXHkiRpZ22e08o5rwHqy98+m1IaUPmRJElqXY82bisA5Jy3AlNSSn8D/LIqU0mS1Iq2jrS+9sFvcs7/AHwZWFLJgSRJ2pVCsRj3FFVKqZhzrvUYkqQ9V2jPnf2cliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksLoUc1fllL6ATAf+DhwBNAHuAbYF7gLWAc05ZzvreZckqQYqnaklVK6Flhf/nZsznkq8BAwEZgEzMs5TwHOSSntU625JElxVCVaKaXzgD8BjeXfuaZ80wpgEHAosLy87U1KR2CSJLVQrbcHL6EUo1T+fscR12BgFaWQDaYUrn6UAidJUguFYrFYtV+WUroU2AwcQilgBwOTgV7AXOAd4MWc8wN7+POKOefKDCtJqoRCu+5czWh1NKMlSeG0K1oueZckhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUnqQhoaGhh10ggO6FmgW7cCB/QsMOqkETQ0NNR6tA5htCSpi5g9ezbnnTOeccOaePlBaH4SXn4Qxg1r4rxzxjN79uxaj9huXjBXkrqAhoYGzjtnPE/PgdF1O9/e2AQTpsP8nz5DfX199Qd8nxfMlaS93Y3Tr2bqBa0HC0rbp1wAM26YVt3BOpjRkqQu4DdNTUw6t+19Jp8LTYsWVWegCjFaktQFbGqGoYe2vc+QQ0r7RWa0JKkL6LUvLF3d9j7LXi/tF5nRkqQu4Li6Ou5f0PY+9y2AuhEjqjNQhRgtSeoC7pgzj3sfL60SbE1jE3zncbj9znuqO1gHM1qS1AXU19dz48xZTJgO0++HxSthy9bS1+n3l5a73zhzVq2Xu7ebn9OSpC6koaGBGTdMo2nRIjY1l85h1Y0Ywe133tNZgtWuz2kZLUlSNfnhYknS3sFoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCkMoyVJCsNoSZLCMFqSpDCMliQpDKMlSQrDaEmSwjBakqQwjJYkKQyjJUkKw2hJksIwWpKkMIyWJCmMHtX4JSmlU4FpwHpgKXAgsB/QB5gEJGAmsBFYkHN+tBpzSZJiqdaR1sHA5Tnny4AxQJ+c8zTgF8CFwPXAdeXbr6jSTJKkYKoSrZzzT4H1KaWbgF8CK8o3rQAGAQNyzqvK24rVmEmSFE9VopVSOgh4EPgP4HvA4eWbBgOrgJUppcPK2wrVmEmSFE9VzmkB9wDDgL8BtlGK1D2UzmlNBF4C7kopNQPzqjSTJCmYQrEY9924lFIx51zrMSRJe65d76a55F2SFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYRgtSVIYRkuSFIbRkiSFYbQkSWEYLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYfSo9QA7pJQOB+4C1gFNOed7azySJKmT6UxHWpOAeTnnKcA5KaV9aj2QJKlz6UzROhRYXv73m0CfGs4iSeqEOs3bg8AyYDClcPUD/rQnd0opVXImSVLHKuacCx/1zoVisdiRw3xkKaVDgbnAO8CLOecHajySJKmT6TTRkiRpdzrTOS1JktpktCRJYRgtSVIYRkuSFIbRkiSFYbQkSWF0pg8X7zGvU7izlNKpwDRgPbAUOBDYj9KVRSYBCZgJbAQW5JwfrdGoNZdS+gEwH/g4cASlx+gaYF/28udVSukI4BZgDaXPTPbF51ELKaXBwGxKV+4pAH/A51ELKaWjgR/lnE9IKU1nN4/Ph/fJOf9xVz876pGW1ync2cHA5Tnny4AxQJ+c8zTgF8CFwPXAdeXbr6jdmLWVUrqWUtgBxuacpwIPARPxeQVwHbACOAx4HZ9HrTkWqAeGAGvxedRC+UIRlwMbUko92f3jc1Ar++xS1Gh5ncIPyTn/FFifUroJ+CWlFx7KXwcBA3LOq8rb9spPlKeUzqN0ebBGSs/9NeWbdjxGPq/gaOAxSi8cl+DzqDXLKUXrS8A4SuECn0cA5JxX55xnUPqfw37s/r+zg1vZZ5eiRmvHdQrhz7hOYVdW/r+VB4H/AL4HHF6+aTCwCliZUjqsvO0jX/cruEuAk4GvUHpRHljevuMx8nkFq4G3c85byt/7PNrZlZSOQIvAW8DQ8nafRztbA/Qv/3tXj8+qVvbZpZCXcfI6hTtLKX0PGEbpfNY2YCVwAKX/y5tYvu1moBn4l5zzghqNWnMppUuBzcAhlM7RHAxMBnqxlz+vUkrDgVmUXmwagTp8HrWQUjqR0mO0DFhC6bHwefQhKaWf5ZzPTilNYzePz4f3yTm/taufGzJakqS9U9S3ByVJeyGjJUkKw2hJksIwWpKkMEJeEUPqilJKZwLfAj6dc96UUhoEPAGcnXNeWdvppM7B1YNSJ5JSupPSJbimAQ3A7JzzU7WdSuo8jJbUiZQv+fMc8AbwnznnWTUeSepUPKcldSLlK1E8AJwJ/EONx5E6HY+0pE4kpTQUeBK4D/g8MC7nvK22U0mdh0daUieRUtoXeITSn2b4NqXLBH2ttlNJnYvRkjqPu4Dncs7/Vv5+CnBxSumM2o0kdS6+PShJCsMjLUlSGEZLkhSG0ZIkhWG0JElhGC1JUhhGS5IUhtGSJIVhtCRJYfx/Jv4n7RGJxGcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gp.plot.plot_data(geo_model, direction='y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.plot.plot_data(geo_model, direction='z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.set_interpolation_data(geo_model,\n",
    "                          compile_theano=True,\n",
    "                          theano_optimizer='fast_compile',\n",
    "                          verbose=[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = gp.compute_model(geo_model, compute_mesh=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_=gp.plot.plot_section_by_name(geo_model, 's1', show_data=True, radius=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "Geological map not computed. Activate the topography grid.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-4-a610e346e5fa>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mgp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_map\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgeo_model\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshow_data\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:/users/elisa/documents/git/gempy\\gempy\\plot\\plot.py\u001b[0m in \u001b[0;36mplot_map\u001b[1;34m(model, contour_lines, show_data, figsize)\u001b[0m\n\u001b[0;32m    169\u001b[0m     \u001b[0mplot\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mPlotSolution\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    170\u001b[0m     plot.plot_map(contour_lines=contour_lines, show_data=show_data,\n\u001b[1;32m--> 171\u001b[1;33m                   figsize=figsize)\n\u001b[0m\u001b[0;32m    172\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    173\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:/users/elisa/documents/git/gempy\\gempy\\plot\\visualization_2d.py\u001b[0m in \u001b[0;36mplot_map\u001b[1;34m(self, solution, contour_lines, show_data, show_all_data, figsize)\u001b[0m\n\u001b[0;32m    493\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    494\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0msolution\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgeological_map\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 495\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Geological map not computed. Activate the topography grid.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    496\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    497\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: Geological map not computed. Activate the topography grid."
     ]
    }
   ],
   "source": [
    "gp.plot.plot_map(geo_model, show_data=False, figsize=(7,7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Stochastic simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### as an example: code to perturb location of fault once. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000\n",
      "0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:/users/elisa/documents/git/gempy\\gempy\\plot\\visualization_2d.py:208: UserWarning: the orientations are not converted to apparent dip.\n",
      "  warnings.warn('the orientations are not converted to apparent dip.')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAGaCAYAAABueRZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VPXd/vH3TPZ9X4AACYQcAoRF1LqAUpeqRcW1LnWrQrVitRR9rK3WX7V7a1t9RKvW9qlr677hrkVFkdYqWwgnIRBCEiD7RvbM/P44CURkz8ycWe7XdXEFJpPkkyHJnfOd77mPw+12IyIiEmicdg8gIiJyOBRgIiISkBRgIiISkBRgIiISkBRgIiISkBRgIiISkMK99Y4Nw8gHnjVNc4ZhGLcAuUASsAiIBO4BGoFi0zSX7Hkf0zTrvDWbiIgEPq8cgRmGkQ3MB3YahhENnGCa5kLgUWABcC1wn2ma1wNzDcNI2Mt9RERE9skrR2CmaW4HfmQYxptAKlA78KoqYCTWEdjWgduagJS93OeA6rZvdkfHpXlqbPEGnScvIgP6+vpJSU1xeOr9eW0JcYhaYDBlcoAarCO/HKwQSx24bc/7HFBPTw8Z2YkeHVZERLyjrbXVo+/P6wFmmmafYRj/MgzjAawjreuAGOAPhmFcBbywj/uIiIjskyOQuxCrK033qDGG3WOIiMhBaGttJSEx0WNLiNpGLyIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAUkBJiIiAckXl1MJaPX19RQXF1NeXs6mTZsoKyujoaGBjo6ddHR00NXVhdvtxuGw+imjIqPIyMxk5MiRZGdnk5OTQ15eHnl5eeTm5hIfH2/zZyQiEhwUYEO0tbXx73//mxUrVrBy5Uq++OJzqqtriIqKZGxOGuNHRVAwspuicd3ERvURG9lLdEQfDty4sQKssyeCurZqdrQUU/p5JO++4WJrTSv1jdZ1cLIyM8ifMIFJkyZTVFTEnDlzmDJlyq4AFBGRgxPSl1Pp6Ohg2bJlvPHGG3ywbBnF69cTFRXBEZMyObawiyNyG5g5wUX+KHAOc7F1Zyds3g5lVWBuS6KkKp7/lrooLt1Genoac+Z8nblz53LWWWeRlqarTItI8PH05VRCKsBcLhfr1q3j3XffZenS1/h4+cdERobx9aPTOaWogeMLO5g6HsLDvDj0HupbYNlqJ++szeTV5d3U1rdy/HHH8a2LLuKiiy4iPT3dd8OIX3G73Vx11VUUFRVx88032z2OyLApwIY4UIB1d3ezZs0aVq5cyYcffsi//vU+DQ2NTCscyelHdnLGjEaOnQwRfrKQ6nbDZya88Gkaz3zgpGpbM+eccy4LFy5k9uzZWmYMISUlJSxcuJCVK1fys5/9TAEmQcHTAeYnP7qHz+12U1lZyYoVK1ixYgUfL1/O2rVr6Xf1M2nCCI6d3M+DNzRx0nQ3qYnVdo+7Vw4HHDURjprYwC+vhA/XwP1vfMwpp7zAuHHjuf3227n44osJDw+a/7aQ197ezne+8x3KyspwOp3MnDmThx56iCVLljB//nzGjBlj94gifiugj8AcDod76tQiGhoaaGhopKurixFZqXxtShyzJjZxjNHOjHyIjbZ70uGpbYIHX0/ij8/2kZKazh13/JTLL7+ciIgIu0eTYXr88cd58sknefPNN+nv7+e6667j1ltvJT8/H4CrrrqKKVOm6AhM/IrL5aKuro7q6mqqq6tpamqivb2d9vZ22tradr1saWmhp6eH6OhoYmNjiQgP5+FHHtESIlgB9sCiaDISu8hIhvEjISfD7qm8p3UnPLA0kd/90010TDyLFv2QBQsWkJSUZPdocpg2b97MCSecwPjx4zn11FM555xzmDx58q7XK8DETh0dHZSXl1NaWkpxcTGrVq1i7Zq1bNlSQW9fHw6Hg9TEZBJj44kJjyTaGUGUO5xoVxhRfU6iXE7C3A76HW56nW4+ay5nR2ezAgxgZLrDXfOs3VN8WUMLpHk5T3Z2wt/eieWe52Opb+rkmmvmc/XVV1NUVKTnyQJQd3c3y5Yt4/333+exxx7j4Ycf5qyzzgIUYOI7fX19fPbZZ7z//vu88847bCjZwPYd2wFIT05lbGIWOV2xZDSGkdkXQyoxJBFN+CH0Yfwvn/Jfd7WeA/NHr34C5/wUXr4bzjzWex8nLgZuOLuD783t4JUV4fz+pee47777GDVqJGefPY958+Yxe/ZsYmJivDeEeMSDDz7IRx99xJNPPslpp53Gjh07+Pzzz3cFmIg3uVwu/vWvf/HQQw/xxuuv09XVReGIcRgdSZzfmEsWk8kgjujmcGi2e9qvUoB5iFkJ1/we/nQ9XP07+OhPYHj5+fewMDh3Vh/nzqqmtgmW/ruOF1a+yjl/exRXv5ujjj6a008/ndNOO42ZM2fq6MwPXXHFFSxbtoxJkyYRGxvLmDFjuPHGG+0eS4JcRUUFjz32GH955BFqa2s5Ln0i83dOwyCd6KrAiQUtIXpA60742kJYfCHMnwuPvAZ/eA5WLoHEON/P090DK0vgvbUpvPV5HJ+trSErK4vzzjuf888/n1mzZmkno0iIaWtr47nnnuOvj/6V5R8vxxiRy7HNmRzVmUkskT6ZwdNLiAqwYXK54Nyfwsg0eHDR7tuv+yNsb4QXfjb8Fo/hamiBV1dG8swnWbz/6Xbi4uM5++x5XHDBBZxyyilERUXZO6CIeM3nn3/OkiVLePrpp4mPiuH4sDyObEhhFIk+n8XTAaY2+mH6+RNWm8a9N3z59vtugLpm6/V2S0uCq77Rw+v/byt1L/Ty5++30V7zLhdfdCHp6Wl873vfo7i42O4xRcRDent7efzxxzlixgyOOuooit/+lLTuSI5rzmRew1hbwssbFGDD8Oon8PBr8Nz/g8g9TsmKjLBuf/g1eG2FHdPtXUIsXHhiH8/+qIr65zt54tYuyla/RlFREbNmHc9zzz1Hf3+/3WOKyGHo6uriwQcfZPy48Xz/+oXkVTj5resbXFc1kcWu4/iQClaxze4xPUZLiIfJrITZP4BXfg7HTNr3/VYUw7w7fLOpYzg2VsOSpek8+loHGZlZ/PjHP+Hyyy8nMtI3a+Micvh6enp4+OGHueuuu3D39HNGfz7HtWcTvcc+vY00cB+fchsnMIIEn8+p58CGsCvA9ty0UV4N978ET71vLSemJ8GlJ8EN58D4UfZv6jgULe3w4OuJ3PMMRETG8qPbbuO73/0u0dEBXmciEoRcLhfPPPMMt956K23NrZzTV0B+RyIfUMGnbKWNbhKI4hhGczLjyCSeD9jMW2zkDuYQg2/bfPQcmM1cLrj8VzBnmhVeb6yEY26AmCj45D7ofst6GRNl3f7GSlhwJpw4Da74tfX2/iwpHn70rVYqn2zlxxe38qtf3Elebi73338/XV1ddo8nIgM++ugjjjjiCK65+moSW+BXrSeQ3BHOr/iQCJz8hBN5hHP4CScSgZO7WcYatnMieRik8wif4SJwD2BAR2CH7K7H4K3/wL/+AFtrrZB65edw7OSv3ndFMZx9O3x6v1Vx9fUfwulHw0+v8O3Mw9HVA395M5ZfPBWJwxnNT++8k6uvvlpLiyI22bJlCz/84Q95+eWXmZI6lta6Rn7MHBrp4G6WcRPHks9Xrym4kQbuZQV3MIcUYvgNH1FEFvMo9NnsOgKz0Z6bNu5/CRbM3Xt4gXX7/G/CkpchKhKe/5n/beo4kOhIq/Vj82PN/OiiNn56+61MmJDPY489ps0eIj7U3d3NnXfeiWEYVHy6jkv6p1BZV8ONHEs4Tt5jEyeSu9fwAsgnjRPI5T02EUEYN/A1PmBzQG/q0BHYQTIr4eiFcOvF1nNcZVXw0KvwxcPW81z7Ul4NM66Fa8+ECTnWc2S//Yf1fJg/b+rYl52d8L+vJPDrp11kZo3gF7/4Jeeffz5Ou092EwliK1eu5Morr6S5toFL2iaS0RfNz1nGNykgnih20M4yNvMzTiKT+H2+n1rauZP3mUMeWcTTTjevU8YdzPHJpg4dgQ1R2xzG4tW++RTuexHioq3lw5UlkJYIO7tgbPb+325MFnR0WfdfWWK9fWy09f4CUVwM/OiiNiqe2MlFx9fznauuYPr06SxdupRA/mVIxB91dHSwaNEijj/+ePIao7mz6Vim9WXyHpuIIpx11LKJRuKJpJs+0ojd7/tLJZYe+oknkk00so5aogjjXcq9/rlc9f4LjJ1V49H3qT6hg7TkJuvPUH98HrZs3/8RWOUO64jtR5d6dz5fS46Hu69o5saz4VfP1nD++ecxefIUbr/9dubNmzfsI7Ly8nLu+eOfeOLJp2hvbSY+MZnLvn0pixf9gPHjx3vosxDxX2+//TYL5i/A1dHD/7hnUVC3e2nwcqZzOdO/fH820kDHfo/AGukgnkjmsu8r2QeSgF5CDA8Pd9/0e/vm/+glF3NGwG+/u+/73PIwfLAdZs8L6IPdA2pvcbPh0zQ++7CdESNGcdttt3HppZceViP+G2+8wQUXX0Lv6Bn0jp4OsUnQ0ULE1lVEbP2C5/7xNGeccYYXPgsR+9XV1XHjjTfy/PPPcfMlCdxxSRMxB9H2duOfHFQszee8/qJ93uf5sDWMO7Oce2+y5+fmeT+FFz5yawnRHxTNcvDQUmu34d6sKIaHl0LR8cHfAh+f5ODI0xq56vZu8qbWc8v/LCY7O5uFCxeyatWqg34/5eXlXHDxJXTMvJC+wpNwxKficIbhiE+lr/AkOmZeyAUXX0J5ufeXPER8qaenh3vvvZeCggLK13/A53/u45dXHVx4Adx0oZvlEZvYSMNeX7+RBj6O2MyNFwTuQcueFGDDkJzu4ORLHJz+Y+tIq7waevusl7c8DKf/GE6+xEFyevAH2KDoGAdTT2jjittamXNBJ8s/fYUjjjiCadOn8ec//5mWlpb9vv09f/wTvaNn4Ejb+w4XR9oYenOm88c/3eeN8UV8zu128+yzz2IYBfzqF3fy++928+kftjEl79Dez/hR8OSd/TwQvZznw9ZQSzt9uKilnefD1vBA9HKevLN/v095BBotIXpAc72btR+7Kf0C2togIQEKZlhHXqEUXvvS1uSmfHUSxSudtDR3cP7557NgwQJmz579lefKElPSaDvmChzxqft8f+72RhJXPk5LY723Rxfxqk8++YQf/OAHFK9by/9cGsXN57UQN8zr0JZXw33POXjqHQdNnS5SYpxceqqbGy9w2x5enl5CVICJz7hdbio3OtiyegSr/7ODrKwsrrzyKq644goKCgoAcIaF4T7nDhzOsP28n36cL/2c/v4+X40u4lFlZWXcfPPNLF26lO+clcbdl9WSve/f2YKGngOTgOVwOhhbACdcuI0F/6+fKbObeP6lv2MYBjNmzODee+8lNj4ROva/zEhHC/FJyb4ZWsSD2traWLRoEZMmFdLT+BlrHunnkRtDI7y8Qdvoxevcbjc9XRAVs/sXr+hYB5OO7mLS0TUcV++gYl0Fv7vnF3S0teBY/nfcRafBCGOvR2IRlV9w+be/7ctPQWRY3G43L7zwAt+/4QYSYvp56zduTprh2XOiQpECTLxueyU8fa+brNxIpszsYcI0iI3fHWbJ6Q6mz2ll2oluyte5ee2xVtz/fh4ionHnzYRxR+KIsS7A526oJKJqFYt+8Khdn47IISkrK2PhwoV89NGH3HFlFDef1/qV6wfK4VGAiddlj4G0UZHs6Cug8T89LHulnIzRYVaYTYWYOCvMHA4H+UUO5l3t4pXHHLhSMnFXr4cNH+DOHE9YdCyRjRU8989/6GRm8XttbW3cfffd/OlPf+TUYzMp+Ws3udnddo8VVBRg4nUOh4OT5vXy4mOb6D11MbhcbN9eRsOnn/OvlzaTNdYKs/wia2kxr9DBlYv7+e+HFWz6Atrcbhz1FbjdLrLHjuGLL74gMzOTGTNmqINR/I7L5eLJJ5/klptvJiHWzas/d3PaUVou9AbtQhSf+ecD4dTEnQz5x+66zd3XDdtKiaj5HNeOCkbkhTN5IMyioh3cM81FZ184I/7vf2hp6+e8iFep3tTGf9ZUk5KSzMknn8LMmTMpLCyksLCQvLw8wsL2vYNRxJs++eQTbrjhBjaWmdx5ZRjfP7tNy4VDeHoXoo7AxGe+Pq+XfzzwPn25R+AIt+oFHOFRMLqIvtFFuHu7qNpmsuPDz3n3uUpGjg9j5nk9nDe7j/PHFfO3DTMpT5vDqhseoKkN3vuimXfXfMiL/1jBbyqaaWxqIzo6msmTJnHkUUdx5JFHcvzxxzNx4kQcDp2PJ96zadMmbr31Vl588UWuPjuNN2/vIDPF7qmCnwJMfCZzlIOxBS42ly3HXXjyV17viIiGMdPoGzMNd08nW2s2cMVvXiEirJ8rjNX8bcNM1jRks6Y+i6npO7jgBBcXnLB919vXt8C6zV18samU/5Tv4LfvvEzZpu2MGjWS0047nblz53LSSSeRnKwt+OIZX3zxBb/85S958cUXOfGoLL54qJ+icbV2jxUyFGDiUyfO7WPLPZ/QN+4YHFFx+7yfIzIGtzOMvJFOzp3dj9O5hTHxzVS2J/N46TR+l/72V94mPQnmTIc509uBdgC2N8Lb/61j6edvMf+aZ2ht6+SIGTP45ty5nHHGGRx99NE6OpNDUl9fz+uvv87f//53li1bxllzRvDJff0cXajnuXxNz4CLTyWmQvbofsLKPtjv/dy9XYQXv8aR5/YRHgZOh5tvF6wG4OmyqfS7Di50slPhilN7+eet1dQ/187KJf2cO7OU95b+lVmzjmf8uHHcfffdbNmyZdifmwSn/v5+PvvsM37zm99w7DHHkJWVxa23/ICJqeso+ZuLl26v5uhCu6cMTQow8amG7VBV5qI/o2C/9wtb/zb5k/oZmevYddHSyw0rwGp2JvJ+9bhD/thOJ8wsgNsuauWj325l27P9LJpXy/NP/5nc3Fy+/vU5PP/88/T29h76JyZBw+12s2HDBpYsWcLZZ59NWloaX/va13j2yfs5o8jkswdd1DzdxJLraykYbfe0oU1LiOJTCckQFhlGf3QC7q52HNFfvfieu2kbzupVfP3KfsA60lq82sk90+o5KrOK/9Tm8IQ5jVNHD++SKulJ8P15HXx/XgfFm+Ght9dzzdVXEhUdy4IF32X+/Pnk5uYO62OIZ3V2drJ9+3a6u7tJTU0lNTWV8PDh/RhrbW1lzZo1rFq1ihUrVvDee++xY8cOJuZn842ZLhbc0soJU90kxVd56LMQT9E2evG5F/8aTsW6HsJyJ9E/8+Ivvc7tdhHxwQOceFIdU4/96jJhruMoblo+l7jwbrZ/57fERXj2aKmjC575MIr7l6bx37U1nHDCbBYs+C7nnXcesbH7v1y7eEZbWxumaVJSUkJxcTFr1qyhtNSkdkctbe3tX7l/YkICCYkJxMfHExcXR3x8/JfOD3S73bhcrl1/ent76e7qpqu7i/b2drZt205kRDgT87M4ssDFKUW1nDS9nyz1E3qcttFLwDvzsl7qtzl47v/28hvt5v+SFNNM0df2/rbre4oJd57Ozr4oXtxUyGXGGo/OFhsNV32jm6u+UUPxZvjru+tZdNNCvve96zj77HlceumlnHrqqURGRnr044aa3t5eKioqKC0t3fWnuLiY0lKTHTtqcTqdjM1JZ1JeNNNzd3LR1GZGpfeTnWo9rxkVAY2t0NAKDa1ttHa00dYZRntXJO3dkez5e3lYmBsnbhwOiAzvJzqij+iIXmKj3BijwRjdR0R4tT0Phhw2BZj4XESkg6wcN/0723H3dlnb5wF3907CS97i9IW9OJx7/yUtNqKD00dv5LUtBo+b0z0eYENNzoN7FjTw6+/AW585eeqjZVx80YuEhUdw8smnUFRUxMSJE5k4cSJRUVG0t7fT3t5OW1vbl/64XC5SU1NJS0sjNTWVvLw8xowZE7S7H7u6umhtbaW1tZW6ujqqq6uprq6mqqqKkpISTHMDFRVb6OvrIyU5gfFjkpk42sVphTu56eRmJo6G/FEuoiL3vx09IRbGZg+9pR/oHPgjoUABJrZwhjlIyY6goakGMq0NGWHFbzBppovMUfv/wX6ZsYrXthi8Vz2Omp0JjIxr8+qsEeFw5jEuzjymhs4b4c3/dPNhyXKWv/MZjz6yk6oa6xLuDoeD2Jgo4mKjiY+LIj42gvhYJ04HNLX109TSRWNzG11dPcTFxTJhwgSmTp3GhAkTyM3NJTc3l3HjxjFixAi/DLfe3l6qqqqorKxk69atVFZWWn+vrKRmWw07duygob6BniGbYMLDw8jKSCY7PZacTAdTRnVwwfnNGDl9GKMhNbEN8O7/nwQvBZjYJie3j4bGKsgch7u+kvDaEmYv2L1xY1/OyjVJiuykpSeGp0qncvOMj30zMBATBefOcnHurN1HB53d4HZDTJQbh6ML6Nrn27vd1rlpxRUdFG8to3hrLe8tXUbltm6qahro6e0jMSGBAqOAoqKpjB07lpSUFFJSUkhOTqavr4/Ozk46Ozvp7+8nJSVl12aG7OxssrKyDrsfsq+vjx07dlBdXU15eTnl5eWUlZVRVlZKZWUl27Ztx+VyERUVyagRKeRkRpGb1ceMrG6+aXQwIrmTrFRIioPEWEiMg8TYfpzOBqDhsGYS2R8FmNhmVK6LkopN9LiOJ2L1C5x8bh9R0Qc+8ogJ7+OC8cU8WnIkT5RO82mA7XWeqIO/r8MBI9KsP6fM7AA6dr3O5YKttVBS2UZJ1QaKq2pYXhpJU6uL5rZeWts6CQt3EhMVQXRUOGFhDlraumlq3kn7TmvZLDIiglGjRjJmbC75+flMmDCBcePGkZOTQ2dnJy0tLbuW9iorK9m0qZzKLZVs276NhoZG3G43TqeTnBGp5OXEYozq4ZwZHeSd2kpuNozNgozkHhyOHR5+FEUOnQJMbJM9BlyN23CUf0paSjvGjIN/28uN1TxaciRrGrJZXZ/FtPTA/4HqdFrP6YzNhtOP3gnsPOi37e2zjuwqtveypXYrFXUtlO8o49XV77C5qp0ddc3EREeREB9DQnwUqYkR5Ga7mJrZxdyCdnLSuhiZDiNSITPFRXhYvfc+UREPUYCJbZLSwOHqxbnhPU7/Yf8hPe8za0QlYxOa2NKWwuPmdKalv+XFSf1fRDiMzrT+zMYFNA/8sbjd4HBog4MEFzVxiG0cDgcjch3MmOUmNfPQNi1Y1VLWDsSny4oOuloqVPnhnhCRYVOAia3mXdXH7Lmuw3rbywtWAbCtI5H3qg69WkpEApsCTGwVEek47C3jRkoDR2daJ0M/UTrNk2OJSABQgEnAGSz3BeucMIAXNk2ivVftGCKhRAEmAWkwxC7OX0e4s5+OvkheKNc1LURCiQJMAtbi1U7SYzo4Y0wZAE+UTrd5IhHxJQWYBLTFq51cNnChy/eq8qhuT7B5IhHxFZ+cB2YYRg5wF9CE1RO0DcgFkoBFQCRwD9AIFJumucQXc0lwWN5Wtrtaqmwqt9jczCEivuGrI7CJwEnAGKAeOME0zYXAo8AC4FrgPtM0rwfmGoYR4aO5JAiEO61qKYAnTO1GFAkVvgqwrVgB9i3g61ghBlAFjASyB+4D1lFako/mkiBxhWEtI65ttKqlRCT4+SrAbgCSTNN0Ay3A2IHbc4AaoHLg7wCpDO3AETkIx4+oJDehCYDHTG3mEAkFvgqwvwF3GYaxBFgJvGwYxgPAd4ElwF+A7xuG8RDwgmmafT6aS4KEVS1lHYU9XTaVPpf2J4kEO59s4jBN83PgrP3cpQW41BezSHBavNrJtcZqfvHfOWzvSOC9qnGcNmaj3WOJiBfp11QJGg9tadpdLaXNHCJBTwEmQSUu3lpGfHFzoaqlRIKcAkyCipFejNOhaimRUKAAk6ASE9FBboqqpURCgQJMgk5hhnWhy/eq8qhqT7R5GhHxFgWYBJ1xqaVEhXXhxslTZVPtHkdEvEQBJkEn3NlHQfo6wNqN6HbbPJCIeIUCTIJSYaa1jLiuMYvVDdk2TyMi3qAAk6A0MqGSvMRGAB7XOWEiQUkBJkHJ4YBvF1hHYU+VqlpKJBjpu1qC1uUFqwDY0ZnAu1XjbJ5GRDxNASZBa0JyI1/Lsq7So2opkeCjAJOgdlnB7mqpth5VS4kEEwWYBK3Fq51clL+OcGc/nX2RvLBpkt0jiYgHKcAkqP2qtIu5Y0sBeKJUy4giwUQBJkGvP9Lajfi+qqVEgooCTIJeXmopUWGduHHyZKmqpUSChQJMgl64s5+C9GLAWkZUtZRIcFCASUgYrJYqbsxiVb2qpUSCgQJMQsLIhEoSo5oAeFzXCRMJCgowCQkOx+7rhKlaSiQ46LtYQkZhpnVSc21nPO9sHW/zNCIyXAowCRkpMY0cM1gtpXPCRAKeAkxCymWGVfD70uaJtPZE2TyNiAyHAkxCykX564hw9tHZF8nz5aqWEglkCjAJKWnRnXxzbBmgZUSRQKcAk5AzeJ2wZdW5bG1TtZRIoFKASUhZvNrJ3NxSUqI6rGqpMh2FiQQqBZiEnB+vc3Ph+IFqKVPVUiKBSgEmIak5zDqpeX1TJl/Uj7B5GhE5HAowCUkjE7aSFN0IwOOmlhFFApECTELS0Gqpp8tULSUSiPRdKyFrMMBqO+N5W9VSIgFHASYhKzmmkWOzKgFrM4eIBBYFmIS0ywyr4PelzYWqlhIJMAowCWmD1VJd/RGqlhIJMAowCWmp0Z2cObYUULWUSKBRgElIW7zauWsZcVl1LpVtSTZPJCIHSwEmIe+bY4dUS5VOtXscETlICjAJeVFh/Xwrfx0Aj5vTVS0lEiAUYCLA5QXWMuKG5gz+WzfS5mlE5GAowCTkLV7t5NjsrYxPbABULSUSKBRgIsDNa3Zv5vjHxiJ6+/WtIeLv9F0qMqCyfy0AdZ3xvL013+ZpRORAFGAiA5JjmhiRMFAtpXPCRPyeAkxkiMGC35c3T6SlW9VSIv5MASYyhJFeTORgtdQmVUuJ+DMFmMgQ0RGdnJlrAvCEOd3maURkfxRgInu4bOC2zioyAAAYQklEQVScsGU1eaqWEvFjCjCRPSxrKSc1qgPQZg4Rf6YAE9lDmLOfnJRiwLrQpaqlRPyTAkxkLyZlqlpKxN8pwET2Iju+iuRoq1rqMW3mEPFLCjCRvXA4dp8T9ugGVUuJ+CN9V4rsw8RMK8A6e+N4S9VSIn5HASayD8nRTYwcrJZSQ72I31GAiexH4cBmjpcrVC0l4m8UYCL7UZBeTJijj+7+CJ4rn2z3OCIyhAJMZD+iw7vISy0FdFKziL9RgIkcwOBuxA9q8tiiaikRv6EAEzmAvJSy3dVS2swh4jcUYCIHEObs56IJ1tWanyhVtZSIv1CAiRyEywca6s3mDD6rVbWUiD9QgIkchGe21zAhqR6Ax0tVLSXiDxRgIgfB4YD0RGsZ8R9lqpYS8Qf6LhQ5SIPVUvVdcby5dYLN04iIAkzkICVHNzEycQug3Ygi/iDcFx/EMIxc4A6gFmgDkoEoIAm4FjCA24AO4DXTNF/0xVwih6owYw01rWN5YdNEmrujSY7qsnskkZDlqyOwxUAVMALYASSZpnkT8AFwPnAzsNg0zWuA7/loJpFDNlgt1e8O57nySXaPIxLSfBVg+cBLwALgMqwwY+DlSCDDNM2agdt0lo34rejwLsalmgA8od2IIrbyVYBtB1pN0+wd+PeogZc5QA1QbRjGiIHbHD6aSeSwDFZLfViTS0Vrss3TiIQuXwXYb4FfGYZxP/AXoNEwjHuBE4FngT8A9xiG8X/AfT6aSeSw5KZsJDp8oFpKBb8itvHJJg7TNEuAb+3nLuuBS30xi8hwhTn7MdLXsXr70TxROo2fzPwAh9YNRHxO2+hFDsOkgQtdljan85/aUQe4t4h4gwJM5DBkxVdTkDxYLaVlRBE7KMBEDoPDAZcNFPz+o6yInv4wmycSCT0KMJHDNBhgDV1xvFmZb/M0IqFHASZymP53cyuzR1QA2o0oYgcFmMgwRMVa54S9WmHQ3B1t8zQioUUBJjIME9LXE+boo7s/gmfLJ9s9jkhIUYCJDMPQaqnHTVVLifiSAkxkmAoHzglbvm0sm1UtJeIzCjCRYcpN3khM+E5AmzlEfEkBJjJMYU4XRsY6wLrQpVvXUxDxCQWYiAcMNtSXtaSzckeOzdOIhAYFmIgHZMVXkxKjaikRX1KAiXiAw7H7KOyfqpYS8QkFmIiHTBwIsMbuWN6onGDzNCLBTwEm4iFJ0c2cMLICsDZziIh3KcBEPOiyglWAVS3V1KVqKRFvUoCJeNCF44uJCuulxxWuaikRL1OAiXjQXRt6mZe7AYAnSlUtJeJNCjARD+uOXAtY1VKbWlJsnkYkeCnARDxsbPJGYiJULSXibQowEQ8Lc7ow0geqpUpVLSXiLQowES8ozLAa6je2pKlaSsRLFGAiXpAVX0NqTB2gaikRb1GAiXiBqqVEvE8BJuIlQ6ulXt+iaikRT1OAiXhJYnQLOYkVgHYjiniDAkzEiwozrc0cr6laSsTjFGAiXjQhbT1hTqta6pnyKXaPIxJUFGAiXhQV3s34VBOAuz5XtZSIJynARLxs8JywbW1juObTNJunEQkeCjARLxubXL6rWqqkbqrN04gEj30GmGEYKwzDOGYvt7/i3ZFEgotVLWUV/JbUTeOHq/R7o4gn7O87KR34u2EY1+xxe4IX5xEJSpMGdiO2dKWyrU3VUiKesL8A2wrMBq40DOMBwzAG76tqUpFDlBm3bVe1VEmdzgkT8YT9rmWYplkLnDTwz2WGYWR4fySR4ONw7D4nrLR+sqqlRDxgfwHmADBNs880zeuBvwEfA1r/EDkMEzOs58G6+lQtJeIJ+wuwO4f+wzTNvwHfBiq8OZBIsEqMaiEncTOgaikRTwjf1ytM0/xwL7f9B/iGVycSCWKFmWuoas3jtQqDxq4YUqM77R5JJGBpP6+ID32pWmrjZLvHEQloCjARHxpaLfVEqaqlRIZDASbiY5MGqqU+2T6G8pYUm6cRCVwKMBEfG5tSTmxEOwCPmzoKEzlcCjARH3M6XBjp6wC4d9003KoGEDksCjARGxQOqZa6dPkYm6cRCUwKMBEbqFpKZPgUYCI2GFotZdZPoVvVUiKHTAEmYpPCjDWAm+6+GF7fUmD3OCIBRwEmYpOEqFZGJ1UA8ISpZUSRQ6UAE7FR4cA5Ya9tKaCxK8bmaUQCiwJMxEb5aSWEOXvpdYXzz41T7B5HJKAowERsFBXeTX7qBkDLiCKHSgEmYrPB3YgrdoxhY0uqzdOIBA4FmIjNxiZvGlItpaMwkYOlABOxmVUtZV2t+clSVUuJHCwFmIgfmDSwjLipNZVPtqtaSuRgKMBE/EBG3HbSYmsBLSOKHCwFmIgfcDh2nxP2d1PVUiIHQwEm4icmZqwF3HT3x3DhBxPtHkfE7ynARPyEVS21GYCSWi0jihyIAkzEj1gFv7C5aQINqpYS2S8FmIgfmZC2nnBnLy53mKqlRA5AASbiRyLDexi/q1pqus3TiPg3BZiInxk8J+zTHaMpa1a1lMi+KMBE/MyYodVSpToKE9kXBZiIn3E6XANb6uHJ0qmqlhLZBwWYiB8aPKl5c2sqH6taSmSvFGAifigjbjuTU3cAqpYS2ZdwX34wwzCeBF4FRgO5QBKwCIgE7gEagWLTNJf4ci4Rf+NwQErCGmg8lWfLp3DvrDeIDu+zeywRv+KzIzDDMH4ItA/88wTTNBcCjwILgGuB+0zTvB6YaxhGhK/mEvFXg9VSzd0xLN1SYPc4In7HJwFmGMZZQDOwYuBj1g68qgoYCWQDWwdua8I6MhMJaUOrpZ4o1TKiyJ58tYR4GVYwGQP/HjwSywFqsEItByvEUrHCTiTkFWauZmvLOF6tKKC+M5b0mA67RxLxGz4JMNM0LwIwDOMqoAvIMgzjASAFuA6IAf4w8PoXTNPUYr8IMCG1hPedZ9LniuDiD4p49/SVdo8k4jd8uonDNM3/28erWoBLfTiKSECIDO8hP7WEDfVTWV87DVCAiQzSNnoRP1eYaTXUb2/PobQ5zeZpRPyHAkzEz1nVUm2ANnOIDKUAE/FzqpYS2TsFmEgA2HWhS1VLieyiABMJABlx20mLtaqlrv9khs3TiPgHBZhIAHA4dh+FldVP5qYvIm2eSMR+CjCRAGEFmJvu/mg2NRaweLW+fSW06TtAJEDER7UxJmkTACV1UwEUYhLS9NUvEkAGzwmraJpAR2+szdOI2EsBJhJA8tNKCHf24HKHUVo/xe5xxA+VV8OiJZB1PoSdYr1ctMS6PdgowEQCSGRYD/lpGwAoqZ1q8zTib95YCcfcADFR8Ml90P2W9TImyrr9jSBrIlOAiQSYSRmrAataqrFD1VJiKa+GK34Nr/wcfjkfxo+C8DDr5S/nW7df8evgOhJTgIkEmNHJm4gbqJbaUKejMLHc/xIsmAvHTt7764+dDPO/CUte9u1c3qQAEwkwTod7V7VUSd1UXG6HzROJP3jqfbjmjP3fZ/434an3fDOPL/j0cioi4hmFmav5b81xtHancMlHY8lJqrR7JLFZfbOLsdn7v8+YLKhrse/0i7UtLo++Px2BiQSgjLgdpMduB6CkTg31AgkJsGX7/u9TucO6X7BQgIkEqMFzwsrqJ9Pn0mJKqCuYAQ8t3f99/rzUul+wUICJBKiJ6WsZWi0loa1oloOHlsKK4r2/fkUxPLwUio4PnudMFWAiAWpotdT6Wi0jhrrkdAcnX+Lg9B/DLQ9b2+V7+6yXtzwMp/8YTr7EQXJ68ASY1h1EAlhh5moqW8azpTmfjt5YYiM67B5JbJRX6OCCG+GDj908cgO0tVnPeRXMgAtuDK7wAgWYSEDLT9vA++U99LoiMeumMGPkv+0eSWyWnO5g9jwHs+fZPYn3aQlRJIBFhvUwPq0E0G5ECT0KMJEAN2ngQpc72kepWkpCigJMJMANrZbSUZiEEgWYSICzqqWso7ANdVNxq1pKQoQCTCQIDJ7U3NqdTHXrGJunEfENBZhIEPhytZQa6iU0KMBEgsSXqqX6dYaMBD8FmEiQmJi+FgcuuvujKW8y7B5HxOsUYCJBIj6qjTHJVrVUSa2WESX4KcBEgkjhwG7ELc35dPTE2TyNiHcpwESCyPi0DUQ4e3C5wzDrp9g9johXKcBEgkhkWA/5u6qltIwowU0BJhJkCjNXA4PVUuk2TyPiPQowkSAzOmkzcZGtgI7CJLgpwESCjFUttRawAkzVUhKsFGAiQWhShrWM2KZqKQliCjCRIJQeV0tGnFUttb5WDfUSnBRgIkGqcOAorKxB1VISnBRgIkHKyLCqpXr6oyhvVLWUBB8FmEiQio9s310tpQtdShBSgIkEscFlxIomVUtJ8FGAiQSx/IFqKTdOVUtJ0FGAiQSxiLBe8tPWA9qNKMFHASYS5AYvdFm7cyQNqpaSIKIAEwlyX66W0lGYBA8FmEiQG1ottaFW1VISPBRgIiFgV7VUTxJVLWNtnkbEMxRgIiHAqpbaBmgZUYKHAkwkRBRmWJs5yhom0dsfYfM0IsOnABMJEUOrpTapWkqCgAJMJER8uVpKF7qUwKcAEwkhQ6uldqpaSgKcAkwkhFjVUt2qlpKgoAATCSERYb1MSC8BoETVUhLgFGAiIWZwGdGqlsqweRqRw6cAEwkxOUkVxA9WS9VqM4cELgWYSIixqqWsc8I21KlaSgKXAkwkBA2e1KxqKQlkCjCREJQeV0umqqUkwCnARELU4GYOVUtJoFKAiYQoI2PdrmqpclVLSQBSgImEqLjIdsamlANaRpTApAATCWGDy4hbmsarWkoCjgJMJISNTzWJDBusliqyexyRQ6IAEwlhEWG95KetB1QtJYFHASYS4gbPCavdOYJ6VUtJAAn3xQcxDOM44CagHdgCxANRQBJwLWAAtwEdwGumab7oi7lEBEYnVRAf2UJ7TxIltdOYnfuu3SOJHBRfHYGlAPNN07wGmAUkmaZ5E/ABcD5wM7B44PXf89FMIgI4HG4mZqwFrGopl6qlJED4JMBM01wKtBuG8RPgY6Bq4FVVwEggwzTNmoHb3L6YSUR2K8y0diO29yRS1ZJr7zAiB8knAWYYRgLwF+BT4K/AqIFX5QA1QLVhGCMGbtOvfyI+lh5bR2ac9TukzgmTQOGT58CAe4EJwHeAfqzAuhfrObAFwCrgHsMweoD7fDSTiAxRmLmG2s0jKasv5KRxS4kI67V7JJH98kmAmaZ59QHush641BeziMjeGelr+XDzN+h1RVHeOHHX82Ii/krb6EUEgLjInburpXShSwkACjAR2WVXtVTzeHb2xNs8jcj+KcBEZJeh1VIb6lQtJf5NASYiu3ypWqpOy4ji3xRgIvIlkwaWEet2jqB+Z6bN04jsmwJMRL4kJ2kLCZEtgI7CxL8pwETkS6xqKavgV9VS4s8UYCLyFYWZVoC19ySytSXP5mlE9k4BJiJfkTa0WkrnhImfUoCJyF4NFvxubCiktz/C5mlEvkoBJiJ7NTF9HQ5cu6qlRPyNAkxE9io2cie5KRsBWF+rhnrxPwowEdmnwoHdiJXN42hXtZT4GQWYiOzT+NQNu6qlTFVLiZ9RgInIPoWH9TEhrRjQhS7F/yjARGS/Bs8Jq9uZrWop8SsKMBHZr5zELSRENQOwXkdh4kcUYCKyX1a1lHV15g11RaqWEr+hABORAxq80OVOVUuJH1GAicgBpcXWkxVfDahaSvyHAkxEDsrgOWEbGybR0x9p8zQiCjAROUjGrmqpSMobVC0l9lOAichBGVotpQtdij9QgInIQRtsqFe1lPgDBZiIHLTxKSaRYV2qlhK/oAATkYMWHtZHQdp6QNVSYj8FmIgcksFlxLqd2dTtzLJ5GgllCjAROSSjEit3VUvpnDCxkwJMRA6Jw+HedU7YhrqpqpYS2yjAROSQDQbYzt4EVUuJbRRgInLIUr9ULaXNHGIPBZiIHJbBgt+NDYWqlhJbKMBE5LAYGetwOvrpdUWyUdVSYgMFmIgcltiIDsYmD1ZLaRlRfE8BJiKHrTDT2sxR2TyO9u4Em6eRUKMAE5HDNlgtBQ421KtaSnxLASYih+1L1VI6qVl8TAEmIsMyWC1V36FqKfEtBZiIDMuoxEoSVS0lNlCAiciwOBxuJg5WS9UXqVpKfEYBJiLDNnhS886eRFVLic8owERk2FJjG3ZVS61XtZT4iAJMRDxiUqaqpcS3FGAi4hFGulUt1eeKZGNDod3jSAhQgImIR8REdJA7WC2l3YjiAwowEfGYwXPCKltULSXepwATEY8Zl1pK1GC1VJ2qpcS7FGAi4jHhzj4mpBcDaqgX71OAiYhHDZ4TVt+RpWop8SoFmIh41KjErSRGNQFQonPCxIsUYCLiUQ6Hm8KBaqmSuiJcbv2YEe/QV5aIeNzghS47ehOobFa1lHiHAkxEPC4lpoHs+CpAmznEexRgIuIVhbuqpSaqWkq8QgEmIl5hpBerWkq8SgEmIl4RE9FBbkoZoIZ68Q4FmIh4zeBuxK0tebR1J9o8jQQbBZiIeI2qpcSbFGAi4jV7Vku53TYPJEFFASYiXjV4ocuGjkzqdmbbPI0EEwWYiHjVyITK3dVSdbpOmHiOAkxEvMrh2L2ZY4OqpcSD9JUkIl43eFKzqqXEkxRgIuJ1KTGNu6uldE6YeIgCTER8Yle1VGMhPX2qlpLhU4CJiE/srpaKoKxR1VIyfOF2DyAioSEmooMjR31MfGQb41JK7R5HgoDfBJhhGKOAe4BGoNg0zSU2jyQiHnb82PftHkGCiMPtJ6fGG4ZxF/CmaZqfGIbxOjDPNM3e/b2Nw+Fwx6leTUQkIHR1QF+v2+Gp9+c3R2BANrB14O9NQBJQv783KCgo8PZMIiLiQYZhuE3T9EiI+VOAVQI5WCGWCjQf6A089SCIiEjg8aclxGzgD0Ab8Jlpmo/YPJKIiPgxvwkwERGRQ6HzwEREJCApwEREJCApwEREJCApwEREJCApwEREJCD503lg+7WvqinDME4BrgAcwIOmaX5i35SByTCM44CbgHZgC3AU0DDw6ruBPPbxGBuGcQuQi3Xi+SLTNOt8N3ngMQxjDPAKsArYBoQBUViP37WAAdwGdACvmab54pC3vRw4AYgFfmaapgoFD8AwjIVYX88RwCxgDfra9jjDMPKBZ03TnLHn4wZEssfP7n09toZhOIGHsU6nijJN8/r9fdxAOgK7Frhv4BOaaxhGxMDtPwTmA98FbrVruACXAsw3TfMaYA4wHugCaoEK9vEYG4YRDZxgmuZC4FFggU+nDkwnAjsG/r4dSDJN8ybgA+B84GZg8cD/xff2eNsrTdNcANw+cD85ANM0l5imeRVQBVyIvrY9buAc3vnAzn08bnv+7E7Yy30GzQE2maa5CKgb+OV6nwIpwPZWNQXgME2zxzTNTiDalskCnGmaS4F2wzB+Avwd+I5pmtdhtaNcxL4f41SsHwRg/YAY6cOxA9W/gauAq4F5WI8b7H78MkzTrBm4bc+TNHv3uK8cBMMwJmL9vFiHvrY9zjTN7aZp/ghrBWdvj9ueP7tT9nKfQUPve8DHPZACbLBqCr5cNdVlGEakYRgxWL9ZySEa+I3oL8CnWEcCg9d8b8A6/N/XY1wLpA38PQeoQQ5kBhBpmqYLa9lq3MDtg49ftWEYIwZu27MqzbXHfeXgXA/8DshEX9vetrfHbc+f3TV7uc+gofc94OMeME0ce1ZNAVOBxcCxwDVYa9z3mKb5mW1DBijDMP4KTMB6/ssB9GF9gycCN2I9h7DrMQa+AP7XNM3rDcO4Cet5mxTgOtM0W3z/GQQOwzBmYi1V1QFrgTFAHNYRwgKs/4fbgR7gOdM0XzMM4yHTNK81DOMi4FSs/5dbTNPcYsfnEGgMw3jdNM1vDjzt8Bf0te0VhmG8aZrm6Xs+bkAMe9QE7uU++cBM4BHgzwz8MjGwvL5PARNgIiIiQwXSEqKIiMguCjAREQlICjAREQlICjAREQlIAdPEIRLsDMM4Ffg9cIxpmp2GYYwE3gJON02z2t7pRPyPdiGK+BHDMH4LxGNVe70P3GWa5jv2TiXinxRgIn5k4Fyl5VjnKq00TfNnNo8k4rf0HJiIHzFNsxfrZM5Tgb/ZPI6IX9MRmIgfMQxjLPA2VhvBucDXTdPst3cqEf+kIzARP2EYRiTwDNblJf6I1Qt3p71TifgvBZiI/7gHWG6a5usD/74euMQwjDn2jSTiv7SEKCIiAUlHYCIiEpAUYCIiEpAUYCIiEpAUYCIiEpAUYCIiEpAUYCIiEpAUYCIiEpAUYCIiEpD+P/GU1y7PmrvzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# execute this cell a couple of times and see how location of fault in section changes\n",
    "\n",
    "mask_surfpoints = geo_model.surface_points.df.surface.isin(['fault']) # perturb only fault\n",
    "indexes = geo_model.surface_points.df[mask_surfpoints].index\n",
    "\n",
    "rand_val = np.random.uniform(-100, 100)\n",
    "we = geo_model.surface_points.df['X'].values[mask_surfpoints] + rand_val\n",
    "ns = geo_model.surface_points.df['Y'].values[mask_surfpoints] + rand_val\n",
    "d = geo_model.surface_points.df['Z'][mask_surfpoints] + rand_val\n",
    "geo_model.modify_surface_points(indexes, X=we, Y=ns, Z=d)\n",
    "\n",
    "geo_model.update_to_interpolator()\n",
    "\n",
    "_=gp.compute_model(geo_model, compute_mesh=False)\n",
    "\n",
    "_=gp.plot.plot_section_by_name(geo_model, 's1', show_data=True, radius=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full simulation in loop with probability calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm_notebook as tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Init probability blocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unique lith ids\n",
    "lith_id = np.unique(np.round(geo_model.solutions.sections[0]).astype(int))\n",
    "\n",
    "# setup solution arrays\n",
    "geomap  = np.round(geo_model.solutions.geological_map[0]).astype(int).ravel()\n",
    "section1  = np.round(geo_model.solutions.sections[0]).astype(int).ravel()\n",
    "block = np.round(geo_model.solutions.lith_block).astype(int)\n",
    "\n",
    "#init counters\n",
    "count_map = np.zeros((len(lith_id), geomap.shape[0]))\n",
    "count_section = np.zeros((len(lith_id), section1.shape[0]))\n",
    "count_block = np.zeros((len(lith_id), block.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_surfpoints = geo_model.surface_points.df.surface.isin(['fault'])\n",
    "indexes = geo_model.surface_points.df[mask_surfpoints].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy initial dataframes to reset data before every iteration\n",
    "import copy\n",
    "\n",
    "west_east = copy.copy(geo_model.surface_points.df['X'])\n",
    "north_south = copy.copy(geo_model.surface_points.df['Y'])\n",
    "depth = copy.copy(geo_model.surface_points.df['Z'])\n",
    "\n",
    "surfindexes = list(geo_model.surface_points.df.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f11e9dd42f7c4d589c050f87c8917982",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=13), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_iter=13 #number of iterations\n",
    "\n",
    "for i in tqdm(range(n_iter)):\n",
    "    plt.cla()\n",
    "    plt.clf()\n",
    "    \n",
    "    # make virgin again\n",
    "    geo_model.modify_surface_points(surfindexes, X=west_east, Y=north_south, Z=depth)\n",
    "    geo_model.update_to_interpolator()\n",
    "    \n",
    "    # perturb data\n",
    "    rand_val = np.random.uniform(-100, 100)\n",
    "    we = west_east[mask_surfpoints] + rand_val\n",
    "    geo_model.modify_surface_points(indexes, X=we)\n",
    "\n",
    "    geo_model.update_to_interpolator()\n",
    "\n",
    "    _=gp.compute_model(geo_model, compute_mesh=False)\n",
    "    \n",
    "    \n",
    "    ##### calculate and update probability fields #####\n",
    "    geomap  = np.round(geo_model.solutions.geological_map[0]).astype(int)[0]\n",
    "    section1  = np.round(geo_model.solutions.sections[0]).astype(int)[0]\n",
    "    block = np.round(geo_model.solutions.lith_block).astype(int)\n",
    "    \n",
    "    \n",
    "    for i, l_id in enumerate(lith_id): #enumerate through all liths\n",
    "        count_map[i][geomap == l_id] += 1 #sum up frequency\n",
    "        count_section[i][section1 == l_id] += 1 \n",
    "        count_block[i][block == l_id] +=1 #block is raveled so no need for indexing\n",
    "\n",
    "#### finish probability calculation and save as numpy arrays ####\n",
    "prob_map = count_map/n_iter\n",
    "prob_section = count_section/n_iter\n",
    "prob_block = count_block/n_iter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gempy.plot import visualization_2d\n",
    "plot = visualization_2d.PlotSolution(geo_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import FixedFormatter, FixedLocator\n",
    "import matplotlib.colors as mcolors\n",
    "from gempy.plot import helpers\n",
    "\n",
    "def plot_probabilities(plot, block, shape, T = False, extent=None):\n",
    "    fig, ax = plt.subplots(1,3, figsize=(12,12))\n",
    "    ax = ax.flatten()\n",
    "    print(int(np.ceil(block.shape[0]/2)))\n",
    "    print(block.shape)\n",
    "    for i in range(block.shape[0]):\n",
    "        #allcolors = list(plot._color_lot.values())[plot.model.faults.n_faults:]\n",
    "        allcolors=['#9f0052', '#ffbe00', '#728f02']\n",
    "        formnames = list(plot._color_lot.keys())[plot.model.faults.n_faults:]\n",
    "        c2 = allcolors[i]\n",
    "        title = formnames[i]\n",
    "        cmap = mcolors.LinearSegmentedColormap.from_list('c',['#FFFFFF',c2])\n",
    "        if T:\n",
    "            plotblock = block[i].reshape(shape).T\n",
    "        else:\n",
    "            plotblock = block[i].reshape(shape)\n",
    "        im = ax[i].imshow(plotblock, origin='lower', cmap=cmap,\n",
    "                    extent = extent)\n",
    "        ax[i].set_title(title)\n",
    "        helpers.add_colorbar(im, label='probability')\n",
    "    fig.tight_layout()\n",
    "    \n",
    "def calculate_ie_masked(prob):\n",
    "    ie = np.zeros_like(prob[0])\n",
    "    for l in prob:\n",
    "        pm = np.ma.masked_equal(l, 0)  # mask where prob is 0\n",
    "        ie -= (pm * np.ma.log2(pm)).filled(0)\n",
    "    return ie\n",
    "\n",
    "def plot_section_by_name(self, section_name, image=None, show_data=True, show_faults=True, show_topo=True,\n",
    "                         show_all_data=False, cmap='viridis', label='cell entropy', title=None, **kwargs):\n",
    "    assert type(section_name) == str, 'section name must be a string of the name of the section'\n",
    "\n",
    "    j = np.where(self.model.grid.sections.names == section_name)[0][0]\n",
    "    l0, l1 = self.model.grid.sections.get_section_args(section_name)\n",
    "    shape = self.model.grid.sections.resolution[j]\n",
    "    if image is None:\n",
    "        image = self.model.solutions.sections[0][0][l0:l1].reshape(shape[0], shape[1]).T\n",
    "    else: image = image.reshape(shape[0], shape[1]).T\n",
    "    extent = [0, self.model.grid.sections.dist[j][0],\n",
    "              self.model.grid.regular_grid.extent[4], self.model.grid.regular_grid.extent[5]]\n",
    "\n",
    "\n",
    "    if show_data:\n",
    "        self.plot_section_data(section_name=section_name, show_all_data=show_all_data, **kwargs)\n",
    "\n",
    "    axes = plt.gca()\n",
    "    im=axes.imshow(image, origin='lower', zorder=-100,\n",
    "                cmap=cmap, extent=extent)\n",
    "    if show_faults:\n",
    "        self.extract_section_fault_lines(section_name, axes)\n",
    "\n",
    "    if show_topo:\n",
    "        if self.model.grid.topography is not None:\n",
    "            alpha = kwargs.get('alpha', 1)\n",
    "            xy = self.make_topography_overlay_4_sections(j)\n",
    "            axes.fill(xy[:, 0], xy[:, 1], 'k', zorder=10, alpha=alpha)\n",
    "\n",
    "    labels, axname = self._make_section_xylabels(section_name, len(axes.get_xticklabels()) - 2)\n",
    "    pos_list = np.linspace(0, self.model.grid.sections.dist[j], len(labels))\n",
    "    axes.xaxis.set_major_locator(FixedLocator(nbins=len(labels), locs=pos_list))\n",
    "    axes.xaxis.set_major_formatter(FixedFormatter((labels)))\n",
    "    if title is None:\n",
    "        title=self.model.grid.sections.names[j]\n",
    "    axes.set(title=title, xlabel=axname, ylabel='Z')\n",
    "    helpers.add_colorbar(im, label=label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "(3, 6400)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA10AAADnCAYAAADlwjzRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+wZGV95/F3X+YCw6AjM4L8mEQkwIOFm5+rtQZjKZutuKHE3Vj5scZkSYSAmIRSTKJR88NocE2BkSz+WNFoNom7G4uQktSaVCqJu4u6ibVqVpSHxI3yYyQzMMMgM4PD6N0/zjnTPff27XtOd5/Tz3nu+0V13b59T/d57oX+0M95vud7BisrK0iSJEmS2rG06AFIkiRJUs6cdEmSJElSi5x0SZIkSVKLnHRJkiRJUoucdEmSJElSi5x0SZIkSesIIZwRQjhp0eNQv21Z9ADUnhDCrwEfizF+asrnPxW4BdgPbANeEWM8Mr8RSsrJrJlTvsZ24APAb83yOpLyMY9smdF/AF4PPLig/SsDTrryd2kI4SeBs4E/AHYAj8YYPxxC+G/AlcANwDeAk4CfB/4SuBu4GXhrjPHvQgi/AzwDiAv4HST1xyyZ807g5cChRQxcUtJeG0J4ELgPuAg4CHwbRWb8GHABcBbwBuBc4N8CJwCfBv6hfPyTwMXAZ4DvAH4L+DpwDbAC3Av8F+DDwB3A91Bk1ncDPwm8vfXfUtmyvDB/dwL/FfgU8IPA7wE/FkK4GPgS8OPATooPOacCzwJOiDG+Isb4uXLC9YPAN2OMTrgkbWSWzPm7GOMvlttJ0qj3xhh/FngRRcZ8rHz8WRSTrMeA3wW+SrEq9TXgEeCF5XZ/G2P8FeCpwNuA9wIvKLd9DHgU+F6KBYl/iDG+Ffhs+fqfocgyaWpOuvL3K8CTgb8BBjHGwxQhciPwLor/Bv48xvg6ihB7gCKkAAghvAm4IMZ4XdcDl9RLM2WOJK3j4Mj9qyhOfbgHGAAfAf4Q+JfAT1FMnN5Csbr1l+VzHi2/Ho4xHgWOUuTRCcC7y0y6o3y82teRcpuVdn4lbSZOuvL3BEUI/TDFByEols4fiTHeB/w+cFkI4V3ADwF7qyeWJUI/ATw7hPD7IYQLOx25pD6aOnMkaYIrQwi/AtxGsUr+QxSlgjuB5wK/ADyNopzwN4EPUpQ4P77B694A3BhC+CCwfZ1t7qFYEZOmNlhZcfK+mYQQnkVxQuhrY4xfXPR4JOXNzJEkyUmXJEmSJLXK8kJJkiRJapGTLkmSJElqUa3rdIUQzgf+KMb4XSGEX6BozbkdeDVwIkVXqn3AXTHGW1ZvE2P0RGmpJ0bf7yOPfT/FNUoGFF2ePtHVGMwcKW9mjqSuLSJ3NlzpCiGcSXFhuIMhhJOB58cYXwW8n6Jl59XAzTHGayk6Uj1pzDaSemD0/b7qR68pH/8Z4Je6GoOZI+XNzJHUtUXlzoYrXTHGB4HXhRA+BuwA9pQ/uh84m+II0H3lY/uB08Zss6EQgh09tGnFGAdNtn/Oec9aObD8RJOnfCXGeG6NcYy+30cNYoxHAMoPJa0xc6T2Nc2cS5973soD+5abPMXMGcPc0WbVNHO+85+ft3L4a40yBxLPnVrlhSP2UFwPAWAXsJtitWwXRSDtKB9bvU0tMcaGw2nfgfhAre0e3/vomse+vu+x476/9/a/WbPNo/es/+fZc+fdtfatfrvhws83fs6B5Sd4/d//s/r7uOD/Pr3xTo73eAjhRIqLSG50zZN5ajdz3nPP3AY6N1tDve1OfOrax7bsPP77018y5vUnXG7vKc+rt2/1Wgg1/xsb8cC+ZeJ7/r7+Pq65wMxZR4qfdaQ2TZM5h7+2zItf+aVGz/nou78t6dxpNOmKMR4NIfxVeVHL04BrgK3ATSGEK4Db1tmmt0566pOP3f/6Q2snVk2cfPrwtcZN0tbse8epw32vmsBJg0bHjKYTQvgd4Hrgt4FbgWXgN9rfc2EzZg5PjJwasnz6bK915KHh/XGTtDX7fnhk3zvX306bUwehY+ZIOqaLDzp0lzu1J10xxheVX9+56kcHgJet2nb1NpLmrcUwGnm//1z50MfLW2fMHCk1Zs7ItmaO1LJBi5kD3edO0/JCSYkYLHVzBEiSCl5lRlJ3BoO8MsdJ15xUpYOTygbPuOSiY/er87tO2lmUEH794bXlg9ufuevYfc/v0hodLbsrUVXp4KSywQN3Du9X53cd3Vd83bJj7fYHvzi87/ldWs3MkdSlzDLHSZfUU4PMwkhS6swcSd1pu7ywa066NnDyzicdu9+0kUbVCGNcE4zRphowfqVLmsjywjxVK1HQvJHG0bIRxuouhnB8Uw0Yv9IlTZRXqY+ktFleKCkJrnRJ6pSZI6lDrnRJSkNeWSQpeYaOpA5ldqDHSVcD1TW7Zr1e11T7nlCqqM1psJTXsrvGqK7ZNev1uqbad1mq6PW6dIyZI6k7lhdKSkNeB4AkpS6zo86S0mZ5oSYabZBRp3181Tr+yReefexnj96ze832Vft4W8er4kqXgOMbZNRpH1+1jj98z/BnWy9cu33VPt7W8TrGzJHUHVe6JKXBo86SOmXmSOpSXpnjpKuB0fbxUP/crjqt4yetitV9LW0uzrk2gdH28VD/3K46reMnrYqN8twuVQwdSR3KrUuzky6ppywvlNQtM0dSdywvlJSGvA4ASUqeoSOpO7k10shrCtlDJ+08lZN2nlp7++3P3HWsqYY2ucGg/k2qHH14WH5Yx8EvDptqaHMzcyR1qUnm9CB3XOmSemqwlH7ASMqJx2kldcfyQh1rqLG6scZq28M5G7/Y5c+ex5DUYzeEMN0Te3BUJ3kvXFn0COqpmlk88dDk7SpH9q59bHVTjr23D+8/8r+Kr6Nt5FerWs4rA2MuEVCLmTMP1//2CYseQi2nnLwDgG1b6zXd2Xby2u22bT2+Ac/F512+Zpszdqz//8BnnH1JrX0rT7mVFzrpknoqt64+khJn5kjqVF6Z46RL6inLCzeRce3a6656refEkfbz41bGVtuyY3h/9aqZNom8Sn002aHH177P6656refg48PcGrcytmb7w8PzT1evmil/lhdKSoNHnSV1ysyR1J3cKnqcdEl9lVcWSUpdZh+AJKXNc7okJcGLI2smp4ycvF6VF24ZKd9Z3VZ+2zOH922qsUmZOZrN3v3DZj3bzirKCyeVEO7Zf/ex+8/YalONzcbyQklpyOsAkKTkGTqSupRX5jjpknoqt1pn1TDaUGNSI42qScZog4yqEUbVBGN5zEnsyxNWuiQzZ1MabagxqZFG1SRjtEFGtYpVrWAdPLw2t0Yfs1mGRuX2OcdJl9RTdi+U1K28Sn0kpc3yQklpyOwIkKTUmTmSutQ0c1ZaGcW8OOmS+spJ1+ZWlQfOer0uGF+OOMnqUkVtDmbOpleVAs56vS4YX444ed/Hlyoqf83LC510SWpBo/LCb7Y3DkmbRZNSH0NH0myalxemnTtOuqS+8qizNlKtYMHkVaztZSvmvbcPH6vax49rqFG1j7d1/CZj5mhj1QoWTF7F+vLuTwBw8XmXD587YTWrah9v6/jNw+t0SUqCcy5JnTJ0JHUqr8xx0iX1lBdH3uSq9u5Nz+madD7W6MpYxdbxOsbM2eyq9vFNz+matII1ujJW8bwtgd0LJaXCo86SOmXmSOqO1+mSlIa8skhS6jL7ACQpdXllTuNJVwhhF/BmYD/FX+OrwLnAduDVwInAjcA+4K4Y4y3zGqykoTbKC0MI5zDm/RtC+AHgJRSZcWeM8UNz3/n6YzJzJlkeKfOZR/v4JqpSRbB9/KZg5mDmAMPW8TCf9vHN9j0sebYMMW9LLZQXLjJzpvltLgIuBb4VeAh4fozxVcD7gauAq4GbY4zXApeFEJbnNVhJIwYNbvWt9/79PuA7gPOBe+Yw+ibMHCkJrYSOmSNpvMGg2a2ehWXONOWF91GE0T8Cfw7cXz5+P3A2xRGg+8rH9lMcGer48KuUv5Zqnc9k/Pv3L4C3AqcAvwtcPvbZ7TBz5qHOBZC3j7RirtrHb72w+Hp4zP+DqtbxYPv4zcDMMXMaqnMB5Kp1PAzbx+/ZFwE4Y0dYs33VOh5sH5+7llrGLyxzplnp+llge4xxBTgAPL18fBewG7i3vA+wA3hk1kFKWmuwNKh9a2C99++vAU+U33d9LqiZIyVhqcGtNjNH0liDwVKjW00Ly5xpXvR3gTeHEO4F/jdwJITwLuA04BpgK3BTCOEK4LYY49F5DVbSiHaOOt/KyPsXeEcI4XqK9/0fAoeAd7Sx4wnMHCkFZo6ZI3Uqr8xpPOmKMf4f4MUTNjkAvGzqEUmqp4UPQDHGBxn//v1QeeucmbOB5TEnktdpqLFRE4w65YjjXs+GGhkzc0qbO3MYXq9rVJ2GGhs1wahTjjju9Wyokac2TqNYZObYMl7qqYZlg5I0o7wuVCopbV4cWVIavGaO2rSlPHJcd6WraqphQ418mTlq2aHDxQpa3ZWuqqmGDTXy1FIjjYVx0iX1lJ9/JHXL0JHUpbwyx0mX1FNtXBxZm8C4c7XGnY+17aLi66GyPXPVOh7Gt4/XJmDmaDrjztUadz5WtXJ1+mlF3lSt42F8+3jlzfJCSWnI6wCQpNS5vC6pQy1dj3RhnHRJPeVKl6RumTmSuuNKl6Q05HUASPOwun18ndbxG7F1vI4xdLTW6vbxdVrHb8TW8cqRky6przJbdpeUODNHUocsL5SUBK/TpU6cUp68flzjjfJo8tGH125v6/iM5VXqo3Tt3V8069l21uTGGxVbx+dpkFnmOOmSeiq3I0CSUmfmSOpObp9znHRJfZVZGElKnJkjqVN5ZY6TLqmnLC/UuqqGGpMaaVQNMmBYOlg1wYBhI4zlMSeyL08oL1TG8ir10XxVDTUmNdKoGmTAsElGVTZYPHdn+dja7Koes1nG5mH3Qklp8KizpE6ZOZK6Y3mhpDRkFkZqwegq1azt48etjE0ybtVM/WbmqIbRVapZ28ePWxmbvO+1q2bqs7wyx0mX1FOWF0rqVl6lPpLSZnmhpCTktuyuBapzAeTtI62Y995efN0yciR59fldVet4sH18NswczU+dCyB/efcnjt2/+LzLi+dNWM2qWseD7eNzMMgsc5x0SX2VVxZJSp0HeiR1KLeDy066pJ4aLOW17C4pdWaOpO5YXigpDXkdAFIblkdKb5o20qgaYYxrgjHaVKNi+/hNwNDRxqrW8dC8kUZVOjiuCcZoU43h69ssI295ZY6TLqmvMlt2l5Q4M0dShywvlJQEuxeqkap9/Kyt46cxadVMPZJXqY/aN7yg8Wyt46fb9/qrZuqHQWaZ46RL6qncjgBJSp2ZI6lDmX3OcdIl9VVmYaQE1L0ActU+vmodD7D1wuLr4XvWbl+1j7d1fL+ZOWpB3QsgV+3jq9bxAHv2RQDO2BHWbF+1j7d1fH/ZMl5SEiwvlNStvEp9JKXN7oWS0uBRZ0mdMnMkdSe30yicdEk9lVkWqW3Lq04mr9tQo27r+EnliHVeS+kzdNTQaPt4qN9Qo27r+EnliHVeS6nLK3OcdEl95cWRJXXKzJHUHcsLJSXBg85auC0jR47rrHTZUKPnDB0t3qHDIxdfrrHSZUON/rKRhqQ0OOuS1CUzR1KXMsscJ11ST9m9UFK38ir1kZS2JcsLJSUhsyNASkzVJKMqG6yaYMCwEca2i4aPHSpKeCZer0s9Z+aoXVWTjKpssGqCAcNGGFW5IMDppxV5M+l6XeqzvDLHSZfUU7m1UpWUODNHUody+5zTeNIVQjgXeBOwB/ga8BTgJGA7cDUQgNcDh4A7Yox/PK/BShrRQnlhCOEc4EZgH3BXjPGW8vEXAS8Bvg78VYzxT9Z5/vXA78UYa3RVqD2mczFz5md163io3z5+ktUrY5OMWzVTD8y/1MfM2RxWt46H+u3jJ1m9MjZx2zGrZkrbIMHMKbedKnem+W2uB+4HzgL+CdgeY7wO+DjwUuC1wPUxxlcAr5zi9SXVMBgMat8auBq4OcZ4LXBZCGG5fPxVFAH1FOAzE55/ELg9hPCREMK/DiHMY2Zo5khJGDS41WbmSBpvMGh2q2fWzIEpc2eaSdf5wO3AVcDLKYKJ8uvZwOkxxt3lYytTvL6kOuYfRABnAveV9/dTHNkFeBbwVuDNwFvWe3KM8T0xxkuAXwV+AvhKCOHXQghPafrrjTBz+uCUUNwqW3YOb6tte+bwpv4wc8ycxOzdfw979w/PHz14+OFjt9X27L/72E39MGj4T00zZQ5MnzvTTLoeBB6NMT5Rfn9O+XUXsBt4IIRwVvlYXsWYUkIGS4PatwbupXgvA+wAHinvfwV4HFj7f7IRIYSnhBCuAd5HcbToOuAuYN1l+hrMHCkJSw1utZk5ksYaDJYa3WqaKXNg+tyZppHG24EbQgh7gFuBi0MI76SYKV4FfBa4MYRwBLh5iteXVEc7J5jeCtwUQrgCuA14R1m7/B+BDwHfoMiA9fwt8PvAj8YYqyNJhBC+c4YxmTlSCswcM0fqUEuNNGbNHJgydxpPumKMXwR+ZMImXwBe1vR1JTXTRhbFGB9k/Pv3I+VtI2+PMb6v+iaE8PMxxptjjG+YYUxmThtGG2pMaqQxrkFG1QhjtAnG8qoT2Udf/+iGBw7VC/MPHTNncxltqDGpkca4BhlVyeBoE4yDh4/PrtHvbZaRgyQzB6bMHVvGS321lM5FA0MI/w64HHhhCOGF5cMnUNRIeyRYyoKZI6k7DUoGOzFr7jjpknoqsctXfAz4KrATeG/52DeBLy1sRKqvWqXqunU8jF81U5rSCh0zp+eqVamuW8cX+167aqb0NGiO0ZWZcsdJl9RXCa10AdtijH8dQti96vFTFzIaSS0wcyR1J7WVLmbMHSddUk+lddCZ64FXMzzyU1kBLu1+OJLmL6nQMXOk3CX2QYcZc8dJl9RXCYVRjPHV5dcXbrSteqoqG4TJpYPbLym+7r19+Fh1ra5xDTWqa3UduHO28al9Zo46VpUNwuTSwS/v/gQAF593+fC5E0oIq2t1PWPrJXMZp9qRWnnhrLnjpEvqqYbX32pVCOGrrL1I6ABYiTGevYAhSZq7dEp9zBwpf6mVF86aO066pL5K66jzWRtvpWRV7d2bNtKommDA2kYYoytjx7axdXy/mTman6p9fNNGGtUKVvHc41exRlfG1ttG/dHSdbqmNmvuOOmSeiqlMAohvDHG+JYQwodZdRQoxuj1bKQcmDmSOpVO5sDsueOkS+qrhMoLgY+WX9+z0FFoNqMXOJ5H+/gmJq2aKRFJlfqYOZk4/oLGs7ePb7bv9VfNtHiplRcyY+446ZJ6KqWVrhjj58q7nwd+GbiwvP+2hQ1K0pyZOZK6k2AjjZlyJ7kppKSaBg1u3fk94B+ANwIPAB/qdO+S2jMY1L91x8yRctUkc3qQO650ST01SOviyJWTY4zvLu9/LoTw0oWORu2ommSMto6vygOr0sDtI62Yq/bxWy8svh6+Z+1rVq3jwfbxyTJztDhVk4zR1vGr28JXreNh2D5+z74IwBk7wprXrFrHg+3jU5RgeWFlqtxx0iX1VUKr7iGE8tM0D4UQfhj4n8BzgH9c3KgkzVc6oWPmSPlLrbxw1txx0iX1VErndHH81dmvLW+w9noWStnymBPJ59FQY9zK2CSrV82UBjNHLahax4+aR0ONcStjE7efcDFlLUZqky5mzB0nXVJfJdS9cL2rs4cQlrsei6S2pFPqY+ZI+UutvHDW3HHSJfVUYitdAIQQrgZeAyxT1CI9QdHdR1LvmTmSOpTg5xyYPnfSmkJKqi/Njj5XAS8A/jvwU8AXuty5ErVlZ3Fb/f2Wdcp4tj3z+MYaSoOZox45dHgfhw4PSxcPHn742G2cPfvvPq6xhhZv0PCfDk2VO066pJ4aLA1q3zr0UIzxq8CTYox/DezYYHtJvbHU4NYZM0fK1GCw1OjWoalyx/JCqa/SXHY/EEL4N8BKufx++qIHpBadOPKvt2qSMa4JxraLiq+HyqPIow07jo4/6qwUmTlavKpBBgybZIxrglGtWp1+2oXlNiPPs1lGP6T5OQemzB1XuqSeGgwGtW8duhL4MvA6ivrmV3a5c0ktSrO80MyRMpVweeFUueNKl9RXCXUvHHEE+D6KEPo8xTUs1Eer28fbOl5pHqc1czKyun28reM3t9S6F46YKneS/W0kTZboSteHgHOAvwAuAD7Q5c4ltWnQ4NYZM0fKVMIrXVPljitdUl8ludDFmTHGHyvv/0kI4eMLHY2k+Unz/AozR8pVmpkDU+aOK11STw2Wlmrf2hZCODGEcCLwjyGEZ5ePfTtwT+s7VxpOPP34xhrjnBKK2yhbx/dIOt0LzRxBUTo42lhjnL3772Hv/uP/s7B1fD+k1r1w1txxpUvqq7QOAEVghWJULwghfB04CXh8oaOSNEdJhY6ZI2Wu45LBOmbKHSddUk91fK7WRDHGZ1T3QwgDivapD8UYv7m4UWkuqoYaTRtpbBm5bEnVCGN5zIns1evbOj59Zo46UjXUaNpIY3T1qmqEMdoqfrjdQ8dto1Slkzkwe+5YXij11dKg/q0jIYQXAF8C/gz4fyGEf9XZziW1LJ3ywoqZI+VrabDU6NaVaXPHSZfUU4l2L3wL8LwY43cBl5TfKwfLTx3eZlWd/7XROWCVLTuGNy1Qkt0LzZyMHTz80LHbzK9Vnv+10Tlgw30/PPHcL3WgybUBe/BZx0mX1FdpBtE3Yoy7AWKMD+D5FVI+zBxJHUq4ZfxUueM5XVJPDdK8OPKjIYSfA/4H8HzAq9pK2UjyOK2ZI2Uq4YsjT5U7TrqkvkropPYRLwfeCLwV+ALw04sdjjpXlQwe2bv+NtsvGd7fe3vxdeuFw8cOr+q+O9o2/sCds41PMzBzlKaqZHDbyeuXP3959yeO3b/4vMsB2LMvHnvsjB3HX85itG38M7ZeghYgzc85MGXuOOmSeirRLHpXjPHHFz0ISS1IM3TMHClTCbaMr0yVO1NPukIIfwB8FPgW4FxgO/Bq4ETgRoqltrtijLdMuw9JE7Rw0eMQwjms8/4NITwZ+BRwaYzxwXVe4uSRCwV+EyDGeGROYzNzFml5pLXytO3jj46pwKizMqZEmDmYOZ2qWsfD9O3jx7WFr7MypsVro7xwDpkDU+bOVL9NCOE1wGPlt8+PMb4KeD9wFXA1cHOM8VrgshDC8jT7kDRZS+e0j33/hhCWgBsoWqROciFwB7Af+Hvg7smb12PmSClopXuhmSNpHU0yp3buzJo5MGXuNJ50hRBeDDwCfLJ8/p7yR/cDZwNnAveVj+2nODIkad7amXWt9/79VeA/ARv1zv1ViqM+dwNHgWua7HwcM0dKhJlj5kgdanJpnAaXx5k1c6ptG+fONOWFL6cYZHXGYXUkaBewmyKgdlH8QjsogkvSnLXUvfBeVr1/QwinA/8CeBrwXOAXgdes8/w3Ac+JMe4JITyNojTnz2cck5mTmupaXU3LDOdhUqmiWtZKJzEzR7VU1+pqWmY4n32vX6qo9rTUvXDWzIEpc6fxpCvG+KMAIYQrKPrSPy2E8C7gNIqZ3lbgpvLnt8UYjzbdh6Qa2jmp/VZG3r/AO4DrY4w/ABBC+CDw9gnPfzjGuAcgxvhPIYRHZx2QmSMlwswxc6QOtdRIY9bMgSlzZ+pGGjHGD67zowPAy6Z9XUn1DFpopFGeOLru+zfGeMUGL/FoCOHPgI8D3wOcEkL4zfK5vzzj2D64zo/MnBRVDTJg2CRj3CpV1T6+ah0Pw/bxq1vHw7B9vK3jF8DMKZk5iaoaZMCwSca4VaqqfXzVOh6G7eNXt46HYft4W8d3q42VrjlkDkyZO7aMl/oqzU6qfzJy/4GFjUJSC5IMHTNHylaSmQNT5o6TLqmnGpw02pkY44cWPQZ1YHnMeQ2znt81bmVskmrVDDy/qytmjhZotH18Zdbzu8atjE3c/vCwx4Lnd7Uvxc85MH3uOOmSeqqlRhqStI5WTmqXpLFaaqSxME66pL5K9AiQpFyZOZK601IjjYVx0iX1lZMu5WTLSKnOscYb5WNHx1w2pWqoATbV6IqZo8wcOjwsWZzUeKNSNdQAm2p0IrPMcdIl9ZTlhZK6lVepj6S0DTLLHCddUl9ldgRIPTTaUGNcc43KKWtbMK/x1Ms33kZzUuPfx1hmjhZvtKHGuOYalb2MueTEKnfx0bmMSRs5f6pnpdpIY1pOuqSeyiyLJKXO0JHUqbwyx0mX1FNtXBxZktZn5kjqjt0LJaUhrwNAkpJn6EjqjuWFktKQWRhJSpyZI6lDtoyXlAS7F0rqVl6lPpLSZnmhpCTktuwuKXVmjqQu5ZU5TrqkvnLSJalLZo6kDuV2cNlJl9RTlhdK6lZepT6S0mZ5oaQ0ZHYESFLqzBxJXcorc5x0SX2VVxZJSp0HeiR1yPJCSUnw4siSumXmSOrOILPMcdIl9VRmB4AkJc/QkdQdV7okpSGzMJKUODNHUqfyyhwnXVJPWV4oqVtmjqTu2L1QUhryOgAkKXmGjqTuWF4oKQ2ZhZGkxJk5kjqVV+Y46ZJ6yosjS+pWXqU+ktK2ZHmhpBTktuwuKXVmjqQu5ZU5TrqkvnLSJalLZo6kDuV2cNlJl9RTlhdK6lZepT6S0uakS1IaMgsjSakzcyR1Ka/McdIl9VVeWSQpdR7okdQhV7okJcGLI0vqlpkjqTteHFlSEjI7ACQpeYaOpC7llTlOuqS+ctYlqUtmjqQODTb7pCuE8L3AdcBjwFeAU4GTgO3A1UAAXg8cAu6IMf7x3EYr6Zg2uheGEM4BbgT2AXfFGG8pH38l8GyK9/t/jjF+dO47X39MZo6UhPmX+pg5ktbTRnnhIjNnmt/mNODKGOMrgOcB22OM1wEfB14KvBa4vvz5K+c2UknHGwzq3+q7Grg5xngtcFkIYbl8/JEY409TvKdfNuffZCNmjpSEQYNbbWaOpLEGg0GjW00Ly5zGk64Y458Cj4UQ3gDcCdxf/uh+4Gzg9Bjj7vKxlbmMUtIag6VB7VsDZwL3lff3UxzZJcb44RDCqRRHh26Y6y+XZSpJAAAHV0lEQVSyATNHSsRgqf6tPjNH0liDwVKjW00Ly5zGk64QwpOAW4FPAR8Azil/tAvYDTwQQjirfCyvYkwpJe2sdN1L8V4G2AE8AhBCuAh4D/DrMca/m+vvsQEzR0pFKytdZo6kdTTJnNpvxYVlzjSNNN4JXAD8FPANivB5J8VM8Srgs8CNIYQjwM3zGqik47V0/YpbgZtCCFcAtwHvCCG8FrgD+BzwGyGEz8cY39bGztdh5kgpMHPMHKlDLTXSWFjmNJ50lfWOk3yB7uuvpc2nhUYaMcYHGf/+PX/uO6vJzJFSMf+T2s0cSetpo5HGIjPHlvFST+V2pXZJqTNzJHUnt885TrqkvsoriySlLrMPQJJSl1fmOOmSemqwNP9ld0lan5kjqTttlBcukpMuqa886iypU2aOpO5YXigpCZllkaTUGTqSOpVX5jjpkvrK8kJJnTJzJHXH8kJJSfCgs6RuGTqSutPSdboWxkmX1FfOuiR1ycyR1CHP6ZKUhEELF0eWpPXlVeojKW2WF0pKQ2ZHgCSlzsyR1KW8MsdJl9RTuS27S0qcmSOpQ7l9znHSJfWV5YWSOpVXqY+ktA0yyxwnXVJP5XYESFLqzBxJHcrsc46TLqmv8soiSanL7AOQpLTZMl5SEgZeHFlSp8wcSd2xe6GkNOR1AEhS8gwdSd3J7TQKJ11ST+UWRpISZ+ZI6lRemeOkS+oruxdK6lRepT6S0mZ5oaQkuNIlqVtmjqTu2EhDUhqcdEnqkpkjqUuZZY6TLqmnBpYXSupUXqU+ktJmeaGkNGR2BEhS6swcSd2xvFBSEpxzSeqUoSOpQ7mdu+6kS+orL44sqVNmjqTuDDLLHCddUk/ldgRIUurMHEkdyuxzjpMuqa/yyiJJqcvsA5CktHlOl6QkDCwvlNQpM0dSd+xeKCkNeR0AkpQ8Q0dShzJbXXfSJfWUK12SOpXZUWdJaVvKbHXdSZfUV3kdAJKUPENHUodc6ZKUArsXSuqWmSOpOzbSkJSGpbzCSFLiLC+U1CEbadQQQjgHuBHYB9wVY7yljf1Im1kbK13rvXdDCN8P/CTFoe53xxg/Mfedz8DMkbpg5lTMHKkDmX3OaWsKeTVwc4zxWuCyEMJyS/uRNq/BoP6tvvXeu68BrgR+Bvil+f4ic2HmSK0bNLjVZuZIGmvQ8J+aFpY5bZUXngncV97fD2wHHmppX9KmNGinvHC99+4gxngEIIRwchs7npGZI7WtnVIfM0fSWC2VFy4sc9qadN0L7KL4pXYAj9R5UgihpeFI2flKCOHpTbavud16793HQwgnAicAjzfYb1fMHKldZs7xpsocgDvefX5bY5Jy0jRzoF7uLCxzBisrK3N/0RDCmcBNwNeAT8cY3zf3nUiau9XvXeDbgeuB5wKvAJaBG2OMn17YIMcwc6R+MnMkdWmRmdPKpEuSJEmSVMirF6MkSZIkJcZJlyRJkiS1yEmXJEmSJLXISZckSZIktchJlyRJkiS1yEmXJEmSJLWorYsj1xZCOAe4EdgH3BVjvGVB4/he4DrgMYqLq50KnERxpeqrgQC8HjgE3BFj/OOOx/cHwEeBbwHOLcf1auBEFvD3CyGcC7wJ2ENxrYOnkMDfK4SwC3gzxVXGB8BXSePvdT7wRzHG7woh/MJGY1q9TYxxbxfj3CxSyB0zp/F4ziXBzCnHllzumDlpMXNqjc/MqT82M6eHUljpuhq4OcZ4LXBZCGF5QeM4DbgyxvgK4HnA9hjjdcDHgZcCrwWuL3/+yi4HFkJ4DUVIAjw/xvgq4P3AVSzu73c9cD9wFvBPpPP3ugi4FPhW4CES+HuVF+K7EjgYQji5xpieNGYbzVcKuWPmNJNq5kBiuWPmJMnMmcDMaczM6aEUJl1nAveV9/dTzHg7F2P8U+CxEMIbgDsp3miUX88GTo8x7i4f6+yK0iGEFwOPAJ+k+Pe1Z9W4FvX3Ox+4neKN8nIS+XtR/C0uBX4EeCFFGI2Oq/O/V4zxwRjj6yj+h7KDjf8dnjZmG83XwnPHzGks1cyBxHLHzEmSmbMOM2cqZk4PpTDpuhfYVd7fQfHG61w5674V+BTwAeCc8ke7gN3AAyGEs8rHBh0O7eXAc4B/T/HGP2PVuBb193sQeDTG+ET5fSp/r5+lOBq1AhwAnr5qXIv+720PsHODMe0es43ma9H/HZg5zaWaOZB27pg5aTBz1mfmNGfm9NBgZaXryfnxyiXJmyjqZT8dY3zfgsbxAeACijrnbwAPANsojg5cVf7sjcAR4CMxxjs6Ht8VwOPA0yjqiE8DrgG2soC/XwjhmcCvU7y5PglcTAJ/rxDCd5fjuhf4crn/hf+9yrF9LMb4ohDCdRuNafU2McYDXY1zM0ghd8ycxuNJMnPKsSWZO2ZOOsycWuO7AjOn7tjMnB5a+KRLkiRJknKWQnmhJEmSJGXLSZckSZIktchJlyRJkiS1yEmXJEmSJLXISZckSZIktchJlyRJkiS1yEmXJEmSJLXo/wPfgq5ilsk8OQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_probabilities(plot, prob_section, geo_model.grid.sections.resolution[0],\n",
    "                  T=True, extent=[0,1000,0,1000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "(3, 10000)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA10AAADvCAYAAAAJkb68AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuwZWdZ5/Hvbvp0utNgk+50CJ0ICULeUDB4hRlIYAg6pQXDZUypDCI2l1QALwhBC4w4XoAgVgLECRcnxsRR0NFQRNBCyzIyQwAHahgcqPAiiCQkxKS707l0+hrO/LHWjn1On8tea+/1rne9+/uhdu1z9tmX9xzoH+tdz/u8a7S4uIgkSZIkqRsb+h6AJEmSJJXMSZckSZIkdchJlyRJkiR1yEmXJEmSJHXISZckSZIkdchJlyRJkgSEEE4LIZzU9zhUno19D0CzEUL4NeDjMcbPtHz9qcBVwN3AVuCVMcYjsxuhpJJMmzn1e2wDrgF+e5r3kVSOWWTLlH4LeDNwR0+fr0I56SrLc0IILwN2AX8EbAfujTF+KITwP4BXAZcBDwInAT8P/C3wZeBK4G0xxn8IIfwOcDYQe/gdJA3HNJnzHuClwAN9DFxS1t4YQrgDuBU4FzgAfBdVZrwYeALwaOBS4CzgPwEPAz4HfLV+/NPAk4DPA98N/DZwGHg1sAjcAvwx8CHgY8D3U2XW9wEvA97Z+W+pueLywrLcBPwJ8BngucAfAC8OITwJ+Brwk8AOqoOchwNPBh4WY3xljPEL9YTrucC3Y4xOuCStZ5rM+YcY4y/Vz5Ok430gxvizwI9QZczH68efTDXJuh/4feBbVFWp+4D9wAX18z4bY/xV4FTgHcAHgGfXz70fuBd4BlXx4asxxrcB/7d+/89TZZk0U066yvKrwHcA/xsYxRgPUoXI5cB7qf77/usY45uoQuw2qpACIITwFuAJMcbXpR64pEGaKnMkaRUHjvv6IqrWh68AI+DPgA8CPwi8nGri9Faq6tbf1q+5t74/GGM8BhyjyqOHAe+rM+lj9ePjzzpSP2exm19J885JV1mOUoXQj1EdCEFVOt8fY7wV+EPgeSGE9wI/Ctw1fmG9ROingKeGEP4whHBO0pFLGqLWmSNJa3hVCOFXgQ9TVcl/lGqp4A7g6cAvAo+iWk74duBaqiXOh9Z538uAy0MI1wLbVnnOV6gqYtJMjRYXndCXKoTwZKqG0DfGGG/uezySymbmSJK0MiddkiRJktQhlxdKkiRJUoecdEmSJElSh7xOl6QlQgiPB/40xvi9xz32Q1TXLRlR7fz0qb7GJ6ksZo6k1PrIHSddkh4SQjid6uKQB5b96A3Ai6i22/1j4IWJhyapQGaOpNT6yp0kk64Qgrt1SKuIMY6aPP9pj3vy4j0LR5t+zDdijGdNMJY7gDeFED6+7EejGOMRgBDC5qYfnpqZI62uaeY85+mPW7xt30LTj5mrzAFzR1pN08z5nh943OLB+xpnDmSeO8kqXTHGVB+lAbt2dGHfQ0hi9+L1AIQQGr/2noWjvPkf/02j11z2hP/32MYftNShEMImqrM/610HJQtFZM7RvfX9nsmef2TZZbCO7Vv5effctPT54885tnft9z9w89rvO0QXzMlx8o3VMU94dfNLMN62b4H4/n9s9Jrw6ifMXeYA/MfXfLXvIUzt5M3bAdi65dSJnr9189Lnbd2yY8XnnbXrGUueP37/1Z4/dtop5070POWpzXHOwfsWeP5rvtb4dR9933dlnTsuL5QSGk+2pjVqdM6ovRDC7wCXAO8GrgYWgN9M8+mSpnbjjMIiUeiYOZKAdAc6pMsdJ13SEHUcRjHGH6nvf65+6BP1TdJcMnMkpTPqOHMgfe446VIW5mVZ4ayMEp4Bkh7issI5ZuaoHy4rnFMFHuc46ZKGaEN5YSQpYyMv6ykpnVGBmeOkS72alwrXrHq5xqx0SS3NS4VrVr1cDzFzJKWTYnlhak66pCEqL4skZc3QkZRQgSeXnXSpF/NS4erKaEN5ZXf1aNKt4odsXipcXSlwqY/6NelW8ZpPGygvc5x0SUNU3gkgSVkzdCQlZKVLUhOz7uUas6dL0opm3ss1ZuZISseeLkl5cNIlKSkzR1JCBR7nOOlSUvZyzcbILeM1jSN39T2CdOzlmg17ujSlrZtP7XsIGhC3jJeUhwLPAEnKmZkjKR2XF0rKg5MuSUmZOZJSKi9znHRJA+TyQklJFbjUR1K+NhSYOU66lIS9XLPl7oVK6ti+vkfQnL1cM2bmKK31ruOlwhV4nOOkSxqi8rJIUtYMHUnp2NMlKQujDeWV3SVlrMClPpLy5e6FkvJQ3gkgSVkzdCQl5PJCqZl57+Ua//67F6+f7RsXGEZK4OievkfQvXnv5Rr//jfOOiPMHLWzdYvX51JzLi+UlAU30pCUlJkjKanyMsdJlzRAbhmvJY7cNeHz5qDSpY6U11+h6WzdPFkFa9LnScezp0tSHjzrLCkpM0dSOiWu6HHSpZma9x6uZMrLIqVwbG+z5x+4uZtxzNK893ClUuABkNI4ecv2Rs8/7ZRzOxqJhsSeLklZcMt4SWmZOZLScXmh5p6VrHZmvotheSeANEvH9k32vKMNK199sJLVzsx3MTR0tLatW3ZM+Dx7vDSJ8jLHSZc0QCWudZaUMTNHUkIlHuc46dISVrIGosAwUg+a9nh1wUrWQJg5mo1JK2Kad+VljpMuaYDcMl5SWuX1V0jKV7uerrxP4jnpEmCFa3CsdGnorHANi5kjKaF2ywvz/v8VJ13SADXOorxzSFL2DB1J6bhlvKQ8NN0y/sFuhqGBOXJX3yPQYDVd6mPoqLJ1s7sVqrl2ywvzzh0nXdIAudJHUlKGjqSkysscJ11zzl6ugfIASE0c3dPudZNe76sJe7kGysxRM22vx+XuhgK3jJeUiRLDSFLOzBxJKZWXOetOukIIZwK/AdxN9Rf4FnAWsA14fYzRJgEptYK3jDdzOnCkZaVLGmvVXzEMZk437OXSNNr1dOVtkkrXucBzgM8CnweeFWN8fgjhAuAi4O0djk/SCrqqdIUQzgAuB/YBX4oxXlU//sPAC6ky46YY43WdDKBi5kjZMXMkpVPicc4k08hbqcLox4ELgPEp028Cu2Y9IEkTGI2a3SZ3MXBljPG1wPNCCAv1488Evht4PPCV2f4yJzBzZu3Y3uomtTZqeJuYmVOok7ds5+Qt2/sehgZq1OI/E+otcyaZdP0ssC3GuAjcAzy2fvxM4PYuBiVpbaMNo0a3Bk6nOgCBaqnNtvrrvwF+CPgx4M0z+0VWZuZIuRltaHabnJkj6QSj0YbGtwn1ljmTLC/8feA3Qgi3AH8PHAkhvBc4BXh1F4OStI7uNtK4hepA41ZgO7C/fvzXqMJoP91vwGPmjLXddVCaOTNnXrTddVCarfIyZ903jTH+H+D5XXy4pHY63LzwauCKEMJu4MPAu0IIl1AdlHwQeAB4V2efjpkj5amz0DFzJJ2gw12ae8sct4yXhqijMIox3gG8ZIUfXVffJM0jM0dSUuVljpMuaYAa9mlJ0pTK275ZUr7mdct4Feja0YV9D0HT8OLI82Gh7q2Ytrdr447q/kiPlxu6YLG/z9YMmDnz4sDBKm+m7e164OC+6n28Xpda6HB5YW+cdElDVGAYScqYmSMpoQZbwA+Gky5pgFxeqEY21WeaH+h3GBqy8pb6qFsHDlUVs52c0/NINEQuL5SUhRLL7pJyZuZISqm8zHHSNWfs5Upr9+L13bxxeVmkUtnLldaNHYWDJ3okJVTiyWUnXdIQFRhGknJm5khKqbzMcdI1B6xupdNZZWsZe7rUyELD3cMO3Dzd51ndSqerytYJyuuvULfGuyBO6rRTzu1oJBoiK12S8lBgGEnKmJkjKSEnXRoUK1zdS1XZWq7EMFIBrHB1L1llazkzR1JK5WWOky5pgFxeKCktlxdKSmeDW8ZrCKxwzU5flax1WelSTqxwzU5vlax1mDmSkiovc5x0SUNUXhZJypqhIymdEtsonHQVZB4rXNlWojo22lBe2V0DNI8VrlwrUZ0zcySlM3J5oaQsFHgGSFLGzBxJSZWXOU66lLV5rWStx+MfqSNzW8laj38XSem4vFBSHgoMI0k5M3MkpVRe5jjpKkCJvVxWuNbmlvFqZdPO6v7IXdO9T4m9XFa41lZgf4XSOHBoDwBbN5/a80g0JG4ZLykPVrokJWXmSEqpvMxx0iUNUIlrnSXlzMyRlE6JxzlOuqQhcnmhpJQKXOojKV9uGa8sjfufSujtspdrMiWeAVIC0/ZyjY37n0ro7bKXa0L+ndSOvVxqp7zMcdIlDVF5WSQpa4aOpHRGBWaOky5pgKx0SUrLzJGUUIHHOU66pCGyp0tSSgX2V0jKl1vGS8qClS5JaZk5klIqL3OcdElD5KRLUlJmjqR0Sjy57KRLWXDXwmZGLi+UpuOuhc0UuNRHUr7cMl5SHgo8AyQpZ2aOpJTKyxwnXdIAOecq3NG9M36/PbN9P80hQ6d0J2/ePtP327rF63OpPbeMl5QHZ12SUjJzJKVUYOY46VIWrh1dCNjbNSl7uubErCpUD8Rl7zvjStoQXbBY3dvbNaHy+iu0sllVqHaeck4n76v54JbxkvJQ4BkgSTkzcySlVF7mrDvpCiGcBbwFuBO4D3gkcBKwDbg4xni4ywFqPljhaqbErVTHzBwlYYWrGTNHUkIlHudMUum6BPgm8Fjgq8BjY4wXhxBeDlwIfLDD8UlaSUfLC0MIZwCXA/uAL8UYr6of/xHghcBh4MYY4w2dDKBi5kjZ6Wapj5kjaSWjAjNnkt/o8cBHgIuAl1IFE/X9rlkPSPNl9+L1VrlaGI1GjW4NXAxcGWN8LfC8EMJC/fjPUAXUI4HPr/biEMIlIYSdbX+vmpkzqSN3VTdN7saRVa5WRg1vEzNzBmbr5lPZutn+LHVsNGp+m8xUmQPtc2eSStcdwL0xxqMhBIAz6sfPBG5v+oGSZqC7svvpwK3113dTLa/ZAzwZ+In6528FXrbK6w8AHwkhfAv4PeDjMcbFhmMwc6TcmDmSEupwy/hpMwda5s4kk653ApeFEO4ErgaeFEJ4Tz3IiyZ4vfQQq1oz0t2J+luoDjRuBbYD++vHvwEcAtbc9i7G+H7g/SGEJwGXAh8IIVwDvDvGuH+t1x7HzGnr2L4JnzdHuxda1ZqRzv6OZs7Abd2yY6bPk4AuT/RMlTnQPnfWnXTFGG8Gfny950lKZ7Shs61UrwauCCHsBj4MvCuEcAnwX4HrgAepDlBWFEJ4JPBiqjNE+4HXUeXMDcC/n2QAZo6Uoe62bzZzJJ2gwy3jp8ocaJ87bhmvJKxwzVZXJ4BijHcAL1nhR39W39bzWeAPgZ+IMY7L94QQvmc2Iyzc+PpZs7o+V1uTVsxyZoVrxrr5e5o5/Tt583ag/+toWQnTUtlmDrTMnfKuPCbNg26aS2fhnTHGXx+HUAjh5wFijJemHISkWetsI41pmTlSgZpuGJZ4i/lWuWOlS52wstWtUUdbxrcVQvjPwAuAC0IIF9QPP4yqMfXK3gam+WFlq1vdLfVpxcyRytbVlvHTmDZ3nHRJQ5TfRQM/DnwL2AF8oH7s28DXehuRpBkycyQllN9xDkyZO066NFNWuBLJL4y2xhj/LoSwfHvlh/cyGi1V8rW8rHAlkt3f2czJnNfy0jQ63DJ+GlPljpMuaYDym3NxCfB6/vXMz9gi8Jz0w5E0W9mFjpkjlSzDAx2mzB0nXZqKla1+dLhlfCsxxtfX9xes91xpKla2+pFZT5eZI5VtlFnmwPS546RLGqLMjjvrq7Ivvxr7CFiMMe7qYUiSZiqv0DFzpLLluLxw2txx0rWGQ3vvA2Dzjkf0PJJ2rEIVLLOye4zx0X2PYdDG1ZuN1fVy+P5Prf38SXu0Ul/vyypUwfL679bMmd7lv/Dgku/fcd0T13z+pD1afV/vS2VIvAX8RKbNHSdd0gBluGX8r8QY3xpC+BDLzgLFGFe6CKGkIclsqY+ZI5Ut0y3jp8odJ11rOLzn3iX3y20LZ6QczrqsbM2R/M4AfbS+f3+voxiK1SpCCzur+3944co/f8oNKz9+bN/Kjx/puNJlZWuOZPfftZnT0PLK1nJv+umbV3x8tQrY1i07Vn7cXQs1C/kd58CUueOkSxqg3MruMcYv1F9+Efhl4Jz663f0NihJM2TmSEonx56uaXPHSdcKxr1c67kn3lY9/66llbDD++5f8v1jXvDU2QxMquW2vPA4fwB8rL5/JnAd8KJeR5ST1SpD416u9YwrYJuWnUneuOyM865XVvfH9k4+NmktmS0vPI6Zs471KlzrWa0Cttw1f1792U/eMmGeSWvIcffC47TKHSdd0hBlVuk6zuYY4/vqr78QQriw19FImhEzR1JC+R7nQMvccdKVwC1//lkA7rzpywD8wG/9VJ/DUQkyy6IQwjn1l3tCCD8G/C/gacDX+xtVRlL3Pt3+e9X9tvOq+0l3O5RWlVfomDnrm7bC1dQrXvCRJd///RevSfr5KkuOywunzR0nXdIA5dbTxdKrs7+2vsGJ17OQNERmjqSEcpx0MWXuOOnqwSdfcdWS7+/9yu1Lvn/uJ9+Wcjgaosx6ula7OnsIYSH1WLKSy+5+jzx/6fePfsXS7//pTenGooHKq7/CzFld6grXav7tk1+x5s+/fvtNiUaiIcqxp2va3HHSJQ1QhpUuAEIIFwNvABao1iMdpdrdR9KgmTmSEsr0OAfa546Triks37WwqcN7p3u95li+YXQR8GzgV4A/BX6h19H0pasK1/JdC5ta8Po5asnMyV4uFS5pFjJdXjjWKnfyq90V7NBd9y6ZqB3eez+H995/wvP+8vxL+cvzL005NA3MaMOo0S2hPTHGbwGPiDH+HeDewX3adOrSidrG7StvT/+4d1Q3aVUbGt6SMXMG6uxd53H2rvP6HoYyNRptaHxLqFXuWOmShijfs873hBBeBCzW5fedfQ9I0iyYOZISyvc4B1rmjpOuFRze027Z3/KLIk/7ftJqcu3pAl4FfBfwJuCNwGv6HU5ibZcVLrQ8Tlx+UeSH3m+Vx6W2zJxsuaxQJcp8eWGr3HHSJQ1Rvll0hOrq7OcAX6S6hoWkwcs2dMwcqUCZT7pa5Y6TruMc2ntfN+/bstI17utyC3ktN9qQbTvmdcA/A39DFUjXAD/d54CSaFvhWqm/ahY2taycjfu6Pv/M2Y1FhTBzclNKhcu+Lq0kxy3jj9Mqd5x0SUOU7wmg02OML66/viGE8IleRyNpNvJdXmjmSCXKN3OgZe446Wqh6VbxK+1QKE0jt56uEMKm+suvhxCeGmP8bAjhKcBX+hzXzPV1seOmW8Wv1usltWbm9KWUipbURI7LC6fNHSdd0hCl3QZ+EhFYpDoye3YI4TBwEnCo11FJmpHslvqYOVLBMl1eOFXuOOmSBii3SleM8ezx1yGEEdX2qXtijN/ub1Qz1FeFS8qFmZOcFS7Nt7wyB6bPnSynkZLWMRo1uyUSQng28DXgr4B/CiH8h2QfLqlDo4a3NMwcqUyj0ajxLZW2ueOkawYO77t/1Wt0SZ3I8/gH4K3A+THG7wXOq7/XrG3cYd+WEss2dMwcqUCjFv9JqFXuOOmSBmi0YUOjW0IPxhhvB4gx3ob9FVIZRhua3dIxc6QCjUYbGt8SapU79nRJA5RZe8Xx7g0h/BzwP4FnAft6Hs907OWSatn+Wygrc7CXSwKyPtChZe5Y6ZKGKNOeLuClwGOAtwHfCbwi5YdL6kq2ywvNHKlAmS8vbJU7VrqAQ3vvA+DwnmbX31rP+P2aXtdLWs8ovy3jx94bY/zJvgeRvY3bq/uFnbN934W6z6vpdb2k9eS5fTOYOVKRMt0yfqxV7kw06Qoh/BHwUarZ3FnANuD1Mca7mn6gpBnIt+y++bgLBX4bIMZ4pOmbmDlSbsrOHDB3pLxkmznQMnfWnXSFEN4AjLfme1aM8fkhhAuAi4C3tx9v+e6peuxOcHivlS9NqaNJVwjhDOByqvXJX4oxXnXcz74D+AzwnBjjHau8xTnAx6iuXXEn8CDwuIZj6D9zhtrLdfK5Kz++YOVL0yo3c+rP6j137OWS/lVXW8DPIHOgZe6sWbsLITwf2A98un7unfWPvgnsWu/NJXWjw5aui4ErY4yvBZ4XQlgACCFsAC6jui7FWv4L1VmfLwPHgFc3+XAzR8pVZz1dvWZO/VnmjpSdppkzce5MmznQMnfWq3S9FLgbCPX347NAZwIrl3EKZm+WctHhNvCnA7fWX99NtbxmD1XA/C7w+nVe/xbgaTHGO0MIj6JaqvPXDT6/38zJrcJlb5Zy0V1/Rd+ZAz3njhUu6UQb8s0caJk7a066Yow/ARBC2E21B/2jQgjvBU6hxdkkSTPS3dzgFqoDjVuB7cD+EMJO4N8BjwKeDvwS8IZVXr83xngnQIzxX0IIjc5UmDlSrjoLnV4zp36duSPlprve9WkzB1rmzkQbacQYr53keZqt537ybX0PQbnqLoyuBq6oDz4+DLwLuCTG+MMAIYRrgXeu8fp7Qwh/BXwC+H7g5BDC2wFijL886SCSZ05uFa6+/NOb+h6BclV45tTPv7bp4Kdhhavy9dtvAuDsXef1PBLlpMMt4KfNHGiZO24ZLw1QV8sL68bRl6zx893rvMUNx3192yzGJCkHZo6kdLraMn4GmQMtc8dJlzREmRZmYozX9T2GRqxwrWx8HbFj+/odhzKS57+VwWUOVrhWc+DgXgC2btnR80iUhzwzB9rnjpMuaYC62kpVklZk5khKqMTjHCdd0hAVGEaScmbmSEqpvMxx0iUN0GhDeWEkKWedbd8sSSfocMv43jjpypC7FmpdVrqmYy/XUu5aqPWYOVOzl2up8a6F0ooKzBwnXdIAFZhFkrJm6EhKp8Mt43vjpKsHJ+34DgAO771/nWdKq+hoy/hizXtl6+ie6n68K6HUmJnTlJUtqb2utozvk5MuaYCsdElKytCRlFR5meOkq4HNO6sK1aG77p3o+dvCrkbPt5dLE/MAaH0lVLeO1BWqTadO9vwHvlw///zJnm8vlyZWwL+nBKxurc9eLk3CLeMlZcHdCyWlVd5SH0n5cnmhGjnp1KoytlqF7KQdDwfs7ZoH144uXPkH57R8wwLPAGkGju6t7lerkB3bV92v19u19YnV/T2ekR6sCxZX+UFo935mjjp2591Vpf7sLef1PBLlwI00JGWhxLK7pJyZOZISKvA4x0nXDJy0va5Y7WtXsbKXa3hWrVylUl4WqYljdUVr4452r7eXa3hWrVylYuhoOvZyqQkrXZKyMHLLeEkpFdhfISlf9nSplfEuhmPj3q4fvOHNfQxHLfRe2VquvBNAmqXxLoacW92Ne7vu+kgvw1ELvVe2ljN01M7ff/EaAE7b3rKfUHOqvMxx0iUNkD1dktIycySlU+JxjpOuBMa7GBJvB+D8a36mx9GoCG4Zr7WMdzEc2//JfsahchS41EfdGle4pDZcXigpCyWeAZKUMzNHUjpupKE1TbuLofKTXS/XmJMuwfS7GCo/2fVyjZk5klIqL3OcdEkD5JxLUlqGjqR0SlzR46Srhc07qx6t8S6Ek/qB3/qpLoajDmRb4Rpzy/j5cmRPdT/ehXBSx/au/Pi4Mrbaz5VethWuWoH9FerG337unQBs3bw0rw4crPJm6xYr81qfPV2SslDiGSBJOTNzJKVUXuY46QI273jEku8P72lWwVrPY17w1Jm+n1RgFs2XY/uWfr+wc7bvv+uV1f2B+npdR+5q9z4bt1f3y8erOWToaG3X/PmLADjtlOr6gMsrXZOyIiYo8+Syky5pgEYuL5SUUoFLfSTly+WFUuGy7+UaK+8EkLow7gF7oN9haA2593I9xNDRZA4cqnpQd3JOzyPRkLllvKQslFh2l5QzM0dSSuVljpOuKay3i6G9XMOze/F6YAAVLydd82m9XQzHvVxjC/XzNtU9Y217u9SdG+t/y7lXvMwcrWLcyzV24GCVU+OKV9veLs23Ek8uO+mSBmi0obwwkpSz8vorJOXLSZcmYoVLnSswjDSFx/xidT/pLoMLXqdLTZk5Wup3/uSZwOS7DI4rYO5KqMmUlzlOuqQhKi+LJOXMEz2SErLSNSdOOrXq1Wp6va5Hnf/ELoYjncAt4wtztO61anq9rse9Y/ZjkVZk5qgyrnBJXXLLeElZKPAEkKSsGTqSUiovc5x0HWfzjkcAcGjvfRM9f1s4Y8m9lIyzrjKMe7A2bp/s+U+5obp3F0KlZubMrXdcV63icRdCpTSX1+kKITwDeB1wP/AN4OHAScA24OIY4+FORyjpRAUfAJk5Uo7MHEnpzGtP1ynAq2KM94UQ/hr4eozx4hDCy4ELgQ92OkJJJ+hqy/gQwhnA5cA+4Esxxqvqx18DPJXqYOS/xxg/2skAKmbOrG2sdwtrWiHbWvep3nPTbMejAeqmv8LMKdcDB6tKftMK2Z13fxmAs7ecN/MxaTi66unqM3PW/Y1ijH8B3B9CuBS4Cfhm/aNvArtmPSBJExiNmt0mdzFwZYzxtcDzQggL9eP7Y4yvAF4DvGTGv80SZo6UITNHUlKjFreJ9JY5kywvfATwbqozPV8BLq1/dCZwexeDyp09XOpbh1X304Fb66/vplpesyfG+KEQwsOpzg5d1tmnY+asaNzLNbap3uVw0srVpvpM8wOzG5LmTWehY+ZkatzLNXbgUH2drQkrV+Pn7+Sc2Q5Mc6HD5YW9Zc4kywvfAzwBeDnwIHBbCOE99SAv6mJQktbW4Zbxt1AdaNwKbAf2A4QQzgV+BXhLjPHrXX14zcyRctPd9s1mjqQTjLq7TEVvmbPupKsutc2l5dfrssKlbHRX6boauCKEsBv4MPCuEMIbgY8BXwB+M4TwxRhjZxeImufMOeF6XcsrXOsZ74I43hVxbGGVM9PjXq9je5t9juZQZ6Fj5mRmeYVrPQcOVvmxdcuOZY/vafR8aYnuKl29ZY5bxktD1FEYxRjvYOW1zI/v5AMlDYSZIymdrraM7zNznHStYHy9rtW+l/pW4laqc21cmTrVnYHrAAAH9UlEQVS/PjN8dOUzxFNr2gs2tloFTXPEzCnVW95f5cLWLd1ch6tpL9hDr7MiNtdKPM5x0iUNUUdbxkvSirrr6ZKkE3S1ZXyfnHRJK9i9eD0A144u7HkkKyvxDJB6sGBPVzZurP9NX7DY7zhWZeZoNsa9XlawtLbyMsdJlzRETrokJWXmSEqnxJPLTrqkNeRa8Rq5vLAs4+rGUStOcy/XileBS33m3SXvfhgAJ2/e3vNIpBN1uGV8b5x0SUNU4BkgSTkzcyQlVOBxjpMuaYAKzKL5NG01o+1uhFJjhk4pxhWuttruRig10dWW8X1y0iUNkbMuSSmZOZISsqdLUhbs6VIrG+vdwqyMqbHy+iuUxgMHq+v7WRlTE24ZLykPBZ4BkpQzM0dSSuVljpMuaQ257Vo4VmLZfa50vTPdxno3smP7lj6+qT7T/EC3H68p5LZr4ZiZM3jT9nKt58DBavfV5dffGveA7eScTj9fZSnxOMdJlzRELi+UlFR5S30k5cst46XC5VrZWq7EM0Bzoe8qxoI9Fdnp+38TEzNzhqrrCtd6Dhzc0+vna6AKPM5x0iUNUXlZJClnBR4AScqXW8ZLysJoQ3ll9ySmrSrc2PL/BJp+7rgidTSTM8Rbn1jd33NTv+NQj8yctqatNF3+Cw8m+dxxRWrrljwq4nfe/WUAzt5yXs8jUR/cvVBSHso7ASQpa4aOpHRKbKNw0iUxnF6usRLDqFOz6psZv8+kFa9U/Tqbdlb3k15/q+nzNXuD6eWqmTmNzaqXavw+k1a8UvVwjXclnPT6W02fr3lXXuY46ZKGyAMgSUmZOZLSsadLGoihVa6aGrllfL+6qlIs1Ne3Obq3m/dfzcb6c48l/tySDK1y1Vh5/RVD01UF64FD1fX8Tt68vZP3X81q1/WSwJ4uSbmw0iUpJTNHUkoFZo6TLhWh9MrWCcrLImlYiq9sLWfoSErH5YWSsuCW8ZLSMnMkpePyQklZKLDq3o25q0Z0bGPd83FsX7/jUHqGzsRS7R44L+z9mk8l7tLspEsaogLDSFLOzBxJKZWXOU66NGhz18s15qRLs7Tg7oUTm9vqqZmj2TpwsL5ulxUsrcCeLklZcMt4SUkV2F8hKV/2dEnKg5Wu9c1tRULqgpkzCfu5pBkp8DjHSZc0QCU2mErKmZkjKR2XF0qZGHov1+7F6wG4LIRWr3d5oZTY0CunN44z45x2ry9wqY+kfLm8UFIerHRJSsrMkZSOlS6pZ0OvcM1MeVmklSycWt0f3dPvOObZ0CtcM2PozIt/3VXw1J5HorlW4MllJ13SANnTJSktM0dSOla6gBDCGcDlwD7gSzHGq2Y+KmkZK1zLzFFPl5mjXljhWqrA/orVmDlS/+zpqlwMXBlj/FQI4S9DCL8bYzw664FJWl1Xla7VDjZCCD8EvIzqdPf7Yoyf6mQAKzNzpN6ZOQk/X1KBxzltJl2nA7fWX98NbANsONBMWdlaR3fLC1c72HgD8CLgYcAfAy/sagArmL/MWdhR3R/dm+bzUn1OzqxsrcPMSfj5vXjg0D4ATt68Pcnn2TOmtXS4vLC3zGkz6boFOJMqkLYD+2c6Iknr6nDL+NUONkYxxiMAIYTNXX34KswcqW/dLfUxcySdoMPlhb1lTptJ19XAFSGE3cCHY4zHJnlRaHk9Is2plpeSGYq21+eqfSOE8Nimr5nweasdbBwKIWyiOgN0qOFnT8vMycapy+5LUvr/XqYKVTNnQh973+M7G1R6357weXeu8/3YzVOMRXOmTebAZLnTW+aMFhddUiGpEkI4HbgCuA/4HPAU4BLg6cArgQXg8hjj53obpKRimDmSUuozc5x0SZIkSVKHytuPUZIkSZIy4qRLkiRJkjrkpEuSJEmSOuSkS5IkSZI65KRLkiRJkjrkpEuSJEmSOtTm4siNhBDOAC4H9gFfijFe1fVnTiKE8AzgdcD9VBdTezhwEtWVqS+OMR7ucXgAhBD+CPgo8J3AWVRje32M8a6ex3UW8BaqKyDeBzySjP52IYQzgd+gutL4CPgWGf391L0cc8fMmWpcZ2HmKGNmTns55o6Zoy6kqHRdDFwZY3wt8LwQwkKCz5zEKcCrYoyvBM4HtsUYXwd8Ariw15EBIYQ3UAUlwLNijD8D/B5wUX+jesglwDeBRwP/QmZ/O+Bc4DnAY4A95Pf3U/dyzB0zpz0zR7kzc1rIOHfMHM1ciknX6cCt9dd3U83Cexdj/Avg/hDCpcBNVP+4qO939TYwIITwfGA/8Gmq/47urH/U+9hqjwc+QvUP+6Vk9Ler3UoVRj8OXEAVSJDP+NS97HLHzJmKmaPcmTkNZZ47Zo5mLsWk6xbgzPrr7VT/wHoXQngEcDXwGeAa4Iz6R2cCt/c1rtpLgacBP031D/60+vEcxgZwB3BvjPFo/X1OfzuAn6U6K7UI3AM8tn48l/Gpe9nljpkzFTNHuTNzmss5d8wczdxocXGx0w8IIZwOXEG1JvZzMcb/1ukHTiiEcA3wBKp1zg8CtwFbqc5OXXTcP7TehBB2A4eARwGBaqnAq2OM9/Q8ricCv051VurTwJPI6G8XQvg+qvHdAvwzcISM/n7qXo65Y+ZMNS4zR1kzc9rLMXfMHHWh80mXJEmSJM0zt4yXJEmSpA456ZIkSZKkDjnpkiRJkqQOOemSJEmSpA456ZIkSZKkDjnpkiRJkqQOOemSJEmSpA79f4fHk1V/NAGjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_probabilities(plot, prob_map, geo_model.grid.topography.resolution,\n",
    "                  T=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 0.9957274520849256 0.0 0.9957274520849256\n"
     ]
    }
   ],
   "source": [
    "ent_map = entropy(prob_map, base=2)\n",
    "ent_sections = entropy(prob_section, base=2)\n",
    "print(ent_sections.min(), ent_sections.max(), ent_map.min(), ent_map.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution=[80, 80]\n",
    "restopo=[100, 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Map entropy')"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAENCAYAAAAlhJ6dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFupJREFUeJzt3X+wXGV5wPFvYiCkJY1cCkqISB3oK1XbMZ1mbOpAxVqpjB0dZixl0EIhhl8tRYgFrM6INrXGYImFtqNYx1FLGwuMYhuk2GZaEaOmjg5knlqphnDFgIHIjxAIpH+cc25OLrv33t275+67u9/PDJOTs2f3vFly3zzned/3eeft378fSZKknMzvdwMkSZImM0CRJEnZMUCRJEnZMUCRJEnZMUCRJEnZMUCRJEnZMUBRWymlo1NKCxv8/OOa+mxJw8X+YvQs6HcD1LyU0s8DHwWeAhYD10fEf87grX8JXJVSehNwd0Tc2+Om/R3wOz3+TElzJKV0DvBh4CURsTeldDzwv8CyiHiwx7ezvxgxBiij4dXAD4D3AYcCv5tSegFFx7IAeCHwLuA04NeBY4ANwHLgHcDPAPemlN4OnFp+xqeAY4E3AtsoOqhV1Q1TSmcCK4HDgVuAI+vXAn8L/GIZ/FwA/BjYCLwNeLK8/kLgc8DW8p5fBt4AfAG4G7gVeEtEWG1Q6p9vAG8B/hE4F/gPgJTSRcArgaXAZ4GHKPqgLwDHAe+OiH3ltb8FvBV4AfBNiiDnT4H/An4FuIqivzgL+G3gGeAO4Nco+qex8vprgAeBJ4D7KPqfn0bEP6SU/gk4JyKebO6rUC85xDMCIuIO4NsUWZQNwCMUP+S/BOwBngVeA/x+RFwCrAK2A/8NfLr2UW+PiHOBc4DLy3P/FhHXAL8w6bbvpgg0dpb3OujaiPgW8L2I+BeKTuRSis7pnoj4Y+CfgbdT/B39AHA1cAnwMeCdwJuAfzU4kfru88Bby4eeFwE7yvPfogha7qb4eQX4SkT8FXA/8Bu1z7gKeAx4FHhdee6bEfHnwOMUAcn3IuJz5Wvrge8Ch0bExcA64I/K1/46Ij5A8bDzaeDMlNIrgO8bnAwWA5QRUGY+tkfEnwAXU3QG84FvRMSVFNmQ7wOHlG85jCLL0e4f//3AvPL4ifLXfZOueY4iqFgL3DXNtU+XHcf88n3V++eV/82nyPQ8FxHjFJ3VOzk4eJLUH3soMqCrgS/Vzq8Ffg7YwoH+ovo3ZxHwdO3aBcAHgfcAXynPVf3F0zz/36pHad1f1O8xPyL2UDycrQdu6PDPpT5ziGc03AlsSCk9RTFU8hngdooni+spnnrOBTamlP4GOAp4P/A/FMHMI+XnfCal9PHy+EPAS6e458fK+xwGXE8R8Ez2eErpHbXf307xJPYRio5tDUXq+MMUWZaPlNfdApzi05CUjRuBm4EEnFGeewh4PbCQ4ucZ4PSU0uEU//Z8vfb+tRQPSvuAL7a5x30ppSuq30TEPSmlZ1NKHwWWUAwffRB4X0rpmbJNADcBKSLun9WfUHNunpsFKmcppU0RcVrt979JMbx0XkTs7FvDJHWk/Nl9TUR8qMF7fAq4spqgm1J6JcVk/ysiYltT91UzDFAkSVJ2nIMiSZKyM6M5KCmlE4CNEfHqlNIa4HiKMb/LKOY0rAd2UazAuH7yNRHxUANtlzRA7EckdWLaDEpK6cXA+cATKaXDgJPLZV03UixHXQ1siIiLKCZALW5xjaQRZj8iqVPTZlDKyUZXppQ2URTDqSYm7qAowHMoxZp2KFZ7HNHimmmllJwMI2UmIuZNf9WMPsd+RBpBs+lDOl1mvJOiwifAMmCcIguzjKJzGSvPTb5mRo773i932Jz+277xVTO+9o0vaz2J/LQl3215fvnCXc87t3bnKdPeZ9uvTi4zInVu+4nfaeqj7UdmYXzNyonjfSsem/H7WvU/7fqeNyzaM3H8k+eK4617xwDYtHv6Pu/Om1ZMHC9dd9cUV2qYzbYP6ShAiYh9KaV/TyndQPGEcwFFwZ1ryz0Zbm5zjSQB9iOSZmbGAUpViyIirpv00m7grEnXTr5maC3YsnjiuJOnmZk4cv6iiePqKWYm6k9YPr0oJ/YjkmbKZcaSJCk7lrqXJPWUmVv1ggHKHLr9vpMmjttNmK27Y8+BIZ7lC4shnvqktnaT1V5/5paJ423rOm6mJEl95xCPJEnKjgGKJEnKjkM8s1Qfa+2kJkpdfaimXV2CakVPvTbKpq7uJkmd66T+yVKcg6LZM4MiSZKyYwalh5qsidIta6JIkgaRGRRJkpQdAxRJkpQdh3j6pNuaKFU9FICrj948cdxuE0FrokjqVCdba0hNMYMiSZKyY4AiSZKy4xBPD9VXyVSrZzpdzTNdTZRudziutwlc0SNJypsBiiSpJ3zwUS85xCNJkrJjBmWEuKJHkjQozKBIkqTsGKBIkqTsOMSTgXZF21qt6KkKtkHnRdskaSr13dKlfjODIkmSsmMGRZKG2Ey20qhUdZYsda8cGKA0pKoHUC+O1mnRtun0omibdQskSTkyQJEkzcqdN60AYCk+8Kh3nIMiSZKyYwZlgLUb1qnv4VNfCVRXFW2zYJskKUcGKJI0QlptQjqVrXvHGmqJNDUDlIbVJ6Fu39g6m1E3XU2UeudS7zjq9Qvqx5s6bK8kSTkwQJEktdVumFhqmpNkJUlSdsygzKEFWxZPHOdSE6Vep8WaKJKkXJhBkSRJ2TFAkSRJ2XGIJ2PtVvS00m6X4+lqolT1UMCaKNKocw8e5cQMiiRJyo4ZFEkacfXaSd1wgr2aYIAyhzot2tZKfZimXUXI+ooei7ZJkgaRQzySNOKOnL/ooAcbKQcGKJIkKTsO8fRJk0XbumXRNklSLjoOUFJKy4BrgEeAecCPgOOBJcBlwKHAemAXcE9EXN+rxkoaDvYjkqbTzRDPy4FTgeOAh4GTI+Ji4EZgFbAa2BARFwGnp5QO6VVjJQ0N+xFJU+pmiOd+io7l/4AvAzvK8zuApRRPPveX5x6heCJ6eHbNVFW0bbqCbdC+aNvVR2+eOF6785Tnvc+ibZpD9iNdymVIWGpaNxmUS4AlEbEf2A28tDy/DBgHtpfHAGPAo7NtpKShYz8iaUrdZFD+HrgmpbQd+DrwdErpBuAI4AJgEXBtSukc4OaI2NerxkoaGvYjkqbUcYASEVuBN09xyW7grK5bNCLqq2Tqq2c6Sd92WrStk302XNGjJtmPSJqOy4wlaYDkWKJAaoKF2iRJUnbMoKglV/RIkvrJAEWSRlwn89OkuWKAMmCqeihwcE2UdhNmp6uJ0qoeiiRJ/eYcFEmSlB0DFEmSlB2HeDLQi5ooM2FNFEmtbN07BsDyhbs6et+dN60AYCn2C+o9MyiSJCk7BiiSJCk7DvGMkFbDOvUVP/WVQHXWRJEkzTUzKJIkKTsGKJIkKTsO8WSmvkpm+8bWQy6VTou2VTP14cBs/fqs/U1dtFeSpCaYQZEkSdkxQJEkSdlxiCdjC7Ysnjhuqmhbp5uEWbRNkjQXzKBIkqTsGKBIkqTsOMQzJNqt6Gnnjj3FEM/yhQeGeCzaJo22bod+pSaYQZEkSdkxQJEkSdlxiCdjnRRta6dd0bZKldIFi7ZJkvJhgCJJApx7orwYoAyIqiZKr+uhzIY1USRJTXEOiiRJyo4BiiRJyo5DPEOok5ooVT0UOLgmytVHb544XrvzlJbvtSaKJKkpZlAkSVJ2DFAkSVJ2HOIZENUqmfrKmU5X9HRSE6XbXY5dzSNJ6gUzKJIkKTtmUCRpANSzp71WryIt5cIART1RrehxNY8kqRcc4pEkSdkxgyJJA6SJ7S6qCfLV5Pite8d6fg+pUwYoQ65d0bZqRU99Nc9sirZJktRLDvFIkqTsGKBIkqTsdDzEk1I6HngvsBN4DHghsBBYAqwGEnAV8CRwW0Tc0qvG6uBCaLMp2jadbou21dtk0Ta1Yz/SrOn24JIGQTcZlMuBHcAxwI+BJRFxKbAZOAO4Arg8Is4DLuxVQyUNFfsRSVPqJkA5AbgVWAWcTdHJUP66FDgqIsbLc/tn3UJJw8h+RNKUulnF8yDw04h4JqUEcGx5fhkwDjyQUjomIn4EzOtNMzXX2g3r1Ff91Pf2qVQF28CibZqS/cgca7X/VjutfraludZNgPJh4C9SSjuBTwCvSCldRzF2vAr4NrA+pfQ0sKFnLZU0TOxHJE2p4wAlIrYBb5viknuBs7pukWasPgl1+8bpn3ha1URpt8NxvVBTfZ+O+vGmDtsrVexHZq6aeN5EgbbZciK8muQyY0mSlB0DFEmSlB1L3Q+JBVsWTxxbE0WSNOgMUCRpxHXy4CHNFYd4JElSdsygjKhqRc9MSmK32+XYmiiSpKaYQZEkSdkxQJEkSdlxiGdIdFq0rZV2Rdvq6it6LNom5aXTXYzrP8NSbsygSJKk7JhBkaSMzbau0VSbBFYZUZcZK0cGKEOoyaJt3bJomySpEw7xSJKk7JhBkaTM1DOOkEcWVJprBigjrirYBgevAGi3oqdd0barj94MwNqdp7S8j0XbJEmdcIhHkiRlxwBFkiRlxyGeIVRfJVMfy+71ip560bZOlim6okeSNB0zKJIkKTsGKJIkKTsO8ain6it+6iuB6lzRI/VWp3vwtLN171hPPkfqBTMokiQpO2ZQNKEXNVHc4VgaHJMnt7fLekr9YAZFkiRlxwyKJGWmyR2M61lOKWcGKENurmuidLptuzVRJEmtGKBIUiYOBOzNbRDYbYFFaa45B0WSJGXHDIpaareip51WT2LWRJEkdcsARZIGVK8KtEk5cohHkiRlxwzKCKmvktm+sbuCTO2KtlUlsutLGC3aJknqlhkUSZKUHTMokjQkpirQJg0aA5QRtWDLYqC5gm1g0TZJUvcc4pEkSdkxQJEkSdlxiEfT6qRo2x17DgzxLF94YIjn6qM3Txyv3XlKy/datE2SVDFAkaQ+q+Zf9XpOWF217P8nzzV2C6mnHOKRJEnZ6TqDklL6LPBF4CXA8cAS4DLgUGA9sAu4JyKun30z1WvVKpluC7ZB+6JtFXdN1XTsR+ZO9fPY6mexKrQo5aSrDEpK6V3A4+VvT46Ii4EbgVXAamBDRFwEnJ5SOqQnLZU0VOxHJE2l4wxKSunNwKPA1ygCnJ3lSzuApRRPPveX5x6heCJ6eNYtVSOqeijQ7Ph3p6yJMtzsR7rnBoEaFd0M8ZxN0WGk8vfVE9AyYJyis1lG0bmMUXRCklRnPyJpSh0HKBHxewAppXOAp4AXpZRuAI4ALgAWAdeWr98cEft61lpJQ8F+RNJ0up4kGxGfavPSbuCsbj9XeWtXE6WaMFufLGtNFE3HfkRSOy4zliRJ2TFAkSR1ZXzNyoMmtEu9ZCXZEVdfIVPvaHLc5djVPJI0OsygSJKk7BigSJKk7DjEo76qr/qpl86vq1b0uJpHkkaHGRRJkpQdAxRJkpQdh3g0odNVMq12BNnGSS3Owkdn9IkWC5X6YfnCXQBs6nM7pDozKJIkKTtmUCRpgLTbbmI6ndYfkvrNDIokScqOGRRJGgFb944B8IZFBzIpk7Mq1bL/dkv+J7MEgJpkBkWSJGXHAEWSJGXHAEWShsSm3a+a8fCMlDsDFEmSlB0nyUrSCKlPjD1y/qKDzlmwTTkxgyJJkrJjBkWSRki13BgOXnIs5cYMiiRJyo4BiiRJyo5DPJI0ZOpLjavqsNKgMYMiSZKyYwZFkgZUtbNxJ7saz0Q962LhN/WLGRRJkpQdMyiS1GdL190FwPialQDsW/FYP5vTsard1Z9D6gUzKJIkKTtmUCRpRNXL3ku5MYMiSZKyY4AiSZKy4xCPJA24bpcbV/vyVHvyOOSjnJhBkSRJ2TGDIklDrCq01knJ++ULdx14f89bJM2MGRRJkpQdMyiSlJkFWxZPHA9a0TapV8ygSJKk7BigSJKk7BigSFImlq67y/1spJIBiiRJyk7Hk2RTSiuBS4HHgR8ChwMLgSXAaiABVwFPArdFxC09a62koWA/Mlxef+YWALat63NDNFS6yaAcAZwfEecBrwWWRMSlwGbgDOAK4PLy9Qt71lJJw8R+RNKUOg5QIuJLwOMppfcAXwV2lC/tAJYCR0XEeHluf09aKWmo2I804/b7Tpooe9+Jnzy3xzL3yk7HAUpKaTHwCeBu4JPAseVLy4Bx4IGU0jHluXm9aKSk4WI/Imk63RRquw44ETgXeJaiI7mOYux4FfBtYH1K6WlgQ68aKmmo2I/MUFW0bbYF26YqeV9tGlgvcV+prq/eL82VjgOUiPjDaS65Fziru+ZIGgX2I5Km4zJjSZKUHffikaTM1Iu1ja9ZOWf3PXL+IoCuJ8xWbbXYnHrBDIokScqOGRRJUlvVxNlNfW6HRo8ZFEmSlB0DFEmSlB0DFEmSlB0DFEmSlB0nyUrSkKnvx/PGl23rY0uk7plBkSRJ2TFAkSRJ2TFAkSRJ2XEOiiSNkPquxJN3Nu62xL3UBDMokiQpOwYokiQpOwYokiQAtu4dY+vesX43QwIMUCRJUoacJCtJA2DBlsUA7FvxWJ9b0t7rz9wCwLZ1fW6IhoIZFEmSlB0zKJI0xKqy952UvD9y/qKJY5ceq1/MoEiSpOwYoEhSxpauu4ul6+7qdzM4bcl3n1fYTWqSAYokScqOAYokScqOAYokScqOAYokScqOy4wlST01vmblxHEOE3w1mMygSJKk7JhBkSRNa/nCXQBs6nM7NDrMoEiSpOwYoEiSpOwYoEiSpOwYoEiSpOwYoEiSpOwYoEiSpOy4zFiSBsiCLYsnjveteKyPLZGaZQZFkiRlxwBFkiRlxwBFkiRlp5E5KCmlY4H1wC7gnoi4von7SBpe9iPSaGsqg7Ia2BARFwGnp5QOaeg+koaX/UjN0nV39XVn4K17x9i6d6xv99foaSpAeTFwf3n8CLCkoftIGl72I9IIa2qZ8XZgGUXnMgY8OqM3nfidhpojaQDZj7Rya+3Pd+vM33b7DM8V6v80VMuanyl/3TqDux24ZvuJM7hcamHe/v37e/6hKaUXA9cCjwHfjIiP9/wmkoaa/Yg02hoJUCRJkmbDZcaSJCk7BiiSJCk7BiiSJCk7BiiSJCk7BiiSJCk7BiiSJCk7TRVqm7FB2G8jpbQSuBR4HPghcDiwkKKy5WogAVcBTwK3RcQtfWpqSymlzwJfBF4CHE/R7suAQ8nwu08pHQ+8F9hJUQPjhQzA951SWgZcQ1H1dB7wIzL+vlNKJwAbI+LVKaU1TNPWyddExEN9aXgL9iPNsx9p3qD1IdBsP5JDBmUQ9ts4Ajg/Is4DXgssiYhLgc3AGcAVwOXl6xf2r5nPl1J6F0WHCHByRFwM3AisIt/v/nJgB3AM8GMG5/t+OXAqcBzwMBl/32URtPOBJ1JKh82grYtbXJOTbL7bKdiPzK1B7EcGpg+B5vuRHAKU7PfbiIgvAY+nlN4DfJXiLz3lr0uBoyJivDyXTeW7lNKbKcqDf43i//XO8qWq3bl+9ydQFPFeBZzNgHzfFN/lqcDbgNdRdDCQ4fcdEQ9GxJUU/+iMMf3fjSNaXJOTbL7bduxH5twg9iMD04dA8/1IDgFKtd8GdLDfxlwqo75PAHcDnwSOLV9aBowDD6SUjinPzZv7FrZ1NrAC+AOKH9Kjy/NVu3P97h8EfhoR1eYfg/J9X0LxlLYf2A28tDyf+/e9EziyPG7X1vEW1+Qk1+92gv3InBvEfmRQ+xBooB/pe6n7QdhvI6X0SeBEinHjZ4EHgJ+liF5Xla/9GfA08PmIuK1PTW0ppXQO8BTwIopx1yOAC4BFZPjdp5ROAt5P8Rf+a8ArGIDvO6W0nKLd24EfULQv6+87pbQpIk5LKV3KNG2dfE1E7O5XuyezH2me/UjzBrEPgeb6kb4HKJIkSZPlMMQjSZJ0EAMUSZKUHQMUSZKUHQMUSZKUHQMUSZKUHQMUSZKUHQMUSZKUnf8Hs+jWJjCD1r4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(9,9))\n",
    "\n",
    "vmin = 0\n",
    "vmax = 1\n",
    "\n",
    "section_ie1 = entropy(prob_section,base=2)[:80*80] # indexing in case there is more than one section\n",
    "map_ie = entropy(prob_map,base=2)\n",
    "\n",
    "extent=[0,1000,0,1000]\n",
    "im = ax[0].imshow(section_ie1.reshape(resolution).T, origin='lower', cmap='viridis',extent=extent, vmin=vmin, vmax=vmax)\n",
    "ax[0].set_title('Section entropy')\n",
    "#add_colorbar(im, ax[0], label='cell entropy')\n",
    "im = ax[1].imshow(map_ie.reshape(restopo), origin='lower', cmap='viridis',extent=extent, vmin=vmin, vmax=vmax)\n",
    "ax[1].set_title('Map entropy')"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}