From d0292ce28d9c74d7b6c092491d49b4d2520e8a2b Mon Sep 17 00:00:00 2001 From: "Kasper D. Fischer" Date: Thu, 20 Nov 2025 12:35:20 +0100 Subject: [PATCH] Updates model uncertainty visualization in Jupyter notebook Refines presentation of model uncertainty by adding error bars to the resistivity and thickness values in the output plot. Integrates a new table displaying layer resistivity, uncertainty, and thickness with corresponding uncertainties to enhance clarity for users. Updates environment requirements to include `pandas` for improved data manipulation capabilities. --- README.md | 64 +++++++++- VES.ipynb | 304 ++++++++++++++++++++++++++++++------------------ environment.yml | 1 + 3 files changed, 253 insertions(+), 116 deletions(-) diff --git a/README.md b/README.md index 1b61bce..fe37113 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,63 @@ -[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/git/https%3A%2F%2Fgit.geophysik.ruhr-uni-bochum.de%2Fkasper%2FresistivityVES.git/HEAD?urlpath=%2Fdoc%2Ftree%2FVES.ipynb) - # resistivityVES -Binder environment using pyGIMLI (https://www.pygimli.org/) to do an VES inversion. \ No newline at end of file +[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/git/https%3A%2F%2Fgit.geophysik.ruhr-uni-bochum.de%2Fkasper%2FresistivityVES.git/HEAD?urlpath=%2Fdoc%2Ftree%2FVES.ipynb) + +Binder environment using [pyGIMLI](https://www.pygimli.org/) to do a VES (Vertical Electrical Sounding) inversion. + +## Overview + +This repository contains a Jupyter notebook (`VES.ipynb`) that demonstrates how to perform 1D DC resistivity inversion using pyGIMLI's built-in VES forward operator. The notebook uses real field data from Bausenberg to invert for a layered earth model. + +## Notebook Contents + +The `VES.ipynb` notebook includes the following workflow: + +### 1. Setup and Imports + +- Imports necessary libraries: `numpy`, `matplotlib`, `pygimli`, and the `VESManager` from `pygimli.physics` + +### 2. Field Data from Bausenberg + +- Uses real VES measurements with AB/2 distances ranging from 1.0 to 100.0 meters +- Apparent resistivity values (`rhoa`) showing variations from ~64 to ~672 Ωm +- Error estimates set at 2% for most measurements, increasing to 5% for deeper soundings +- MN/2 spacing fixed at 0.5 meters + +### 3. Inversion Setup + +- Configures a 3-layer earth model (`nlay=3`) +- Uses regularization parameter `lam=1000` with a reduction factor of 0.8 +- Inverts the apparent resistivity data to determine layer thicknesses and resistivities + +### 4. Visualization + +The notebook provides comprehensive visualization including: + +- **Model plot**: Displays the inverted resistivity model as a function of depth (up to 50m) +- **Data fit plot**: Compares measured apparent resistivity data with the model response +- Both plots are displayed side-by-side for easy comparison + +### 5. Uncertainty Analysis + +- Computes model covariance matrix to assess parameter uncertainties +- Displays correlation matrix showing interdependencies between layer parameters +- Generates error bars for both resistivities and layer thicknesses +- Visualizes uncertainties at layer midpoints and boundaries + +## Requirements + +See `environment.yml` for the complete list of dependencies. Main requirements: + +- pyGIMLI>=1.5.0, which requires at least: + - numpy + - matplotlib + - suitesparse=5.10.1 +- jupyterlab + +## Usage + +Click the Binder badge above to launch an interactive session, or run locally with: + +```bash +jupyter notebook VES.ipynb +``` diff --git a/VES.ipynb b/VES.ipynb index 006e08b..61a1a4e 100644 --- a/VES.ipynb +++ b/VES.ipynb @@ -1,19 +1,5 @@ { "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - } - }, - "outputs": [], - "source": [ - "# Checkout www.pygimli.org for more examples" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -21,21 +7,29 @@ "\n", "# VES inversion for a blocky model\n", "\n", + "## PyGIMLI\n", + "\n", + "This notebook uses the [pyGIMLI](https://www.pygimli.org/) library for geophysical modelling and inversion. Check out the [pyGIMLI documentation](https://www.pygimli.org/) for installation instructions and further tutorials and examples.\n", + "\n", + "\n", + "## Introduction\n", + "\n", "This tutorial shows how an built-in forward operator is used for inversion.\n", - "A DC 1D (VES) modelling is used to generate data, noisify and invert them.\n" + "A DC 1D (VES) modelling is used to invert meassured data them.\n", + "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We import numpy, matplotlib and the 1D plotting function\n", + "We import numpy, matplotlib, pandas, pygimli and the VESmanager\n", "\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": { "collapsed": false, "jupyter": { @@ -46,6 +40,7 @@ "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", "import pygimli as pg\n", "from pygimli.physics import VESManager" ] @@ -54,70 +49,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "some definitions before (model, data and error)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - } - }, - "outputs": [], - "source": [ - "#ab2 = np.logspace(-0.5, 2.5, 40) # AB/2 distance (current electrodes)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "define a synthetic model and do a forward simulatin including noise\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - } - }, - "outputs": [], - "source": [ - "#synres = [100., 500., 30.] # synthetic resistivity\n", - "#synthk = [0.5, 3.5] # synthetic thickness (nlay-th layer is infinite)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "the forward operator can be called by f.response(model) or simply f(model)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - } - }, - "outputs": [], - "source": [ - "#synthModel = synthk + synres # concatenate thickness and resistivity\n", - "#ves = VESManager()\n", - "#rhoa, err = ves.simulate(synthModel, ab2=ab2, mn2=0.5,\n", - "# noiseLevel=0.03, seed=1337)" + "Here we use data measured at a real site (at Bausenberg in the Eifel, Germany).\n", + "\n", + "- ab2 is the half current electrode spacing $ab/2$.\n", + "- rhoa is the apparent resistivity $\\rho_a$ (Ohm.m), calculated from measured resistances $R$ and geometry factor $K$ as $\\rho_a = R \\times K$.\n", + "- err is the estimated error of apparent resistivity values (here 2% of $\\rho_a$ for most values and 5% for the values at larger $ab/2$).\n", + "- mn2 is the potential electrode spacing $mn/2$ (here fixed to 0.5 m).\n", + "\n", + "Replace these data with your own measurements to perform a VES inversion for your site of interest." ] }, { @@ -134,37 +73,103 @@ "# use data from Bausenberg\n", "ab2 = [1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0, 12.5, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 75.0, 100.]\n", "rhoa = [64.0, 70., 77., 89., 103., 115., 152., 179., 206., 244., 305., 360., 410., 489., 540., 588., 672., 605., 535., 433., 283.]\n", - "err = np.ones(len(rhoa))*0.02\n", - "err[-4:] = 0.05\n", + "err = np.ones(len(rhoa))*0.03\n", + "err[-4:] = 0.07\n", "mn2 = np.ones(len(rhoa))*0.5" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " AB/2 (m) ρₐ (Ωm) δρₐ (Ωm) Error (%)\n", + " 1.00 64.0 1.92 3.0\n", + " 1.25 70.0 2.10 3.0\n", + " 1.50 77.0 2.31 3.0\n", + " 2.00 89.0 2.67 3.0\n", + " 2.50 103.0 3.09 3.0\n", + " 3.00 115.0 3.45 3.0\n", + " 4.00 152.0 4.56 3.0\n", + " 5.00 179.0 5.37 3.0\n", + " 6.00 206.0 6.18 3.0\n", + " 7.50 244.0 7.32 3.0\n", + " 10.00 305.0 9.15 3.0\n", + " 12.50 360.0 10.80 3.0\n", + " 15.00 410.0 12.30 3.0\n", + " 20.00 489.0 14.67 3.0\n", + " 25.00 540.0 16.20 3.0\n", + " 30.00 588.0 17.64 3.0\n", + " 40.00 672.0 20.16 3.0\n", + " 50.00 605.0 42.35 7.0\n", + " 60.00 535.0 37.45 7.0\n", + " 75.00 433.0 30.31 7.0\n", + " 100.00 283.0 19.81 7.0\n" + ] + } + ], "source": [ - "print(ab2)\n", - "print(rhoa)\n", - "print(err)" + "# Create a table with ab2, rhoa, and err\n", + "data_table = pd.DataFrame({\n", + " 'AB/2 (m)': ab2,\n", + " 'ρₐ (Ωm)': rhoa,\n", + " 'δρₐ (Ωm)': rhoa*err,\n", + " 'Error (%)': err*100\n", + "})\n", + "\n", + "print(data_table.round(2).to_string(index=False))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "set up VES manager and perform inversion, you can change the number of layers `nlay` as desired.\n", + "\n", + "Finally, we print a table with the data used for inversion." ] }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - } - }, + "execution_count": 40, + "metadata": {}, "outputs": [], "source": [ - "nlay = 3\n", + "nlay=3\n", "ves = VESManager()\n", - "ves.invert(rhoa, err, ab2=ab2, mn2=mn2,\n", - " nLayers=nlay, lam=1000, lambdaFactor=0.8, cTpye=2, verbose=False)" + "result =ves.invert(rhoa, err, ab2=ab2, mn2=mn2, nLayers=nlay,\n", + " lam=1000, lambdaFactor=0.8, cType=2, verbose=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Layer Thickness (m) Resistivity (Ωm)\n", + " 1 1.4 62\n", + " 2 31.49 945\n", + " 3 ∞ 61\n" + ] + } + ], + "source": [ + "data_table_result = pd.DataFrame({\n", + " 'Layer': np.arange(1, nlay+1),\n", + " 'Thickness (m)': list(result[0:nlay-1]) + [np.inf],\n", + " 'Resistivity (Ωm)': result[nlay-1:nlay*2-1]\n", + "})\n", + "data_table_result['Resistivity (Ωm)'] = data_table_result['Resistivity (Ωm)'].round(0).astype(int)\n", + "data_table_result['Thickness (m)'] = data_table_result['Thickness (m)'].round(2).replace({np.inf: '∞'})\n", + "print(data_table_result.to_string(index=False))" ] }, { @@ -177,18 +182,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 51, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArYAAAIdCAYAAADI2D7BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB8CElEQVR4nO3deVxU1f/H8dedYXMBFBTEnaxUxF0szVxyt1zavm1ulZVGltlmK9m+72hqmpWV/iqXLHOpNDU31NAMLTXcMRSUVbaZ+/uDpAhUQOAC834+HvNw5t4zd96DzOHD4dxzDdM0TUREREREKjmb1QFEREREREqDClsRERERqRJU2IqIiIhIlaDCVkRERESqBBW2IiIiIlIlqLAVERERkSpBha2IiIiIVAkqbEVERESkSlBhKyIiIiJVggpbEREREakSVNiKiIiISJWgwlZERPjmm29o3rw5F110ER988IHVcURESsQwTdO0OoSIiFgnJyeHkJAQVq5ciY+PDx06dGDjxo34+flZHU1EpFg0Yisi4uI2bdpEq1ataNCgAd7e3gwaNIhly5ZZHUtEpNhU2IqIVHKrV69m8ODB1K9fH8MwWLhwYYE2U6ZMITg4GC8vLzp27MiaNWvy9h05coQGDRrkPW7YsCGHDx8uj+giIqXKzeoAVnM6nRw5cgRvb28Mw7A6johUMaZpkpKSQv369bHZymYsIS0tjbZt23Lrrbdy7bXXFtg/b948JkyYwJQpU7jsssuYNm0aAwcOJCYmhsaNG1PYjLSz9YeZmZlkZmbmPXY6nSQmJuLv769+VETKRJH7UtPFHTx40AR000033cr0dvDgwXLp0wBzwYIF+bZ17tzZHDt2bL5tLVq0MCdNmmSapmn+/PPP5rBhw/L23Xvvveann356xteIiIiw/Oupm266uebtXH2py4/Yent7A3Dw4EF8fHwsTlO2srOzWb58Of369cPd3d3qOCJVxtk+W8nJyTRq1CivrylvWVlZbNmyhUmTJuXb3q9fP9atWwdA586d2bFjB4cPH8bHx4clS5bw1FNPnfGYjz76KBMnTsx7nJSUROPGjYmNjT3j+9y8/wR3fPJLge3dbdFM8XiX/c66XJn1EjNGtKdTk9oleavlKjs7m5UrV9KrVy/1pyLlICUlheDg4HP2pS5f2J7+s5mPj49LFLbVq1fHx8dHHbFIKSrKZ8uqP9EfP34ch8NBYGBgvu2BgYEcPXoUADc3N15//XV69eqF0+nk4Ycfxt/f/4zH9PT0xNPTs8B2Pz+/M/ajfWr70WDZfo4mZWD+a3s07ajuYae1cZyOvjn0adcMu63iT2c4/X/u7++v/lSkHJz+nJ2rL9XJYyIiLuC/PwxM08y3bciQIfzxxx/s2bOHO++8s9Rf324ziBgckpvlX9tTqU60eSEAz7ROqBRFrYhUXCpsRUSqsDp16mC32/NGZ0+Lj48vMIpb1gaEBjF1eAfq+Xrl276J1gC0ziw4VUFEpDhU2IqIVGEeHh507NiRFStW5Nu+YsUKunbtWu55BoQGsfaRK/j8jkuZ2PciAH525I7kErsadM0gETkPLjvHNjIyksjISBwOh9VRRETOS2pqKnv27Ml7HBsbS3R0NH5+fjRu3JiJEycyYsQIOnXqRJcuXZg+fToHDhxg7NixluS12wy6NPOnSzN/Vv5+jKgDF5Lt6YV7WjzE74TAEEtyiRTG4XCQnZ1tdYwqz93dHbvdft7HcdnCNjw8nPDwcJKTk/H19bU6johIiW3evJlevXrlPT69YsGoUaOYPXs2N9xwAwkJCTzzzDPExcURGhrKkiVLaNKkiVWR89x2WTDjD5wkymxBV6Lhz1UqbKVCME2To0ePcvLkSaujuIxatWpRr1698zrZ1mULWxGRqqJnz56FXmTh3+6++27uvvvuckpUdANC61HPx4tVaS3p6h4NsT9Bl4qXU1zP6aI2ICCA6tWr6+IjZcg0TdLT04mPjwcgKCioxMdSYSsiIpZxt9sY2bUJ3y5rDXyOue9nDEc22LWElljH4XDkFbVnW/pOSk+1atWA3BNbAwICSjwtQSePiYiIpW4Ka8yfbk1JNGtiZKXA4a1WRxIXd3pObfXq1S1O4lpOf73PZ06zClsREbFU7RoeDGvfmPXO06sj/GRtIJG/afpB+SqNr7cKWxERsdxtlzXlZ2fuerYZf/xocRoRqaxU2IqIiOUuCvQms1E3ANyPREFWmsWJRKqWnj17MmHCBKtjlDkVtiIiUiFc1eMyDpl1sJs5nNr7s9VxREqFw2myfm8Ci6IPs35vAg5nxb8IyapVqzAMo1IudaZVEUREpELo0TyApe7taJjzPXs3fktoy35WRxI5L0t3xDF5cQxxSRl524J8vYgYHMKA0JIvaSVnphFbERGpEGw2A++WvQHwOLAaZyUY2RI5k6U74hg3Z2u+ohbgaFIG4+ZsZemOuDJ77bS0NEaOHEnNmjUJCgri9ddfz7d/zpw5dOrUCW9vb+rVq8fNN9+ct4bsvn378i74Urt2bQzDYPTo0bnvaelSunXrRq1atfD39+eqq65i7969ZfY+SkKFrYiIVBgdew4F4EJHLGt//d3iNCL/ME2T9KycIt1SMrKJ+Po3CvvV7PS2p7+OISUju0jHO9cFWP7roYceYuXKlSxYsIDly5ezatUqtmzZkrc/KyuLZ599lm3btrFw4UJiY2PzitdGjRrx1VdfAfD7778TFxfH22+/DeQWzBMnTiQqKooffvgBm83G1VdfjdPpLO6Xs8xoKoKIiFQY1f0bcKzaBdQ99Sc/r1jACXMkAd5edA72w27T0ktinVPZDkKeWlYqxzKBo8kZtH56eZHaxzzTn+oeRSvZUlNTmTlzJh9//DF9+/YF4KOPPqJhw4Z5bW677ba8+xdccAHvvPMOnTt3JjU1lZo1a+Ln5wdAQEAAtWrVymt77bXX5nutmTNnEhAQQExMDKGhoUXKV9Y0YltSppl71m5WWu59EREpFQkBXQBodDKK++ZGc9OMDXR7+ccy/dOtSFWxd+9esrKy6NKlS942Pz8/mjdvnvf4l19+YejQoTRp0gRvb2969uwJwIEDB8557JtvvpkLLrgAHx8fgoODi/S88qQR25LKTocX6ufef+wIeNSwNo+ISBWwdEccX+6uxwcecJltR9720/MSpw7voJNuxBLV3O3EPNO/SG03xSYy+sOoc7abfWsYnYP9ivTaRXWuaQtpaWn069ePfv36MWfOHOrWrcuBAwfo378/WVlZZ33u4MGDadSoETNmzKB+/fo4nU5CQ0PP+bzy5LKFbWRkJJGRkTgcDqujiIgIucsiTV4cQ4qzJTmmjWDbX9TnOEeogwkYwOTFMfQNqadpCVLuDMMo8nSAyy+qS5CvF0eTMgqdZ2sA9Xy9uPyiuqX+vXzhhRfi7u7Ohg0baNy4MQAnTpzgjz/+oEePHuzatYvjx4/z0ksv0ahRIwA2b96c7xgeHh4A+WqkhIQEdu7cybRp07j88ssBWLt2balmLw0uOxUhPDycmJgYoqLO/RuViIiUvU2xicQlZZBKdbaZzQDoZ//nB64JxCVlsCk20aKEIkVjtxlEDM69RPR/y9bTjyMGh5TJL2g1a9bk9ttv56GHHuKHH35gx44djB49Gpstt+Rr3LgxHh4evPvuu/z55598/fXXPPvss/mO0aRJEwzD4JtvvuHYsWOkpqZSu3Zt/P39mT59Onv27OHHH39k4sSJpZ7/fLlsYSsiIhVLfMo/yyItcnQF4Cb7j/CfMa9/txOpqAaEBjF1eAfq+Xrl217P16vMp9S8+uqrdO/enSFDhtCnTx+6detGx44dAahbty6zZ8/miy++ICQkhJdeeonXXnst3/MbNGjA5MmTmTRpEoGBgdxzzz3YbDbmzp3Lli1bCA0N5f777+fVV18ts/dQUi47FUFERCqWAO9/CoAFjsuZ5DaX5rZDdDJ+Z7PZotB2IhXZgNAg+obUY1NsIvEpGeW2wkfNmjX55JNP+OSTT/K2PfTQQ3n3b7rpJm666aZ8z/nv3Nwnn3ySJ598Mt+2Pn36EBMTc9bnWU0jtiIiUiF0DvYjyNcLA0ihOosduWd13+z2Y16bIF+vIp1sI1JR2G0GXZr5M7RdA7o089f88DKmwlZERCqE/85L/NSRexWyK20bqUUKAE9eWTbzEkWkalBhKyIiFca/5yVuNy9gh7MpnkY219lXA5CWlWNxQhGpyFTYiohIhTIgNIi1j1zB53d0IafDaADu8VkLmLy67HfSMlXcikjhVNiKiEiFc3peYruBY8CjJrVO7WdYrT+JT8nk/Z/2Wh1PRCooFbYiIlJxeXpD6+sBeKTOegCmr/6TwydPWZlKRCooFbYiIlKxdboVgHpHVtC3sY3MHCcvf7fL4lAiUhGpsBURkYotqC006IjhzObZJtEYBny97QhbD5ywOpmIVDAqbEVEpOLrdBsA9fbM5X8d6gPwzOKYCrc4vIhYy2UL28jISEJCQggLC7M6ioiInEura8DTF07s49Hmf1Hdw070wZN8ve2I1clEXFLPnj2ZMGFCkdvPnj2bWrVqlVme01y2sA0PDycmJoaoqCiro4iIyLl4VIe2NwJQK+YTwntdCMDL3+3iVJbDymQiUoG4bGErIiKVzN8nkfH7d9ze1osGtapxJCmD6av3sn5vAouiD7N+bwIOp6YniLgqFbYiIlI5BLSERpeC6cDr1895ZGALAN78fjc3zdjAfXOjuWnGBrq9/CNLd8RZHFbEGj179mT8+PFMmDCB2rVrExgYyPTp00lLS+PWW2/F29ubZs2a8d133+U956effqJz5854enoSFBTEpEmTyMn550IoaWlpjBw5kpo1axIUFMTrr79e4HWzsrJ4+OGHadCgATVq1OCSSy5h1apV5fGW81FhKyIilcffJ5GxZTbuhrPQJkeTMhg3Z6uKWyldpglZadbcinmS5EcffUSdOnXYtGkT48ePZ9y4cVx//fV07dqVrVu30r9/f0aMGEF6ejqHDx9m0KBBhIWFsW3bNqZOncrMmTN57rnn8o730EMPsXLlShYsWMDy5ctZtWoVW7Zsyfeat956Kz///DNz585l+/btXH/99QwYMIDdu3eXype/qNzK9dVERETOR8hQWPoIJB/i+68/BUILNDEBA5i8OIa+IfWw24zyTilVUXY6vFDfmtd+7Ah41Chy87Zt2/LEE08A8Oijj/LSSy9Rp04d7rjjDgCeeuoppk6dyvbt21m8eDGNGjXivffewzAMWrRowZEjR3jkkUd46qmnSE9PZ+bMmXz88cf07dsXyC2cGzZsmPd6e/fu5fPPP+fQoUPUr5/7NXrwwQdZunQpH374IS+88EJpfSXOSSO2IiJSebh7QbtbALg2cyG5ZWxBJhCXlMGm2MRyiyZSUbRp0ybvvt1ux9/fn9atW+dtCwwMBCA+Pp6dO3fSpUsXDOOfXwAvu+wyUlNTOXToEHv37iUrK4suXbrk7ffz86N58+Z5j7du3Yppmlx88cXUrFkz7/bTTz+xd2/5XgJbI7YiIlK5dL4Dx8bpdCWGgY5NfOe85IxN41MyyjGYVGnu1XNHTq167eI0d3fP99gwjHzbThexTqcT0zTzFbVA3vrQhmEUaa1op9OJ3W5ny5Yt2O32fPtq1qxZrOznS4WtiIhULrWbEtfqLhr++i5PuM9hVWZbTuFVaNMA78K3ixSbYRRrOkBlERISwldffZWvwF23bh3e3t40aNCA2rVr4+7uzoYNG2jcuDEAJ06c4I8//qBHjx4AtG/fHofDQXx8PJdffrll7wU0FUFERCqhoKseJY66NDASCHdbVGC/AQT5etE52K/8w4lUInfffTcHDx5k/Pjx7Nq1i0WLFhEREcHEiROx2WzUrFmT22+/nYceeogffviBHTt2MHr0aGy2f0rIiy++mFtuuYWRI0cyf/58YmNjiYqK4uWXX2bJkiXl+n5U2IqISKVj96xBXJcIAO6wf0tTo+AKCBGDQ3TimMg5NGjQgCVLlrBp0ybatm3L2LFjuf322/NOPgN49dVX6d69O0OGDKFPnz5069aNjh075jvOhx9+yMiRI3nggQdo3rw5Q4YMYePGjTRq1Khc349huviFtpOTk/H19SUpKQkfH5+iPzEr7Z+zI4t5tqJVsrOzWbJkCYMGDSow/0ZESu5sn60S9zGViGXv0TQ59v5V1P1rLSsdbbk1+2Fyx2phSJsg3rm5Q5m9tPrTqi0jI4PY2FiCg4Px8tJ0lvJytq97UfsZlx2xjYyMJCQkhLCwMKujiIhISRgGda9/G9PmTi/7Nub1SOTeK3IvtfvDrniOp2ZaHFBEypvLFrbh4eHExMQQFRVldRQRESmpOhdidB0PwCW/v8qEHo1o29CXtCwH7/xQvgvDi4j1XLawFRGRKqL7g+DTAE4ewLbubSYNbAnAZxsP8OexVIvDiUh5UmErIiKVm0cN6P987v21b9KldjK9WwSQ4zR5ddnv1mYTkXJVJQrbb775hubNm3PRRRfxwQcfWB1HRETKW8gwCO4BjkxY9hiPDGyBzYDvdhxly/4TVqeTSsrFz68vd6Xx9a70hW1OTg4TJ07kxx9/ZOvWrbz88sskJuoSiiIiLsUwYNCrYHOD35dwcdJ6/tcpd5mhF5fsVIEixXJ6pYv09HSLk7iW01/v81lppNJfeWzTpk20atWKBg0aADBo0CCWLVvGTTfdZHEyEREpV3Wbw6XjYN278N3D3D9iFQujD7N5/wmWx/xF/1b1rE4olYTdbqdWrVrEx8cDUL169QKXnZXSY5om6enpxMfHU6tWrQKX5S0Oywvb1atX8+qrr7Jlyxbi4uJYsGABw4YNy9dmypQpvPrqq8TFxdGqVSveeuutvEu2HTlyJK+oBWjYsCGHDx8uz7cgIiIVRY9HYPsXcCKWwB0zuOPyIbz74x5e/m4XV7QIwN1e6f9QKeWkXr3cX4ROF7dS9mrVqpX3dS8pywvbtLQ02rZty6233sq1115bYP+8efOYMGECU6ZM4bLLLmPatGkMHDiQmJgYGjduXOifl/RblYiIi/L0hn7PwfwxsPp17rrjWj7d6MGfx9OYF3WQ4Zc2sTqhVBKGYRAUFERAQADZ2dlWx6ny3N3dz2uk9jTLC9uBAwcycODAM+5/4403uP322xkzZgwAb731FsuWLWPq1Km8+OKLNGjQIN8I7aFDh7jkkkvOeLzMzEwyM/9ZtDs5ORnIvYpMsb5xs7Nxz7ubDUbF/6Y//f70ARUpXWf7bOnzZoHW18GW2bB/LTVXRXBf72eJ+Po33vr+D4a1b0BNT8t/9EklYrfbS6XgkvJRoT/dWVlZbNmyhUmTJuXb3q9fP9atWwdA586d2bFjB4cPH8bHx4clS5bw1FNPnfGYL774IpMnTy6wffny5VSvXr3I2eyOTK76+/6yZctx2D2L/FyrrVixwuoIIlVSYZ8tnXxiAcOAQa/A+5fDzq+5uf1oPvSvzr6EdGas/pP7+15sdUIRKSMVurA9fvw4DoeDwMDAfNsDAwM5evQoAG5ubrz++uv06tULp9PJww8/jL+//xmP+eijjzJx4sS8x8nJyTRq1Ih+/foV7xrnWWmwPfdu//79ctdRrOCys7NZsWIFffv21bXNRUrR2T5bp/8qJOUssBV0vhM2TsV92SM80u8rxn2+gxlr/uSWSxoT4ON17mOISKVToQvb0/47Z9Y0zXzbhgwZwpAhQ4p0LE9PTzw9C46uuru7F6/YM/9p6+7uDpWoUCz2exWRIinss6XPmoV6PQo7voKE3QxIWUC7RmFEHzzJWz/s5oWrW1udTkTKQIU+PbROnTrY7fa80dnT4uPjC4ziioiI5OPlC32fAcD46RUietQCYF7UQfbEp1gYTETKSoUubD08POjYsWOBeWsrVqyga9euFqUSEZFKo80N0OgSyE6j/c7X6BsSiMNp8vJSXWpXpCqyvLBNTU0lOjqa6OhoAGJjY4mOjubAgQMATJw4kQ8++IBZs2axc+dO7r//fg4cOMDYsWMtTC0iIpWCzQaDXgPDBr/N5+nWCdhtBiti/mJTrK5SKVLVWF7Ybt68mfbt29O+fXsgt5Bt37593soGN9xwA2+99RbPPPMM7dq1Y/Xq1SxZsoQmTbQWoYiIFEFQG+h0GwAN1j3FTZ2CAHhBl9oVqXIsP3msZ8+e5+xY7r77bu6+++5ySiQiIlXOFU/Abwvg2C4eCfmJ+R4tiD54ku92HGVQ6yCr04lIKbF8xFZERKTMVasNfZ4GwHv9a0y4JHd5x1eW7iIrx2lhMBEpTSpsRUTENbQbDg06QlYKt6bPok5NT/YlpPPpxv2s35vAoujDrN+bgMOp6QkilZXlUxFERETKxekTyWZcgftvX/B850HctdqTZ76J4d8z4oJ8vYgYHMKAUE1REKlsNGIrIiKuo0EH6DgKgC6/v4QdB/89zeNoUgbj5mxl6Y44CwKKyPlQYSsiIq7liqcwq9XGJ+l3htu/L7D7dJ07eXGMpiWIVDIqbEVExLXU8Ce2zUQAHnKbRzPjcIEmJhCXlKG1bkUqGRW2IiJSIpGRkYSEhBAWFmZ1lGL7td4w1jtCqGlkMM39TWpwqtB28SkZ5ZxMRM6HClsRESmR8PBwYmJiiIqKsjpKsQX41OCe7PHEmX5caDvCa+7v888khH+18/Yq/3AiUmIqbEVExOV0DvbDwzeQu7PuI8u0M9AexVj74rz9BrmrI3QO9rMupIgUmwpbERFxOXabQcTgEKLNi3g6ZzSQO9/2MtuvQO7YbcTgEOw2w7qQIlJsKmxFRMQlDQgNYurwDqysMYj/y+mB3TB51/1dGnAMgNrVPSxOKCLFpcJWRERc1oDQINZO6k2TEVM5WasVfkYqc2tNwZMs7p8XTdKpbKsjikgxqLAVERGXZrcZXNK8AbVGz4VqfjTK+J03an7CkaRTPL7gV8z/XsFBRCosFbYiIiIAtRrDdTPBsHFlzg/c4raSb7bHseCXguvcikjFpMJWRETktGZXwBVPAvCM+0e0M/bw1KLfOJCQbnEwESkKFbYiIiL/1u1+aHEVdjObmdXexjMzgQnzfiHH4bQ6mYicgwpbERGRfzMMGDYV/C/C35nAVM/32HYggfdW7rE6mYicgwpbERGR//LygRs/BY+adDZ+4xG3ubzzw2627E+0OpmInIUKWxERkcLUbQ7DpgBwp9u3DDQ2MGFeNCkZWgJMpKJSYSsiInImIUPhsvsAeNVjOl4ndhPx9W8WhxKRM1FhKyIicjZXPAXBPahOBtPd32DF1t18sz3O6lQiUggVtiIiImdjd4PrZoFPQ4JtR3ndfSoRi38jMdPqYCLyXypsRUREzqVGHbjhY0y7B/3sWxiRPZ85u+04nLoqmUhF4mZ1gCohq5Is3J2djd2RCVlpYLpbnUak8nKvnrsklLiWBh0xrnwdvh7PA25f8GtaMNNWX0jnC+oSn5JBgLcXnYP9sNv0vSFiFRW2peG1C61OUCTuwFUA2y0OIlLZNboUbluq4tYVdRgJhzZj2/oR77i/x1U/BvHmDwF5u4N8vYgYHMKA0CALQ4q4Lk1FKCn36rk/3ETE9RzcANmV5C81UvoGvcqJ2q2pZaQxzf1NvPhnsu3RpAzGzdnK0h06uUzEChqxLSnDyB2xqUQ/3LKzs1m2bDn9+/fD3V1TEUSKLSu90vyFRsqOw+bBqNR7mGU+RCvbfp53n8UD2WMBAxMwgMmLY+gbUk/TEkTKmQrb82EY4FHD6hRFZ2TjsHvmZlZhKyJSIptiE9me4s1423jmuL/AtfY17HA25UPHQABMIC4pg02xiXRp5m9tWBEXo6kIIiIixRCfkgHAemcrXsy5GYAn3eYwwLap0HYiUn5U2IqIiBRDgLdX3v0PHIOYk9Mbm2HytnsknYxdhbYTkfKhwlZERKQYOgf7Uc/HE/6eUftUzq2scHTE08jmA4/XaWYcJsg3d+kvESlfKmxFRESKwW4zeGJQCyD3RDEnNsZn38NW54XUMtL4yONlnuyh9WxFrKDCVkREpJj6twrktoudBPp4ApCBJ7dnPUisWY+GxnHarb4DMlMsTinierQqgoiISAm09Td5+Jbu/HIoJe/KY4nJzam5cCj1M3ZzYvZN1B6zAOxahUakvGjEVkREpITsNoMuzfwZ2q4BXZr507F9B75q/gbppie149aQvSAcTNPqmCIuQ4WtiIhIKRp+7dU85fEQOaYN9x3z4MfnrI4k4jJU2IqIiJSimp5uDLl+NI/l3J67Yc1rsHmWtaFEXIQKWxERkVLW/eK60H4Eb+VcA4D57QPw+3cWpxKp+lTYioiIlIHHrwzhM6+bmZvTE8N0whe3wqHNVscSqdJU2IqIiJQB32ruPH9NG57IuY1VjraQcwo++x8k7LU6mkiVpcJWRESkjPQNCeTKdo25O/s+/rBfCOkJMOdaSD1mdTSRKkmFrYiISBmKGNyKajV8uDntAU561ocTsZz44Go27jqIw6mlwERKkwpbERGRMuRXw4NnhoZyHF+uTn6ARLMmtU/+SuqnI+jx0gqW7oizOqJIlaHCVkREpIzZ//5pG2sGMSbrQTJMd3rbfyE8fQrj5mxRcStSSlTYioiIlCGH02Ty4pi8x1vNixmfPR6HaXCT20rG2xcweXGMpiWIlAIVtiIiImVoU2wicUkZ+batcHYiImc0ABPdv6Rb6lI2xSZakE6kalFhKyIiJRIZGUlISAhhYWFWR6nQ4lMyCt0+x9GXyJwhALzkNgPb79+UZyyRKkmFrYiIlEh4eDgxMTFERUVZHaVCC/D2OuO+V3NuYF5OT+yGSdjmB+HPn8oxmUjVo8JWRESkDHUO9iPI1wuj0L0Gj+Xczo90xubMgrk3w+Et5ZxQpOpQYSsiIlKG7DaDiMEhAIUWtw7s3JMVTkJAF8hKhTnXwbHfyzekSBWhwlZERKSMDQgNYurwDtTzzT8toZ6vF+0b1SLd6U6vQ3dwolZrOJUIHw+DkwesCStSiblZHUBERMQVDAgNom9IPTbFJhKfkkGAtxedg/0wTZOHv9rO/K2HueKve1jp/wq1UvbmFre3LYWaAVZHF6k0NGIrIiJSTuw2gy7N/BnargFdmvljtxm42W28dl1bhl/amBOmN/2P30+KVxAk7oU510BGktWxRSoNFbYiIiIWs9kMnh0ayl09LuAv/Bic9CBp7n5w9Ff47EbISrc6okiloMJWRESkAjAMg0kDWvBA34vZZwZxfeqDZNhrwIF18MVocGRbHVGkwlNhKyIiUkEYhsH43hfx5FUhxJhNGZ7+AFmGJ+xeBgvHgdNpdUSRCk2FrYiISAVze7dgXr62NVtowZ2Z9+LADr9+AUsfAdO0Op5IhaXCVkREpAK6Iawxb9/YnrV04P6scTgxYNN0WPWi1dFEKiwVtiIiIhXUkLb1eX94R5bauvFU9ujcjT+9DBumWppLpKJSYSsiIlKB9QkJ5MPRYXxlG8Br2dfnblw6CaI/tzaYSAXksoVtZGQkISEhhIWFWR1FRETkrC67sA5zxnTmI/fr+CBnIADmonDYtcTiZCIVi8sWtuHh4cTExBAVFWV1FBERkXPq2MSPz+/owhT30Xzp6I5hOjC/GA2xa6yOJlJhuGxhKyIiUtmENvBl3tjLeN3zHpY7OmI4MnF+fiMc+cXqaCIVggpbERGRSuSiQG/mjbucl2o8zDpHCLasVBwfXwPH/rA6mojlVNiKiIhUMo39q/PpuB487/sk25wXYM9IJHv2EDh50OpoIpZSYSsiIlIJBflW46OxvXmh9rPscdbHPS2OjA+HQOoxq6OJWEaFrYiISCVVp6Yn0+/qz0t1XuSQWQevpD9JmzUUMpKsjiZiCRW2IiIilZhvdXfeuusqXg98ieOmDzUSf+PkrGsh+5TV0UTKnQpbERGRSq6mpxsv3nEN7wW9RLJZjVrxUcTPugkc2VZHEylXKmxFRESqAC93O4+NuYmZjV4kw3QnIG4lBz68FZxOq6OJlBsVtiIiIlWEh5uN8beO4tMmz5Jj2mh8aDG7ZoeDaeJwmqzfm8Ci6MOs35uAw2laHVek1LlZHUBERERKj5vdxq2jxzL/o3Su2/8MLQ58xtIpnkxOHkxcUkZeuyBfLyIGhzAgNMjCtCKlSyO2IiIiVYzNZnDt6ImsaPogAAOOfUi/1EX52hxNymDcnK0s3RFnRUSRMqHCVkREpAoyDIMrRj7Be+b/AJjs/hHDbGvz9p+eiDB5cYymJUiVocJWRESkitoUm8hrmUP5MKc/AK+5v09v25a8/SYQl5TBpthEixKKlC4VtiIiUiKRkZGEhIQQFhZmdRQ5g/iUDMDgmZwRfOXohpvhJNL9HcKMXYW0E6n8VNiKiEiJhIeHExMTQ1RUlNVR5AwCvL0AMLHxSPadfO9oj5eRzUyP12hp7C/QTqSyU2ErIiJSRXUO9iPI1wsDyMGN8Oz72ORsjo+RzsceL9HY+IsgXy86B/tZHVWkVKiwFRERqaLsNoOIwSEAGEAmHozJepCdzsbUNZL4xP1FJl/hj91mWBtUpJSosBUREanCBoQGMXV4B+r55k43SKYGI7MeYb8zgCa2eNqsvB1OnbQ2pEgp0QUaREREqrgBoUH0DanHpthE4lMyCPD24lDCxVT/9jrqZezh2PSrqTvuW/CobnVUkfOiEVsREREXYLcZdGnmz9B2DejSzJ/LOndmRYepJJvVqXtiK0mfDAdHttUxRc6LClsREREXdePgQUyt/zwZpju+B38g46tx4HRaHUukxFTYioiIuCibzWDcqBE8U/0RckwbXjFf4Fj2OJi6EplUTipsRUREXJiPlzu33TqOpxgHgH3jFFj7hsWpREpGha2IiIiLuzCgJj3+dy/PZg/P3fDDM7BltqWZREpCha2IiIjQv1U9qve4l8icIQCY39wPMYssTiVSPCpsRUREBIAJfS5m8wX38FlOLwzTifnVGPjzJ6tjiRSZClsREREBcpcEe+umDnzgcw/fOcIwHFmYc2+Gw1utjiZSJCpsRUREJI9vNXemjryEScZ9/OxohZGVCp9eB8f+sDqayDmpsBUREZF8mtfz5oXrOnFn9kS2O4MhPQE+uRqSDlsdTeSsVNiKiIhIAVe2CWJEj1BGZz1CrBkEyYdyi9v0RKujiZyRClsREREp1EP9m9Pqogu4JfNR4g1/OP577rSEzFSro4kUSoWtiIiIFMpuM3jnxvbYajfi5oxHSLV5w+EtnJx9I19vjWX93gQcTl2lTCoOFbYiIiJyRrVreDBtREcOuTVmxKkHSTc9qRW3BmPBWG6ZsY5uL//I0h1xVscUAVTYioiIyDm0qu/LzZ0b84t5EWOzJ5Bl2hls38Bkt9kcTTrFuDlbVdxKhaDCVkRERM7K4TT5bsdRAFY72zIx+26cpsEIt++Z4PYVAJMXx2haglhOha2IiIic1abYROKSMvIef+PswlM5owG4z20+N9u/Jy4pg02xWjFBrKXCVkRESiQyMpKQkBDCwsKsjiJlLD4lo8C2OY6+vJVzDQDPuH1IP1tUoe1EypMKWxERKZHw8HBiYmKIioqyOoqUsQBvr0K3v5VzLZ/n9MJumLzj/h4XnPqtnJOJ5Oeyha1GGkRERIqmc7AfQb5eGAX2GDyRcxs/ONrjZWQTuvpOXXpXLOWyha1GGkRERIrGbjOIGBwCUKC4dWDnnuzxHKnZCuPUCZhzLaQcLf+QIrhwYSsiIiJFNyA0iKnDO1DPN/+0hOoedk7hxQ0p95Pl0xSSDuRenSwj2Zqg4tLcrA4gIiIilcOA0CD6htRjU2wi8SkZBHh70a5RLW6dvYkNf8LwzIeZW/0pbEd/hf8bATd/AW4eVscWF6IRWxERESkyu82gSzN/hrZrQJdm/lTzsDP1lo409a/OpqRaPOr1JKZ7DfhzFSwKB6fT6sjiQlTYioiIyHmpXcODD0aF4e3lxrwjdZkRFIFpc4Nf/w9+mGx1PHEhKmxFRETkvF0YUJMpt3TAbjN44Y+GrLz4ydwdP78FG6dZmk1chwpbERERKRWXX1Q3b/WE27ddxB+h9+fu+O4RiFlkYTJxFSpsRUREpNSM7NKUkV2aYJowbPslJIaMAEz46g7Yv87qeFLFqbAVERGRUvXUVSF0u7AO6VlOhuwZQuaFA8GRCZ/fCPG7rI4nVZgKWxERESlVbnYbkTd34IK6NTiUnM2Ik3fibNgZMpJyL+CQdNjqiFJFqbAVERGRUudb3Z2Zo8LwrebOpkOneMzrcUz/iyD5EMkzh7JpZywOp2l1TKliVNiKiIhImQiuU4OpwzvgZjOYuyON/sfvI96shU/ybhyf3Uyvl5axdEec1TGlClFhKyIiImWma7M63BDWCIA/Mv0YnfUwKWY1uthjeOjUm9w9Z7OKWyk1KmxFRESkzDicJj/uis97HGM25a7s+8ky7Qy2b+Axt8+YvDhG0xKkVKiwFRERkTKzKTaRuKSMfNvWOUN5KPsuAMa4LaF/6kI2xSZaEU+qGBW2IiJSIpGRkYSEhBAWFmZ1FKnA4lMyCt2+yNmNF7NvAuBJt08wdi8tz1hSRamwFRGREgkPDycmJoaoqCiro0gFFuDtdcZ90xxX8VlOL+yGSdjmByFuWzkmk6pIha2IiIiUmc7BfgT5emEUutfgqZxb2Wi0xZ6TDp/dAMlHyjmhVCUqbEVERKTM2G0GEYNDAAotbnNwY32H16FuC0iJyy1uM1PLN6RUGSpsRUREpEwNCA1i6vAO1PPNPy2huocdgBlRCfzZbxZUrwNHt8NXY8DpsCKqVHJuVgcQERGRqm9AaBB9Q+qxKTaR+JQMAry9aN+4FqM/3MSGPxMZvSCeb67+BJ95V8Mf38HyJ2DAi1bHlkpGI7YiIiJSLuw2gy7N/BnargFdmvnj5W5nyi0daeRXjQOJ6dy50iBn6NTcxhumwKYZ1gaWSkeFrYiIiFjGr4YHM0eFUcPDzoY/E3l670XQ+6ncnd89DLtXWBtQKhUVtiIiImKpiwO9efvG9hgGzNlwgE/croV2t4DphC9Gw9EdVkeUSkKFrYiIiFiuT0ggD/VvDsDTi2NYH/IENL0cslJzV0pIOWpxQqkMilXY/v777zz99NP07t2bZs2aERQURJs2bRg1ahSfffYZmZmZZZVTREREqrhxPZoxrF19HE6TcXN3cLDvNPC/CJIPwec3Qlaa1RGlgitSYfvLL7/Qt29f2rZty+rVqwkLC2PChAk8++yzDB8+HNM0efzxx6lfvz4vv/yyClwREREpNsMweOnaNrRt6MvJ9Gxum7eH1Os+g2p+cOQXmH8nOJ1Wx5QKrEjLfQ0bNoyHHnqIefPm4efnd8Z269ev58033+T111/nscceK7WQZSEyMpLIyEgcDq2TJyIiUlF4uduZPrITg99dy+74VO5blsT0Gz7F/slQ2PUNfP8U9HvO6phSQRWpsN29ezceHh7nbNelSxe6dOlCVlbWeQcra+Hh4YSHh5OcnIyvr6/VcURERORvgT5eTB/Zif9NW88Pu+J5rV4zHhk6BeaPgXXvgl8z6HSr1TGlAirSVISiFLXn015ERETk39o1qsWr17UBYOqqvSx0dIWef/81+NsHYO+PFqaTiqpEVx7btGkTq1atIj4+Hud/5rq88cYbpRJMREQqNk3pkrI2tF0Ddh1NYeqqvTz81Xaa3nkn7RL3wvZ58H+j4PblENDS6phSgRS7sH3hhRd44oknaN68OYGBgRiGkbfv3/dFRKRq05QuKQ8P9WvO7r9S+H5nPHd+soWvx75KvZMH4cA6+PR/cMcPUDPA6phSQRS7sH377beZNWsWo0ePLoM4IiIiIv+w2QzeurE910z5mT/+SuXOz3/l81s+wjarD9WS9pMy+3qq3/Edds/qVkeVCqDYF2iw2WxcdtllZZFFREREpICanm58MDKM2tXd2X4oiY5vbGHQ8Xs5YdbE+3g0q16+lqW/HrY6plQAxS5s77//fiIjI8sii4iIiEihGvtXZ3TXYAAysp3EmkGMzbqfLNNOb+c69s6bxNIdcRanFKsVeyrCgw8+yJVXXkmzZs0ICQnB3d093/758+eXWjgRERERAIfTZG7UgXzbNpotmZR9B294vE+429c8t7ARfUOex27TOT+uqtgjtuPHj2flypVcfPHF+Pv74+vrm+8mIiIiUto2xSYSl5RRYPt8Z3feyRkGwCPZ77Nr/bflnEwqkmKP2H788cd89dVXXHnllWWRR0TEZWRnZ3P06FHS09OpW7fuWa/sKOLq4lMKFrWnvZFzPU2NvxhiX8+FP90DrdpArcblmE4qimKP2Pr5+dGsWbOyyCIiUuWlpqYybdo0evbsia+vL02bNiUkJIS6devSpEkT7rjjDqKioqyOKVLhBHh7nWWvwUPZd7HdGYxn1gmYewtkpZdbNqk4il3YPv3000RERJCerm8YEZHiePPNN2natCkzZszgiiuuYP78+URHR/P777+zfv16IiIiyMnJoW/fvgwYMIDdu3dbHVmkwugc7EeQrxdnmj2biQdPeU3CrF4Hjm6HxfeBaZZrRrFesacivPPOO+zdu5fAwECaNm1a4OSxrVu3llo4EZGqZN26daxcuZLWrVsXur9z587cdtttvP/++8ycOZOffvqJiy66qJxTilRMdptBxOAQxs3ZigEUVrI2anoxRpfZ8PFQ+PX/oH476BJevkHFUsUubIcNG1YGMUREqr4vvviiSO08PT25++67yziNSOUzIDSIqcM7MHlxTL4TyXyruZN0KpvF2+PoH9qeq/q/AEsfgeVPQmAoXNDDwtRSnopd2EZERJRFDhEREZFzGhAaRN+QemyKTSQ+JYMAby86B/vxwpKdzFwbywP/t42Gd95Au7bRsO1z+GI03LkKajexOLmUh2IXtkVhmiaGoTXkRETOJCMjg3fffZeVK1cSHx+P0+nMt1/TukTOzG4z6NLMP9+2xwa1JPZ4Gj/uimfMx1tYNPYFGsTvhLhomHcL3LYcPHTZ3aquSCePtWzZks8++4ysrKyzttu9ezfjxo3j5ZdfLpVwIiJV1W233cYrr7xCkyZNuOqqqxg6dGi+W2UQGRlJSEgIYWFhVkcRwW4zeOem9rSo583x1EzGfPob6dd8DNXrwNFfYfG9OpnMBRRpxDYyMpJHHnmE8PBw+vXrR6dOnahfvz5eXl6cOHGCmJgY1q5dS0xMDPfcc4/mhomInMO3337LkiVLuOyyy6yOUmLh4eGEh4eTnJysC/RIhVDT040PRnViWOTP7IxL5t4lx5h23Wzsc4bBr19AUDvoeo/VMaUMFamwveKKK4iKimLdunXMmzePzz77jH379nHq1Cnq1KlD+/btGTlyJMOHD6dWrVplHFlEpPJr0KAB3t7eVscQqXIa1q7O9JGduHH6Br7f+Rcv172Ax/q/AN89DCuehHqhcEFPq2NKGSnWHNuuXbvStWvXssoiIuIyXn/9dR555BHef/99mjTRSS0ipalD49q8el0b7psbzfTVf9LsmgHc0DYatn0GX9wKd66E2k2tjilloNgXaBARkfPXqVMnMjIyuOCCC/D29sbPzy/fTUTOz9B2Dbivd+460I8v/I0NrZ6E+u3hVCLMHa4rk1VRZbIqQmUQGRlJZGQkDofD6igi4oJuuukmDh8+zAsvvEBgYKBWkhEpAxP6XMSfx9NYvO0Id32+g8Ujp9P4y0Hw16/w9T1w7UzQZ69KcdnCVic9iIiV1q1bx/r162nbtq3VUUSqLMMwePW6NhxMTCf64ElGf3WEr4fMpOa8a2DHV7knk112r9UxpRRpKoKIiAVatGjBqVOnrI4hUuV5uduZPrIj9X29+PN4Gneu9sTR7/ncnd9HwN4frQ0opUqFrYiIBV566SUeeOABVq1aRUJCAsnJyfluIlJ6Ary9mDk6jBoedtbtTeCJw10x290MphO+vA1O7LM6opSSEk1FcDqd7Nmzp9Cr5XTv3r1UgomIVGUDBgwAoHfv3vm2n75yo+b/i5SulkE+vH1je+74ZDOfRx2k+cD7GF1/FxzZCnNvgduXg0cNq2PKeSp2YbthwwZuvvlm9u/fj/mfK3ioMxYRKZqVK1daHUHE5fQJCeTxQS157tudTF66lwuve4tuSdfCXztg0T1w3SydTFbJFbuwHTt2LJ06deLbb78lKChIZ/KKiJRAjx49rI4g4pJu7xbM3mOpfL7pIHctiuPbYVNp+s2N8Nt8qN8OLrvP6ohyHoo9x3b37t288MILtGzZklq1auHr65vvJiIihTtw4ECx2h8+fLiMkpSOyMhIQkJCCAsLszqKSJEZhsEzQ0PpcoE/aVkOblluJ6XXc7k7v39aJ5NVcsUubC+55BL27NlTFllERKq0sLAw7rjjDjZt2nTGNklJScyYMYPQ0FDmz59fjumKLzw8nJiYGKKioqyOIlIs7nYbU4d3ILhODQ6fPMXI7a3JaXtL7slkX9wKibFWR5QSKvZUhPHjx/PAAw9w9OhRWrdujbu7e779bdq0KbVwIiJVyc6dO3nhhRcYMGAA7u7udOrUifr16+Pl5cWJEyeIiYnht99+o1OnTrz66qsMHDjQ6sgiVVat6h7MHNWJq6es45eDSTxcaySvN9iFcXhL7slkY1boZLJKqNiF7bXXXgvAbbfdlrfNMAydySsicg5+fn689tprPPfccyxZsoQ1a9awb98+Tp06RZ06dbjlllvo378/oaGhVkcVcQkX1K3J1OEdGDlzE/N/TSDk8mcZc3I0xP8Gi8Lhug91MlklU+zCNjZWw/MiIufDy8uLa665hmuuucbqKCIur2uzOjw3LJRJ83/luTUnadH3Dbr9PBp+W5B7ZbJuEyxOKMVR7MK2SZMmZZFDREQsdPXVV7Nq1Sp69+7Nl19+aXUckXJ1Y+fG7D2Wyow1sdy20o0fLo+g0fonMX+YzPr0+hjNetM52A+7TaO3FV2RCtuvv/6agQMH4u7uztdff33WtkOGDCmVYCIiUn7uvfdebrvtNj766COro4hYYtLAlvx5LI0fdsUzaP3FPGv0Zpj5AyE/T2DIyufI9mlCxOAQBoQGWR1VzqJIhe2wYcM4evQoAQEBDBs27IztNMdWRKRy6tWrF6tWrbI6hohl7DaDt29qT/83f+LwyQweZiRNPPbT3raHae5vcG3SZMbN2crU4R1U3FZgRVruy+l0EhAQkHf/TDcVtSIi5W/16tUMHjyY+vXrYxgGCxcuLNBmypQpBAcH4+XlRceOHVmzZk35BxWp4Kq528ly5F5VNQt3xmZNIN6sRUvbQV5xnw6YTF4cg8Npnv1AYpliz7EVEZGKJS0tjbZt23LrrbfmrVzzb/PmzWPChAlMmTKFyy67jGnTpjFw4EBiYmJo3LhxsV8vMzOTzMzMvMfJyckAZGdnk52dXfI3Uomcfp+u8n5dxcbYRI6l/PO9/Rd+jMu6j889nuMq+wa2O4OZnjSY9XviuSTYz8KkrqeonzUVtiIiFnA6ndhsBf9o5nQ6OXToULEKzoEDB551zds33niD22+/nTFjxgDw1ltvsWzZMqZOncqLL75Y7OwvvvgikydPLrB9+fLlVK9evdjHq8xWrFhhdQQpRVuOG4A9/zazOZNzRvG8+ywedpvHZmdzlq9xkLBTo7blKT09vUjtXLawjYyMJDIyUtMnRKRcJScnM2bMGBYvXoyPjw9jx47lqaeewm7P/WF67NgxgoODS61vysrKYsuWLUyaNCnf9n79+rFu3boSHfPRRx9l4sSJeY+Tk5Np1KgR/fr1w8fH57zyVhbZ2dmsWLGCvn37FrhQkVRe/rGJfLx7c4Htnzp6c4ltJ0Ps63nH4z0OXvIdnZoHW5DQdZ3+y9C5uGxhGx4eTnh4OMnJyfj6+lodR0RcxJNPPsm2bdv45JNPOHnyJM899xxbtmxh/vz5eHh4AGCapTcSdPz4cRwOB4GBgfm2BwYGcvTo0bzH/fv3Z+vWraSlpdGwYUMWLFhAWFhYocf09PTE09OzwHZ3d3eXK/Jc8T1XZV0uDCDI14ujSRnk/xQaPJZ9O22MP2lq+4sGOyZjtJqjizeUo6J+zop08piIiJSOhQsXMm3aNK677jrGjBnDli1bOH78OIMHD86bt2qUwQ/L/x7z9NUiT1u2bBnHjh0jPT2dQ4cOnbGoFanK7DaDiMEhAPz3U5hKdcZnj8dpuGPs+gaiPij/gHJOJSpsnU4nf/zxB2vXrmX16tX5biIicmbHjx/Pd6Ebf39/VqxYQUpKCoMGDSryPLKiqlOnDna7Pd/oLEB8fHyBUVwRgQGhQUwd3oF6vl75thvAr+YFLA0al7th2WMQt738A8pZFXsqwoYNG7j55pvZv39/gT+XaR1bEZGza9SoETt37iQ4+J/5ed7e3ixfvpx+/fpx9dVXl+rreXh40LFjR1asWJHv2CtWrGDo0KGl+loiVcWA0CD6htRjU2wi8SkZBHh7kXIqm7s+3cLdf17Cz4170iB+FXx5K9z5E3jWtDqy/K3YI7Zjx46lU6dO7Nixg8TERE6cOJF3S0xMLIuMIiJVRt++ffnwww8LbK9ZsybLli3Dy8urkGedXWpqKtHR0URHRwMQGxtLdHQ0Bw4cAGDixIl88MEHzJo1i507d3L//fdz4MABxo4de17vRaQqs9sMujTzZ2i7BnRp5k+/0HpM7HMxYDDs8M1kVQ+ChD3w7QNWR5V/KfaI7e7du/nyyy+58MILyyKPiEiV9swzz3DkyJFC93l7e/P999+zZcuWYh1z8+bN9OrVK+/x6RULRo0axezZs7nhhhtISEjgmWeeIS4ujtDQUJYsWZJvSoSInFt4rwv57UgyS387Snjm3Uw3IjC2z4ULekK7m6yOJ5RgxPaSSy5hz549ZZFFRKTKCw8P57fffjvj0jU1a9akR48exTpmz549MU2zwG327Nl5be6++2727dtHZmYmW7ZsoXv37ufzNkRcks1m8Pr/2tI80JsVac34rNrNuTu+fQCO77Y2nABFHLHdvv2fydHjx4/ngQce4OjRo7Ru3brA8gtt2rQp3YQiIlVI8+bNefnllxk5ciTdu3dn6NChDBkyhEaNGlkdTUSKoIanG9NHdmTIez/zZOIAOvvHcFHaVvhiNIz5AdyLP51ISk+RCtt27dphGEa+k8Vuu+22vPun9+nkMRGRs4uIiCAiIoJDhw7x9ddfs2jRIh544AFCQkIYMmQIQ4cOpX379lbHFJGzaOJfg3dvas/oDzdxc8LtrPaOpdpfO2D543Dl61bHc2lFKmxjY2PLOoeIiEtp2LAhd999N3fffTcpKSl89913LFq0iN69e+Pt7c3gwYMZN24crVq1sjqqiBSi+8V1mTSwBS8s2cXd6XfyofvLuWvbBneHEK04YpUizbFt0qRJ3m3//v00aNAg37YmTZrQoEED9u/fX9Z5RUSqHG9vb/73v//x6aefcuzYMWbNmoXdbmf9+vVWRxORs7jj8gsY1q4+Kx1tmW38XcwuGg8nVA9ZpdirIvTq1Yu4uDgCAgLybU9KSqJXr16aiiAiUgJZWVlkZWVRs2ZNevfuTe/eva2OJCLnYBgGL13bhj3HUnnu8LVcUmMXLTN/h69uh1u/A7sut1zeir0qwn8vw3haQkICNWrUKJVQIiJV2Ycffsj48eP59NNPAXj00Ufx9vbG19eXvn37kpCQYHFCESkqL3c700Z0wrdGde5Iv5t0W004FAU/Pmt1NJdU5BHba665Bsj97WT06NF4enrm7XM4HGzfvp2uXbuWfkIRkSrk+eef5/nnn6dr16589tlnrF27loULF/LMM89gs9l45513eOKJJ5g6darVUUWkiBrUqsaUWzpwywfZ3J8xhmkeb8HPb0PT7nBRH6vjuZQiF7a+vr5A7oitt7c31apVy9vn4eHBpZdeyh133FH6CUVEqpDZs2czc+ZMbrrpJjZv3swll1zCvHnzuO666wAIDQ3VFcFEKqFLLvDnqcEhPLXI5BNHX0bYV8CCu2DsWvAJsjqeyyhyYXv6EpBNmzblwQcfrPTTDiIjI4mMjNScYBEpVwcOHKBbt24AdOrUCTc3N1q3bp23v02bNsTFxVkVT0TOw4hLm/Db4WSe23wLne1/0Dx9P8y/A0YuApvd6nguodhzbCMiIqhRowbx8fGsWbOGtWvXEh8fXxbZylR4eDgxMTFERUVZHUVEXEh2dna+qVweHh75LnTj5uamX7hFKinDMHhmWCtCGgdwd+Y9nMIT9q3hwKLnWBR9mPV7E3A4zXMfSEqs2KsiJCcnEx4ezty5c/M6X7vdzg033EBkZGTelAURESlcTEwMR48eBXKnd+3atYvU1FQAjh8/bmU0ETlPnm523h/ekcHvnuKJtFt53eN9GkS/yQObahJltiDI14uIwSEMCNX0hLJQ7BHbMWPGsHHjRr755htOnjxJUlIS33zzDZs3b9YcWxGRIujduzft2rWjXbt2pKenc9VVV9GuXTvat29Pnz460USksgv08WJ016Z85ezOV45u2A2Ttz3eoxYpHE3KYNycrSzdoSlHZaHYI7bffvsty5Yty5sjBtC/f39mzJjBgAEDSjWciEhVoys5ilR9DqfJJxtyL9LwZPZttDP20swWx2vu7zMm+0EMDCYvjqFvSD3stoJLqErJFbuw9ff3L3S6ga+vL7Vr1y6VUCIiVVWTJk3O2SY6OrpI7USkYtoUm0hcUgYA6XgxPns8Czwi6GP/hZudP/KZozdxSRlsik2kSzN/i9NWLcWeivDEE08wceLEfGftHj16lIceeognn3yyVMOJiLiKpKQkpkyZQocOHejYsaPVcUTkPMSnZOR7HGM25ZWcGwB43G0OjYy/Cm0n56/YI7ZTp05lz549NGnShMaNGwO5y9d4enpy7Ngxpk2bltd269atpZdURKQK+vHHH5k1axbz58+nSZMmXHvttcycOdPqWCJyHgK8vQpsm+UYQD/7Zi6x7eI192ncmPVEoe3k/BS7sB02bFgZxBARcR2HDh1i9uzZzJo1i7S0NP73v/+RnZ3NV199RUhIiNXxROQ8dQ72I8jXi6NJGZxe3MvExoPZd7HUYxKX2HZxX40f6Bx8laU5q6JiF7YRERFlkUNExCUMGjSItWvXctVVV/Huu+8yYMAA7HY777//vtXRRKSU2G0GEYNDGDdnKwbkFbcHzUBeyLmF591ncY/5GfbEcKhzkZVRq5xiz7EFOHnyJB988AGPPvooiYmJQO60g8OHD5dqOBGRqmb58uWMGTOGyZMnc+WVV2K362pEIlXRgNAgpg7vQD3f/NMNvjL6strRGjdnJjlf3QWOHIsSVk3FHrHdvn07ffr0wdfXl3379nHHHXfg5+fHggUL2L9/Px9//HFZ5BQRqRLWrFnDrFmz6NSpEy1atGDEiBHccMMNVscSkTIwIDSIviH12BSbSHxKBgHeXrQI8uau99Jol3YvPnFbcP78NrbuD1gdtcoo9ojtxIkTGT16NLt378bL65/fQgYOHMjq1atLNZyISFXTpUsXZsyYQVxcHHfddRdz586lQYMGOJ1OVqxYQUpKitURRaQU2W0GXZr5M7RdA7o086d2dQ+eHTmAF8zRAJgrX4CjO6wNWYUUu7CNiorirrvuKrC9QYMGeZeIrAwiIyMJCQkhLCzM6igi4oKqV6/Obbfdxtq1a/n111954IEHeOmllwgICGDIkCFWxxORMtS8njeXDgtnhaMjdjOHlLljICfL6lhVQrELWy8vL5KTkwts//3336lbt26phCoP4eHhxMTEEBUVZXUUEXFxzZs355VXXuHQoUN8/vnnVscRkXIwrENDtraJINGsiffJnSQtf8HqSFVCsQvboUOH8swzz5CdnQ2AYRgcOHCASZMmce2115Z6QBERV2G32xk2bBhff/211VFEpBxMuLobM33uAaDmprfJ3K/BtvNV7ML2tdde49ixYwQEBHDq1Cl69OjBhRdeiLe3N88//3xZZBQRERGpcjzd7Nxy+wSW0RU7Tk5+djtk62pk56PYqyL4+Piwdu1afvzxR7Zu3YrT6aRDhw706dOnLPKJiIiIVFn1a1Vj/3Vvc+yLPgRm7mfX54/QYuTbVseqtIpd2J52xRVXcMUVV5RmFhERERGX0yX0Yr7ZGcFVv03k4r0f8eeWq7igY1+rY1VKxZqK4HQ6mTVrFldddRWhoaG0bt2aIUOG8PHHH2Oa5rkPICIiIiIFDLr2NtbW7I/NMPH45h4SE0+wfm8Ci6IPs35vAg6n6qyiKPKIrWmaDBkyhCVLltC2bVtat26NaZrs3LmT0aNHM3/+fBYuXFiGUUVERESqJpvNoPXtU/nrnc40NI8y5607eSJrdN7+IF8vIgaHMCA0yLqQlUCRR2xnz57N6tWr+eGHH/jll1/4/PPPmTt3Ltu2beP777/nxx9/1FXHRERERErIt7Y/m9s8C8Bw23K62v65cMPRpAzGzdnK0h1xVsWrFIpc2H7++ec89thj9OrVq8C+K664gkmTJvHpp5+WajgRERERV+Fwmjy3K5BPcnJPyH/FfTrepANweiLC5MUxmpZwFkUubLdv386AAQPOuH/gwIFs27atVEKJiIiIuJpNsYnEJWXwYs7N7HcG0NA4ziS3fy7aYgJxSRlsik20LmQFV+TCNjExkcDAwDPuDwwM5MSJE6USqjzokroiIiJSkcSn5K5hm44XD2ffBcBN9h8JNf4stJ0UVOTC1uFw4OZ25nPN7HY7OTk5pRKqPOiSuiIiIlKRBHh75d3faLZkoaMrNsPkGffZGDgLbSf5FWtVhNGjR+Pp6Vno/szMzFILJSIiIuJqOgf7EeTrxdGkDEzghexb6GPbSgfbHq6zr+YLR0+CfL3oHOxnddQKq8gjtqNGjSIgIABfX99CbwEBAYwcObIss4qIiIhUWXabQcTgEAAMIJ7avJ1zDQCPuM3FhzQiBodgtxkWpqzYijxi++GHH5ZlDhERERGXNyA0iKnDOzB5cQxxSRnMdgzgBvsqLrQd4X63LzHNy62OWKEV68pjFdXVV19N7dq1ue6666yOIiIiInJeBoQGsfaRK/j8jkt57cZOZPV7CYCR9uXMXvANiWlZFiesuKpEYXvvvffq4hAiIiJSZdhtBl2a+TO0XQNCug3F0XIIdsNkYs4HRCzace4DuKgqUdj26tULb29vq2OIiIiIlAl7/xdwulXjEtsubDu+0BXIzsDywnb16tUMHjyY+vXrYxgGCxcuLNBmypQpBAcH4+XlRceOHVmzZk35BxURERGxSq1G2Lo/CMDj7p/x4oJNmpJQCMsL27S0NNq2bct7771X6P558+YxYcIEHn/8cX755Rcuv/xyBg4cyIEDB8o5qYiIiIiFuo7HWfsCAoyT3Jw5j6e//s3qRBVOkVdFKCsDBw5k4MCBZ9z/xhtvcPvttzNmzBgA3nrrLZYtW8bUqVN58cUXi/16mZmZ+dbcTU5OBiA7O5vs7OxiH68yOf3+qvr7FCkz2dm4593NBiP/Z6qwz5Y+byJSatw8sQ16BT69jtvsSxm4vQdLWwcxILSe1ckqDMsL27PJyspiy5YtTJo0Kd/2fv36sW7duhId88UXX2Ty5MkFti9fvpzq1auX6JiVzYoVK6yOIFIp2R2ZXPX3/WXLluOw579gTWGfrfT09HJIJiIu46K+0HwQ7r8v4Wm3j5iw4AIuCfajdg0Pq5NVCBW6sD1+/DgOh4PAwMB82wMDAzl69Gje4/79+7N161bS0tJo2LAhCxYsICwsrNBjPvroo0ycODHvcXJyMo0aNaJfv374+PiUzRupILKzs1mxYgV9+/bF3d393E8Qkfyy0mB77t3+/fuBRw3g7J+t038VEhEpNQNexNzzA934jbBTa3h6cR3evrG91akqhApd2J5mGPmvsGGaZr5ty5YtK/KxPD09C70ssLu7u8sUe670XkVKlfnP58bd3R3+8zkq7LOlz5qIlLraTTG63Q8/vcST7nPoHd2WZa2D6N9KUxIsP3nsbOrUqYPdbs83OgsQHx9fYBRXRERExGV0mwC1GlPfSCDcbRGPL9jBCa2SULELWw8PDzp27Fhg3tqKFSvo2rWrRalERERELOZeDQbkXpHsTrdvqZm2j8mLtUqC5VMRUlNT2bNnT97j2NhYoqOj8fPzo3HjxkycOJERI0bQqVMnunTpwvTp0zlw4ABjx461MLWIiIiIxZoPggv74LHne552+5jR0fUY1DqIfi48JcHywnbz5s306tUr7/HpE7tGjRrF7NmzueGGG0hISOCZZ54hLi6O0NBQlixZQpMmTayKLCIiImI9w4CBr8CUS+nJNvo4tvL4Qi86B/tRq7prrpJg+VSEnj17Yppmgdvs2bPz2tx9993s27ePzMxMtmzZQvfu3a0LLCIiIlJR+DeDLvcA8KznJySnpDB5cYzFoaxjeWErIiIiIueh+4Pg04AgM55xbotZ8MthVsT8ZXUqS6iwFREREanMPGpA/+cBCHf/hoZGPI8t+JWT6a63SoIKWxEREZHKLmQYBHfH3czk5Rqfcywlk2dccEqCClsRERGRys4wYOCrYHPjspyN9LJHM/+Xw3zvYlMSVNiKiIiIVAUBLeCS3OVQX6/5OR5k89iCX0lKz7Y4WPlRYSsiIiJSVfR4BGoG4pd5kId8vic+JZOnF+9g/d4EFkUfZv3eBBxO0+qUZcbydWxFREREpJR4+UC/52D+Hdzm/JIP6cyCX2DBL0fymgT5ehExOIQBoUEWBi0bGrEVERERqUpaXw+Nu2DPOcVj7nMK7D6alMG4OVtZuiPOgnBlS4WtiIiISFViGDgGvIIDG1fZN9LVtiPf7tMTESYvjqly0xJU2IqIiIhUMZtONeDjnL4APOH2Kf+Us7lMIC4pg02xieUfrgypsBURkRKJjIwkJCSEsLAwq6OIyH/Ep2TwVs61pJpehNj2c4XtlzO2q0pU2IqISImEh4cTExNDVFSU1VFE5D8CvL1IoiafOHJHbce7LeS/o7an21UlKmxFREREqpjOwX4E+XoxK2cQGaY77W176Gr7LW+/Qe7qCJ2D/awLWQZU2IqIiIhUMXabQcTgEI7jy+eOKwAYb1+Yt98EIgaHYLcZ1gQsIypsRURERKqgAaFBTB3egQXVryXLtNPFHkNH43cgd0RX69iKiIiISKUxIDSIBZP+x4mLrgfg3QY/ALB5XyL7E9KsjFYmVNiKiIiIVGF2m0HgoEfAsFH/+FpGNT2B04T3f9prdbRSp8JWREREpKrzuyD3imTABM/FAHy55RBHTp6yMlWpU2ErIiIi4gq6TQQMau9fynWNksl2mExf/afVqUqVClsRERERVxDQAloOBuDhGt8B8PmmAxxLybQyValSYSsiIiLiKro/CEDd/YsZUD+dzBwnM9fGWhyq9KiwFREREXEVQW3hon4YppMnfZcC8Mn6fZxMz7I4WOlQYSsiIiLiSro/BED9/YvoHphBWpaDD3/eZ22mUqLCVkRERMSVNOoMTS/HcGbzjH/uuraz1+0jJSPb4mDnT4WtiIiIiKv5e9S2yYGv6OSfRdKpbOZsOGBxqPOnwlZERETE1QR3h4ZhGDkZPF/vJwBmrv2TU1kOi4OdHxW2IiIiIq7GMPJGbS8++AWtamdzPDWLuVGVe9RWha2IiIiIK7qoH9RrjZGVyosN1gHw/qq9rP4jnkXRh1m/NwGH07Q4ZPG4WR1ARERERCxgGHD5g/DFKFofmkt9ry4cSYGRs6LymgT5ehExOIQBoUEWBi06jdiKiIiIuKqWQ6DOxRgZSQzNXlpg99GkDMbN2crSHXEWhCs+FbYiIiIirspmw9ltIgC3uy3Bi/yX1z09EWHy4phKMS1Bha2IiIiIC9tY4woOOOtSx0jmRvvKAvtNIC4pg02xieUfrphU2IqIiIi4sPi0HKY7rgLgBvuqM7dLySifQOdBha2IiIiICwvw9uJrRxcyTXda2g4QYuw7Y7uKToWtiIiIiAvrHOxHDd86fO/sAMA19jX59hvkro7QOdjPgnTFo8JWRERExIXZbQYRg0OY77gcgKH2n3EjJ1+biMEh2G2GFfGKRYWtiIiIiIsbEBrE/24cRSI+1DWSudz2a96+4V2aaB1bEREREak8+rdpTK1LbgHguabbGX5JYwBW/R5PtsNpZbQiU2ErIiIiIgDY2t0EQIO/VvH4FfWpU9OTg4mn+GrLIYuTFY0KWxERERHJVa8NBLQCRybVdn/NuJ7NAHj3xz1k5VT8UVsVtiIiIiKSyzCg7Y2597d9zi2XNCbA25PDJ0/xf5sPWputCFTYioiIiMg/2vwPDBsc3IhX8j7Ce10IQOTKPWRkOywOd3YqbEVERETkH971oNkVufe3zeWGsEYE+XoRl5TBvKiKPWqrwlZERERE8mubexIZ2+fiZTcqzaitClsRERERya/5IPDwhpMH4MA6/tepEQ1qVSM+JZNPNx6wOt0ZqbAVERERkfw8qkOrYbn3t32Oh5uNe67IHbWdsnIPP/0ez6Low6zfm4DDaVqX8z9U2IqIiIhIQaenI/y2CLLSua5jQ/xrepCQlsWoD6O4b240N83YQLeXf2Tpjjhrs/5Nha2IiIiIFNS4C9RqAlkpsOtbftj5FwmpWQWaHU3KYNycrRWiuFVhKyIiIiIF2Wx5a9qa2z5n8uKYQpudnogweXGM5dMSVNiKiIiISOHa3JD7758rcSSdeUTWBOKSMtgUm1g+uc5Aha2IiJRIZGQkISEhhIWFWR1FRMqKfzNoGIZhOulh33bO5vEpGeUQ6sxU2IqISImEh4cTExNDVFSU1VFEpCw1ugSAUCP2nE0DvL3KOs1ZqbAVERERkTMLagdAB/cDGGdoYgBBvl50DvYrr1SFUmErIiIiImcW1BaAENt+bDjPWNxGDA7BbjvT3vKhwlZEREREzsy/GbjXwO7IYPaQWtTzLTjd4Mo2QQwIDbIgXH5uVgcQERERkQrMZod6reHgBi6vcYS1j/yPTbGJxKdksPdYKu/8sIef/jhG0qlsfKu5WxvV0lcXERERkYrv7+kIxEVjtxl0aebP0HYNmND7YpoHepOSkcPMtec+uaysqbAVERERkbPLK2zzL/llsxnc1+ciAD5cG0tSenZ5J8ufx9JXFxEREZGKr3673H/jtoPTmW/XgFb1aFHPm5TMHD5Y+2f5Z/sXFbYiIiIicnZ1moObF2SlwIn8Uw5sNoMJp0dtf97HyfQsKxLmZrHslS2mK+aIiIiIFJHdDQJb5d6Piy6wu19IPVoG+ZCamcOMNdaN2rpsYasr5oiIiIgUwxnm2UL+UdvZP+8jMc2aUVuXLWxFREREpBjOUtgC9AsJpFV9H9KyHJaN2qqwFREREZFzO13YHokG0yyw2zAMJvS5GICP1u0jITWzHMPlUmErIiIiIucWEAI2d8g4CScPFNqkT8sAWjfwJT3LwXQLRm1V2IqIiIjIubl5QkDL3PtnmI6QO2qbO9f243X7OV7Oo7YqbEVERESkaM4xzxbgihYBtG3oy6lsB9N+2ltOwXKpsBURERGRoilCYfvvubafbNhPfEpGeSQDVNiKiIiISFEFtcv9Ny660BPITuvZvC7tGtUiI9vJtJ/Kb66tClsRERERKZrAVmDYIO0YpBw9YzPDMLi/b+6o7ZwN+4lPLp9RWxW2IiIiIlI0HtVzL68LhV6B7N+6X1SHDo1rkZnjZGo5zbVVYSsiIiIiRVe/Xe6/Z5lnC/lHbT/deIC/ymHUVoWtiIiIiBRdEU4gO63bhXXo1KQ2WTlOpq4q+1FbFbYiIiIiUnTFKGz/PWr72aYDHE0q21FbFbYiIiIiUnT1Wuf+m3wYUo+ds3nXZv50bupHVo6T91buZv3eBBZFH2b93gQczjOvrFASbqV6NBERERGp2jy9wf9CSNgDR7fBhX3O2twwDCb0vYibZ2xkzoYDzNnwz+V4g3y9iBgcwoDQoFKJphFbERERESmeYkxHAEg+lV3o9qNJGYybs5WlO+JKJZYKWxEREREpntOF7ZHoczZ1OE0mL44pdN/piQiTF8eUyrQEFbYiIiIiUjz/vgLZOWyKTSTuLCeNmUBcUgabYhPPO5YKWxERKZHIyEhCQkIICwuzOoqIlLfaTXL/TY0/Z9P4lKKthFDUdmejwlZEREokPDycmJgYoqKirI4iIuUtKz33X4+a52wa4O1VpEMWtd3ZqLAVERERkeLJTMn919P7nE07B/sR5HvmotUgd3WEzsF+5x1Lha2IiIiIFE/W6cL23CO2dpvBvb0vLHSf8fe/EYNDsNuMQtsUhwpbERERESmevBFbnyI133E4GQAPe/7itZ6vF1OHdyi1dWx1gQYRERERKZ7M1Nx/izDH9mBiOv+3+SAAH992CSa5J4oFeOdOPyiNkdrTVNiKiIiISPEUY47tuz/uJtthcvlFdbi0mX+ZxtJUBBEREREpnqy/R2zPMcd23/E0vtp6GID7+15c1qlU2IqIiIhIMWXmzpk914jt2z/sxuE06dW8Lh0a1y7zWCpsRURERKR48ubYnrmw3ROfwsLo3NHaiX2bl0cqFbYiIiIiUkxFmGP71ve7MU3oFxJI64a+5RJLha2IiIiIFM855tjuOprMN9vjgPKZW3uaClsRERERKZ5zjNi+ueIPAK5sE0TLoKKtdVsaVNiKiIiISPGcpbDdcTiJZb/9hWHAhN4XlWssFbYiIiIiUjynC9tCTh574+/R2qFt63NR4LnXuS1NKmxFREREpHjOMGK79cAJftwVj91mcF+f8ptbe5oKWxEREREpnjOcPHZ6bu017RsQXKdGeadSYSsiIiIixeDIhpyM3Pv/GrHdFJvImt3HcbMZ3FvOc2tPc7PkVUVERESkcjo9DQFwuNdk094E4lMyeP+nvQD8L6wRjfyqWxJNha2IiIiIFN3fha3D7kW3V1cTl5SRb3do/fJb3uu/NBVBRERERIru7/m1J3I8ChS1AI8v2MHSHXHlnQpw4cI2MjKSkJAQwsLCrI4iIiIiUmk4TiUDkGJWO2ObyYtjcDjN8oqUx2UL2/DwcGJiYoiKirI6ioiIiEilsefPPQCkUPg8WhOIS8pgU2xiOabK5bKFrYiInB/95UvENXnt+wGALc6zr1Mbn1JwmkJZU2ErIiIlor98ibggp4P6f/0EwHJnp7M2DfD2Ko9E+WhVBBEREREpmoMbcc9MJImabHK2KLSJAdTz9aJzsF/5ZkMjtiIiIiJSVLu+BWC9vRMO7AV2G3//GzE4BLvNKLC/rKmwFREREZFzM03Y9Q0AC061o4aHnQBvz3xN6vl6MXV4BwaEBlmRUFMRis00ITvd6hQlk52N3ZEJWWlguludRqTyyaqkn30RkdIQvxNO7CPTdGeNsw3v3tSens0D2BSbSHxKBgHeudMPrBipPU2FbXGYJszqDwc3Wp2kRNyBqwC2WxxEREREKp2UbQvxBtY4Q7nhshb0bhkIQJdm/tYG+xdNRSiO7PRKW9SKSClqdCm4W3MddBERKzicJvFR8wH4zbsbkwYWfuKY1TRiW1IP7gGPyvWDLTs7m2XLltO/fz/c3TUVQaTE3KuDYd2f2kREyttHS3/mtuzdOE2DYTfcjqdbwRPHKgIVtiXlUR08alidoniMbBx2z9zcKmxFRESkCDbFJrLv5y/BHU74taVJk2CrI52RpiKIiIiISKFOpmcxYe4v9LFtBsC/0zUWJzo7FbYiIiIiUoBpmjz85XZSkxLoao/J3djiKmtDnYMKWxEREREp4JMN+1ke8xe93bbhhgPqNAf/ZlbHOivNsRURERERHE4zb03aU1kOnv0md5R2fIM/4C+gxZXWBiwCFbYiIiIiLm7pjjgmL44hLikj3/b29asRfGJ97gMVtiIiIiJSkS3dEce4OVsxC9nnc3QDhkcK1KwH9TuUe7bi0hxbERERERflcJpMXhxTaFEL0Pfv1RCczQeCreKXjRU/oYiIiIiUiU2xiQWmH5xWnQz626MA+N23e3nGKjEVtiIiUiKRkZGEhIQQFhZmdRQRKaH4lMKLWoCJbl9Q10jmoLMuu2u2L8dUJafCVkRESiQ8PJyYmBiioqKsjiIiJRTg7VXo9lDjT261LwXgyZxbqevrU56xSkyFrYiIiIiLSsvMLrDNjoOX3D/Abph87ejC796X0jnYz4J0xafCVkRERMQFLd1xlHGfbs17bPz972j7UkJt+0gyq/Ns9kgiBodgtxmFH6SCUWErIiIi4mIWbztC+GdbyXaYDG5bn/dubk89Xy8acIwH3L4EINJ9FM8O78WA0CCL0xad1rEVERERcSHztx7iwS+24TThmg4NePW6tthtBgNb1SN51jVUP5xJckAYj9z1PHa73eq4xaLCVkRERMRFzIs6wKT5v2KacGNYI164ujW2v6cZ2HcupPbhlWBzx+f6SKhkRS2osBURERFxCZ+s38eTi34DYGSXJjw9uFVeUcupk7B0Uu79yx+Aus2tCXmeXHaOrdZfFBEREVfxwZo/84raMd2CmTzkX0UtwPdPQ+pf4H8RXD7RmpClwGULW62/KCIiIq4gcuUenvt2JwDhvZrx+JUtMYx/FbX718OWD3PvD34L3DzLP2Qp0VQEERERkSrINE3e+n43b/+wG4D7+1zMvb0vzF/U5mTBNxNy77cfAU27lX/QUqTCVkRERKSKMU2TV5b9ztRVewF4ZEALxvVsVrDhz2/DsV1Qoy70faacU5Y+FbYiIiIiVYhpmjz37U5mro0F4MmrQri9W3DBhsf3wOpXc+/3fxGqV46ri52NClsRERGRKsLpNIn4+jc+2bAfgGeHhTLi0iYFG5pm7hQERyY0uwJaX1e+QcuIClsRERGRKsDhNHl8wa/MjTqIYcBL17TmhrDGhTeO/gz2rQG3anDlG2BUjkvmnovLroogIiLnR8smilQcOQ4nD32xjblRB7EZ8Mb/2p65qE07Dssfz73fcxL4FTJNoZJSYSsiIiWiZRNFKoZsh5MJ86KZ/8th7DaDd25qz9XtG575Ccseg1MnIDAUuoSXX9ByoKkIIiIiIpWEw2myKTaR+JQMAry9aNvIlwlzo1ke8xfudoP3bu5A/1b1znyAvT/C9nmAAYPfAbt7uWUvDypsRURERCqBpTvimLw4hrikjLxtnm42MnOceLjZmDa8I71aBJz5AFnp8M39ufc73wkNO5Zx4vKnwlZERESkglu6I45xc7Zi/md7Zo4TgLE9Ljh7UQuw+hU4sQ+868MVT5RJTqtpjq2IiIhIBeZwmkxeHFOgqP23LzYfwuE8S4ujO2Ddu7n3r3wNvHxKNWNFocJWREREpALbFJuYb/pBYeKSMtgUm1j4zqw0+Ho8OHOgxVXQ4soySFkxqLAVERERqcDiU85e1J61XVoCfDQYjmwFTx8Y9Gopp6tYNMdWREREpAJLTMsqUrsAb6/8G07sgznXQsIeqFYbbv4/8Klf+gErEBW2IiIiIhXQibQsXvpuF/M2HzxrOwOo5+tF52C/fzYe/TW3qE39C3wbwfD5UPfisg1cAaiwFREREalATNNk/tbDPL9kZ95obbcL6/DznuO5+//V9vSFcCMGh2C3/f0odjXMvQUykyGgFQz/CnyCyu8NWEiFrYiIiEgF8eexVJ5YuIN1exMAaB7ozfNXh9KpqV+h69jW8/UiYnAIA0L/Llx3zIcFd4EjC5p0gxs/hWq1LHgn1lBhKyIiImKxzBwHU1ftZcrKvWQ5nHi527iv98WMuTwYd3vuuf4DQoPoG1Iv35XHOgf7/TNSu3EafPcIYELLIXDNDHD3OvOLVkEqbEVEREQstG7vcZ5YsIM/j6cB0OPiujw3LJRGftULtLXbDLo088+/0TThh8mw9s3cx2F3wMCXwWYv6+gVjgpbEREpkcjISCIjI3E4HFZHEamUElIzeX7JTuZvPQxAXW9PIgaHcGXrIAzDOMez/+bIzl2jdtvnuY+veBIufwCK+vwqxmULW3XIIiLnJzw8nPDwcJKTk/H19bU6jkilYZomX2w+xAvf7eRkejaGAcMvacJDA5rj4+Ve9ANlpcH/jYI9K8Cww+C3ocOIsgteCbhsYVuiDtm9Ojx25J/7IiIiIsWwJz6FxxbsyLtKWMsgH164OpT2jWsX70Bpx+HT63MvvOBWDf73EVzcvwwSVy4uW9iWiGGARw2rU4iIiEglk5Ht4L0f9zBt9V6yHSbV3O1M7Hsxt17WFDd7MS8Ee2IffHINJO79+8ILX0CjsDLJXdmosBUREREpQ2t2H+OJhTvYn5AOQJ+WATw9pBUNa5fgr79x23JHalP/At/GuWvUusCFF4pKha2IiIhIGTiWkslz38awKDp3GmM9Hy+eHtKK/q0Ci35y2L/9+VPuhReyUlzuwgtFpcJWREREpBQ5nSZzow7y0nc7Sc7IwWbAqK5NeaBfc2p6lrD02vEVzL8LnNkueeGFolJhKyIiIlJKdh1N5vEFO9iy/wQAoQ18ePHqNrRueB4rh2yYCksn5d4PGQpXT3e5Cy8UlQpbERERkfN0KsvB2z/s5oM1f5LjNKnhYeeBfs0Z2aVJ8U8OO8004fun4ee3ch+78IUXikqFrYiIiMh5WLkrnicX7eDQiVMA9G8VyNNDWhHkW63kB/3vhRd6PwXdJrrshReKSoWtiIiISAn8lZzBM4tj+PbXOAAa1KrG5CGt6BMSeH4HzkyFL0bBnu9zL7ww5F1of0spJK76VNiKiIiIFIPDafLpxv28uvR3UjJzsNsMbrusKRP6XEyNkp4cdpouvHBeVNiKiEiJ6NLk4op+O5LEY/N/ZduhJADaNqrFC1eH0qp+KVxWOjEW5lwDiX9CNT+45Qto2On8j+tCVNiKiEiJlOjS5CKVVFpmDm+u+IMP1+3D4TTx9nTj4QHNufmSJthtpTDvNW4bzLkO0uJzL7wwYj7Uuej8j+tiVNiKiIiInMWKmL+IWLSDI0kZAFzZJoinrgoh0KeUltz6cxXMHZ574YXAULjlS114oYRU2IqIiIgU4sjJUzz99W8sj/kLgIa1q/HssFB6NQ8ovRf59UtYMDb3wgtNL8+98IKX/gJSUipsRURERP4lx+Hko/X7eWP576RlOXCzGYy5/ALu630R1TxKcQ3Z9VNg2aO590OGwTXTwc2z9I7vglTYioi4uIMHDzJixAji4+Nxc3PjySef5Prrr7c6lkiZcjhNNsUmEp+SQYC3F52D/bDbDLYfOsljC35lx+FkADo2qc3zV4fSop5P6b240wnfR8C6d3Ifd74LBrwEthJeyEHyqLAVEXFxbm5uvPXWW7Rr1474+Hg6dOjAoEGDqFGjhtXRRMrE0h1xTF4cQ9zfc2YBAn08CQny4ac/juE0wcfLjUkDW3JjWCNspXFy2GmObFh0D2yfm/u4dwR0u18XXiglLlvYapkaEZFcQUFBBAXlnqgSEBCAn58fiYmJKmylSlq6I45xc7Zi/mf7X8mZ/JV8DICh7erzxJUh1PUu5WkBmanwfyNh7w+68EIZcdkx7/DwcGJiYoiKirI6iojIeVm9ejWDBw+mfv36GIbBwoULC7SZMmUKwcHBeHl50bFjR9asWVPosTZv3ozT6aRRo0ZlnFqk/DmcJpMXxxQoav/Nr4YHb/yvXekXtanH4KOrcota9+pw01wVtWXAZUdsTzPN3G/v5ORki5OUvezsbNLT00lOTsbd3d3qOCJVxtk+W6f7ltN9TVlIS0ujbdu23HrrrVx77bUF9s+bN48JEyYwZcoULrvsMqZNm8bAgQOJiYmhcePGee0SEhIYOXIkH3zwwVlfLzMzk8zMzLzHSUm5C9UnJiaSnZ1dSu+qYjv9f56QkKD+tBLZvP8Eh+MTz9rmeGY630fvpVOT2qX3wicP4PblSIyT+zG9auMYMgvTvx0kJJTea1RxKSkpQBH6UtPFHTx40AR000033cr0dvDgwXLp0wBzwYIF+bZ17tzZHDt2bL5tLVq0MCdNmpT3OCMjw7z88svNjz/++JyvERERYfnXUzfddHPN27n6Upcfsa1fvz4HDx7E29sbowgTt8PCwko8faG4zy1O+6K0TU5OplGjRhw8eBAfn1I8u7OSOZ//w7JUXrnK4nVK45hV9bNlmiYpKSnUr1+/yPlKU1ZWFlu2bGHSpEn5tvfr149169YBuRlHjx7NFVdcwYgRI855zEcffZSJEyfmPXY6nSQmJuLv71+gHz3f743y/L4oznPUnxZfRe17z8TKvOXx2qX5GuXxM6CofanLF7Y2m42GDRsWub3dbi9xJ1bc5xanfXHa+vj4uHRHfD7/h2WpvHKVxeuUxjGr8mfLysvNHj9+HIfDQWBgYL7tgYGBHD16FICff/6ZefPm0aZNm7z5uZ988gmtW7cu9Jienp54euaff1irVq1C257v90Z5fl+U5Dmu3p8WR0Xte8/Eyrzl8dql+Rrl9TOgKH2pyxe2xRUeHl5uzy1O+/PJ5Woq6teqvHKVxeuUxjH12Spb/x1JNU0zb1u3bt1wOp1l8rrn+/Urz++L8309ObvK9rW1Mm95vHZpvobVPwP+zfh7Tpa4gOTkZHx9fUlKSqpUvzWLVHQV6bNlGAYLFixg2LBhQO5UhOrVq/PFF19w9dVX57W77777iI6O5qeffrIoaeVWkf7PReQfLrvclyvy9PQkIiKiwJ8QReT8VOTPloeHBx07dmTFihX5tq9YsYKuXbtalKryq8j/5yKuTCO2IiKVXGpqKnv27AGgffv2vPHGG/Tq1Qs/Pz8aN27MvHnzGDFiBO+//z5dunRh+vTpzJgxg99++40mTZpYnF5EpPSosBURqeRWrVpFr169CmwfNWoUs2fPBnIv0PDKK68QFxdHaGgob775Jt27dy/npCIiZUuFrYiIiIhUCZpjKyIiIiJVggpbEREREakSVNiKiIiISJWgwlYAOHjwID179iQkJIQ2bdrwxRdfWB1JpNJLSUkhLCyMdu3a0bp1a2bMmGF1JCkH6k9FrKOTxwSAuLg4/vrrL9q1a0d8fDwdOnTg999/p0aNGlZHE6m0HA4HmZmZVK9enfT0dEJDQ4mKisLf39/qaFKG1J+KWEeX1BUAgoKCCAoKAiAgIAA/Pz8SExPVEYucB7vdTvXq1QHIyMjA4XCgsYSqT/2piHU0FaGKWL16NYMHD6Z+/foYhsHChQsLtJkyZQrBwcF4eXnRsWNH1qxZU+ixNm/ejNPppFGjRmWcWqRiK43P1cmTJ2nbti0NGzbk4Ycfpk6dOuWUXkpK/alI5aXCtopIS0ujbdu2vPfee4XunzdvHhMmTODxxx/nl19+4fLLL2fgwIEcOHAgX7uEhARGjhzJ9OnTyyO2SIVWGp+rWrVqsW3bNmJjY/nss8/466+/yiu+lJD6U5FKzJQqBzAXLFiQb1vnzp3NsWPH5tvWokULc9KkSXmPMzIyzMsvv9z8+OOPyyOmSKVS0s/Vv40dO9b8v//7v7KKKGVA/alI5aIRWxeQlZXFli1b6NevX77t/fr1Y926dQCYpsno0aO54oorGDFihBUxRSqVonyu/vrrL5KTkwFITk5m9erVNG/evNyzSulRfypSsamwdQHHjx/H4XAQGBiYb3tgYCBHjx4F4Oeff2bevHksXLiQdu3a0a5dO3799Vcr4opUCkX5XB06dIju3bvTtm1bunXrxj333EObNm2siCulRP2pSMWmVRFciGEY+R6bppm3rVu3bjidTitiiVRqZ/tcdezYkejoaAtSSVlTfypSMWnE1gXUqVMHu92eN5pwWnx8fIFRBxEpGn2uXJP+30UqNhW2LsDDw4OOHTuyYsWKfNtXrFhB165dLUolUrnpc+Wa9P8uUrFpKkIVkZqayp49e/Iex8bGEh0djZ+fH40bN2bixImMGDGCTp060aVLF6ZPn86BAwcYO3ashalFKjZ9rlyT/t9FKjFL12SQUrNy5UoTKHAbNWpUXpvIyEizSZMmpoeHh9mhQwfzp59+si6wSCWgz5Vr0v+7SOVlmKau7ygiIiIilZ/m2IqIiIhIlaDCVkRERESqBBW2IiIiIlIlqLAVERERkSpBha2IiIiIVAkqbEVERESkSlBhKyIiIiJVggpbEREREakSVNiKiIiISJWgwlZEREREqgQVtiIiIiJSJaiwlUqvZ8+eTJgwoVyff76veSYJCQkEBASwb9++Uj/2mVx33XW88cYb5fZ6IiKFKe++vKz6cVBfbilTpAyNGjXKBEzAtNvtZqNGjcyxY8eaiYmJpfYaCQkJZnJycpHa9ujRw7zvvvtK/PzCnlPYMUvqgQceMG+77bYC27dt22ZeffXVpp+fn+np6WmGhISYr7zyipmdnX3er7lt2zbTz8/PTEpKOu9jiVRlP//8s2mz2cz+/ftbHaVMlGZfVhLl3Zf/t7368qpBI7ZS5gYMGEBcXBz79u3jgw8+YPHixdx9992ldnw/Pz+8vb3L9fnn+5qFOXXqFDNnzmTMmDH5tv/0009ceumlVKtWjUWLFrFt2zYefvhhXnvtNa655hqcTud5vW6bNm1o2rQpn3766XkdR6SqmzVrFuPHj2ft2rUcOHDA6jgAZGVlWR2h1DKUd19eFv04qC+3nNWVtVRto0aNMocOHZpv28SJE00/P7+8x06n03z55ZfN4OBg08vLy2zTpo35xRdf5HvOF198YYaGhppeXl6mn5+f2bt3bzM1NdU0zYK/ZZ+p7b9Hj0/fYmNj8z3//fffN+vXr286HI58rz948GBz5MiReY9PP+dMx/zoo49MPz8/MyMjI99xrrnmGnPEiBGFfq2++uors06dOvm25eTkmMHBweYtt9xSoP1vv/1muru7mx988EG+XPfcc4953333mbVq1TIDAgLMadOmmampqebo0aPNmjVrmhdccIG5ZMmSfMd6+umnzcsvv7zQXCJimqmpqaa3t7e5a9cu84YbbjAnT56cb3+PHj3M8PBwMzw83PT19TX9/PzMxx9/3HQ6ncVq891335mXXXZZ3v4rr7zS3LNnT4Fj3H///aa/v7/ZvXt30zTP3Y/26NHDHD9+vPnQQw+ZtWvXNgMDA82IiIi8/WfqywpT0gxn68dPH/d0X3y2tqXVl/+7fWHHnDx5crH7cdM8/768pP24aaovN03TVGErZeq/he3evXvNkJAQMzAwMG/bY489ZrZo0cJcunSpuXfvXvPDDz80PT09zVWrVpmmaZpHjhwx3dzczDfeeMOMjY01t2/fbkZGRpopKSmmaebvnM7W9uTJk2aXLl3MO+64w4yLizPj4uLMnJycfM9PSEgwPTw8zO+//z4vX2Jiounh4WEuW7Ysb9vp55zpmOnp6aavr6/5f//3f3nPOXbsmOnh4WH++OOPhX6t7rvvPnPAgAH5tq1bt84EzOjo6EKfM3ToULNv3775cnl7e5vPPvus+ccff5jPPvusabPZzIEDB5rTp083//jjD3PcuHGmv7+/mZaWlve8JUuWmJ6engU6cBHJNXPmTLNTp06maZrm4sWLzaZNmxYoWmvWrGned9995q5du8w5c+aY1atXN6dPn16sNl9++aX51VdfmX/88Yf5yy+/mIMHDzZbt26dV6CdPsZDDz1k7tq1y9y5c6dpmufuR3v06GH6+PiYTz/9tPnHH3+YH330kWkYhrl8+XLTNM0z9mWFKUmGc/Xjp4973333nbNtafXl/25f2DFTU1OL3Y+b5vn35SXtx01TfblpqrCVMjZq1CjTbrebNWrUML28vPJ+E37jjTdM08wdBfHy8jLXrVuX73m33367edNNN5mmaZpbtmwxAXPfvn2Fvsa/O6fitD3TtiFDhuSbGzVt2jSzXr16+Tr5fz/nTPOyxo0bZw4cODDv8VtvvWVecMEF+X4Y/tvQoUMLzMn69NNPTaBA53XaQw89ZF588cX5cnXr1i3vcU5OjlmjRo18owtxcXEmYK5fvz5v27Zt2876dRNxdV27djXfeust0zRNMzs726xTp465YsWKvP09evQwW7Zsme/z/cgjj5gtW7YsVpv/io+PNwHz119/zTtGu3bt8rUpSj/6377BNE0zLCzMfOSRR/LlK8oc05JkOFff/O/XL07bs207V1/+3/aFHbO4/bhpnn9fXtJ+3DTVl5um5thKOejVqxfR0dFs3LiR8ePH079/f8aPHw9ATEwMGRkZ9O3bl5o1a+bdPv74Y/bu3QtA27Zt6d27N61bt+b6669nxowZnDhxotDXKk7bM7nlllv46quvyMzMBODTTz/lxhtvxG63F+s4d9xxB8uXL+fw4cMAfPjhh4wePRrDMAptf+rUKby8vPJt8/HxASAxMbHQ55w4cSKvzWlt2rTJu2+32/H396d169Z52wIDAwGIj4/P21atWjUA0tPTi/TeRFzJ77//zqZNm7jxxhsBcHNz44YbbmDWrFn52l166aX5Pt9dunRh9+7dOByOIrfZu3cvN998MxdccAE+Pj4EBwcD5JvT26lTp3yvW5R+FPL3DQBBQUH5+oHiKG6G8u7HoXT68uL241A6fXlJ+nFQXw5a7kvKQY0aNbjwwgtp06YN77zzDpmZmUyePBkgb7L8t99+S3R0dN4tJiaGL7/8Esj9UK9YsYLvvvuOkJAQ3n33XZo3b05sbGyB1ypO2zMZPHgwTqeTb7/9loMHD7JmzRqGDx9e7Pfdvn172rZty8cff8zWrVv59ddfGT169Bnb16lTp0Dn3aVLF9zd3Vm8eHGB9g6Hg+XLl9OtW7d8293d3fM9Ngwj37bTHfK/T1Q43dnWrVu3aG9OxIXMnDmTnJwcGjRogJubG25ubkydOpX58+eXqOA6m8GDB5OQkMCMGTPYuHEjGzduBPKfoFWjRo18zylKPwqF9w0lPWGpuBnKux+H0unLi9uPQ+n05SXpx0F9OaiwFQtERETw2muvceTIEUJCQvD09OTAgQNceOGF+W6NGjXKe45hGFx22WVMnjyZX375BQ8PDxYsWFDo8c/W1sPDI9/oSWGqVavGNddcw6effsrnn3/OxRdfTMeOHc/Y/mzHHDNmDB9++CGzZs2iT58++d7Tf7Vv356YmJh82/z9/bn33nt57rnnOHLkSL59b775JgkJCdx///1nfT9FsWPHDho2bEidOnXO+1giVUlOTg4ff/wxr7/+er6Cbdu2bTRp0iTfGegbNmzI99wNGzZw0UUX5RshPFubhIQEdu7cyRNPPEHv3r1p2bJlkQrnovaj51KU/vF8MpRWP17UrKXVlxenHwf15VZzszqAuJ6ePXvSqlUrXnjhBd577z0efPBB7r//fpxOJ926dSM5OZl169ZRs2ZNRo0axcaNG/nhhx/o168fAQEBbNy4kWPHjtGyZcsCxz5X26ZNm7Jx40b27dtHzZo18fPzKzTjLbfcwuDBg/ntt9/O+Rt+Yce02Wx5x3nwwQeZMWMGH3/88VmP079/fx599FFOnDhB7dq1AUhNTeXee+9lw4YN9OrVi88//5wOHTrw6quv8vjjjzNt2rS8zri4UyX+bc2aNfTr16/Ezxepqr755htOnDjB7bffjq+vb7591113HTNnzuSee+4B4ODBg0ycOJG77rqLrVu38u677/L666/ne87Z2tSuXRt/f3+mT59OUFAQBw4cYNKkSefM6O3tfc5+tCjO1pedb4YWLVqUWj9+pqyFKY2+vDj9OJx/X34+1Jej5b6kbBW23Jdp5k6k9/DwMA8cOGA6nU7z7bffNps3b266u7ubdevWNfv372/+9NNPpmmaZkxMjNm/f3+zbt26pqenp3nxxReb7777bt6x/j3h/1xtf//9d/PSSy81q1WrVugSMafl5OSYQUFBJmDu3bu3QP5/P6ewY/7biBEjCl0ypjCXXnqp+f777+c9joiIyLf8zKhRo0zTNM+4LE9h76VJkybmm2++mW8bYC5YsMA0TdM8deqU6ePjU+AkBBExzauuusocNGhQoftOn+S0ZcsWs0ePHubdd99tjh071vTx8TFr165tTpo0qcDKCedqs2LFCrNly5amp6en2aZNG3PVqlX5Pq9nOsHrXP1oYc8bOnRoXp9imufuy/79Poqb4Vx987+PW5S2pdGX/7f92d5/cfpx0zy/vrwk/bhpqi8/zTBN0yy3KlrEBfXt25eWLVvyzjvvnLPtkiVLePDBB9mxY0eRR0rOV2RkJIsWLWL58uXl8noiVVHPnj1p164db7311nm1kYqpOP04qC+3kqYiiJSRxMREli9fzo8//sh7771XpOcMGjSI3bt3c/jw4WLNjTsf7u7uvPvuu+XyWiIilUlJ+nFQX24lFbYiZaRDhw6cOHGCl19+mebNmxf5effdd18ZpirozjvvLNfXExGpLEraj4P6cqtoKoKIiIiIVAla7ktEREREqgQVtiIiIiJSJaiwFREREZEqQYWtiIiIiFQJKmxFREREpEpQYSsiIiIiVYIKWxERERGpElTYioiIiEiVoMJWRERERKoEFbYiIiIiUiX8P/3azSRk0jn2AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "fig, ax = plt.subplots(ncols=2, figsize=(8, 6)) # two-column figure\n", "#ves.showModel(synthModel, ax=ax[0], label=\"synth\", plot=\"semilogy\", zmax=20)\n", "ves.showModel(ves.model, ax=ax[0], label=\"model\", zmax=50, color=\"C1\")\n", + "ax[0].invert_yaxis()\n", "ves.showData(rhoa, ax=ax[1], label=\"data\", color=\"C0\", marker=\"o\")\n", "out = ves.showData(ves.inv.response, ax=ax[1], label=\"model\", color=\"C1\")" ] @@ -197,20 +214,37 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We are interested in the model uncertaincies and through model covariance\n", - "\n" + "Calculate the model covariances and plot it. This can be usesful to assess the uncertainty of the inversion result." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 57, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "20/11/25 - 12:29:37 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - [0.04565188 1.26457602 0.02521144 0.06352483 3.10585657]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAGmCAYAAABycI2BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA160lEQVR4nO3de3RU5b3/8U9IwgQoGYlzSEJJkXrhIgm5CSRWwFuA1IpBJB5qlCMGKFpKU9axqW0NtEdqtXKttCg9AaEaLXLxGNEgN5FAgVxQRKCWmmgT0AAZoJAgmd8fNvNzyHXYE2a2+/1aa6/V2fPMM9+9q37yfWbP7CCXy+USAACwjE7+LgAAAFxehD8AABZD+AMAYDGEPwAAFkP4AwBgMYQ/AAAWQ/gDAGAxhD8AABZD+AMAYDGEPwAAFkP4AwDQDtu2bdP3vvc99erVS0FBQVq7dm2br9m6dauSkpIUFhamb3/72/rDH/7QZMzq1as1cOBA2Ww2DRw4UGvWrOmA6j0R/gAAtMOZM2c0ePBgLV68uF3jjxw5ovT0dN10000qLS3Vz372M82YMUOrV692jykuLlZmZqaysrJUXl6urKwsTZgwQbt27eqow5AkBXFjHwAAvBMUFKQ1a9borrvuanHMo48+qvXr1+vAgQPufdOmTVN5ebmKi4slSZmZmXI6nXrjjTfcY0aPHq0ePXroxRdf7LD6Qzps5kvQ0NCgf/7zn+revbuCgoL8XQ4AwEsul0unTp1Sr1691KlTxywunzt3TvX19T6Zy+VyNckbm80mm81meO7i4mKlpaV57Bs1apSWLVum8+fPKzQ0VMXFxfrxj3/cZMz8+fMNv39rAir8//nPfyomJsbfZQAADKqsrFTv3r19Pu+5c+fUt0sXVftovm984xs6ffq0x77HH39ceXl5hueurq5WZGSkx77IyEh98cUX+vzzzxUdHd3imOpqXx1h8wIq/Lt37y5JqpQU7t9SAt6iX9f6uwRT+OFzcf4uwRwSE/1dgTmMHu3vCgKe8+xZxcyc6f7vua/V19erWr7JCaekmNOnVVlZqfDw/z+bL7r+RhevKjR+0v7V/c2N6ejV74AK/8aDDRfh35awMM5Qe4R30LLj105oqL8rMIcuXfxdgWl0dHiFBwUp3Oh7uFySy6Xw8HCP8PeVqKioJh38sWPHFBISoiuvvLLVMRevBvga/2UEAJhPp06+2TpQSkqKioqKPPa99dZbSk5OVui//+BuaUxqamqH1hZQnT8AAIHq9OnT+tvf/uZ+fOTIEZWVlSkiIkLf+ta3lJubq08//VQrVqyQ9OWV/YsXL1ZOTo6ys7NVXFysZcuWeVzF/6Mf/UjDhw/Xk08+qbFjx2rdunXauHGjtm/f3qHHQucPADAfP3T+e/bsUUJCghISEiRJOTk5SkhI0C9/+UtJUlVVlSoqKtzj+/btq8LCQm3ZskXx8fH61a9+pYULF+ruu+92j0lNTdVLL72k//3f/1VcXJzy8/NVUFCgoUOH+uAktYzOHwBgPp06Sb74zN8LI0eOVGs/jZOfn99k34gRI1RSUtLqvOPHj9f48eO9qsUowh8AYD5+CP+vE5b9AQCwGDp/AID50PkbQvgDAMyH8DeEZX8AACyGzh8AYD50/oYQ/gAA8yH8DWHZHwAAi6HzBwCYT1CQ8d/mb2jwTS0mRPgDAMznMtyY5+uMMwcAgMXQ+QMAzIfO3xDCHwBgPoS/IYQ/AMB8CH9DOHMAAFgMnT8AwHzo/A0h/AEA5kP4G8KZAwDAYuj8AQDmQ+dvCOEPADAfwt8QzhwAABZD5w8AMB9f3NjHwrf0JfwBAObji2V/C4c/y/4AAFgMnT8AwHzo/A3xWef/m9/8RikpKb6aDgCAljWGv9HNonzW+ZeXl2vw4MG+mg4AgJbR+Rvisz97ysvLFR8f76vpAABAB7mk8D9w4IBGjhypLl26KCEhQXv27NGhQ4fo/AEAlwfL/oZ4vez/4YcfaujQoZoyZYqWLVumffv2aezYsXK5XIqLi/Nqrrq6OtXV1bkfO51Ob8sBAFgRy/6GeH3mHn74YY0dO1ZPP/20rr76amVkZGj48OG65ppr1K1bN2VkZKhHjx4aP358m3PNnTtXdrvdvcXExFzSQQAAgPbzKvw//vhjbdq0STk5OR77Q0ND3Uv+M2bM0IoVK9o1X25urmpra91bZWWlN+UAAKyKZX9DvDrysrIyhYSEKDY21mN/SUmJ+2K/m2++Wd27d2/XfDabTeHh4R4bAABtIvwN8erIO3XqpIaGBtXX17v3FRYWav/+/VzpDwCASXgV/klJSQoNDdWsWbP097//XevXr1d2drYkcaU/AODyabyxj5EtKMjfR+E3XoV/r1699Pzzz+u1115TcnKyFixYoAceeEAOh0Pf/OY3O6pGAAA8+WnZ/9lnn1Xfvn0VFhampKQkvfPOOy2OnTRpkoKCgpps119/vXtMfn5+s2POnTt3Saelvbw+8vvuu0+VlZU6fvy43n77bT3xxBP67LPPOqI2AAACRkFBgWbOnKnHHntMpaWluummmzRmzBhVVFQ0O37BggWqqqpyb5WVlYqIiNA999zjMS48PNxjXFVVlcLCwjr0WHx+Y59Ro0appKREZ86cUe/evbVmzRrdcMMNvn4bAICV+eKCPS9f/8wzz2jy5Ml66KGHJEnz58/Xm2++qSVLlmju3LlNxjd+jb3R2rVrdeLECf3Xf/2Xx7igoCBFRUVdwgFcOp+H/5tvvunrKQEA8OTD8L/4B+ZsNptsNpvHvvr6eu3du1c//elPPfanpaVpx44d7Xq7ZcuW6bbbblOfPn089p8+fVp9+vTRhQsXFB8fr1/96ldKSEjw9mi8Yt3vOQAAzMuHn/nHxMR4/OBcc138559/rgsXLigyMtJjf2RkpKqrq9sst6qqSm+88YZ71aBR//79lZ+fr/Xr1+vFF19UWFiYbrzxRh0+fNjAyWmbzzt/AADMpLKy0uN3Zi7u+r8q6KJvCLhcrib7mpOfn68rrrhCd911l8f+YcOGadiwYe7HN954oxITE7Vo0SItXLiwnUfgPcIfAGA+Plz2b8+PzDkcDgUHBzfp8o8dO9ZkNeBiLpdLf/rTn5SVlaXOnTu3UVIn3XDDDR3e+bPsDwAwn8v8Vb/OnTsrKSlJRUVFHvuLioqUmpra6mu3bt2qv/3tb5o8eXKb7+NyuVRWVqbo6Oh213Yp6PwBAGiHnJwcZWVlKTk5WSkpKVq6dKkqKio0bdo0SV/er+bTTz9tcn+bZcuWaejQoRo0aFCTOWfPnq1hw4bp2muvldPp1MKFC1VWVqbf//73HXoshD8AwHz88FW/zMxM1dTUaM6cOaqqqtKgQYNUWFjovnq/qqqqyXf+a2trtXr1ai1YsKDZOU+ePKkpU6aourpadrtdCQkJ2rZtm4YMGXJpx9ROhD8AwHz8EP6SNH36dE2fPr3Z5/Lz85vss9vt+te//tXifPPmzdO8efO8rsMoPvMHAMBi6PwBAObTeGMfo3NYFOEPADAfPy37f11Y98gBALAoOn8AgPnQ+RtC+AMAzIfwN4TwBwCYD+FviHWPHAAAi6LzBwCYD52/IYQ/AMB8CH9DrHvkAABYFJ0/AMB86PwNIfwBAOZD+Bti3SMHAMCi6PwBAObDjX0MIfwBAObDsr8h1j1yAAAsis4fAGA+dP6GEP4AAPMh/A0h/AEA5kP4G2LdIwcAwKICsvNf9OtahYWF+7uMgPaTWdb9ioo39u5x+bsEU3jwQX9XYA4fvOrvCgKfy+WUNLXj34jO35CADH8AAFpF+Bti3SMHAMCi6PwBAOZD528I4Q8AMB/C3xDrHjkAABZF5w8AMB9u7GMI4Q8AMB+W/Q2x7pEDAGBRdP4AAPOh8zeE8AcAmA/hbwjhDwAwH8LfEOseOQAAFkXnDwAwHzp/Q6x75AAA82oMf6Obl5599ln17dtXYWFhSkpK0jvvvNPi2C1btigoKKjJ9uGHH3qMW716tQYOHCibzaaBAwdqzZo1XtflLcIfAIB2KCgo0MyZM/XYY4+ptLRUN910k8aMGaOKiopWX3fw4EFVVVW5t2uvvdb9XHFxsTIzM5WVlaXy8nJlZWVpwoQJ2rVrV4ceC+EPADAfP3T+zzzzjCZPnqyHHnpIAwYM0Pz58xUTE6MlS5a0+rqePXsqKirKvQUHB7ufmz9/vm6//Xbl5uaqf//+ys3N1a233qr58+dfyllpN8IfAGA+lzn86+vrtXfvXqWlpXnsT0tL044dO1p9bUJCgqKjo3Xrrbdq8+bNHs8VFxc3mXPUqFFtzmkUF/wBACzN6XR6PLbZbLLZbB77Pv/8c124cEGRkZEe+yMjI1VdXd3svNHR0Vq6dKmSkpJUV1enF154Qbfeequ2bNmi4cOHS5Kqq6u9mtNXCH8AgPn48MY+MTExHrsff/xx5eXltfASz5sBuVyuJvsa9evXT/369XM/TklJUWVlpZ5++ml3+Hs7p68Q/gAA8/HhV/0qKysVHh7u3n1x1y9JDodDwcHBTTryY8eONencWzNs2DCtXLnS/TgqKsrwnJeCz/wBAJYWHh7usTUX/p07d1ZSUpKKioo89hcVFSk1NbXd71VaWqro6Gj345SUlCZzvvXWW17NeSno/AEA5uOHH/nJyclRVlaWkpOTlZKSoqVLl6qiokLTpk2TJOXm5urTTz/VihUrJH15Jf9VV12l66+/XvX19Vq5cqVWr16t1atXu+f80Y9+pOHDh+vJJ5/U2LFjtW7dOm3cuFHbt283dmxtIPwBAObjh/DPzMxUTU2N5syZo6qqKg0aNEiFhYXq06ePJKmqqsrjO//19fWaNWuWPv30U3Xp0kXXX3+9Xn/9daWnp7vHpKam6qWXXtLPf/5z/eIXv9DVV1+tgoICDR061NixtSHI5XK5OvQdvOB0OmW32/XrX9cqLCy87RdY2E9mdezFIF8Xe/cEzD/eAe3BB/1dgTl88IG/Kwh8LpdTFy7YVVtb6/E5uq805kRtXp7Cw8KMzXXunOx5eR1WayDjM38AACyGZX8AgPlwYx9DCH8AgPkQ/oZY98gBALAoOn8AgPnQ+RtC+AMAzIfwN8S6Rw4AgEX5LPx/85vfKCUlxVfTAQDQssYb+xjZOvjmOYHMZ8v+5eXlGjx4sK+mAwCgZSz7G+KzIy8vL1d8fLyvpgMAAB3kksL/wIEDGjlypLp06aKEhATt2bNHhw4dovMHAFweRpf8fbFyYGJeL/t/+OGHGjp0qKZMmaJly5Zp3759Gjt2rFwul+Li4ryaq66uTnV1de7HTqfT23IAAFbEsr8hXh/5ww8/rLFjx+rpp5/W1VdfrYyMDA0fPlzXXHONjh8/rpEjR2rgwIGKi4vTK6+80upcc+fOld1ud28xMTGXfCAAAAuh8zfEq87/448/1qZNm1RSUuKxPzQ0VIMHD1ZISIjmz5+v+Ph4HTt2TImJiUpPT1e3bt2anS83N1c5OTnux06nkz8AAADoYF6Ff1lZmUJCQhQbG+uxv6SkRBMnTlR0dLSio6MlST179lRERISOHz/eYvjbbDbZbLZLLB0AYFks+xviVfh36tRJDQ0Nqq+vV0jIly8tLCzU/v37m1zpv2fPHjU0NNDJAwB8j/A3xKsjT0pKUmhoqGbNmqW///3vWr9+vbKzsyXJ40r/mpoa3X///Vq6dKlvqwUAAIZ5Ff69evXS888/r9dee03JyclasGCBHnjgATkcDn3zm9+U9OUV/BkZGcrNzVVqamqHFA0AsDgu+DPE66/63Xfffbrvvvs89j3xxBOSJJfLpUmTJumWW25RVlaWbyoEAOBiLPsb4tMjf/fdd1VQUKC1a9cqPj5e8fHxeu+993z5FgAAwCCf3tL3O9/5jhoaGnw5JQAATTXe2MfoHBbl0/AHAOCyYNnfEOseOQAAFkXnDwAwHzp/Qwh/AID5EP6GEP4AAPMh/A2x7pEDAGBRdP4AAPOh8zeE8AcAmA/hb4h1jxwAAIui8wcAmA+dvyGEPwDAfAh/Q6x75AAAWBSdPwDAfOj8DSH8AQDmw139DLHunz0AAFgU4Q8AMJ/GZX+jm5eeffZZ9e3bV2FhYUpKStI777zT4thXX31Vt99+u/7jP/5D4eHhSklJ0ZtvvukxJj8/X0FBQU22c+fOeV2bNwh/AID5+CH8CwoKNHPmTD322GMqLS3VTTfdpDFjxqiioqLZ8du2bdPtt9+uwsJC7d27VzfffLO+973vqbS01GNceHi4qqqqPLawsLBLPjXtwWf+AADz8cMFf88884wmT56shx56SJI0f/58vfnmm1qyZInmzp3bZPz8+fM9Hj/xxBNat26dXnvtNSUkJLj3BwUFKSoqyvv6DaDzBwBYmtPp9Njq6uqajKmvr9fevXuVlpbmsT8tLU07duxo1/s0NDTo1KlTioiI8Nh/+vRp9enTR71799Ydd9zRZGWgIxD+AADz8eGyf0xMjOx2u3trrov//PPPdeHCBUVGRnrsj4yMVHV1dbtK/t3vfqczZ85owoQJ7n39+/dXfn6+1q9frxdffFFhYWG68cYbdfjwYQMnp20s+wMAzMeHy/6VlZUKDw9377bZbC2+JOiirwe6XK4m+5rz4osvKi8vT+vWrVPPnj3d+4cNG6Zhw4a5H994441KTEzUokWLtHDhwnYfircIfwCApYWHh3uEf3McDoeCg4ObdPnHjh1rshpwsYKCAk2ePFmvvPKKbrvttlbHdurUSTfccEOHd/4s+wMAzOcyX+3fuXNnJSUlqaioyGN/UVGRUlNTW3zdiy++qEmTJunPf/6zvvvd77b5Pi6XS2VlZYqOjm53bZeCzh8AYD5+uNo/JydHWVlZSk5OVkpKipYuXaqKigpNmzZNkpSbm6tPP/1UK1askPRl8N9///1asGCBhg0b5l416NKli+x2uyRp9uzZGjZsmK699lo5nU4tXLhQZWVl+v3vf2/s2NpA+AMA0A6ZmZmqqanRnDlzVFVVpUGDBqmwsFB9+vSRJFVVVXl85/+Pf/yjvvjiCz388MN6+OGH3fsfeOAB5efnS5JOnjypKVOmqLq6Wna7XQkJCdq2bZuGDBnSoccS5HK5XB36Dl5wOp2y2+369a9rFRbW+ucvVveTWdb9TWpv7N0TMP94B7QHH/R3BebwwQf+riDwuVxOXbhgV21tbZufo1+KxpyofftthXfrZmyuM2dkv/XWDqs1kNH5AwDMhxv7GMIFfwAAWExAdv4/fC5O4Ra+z3J7sJzdPknJ1v3L3hvlr77q7xLMYeBAf1cQ8JynT8uefBneyA8X/H2dBGT4AwDQKsLfEMIfAGA+hL8h1j1yAAAsis4fAGA+dP6GEP4AAPMh/A2x7pEDAGBRdP4AAPOh8zeE8AcAmA/hb4h1jxwAAIui8wcAmA+dvyGEPwDAfLixjyHW/bMHAACLovMHAJgPy/6GEP4AAPMh/A0h/AEA5kP4G2LdIwcAwKLo/AEA5kPnbwjhDwAwH8LfEOseOQAAFkXnDwAwHzp/Qwh/AID5EP6GWPfIAQCwKDp/AID50PkbQvgDAMyHG/sYYt0/ewAAsCg6fwCA+bDsbwjhDwAwH8LfEMIfAGA+hL8h1j1yAAAsis4fAGA+dP6GEP4AAPMh/A2x7pEDAGBRhD8AwHwaO3+jm5eeffZZ9e3bV2FhYUpKStI777zT6vitW7cqKSlJYWFh+va3v60//OEPTcasXr1aAwcOlM1m08CBA7VmzRqv6/LWJYf/9u3bNWTIEIWFhcnhcGjevHm+rAsAgJb5IfwLCgo0c+ZMPfbYYyotLdVNN92kMWPGqKKiotnxR44cUXp6um666SaVlpbqZz/7mWbMmKHVq1e7xxQXFyszM1NZWVkqLy9XVlaWJkyYoF27dhk6PW0JcrlcLm9fVFhYqAceeEBPPfWUUlNTtXz5cj3xxBP629/+pquvvvqSi3E6nbLb7art00fhFv4spj32vvJ3f5dgCknJ1v35Tq+8+qq/KzCHgQP9XUHAc54+LXtysmpraxUeHu77+Rtz4vhxw/M7nU7ZIyLaXevQoUOVmJioJUuWuPcNGDBAd911l+bOndtk/KOPPqr169frwIED7n3Tpk1TeXm5iouLJUmZmZlyOp1644033GNGjx6tHj166MUXXzRyeK3yOmHPnTunqVOnasGCBZo0aZKuu+46zZkzR927d9fWrVuVkZGhHj16aPz48W3OVVdXJ6fT6bEBANAmH3b+F+dQXV1dk7err6/X3r17lZaW5rE/LS1NO3bsaLbE4uLiJuNHjRqlPXv26Pz5862OaWlOX/E6/Ddt2qSzZ88qMzPTvS84OFghISGy2WyaMWOGVqxY0a655s6dK7vd7t5iYmK8LQcAYEEuBflkk6SYmBiPLGqui//888914cIFRUZGeuyPjIxUdXV1szVWV1c3O/6LL77Q559/3uqYlub0Fa+/6rd582YNHjxYwcHB7n1HjhzRiRMnlJiYqAEDBmjLli3tmis3N1c5OTnux06nkz8AAACXVWVlpceyv81ma3Fs0EV3AnS5XE32tTX+4v3ezukLXod/aWmp6uvrPfYtWrTIHfzesNlsrZ5kAACa09Dw5WZ0DkkKDw9v8zN/h8Oh4ODgJh35sWPHmnTujaKiopodHxISoiuvvLLVMS3N6SteL/uXlZVp//79Wr58uQ4dOqRf//rXWrJkiZ577rmOqA8AgCYaw9/o1l6dO3dWUlKSioqKPPYXFRUpNTW12dekpKQ0Gf/WW28pOTlZoaGhrY5paU5f8arzr6ioUE1NjV5//XU9+uijOnTokOLi4rRhwwYlJiZ2VI0AAHjwZeffXjk5OcrKylJycrJSUlK0dOlSVVRUaNq0aZK+/Cj7008/dV/3Nm3aNC1evFg5OTnKzs5WcXGxli1b5nEV/49+9CMNHz5cTz75pMaOHat169Zp48aN2r59u7GDa4NX4V9WVqaIiAilp6crPT29o2oCACDgZGZmqqamRnPmzFFVVZUGDRqkwsJC9enTR5JUVVXl8Z3/vn37qrCwUD/+8Y/1+9//Xr169dLChQt19913u8ekpqbqpZde0s9//nP94he/0NVXX62CggINHTq0Q4/Fq/AvLS1VbGxsq2NGjRqlkpISnTlzRr1799aaNWt0ww03GCoSAICv8kfnL0nTp0/X9OnTm30uPz+/yb4RI0aopKSk1TnHjx/frq/H+5LX4R8XF9fqmDfffNNQQQAAtMVf4f914VX4r127toPKAAAAlwu39AUAmA6dvzGEPwDAdAh/Y7h7DgAAFkPnDwAwHTp/Ywh/AIDpuFzGw9v7G9p/fbDsDwCAxdD5AwBMh2V/Ywh/AIDpEP7GEP4AANMh/I3hM38AACyGzh8AYDp0/sYQ/gAA0yH8jWHZHwAAi6HzBwCYDp2/MYQ/AMB0CH9jWPYHAMBi6PwBAKZD528M4Q8AMB1u7GMMy/4AAFgMnT8AwHRY9jeG8AcAmA7hbwzhDwAwHcLfGD7zBwDAYuj8AQCmQ+dvDOEPADAdwt8Ylv0BALAYOn8AgOnQ+RtD+AMATIfwNyYwwz8xUQoN9XcVAe3BB/1dgTmUv/qqv0swh3Hj/F2BOYwZ4+8KAt/58/6uAO0QmOEPAEAr6PyNIfwBAKbDjX2M4Wp/AAAshs4fAGA6LPsbQ+cPADCdxvA3unWUEydOKCsrS3a7XXa7XVlZWTp58mSL48+fP69HH31UsbGx6tatm3r16qX7779f//znPz3GjRw5UkFBQR7bvffe63V9hD8AwHQCPfwnTpyosrIybdiwQRs2bFBZWZmysrJaHP+vf/1LJSUl+sUvfqGSkhK9+uqrOnTokO68884mY7Ozs1VVVeXe/vjHP3pdH8v+AAD40IEDB7Rhwwbt3LlTQ4cOlSQ999xzSklJ0cGDB9WvX78mr7Hb7SoqKvLYt2jRIg0ZMkQVFRX61re+5d7ftWtXRUVFGaqRzh8AYDq+7PydTqfHVldXZ6i24uJi2e12d/BL0rBhw2S327Vjx452z1NbW6ugoCBdccUVHvtXrVolh8Oh66+/XrNmzdKpU6e8rpHOHwBgOr684C8mJsZj/+OPP668vLxLnre6ulo9e/Zssr9nz56qrq5u1xznzp3TT3/6U02cOFHh4eHu/d///vfVt29fRUVF6f3331dubq7Ky8ubrBq0hfAHAFhaZWWlR8DabLZmx+Xl5Wn27NmtzrV7925JUlBQUJPnXC5Xs/svdv78ed17771qaGjQs88+6/Fcdna2+38PGjRI1157rZKTk1VSUqLExMQ2525E+AMATMeXnX94eLhH+LfkkUceafPK+quuukr79u3T0aNHmzz32WefKTIystXXnz9/XhMmTNCRI0e0adOmNutKTExUaGioDh8+TPgDAL7e/PE9f4fDIYfD0ea4lJQU1dbW6q9//auGDBkiSdq1a5dqa2uVmpra4usag//w4cPavHmzrrzyyjbfa//+/Tp//ryio6PbfyDigj8AAHxqwIABGj16tLKzs7Vz507t3LlT2dnZuuOOOzyu9O/fv7/WrFkjSfriiy80fvx47dmzR6tWrdKFCxdUXV2t6upq1dfXS5I++ugjzZkzR3v27NE//vEPFRYW6p577lFCQoJuvPFGr2qk8wcAmE6g/8LfqlWrNGPGDKWlpUmS7rzzTi1evNhjzMGDB1VbWytJ+uSTT7R+/XpJUnx8vMe4zZs3a+TIkercubPefvttLViwQKdPn1ZMTIy++93v6vHHH1dwcLBX9RH+AADTCfQb+0RERGjlypVtvP//L+Cqq67yeNycmJgYbd261Sf1Ef4AANMJ9M4/0PGZPwAAFkPnDwAwHTp/Ywh/AIDpEP7GsOwPAIDF0PkDAEyHzt8Ywh8AYDqEvzEs+wMAYDF0/gAA06HzN4bwBwCYDuFvDMv+AABYDJ0/AMB06PyNIfwBAKZD+BtD+AMATCfQ7+oX6PjMHwAAi6HzBwCYDsv+xhD+AADTIfyNYdkfAACLueTw3759u4YMGaKwsDA5HA7NmzfPl3UBANCixs7f6GZVlxT+hYWFysjI0PTp07Vv3z5NnTpVOTk5+uijj3xdHwAATRD+xngd/ufOndPUqVO1YMECTZo0Sdddd53mzJmj7t27a+XKlRo5cqQGDhyouLg4vfLKK63OVVdXJ6fT6bEBAICO5fUFf5s2bdLZs2eVmZnp3hccHKyQkBBdc801Gjt2rOLj43Xs2DElJiYqPT1d3bp1a3auuXPnavbs2ZdePQDAkrjgzxivO//Nmzdr8ODBCg4Odu87cuSITpw4ocTERMXHx0uSevbsqYiICB0/frzFuXJzc1VbW+veKisrvT8CAIDlsOxvjNedf2lpqerr6z32LVq0SImJiRowYIB73549e9TQ0KCYmJgW57LZbLLZbN6WAAAADPA6/MvKytTQ0KDly5crJSVFL7/8spYsWaJ3333XPaampkb333+/nn/+eZ8WCwCAxLK/UV6Ff0VFhWpqavT666/r0Ucf1aFDhxQXF6cNGzYoMTFR0pcX8WVkZCg3N1epqakdUjQAwNoIf2O8Cv+ysjJFREQoPT1d6enpTZ53uVyaNGmSbrnlFmVlZfmsSAAAvoob+xjj1QV/paWlio2NbfH5d999VwUFBVq7dq3i4+MVHx+v9957z3CRAADAd7zq/EtLSxUXF9fi89/5znfUYOV1FADAZcGyvzFehf/atWs7qAwAANqP8DeGG/sAAGAx3NIXAGA6dP7GEP4AANMh/I1h2R8AAIuh8wcAmA6dvzF0/gAA0wn0G/ucOHFCWVlZstvtstvtysrK0smTJ1t9zaRJkxQUFOSxDRs2zGNMXV2dfvjDH8rhcKhbt26688479cknn3hdH+EPAICPTZw4UWVlZdqwYYM2bNigsrKydv3y7ejRo1VVVeXeCgsLPZ6fOXOm1qxZo5deeknbt2/X6dOndccdd+jChQte1ceyPwDAdAJ52f/AgQPasGGDdu7cqaFDh0qSnnvuOaWkpOjgwYPq169fi6+12WyKiopq9rna2lotW7ZML7zwgm677TZJ0sqVKxUTE6ONGzdq1KhR7a6Rzh8AYDq+XPZ3Op0eW11dnaHaiouLZbfb3cEvScOGDZPdbteOHTtafe2WLVvUs2dPXXfddcrOztaxY8fcz+3du1fnz59XWlqae1+vXr00aNCgNue9GOEPADCdxhv7GNkab+wTExPj/mzebrdr7ty5hmqrrq5Wz549m+zv2bOnqqurW3zdmDFjtGrVKm3atEm/+93vtHv3bt1yyy3uP0aqq6vVuXNn9ejRw+N1kZGRrc7bHJb9AQCWVllZqfDwcPdjm83W7Li8vDzNnj271bl2794tSQoKCmrynMvlanZ/o8zMTPf/HjRokJKTk9WnTx+9/vrrGjduXIuva2ve5hD+AADT8eVn/uHh4R7h35JHHnlE9957b6tjrrrqKu3bt09Hjx5t8txnn32myMjIdtcXHR2tPn366PDhw5KkqKgo1dfX68SJEx7d/7Fjx5SamtrueSXCHwBgQv644M/hcMjhcLQ5LiUlRbW1tfrrX/+qIUOGSJJ27dql2tpar0K6pqZGlZWVio6OliQlJSUpNDRURUVFmjBhgiSpqqpK77//vn772996dSx85g8AgA8NGDBAo0ePVnZ2tnbu3KmdO3cqOztbd9xxh8eV/v3799eaNWskSadPn9asWbNUXFysf/zjH9qyZYu+973vyeFwKCMjQ5Jkt9s1efJk/eQnP9Hbb7+t0tJS3XfffYqNjXVf/d9edP4AANMJ5K/6SdKqVas0Y8YM95X5d955pxYvXuwx5uDBg6qtrZUkBQcH67333tOKFSt08uRJRUdH6+abb1ZBQYG6d+/ufs28efMUEhKiCRMm6OzZs7r11luVn5+v4OBgr+oj/AEAphPo4R8REaGVK1e2OsbV+HUDSV26dNGbb77Z5rxhYWFatGiRFi1aZKg+lv0BALAYOn8AgOkEeucf6Ah/AIDpEP7GsOwPAIDF0PkDAEyHzt8Ywh8AYDqEvzGEPwDAdBpv7GN0DqviM38AACyGzh8AYDos+xtD+AMATIfwN4ZlfwAALIbOHwBgOnT+xgRm+I8eLXXp4u8qAtoHr/q7ApMYONDfFZjDmDH+rsAc3njD3xXg3wh/Y1j2BwDAYgKz8wcAoBV0/sYQ/gAA0yH8jWHZHwAAi6HzBwCYDp2/MYQ/AMB0CH9jCH8AgOlwYx9j+MwfAACLofMHAJgOy/7GEP4AANMh/I1h2R8AAIuh8wcAmA6dvzGEPwDAdAh/Y1j2BwDAYuj8AQCmQ+dvDOEPADAdwt8Ylv0BALAYOn8AgOnQ+RtD+AMATIfwN4bwBwCYDjf2MYbP/AEAsBg6fwCA6bDsbwydPwDAdBrD3+jWUU6cOKGsrCzZ7XbZ7XZlZWXp5MmTrb4mKCio2e2pp55yjxk5cmST5++9916v66PzBwDAxyZOnKhPPvlEGzZskCRNmTJFWVlZeu2111p8TVVVlcfjN954Q5MnT9bdd9/tsT87O1tz5sxxP+7SpYvX9RH+AADTCeRl/wMHDmjDhg3auXOnhg4dKkl67rnnlJKSooMHD6pfv37Nvi4qKsrj8bp163TzzTfr29/+tsf+rl27NhnrLZb9AQCmE8jL/sXFxbLb7e7gl6Rhw4bJbrdrx44d7Zrj6NGjev311zV58uQmz61atUoOh0PXX3+9Zs2apVOnTnldI50/AMDSnE6nx2ObzSabzXbJ81VXV6tnz55N9vfs2VPV1dXtmmP58uXq3r27xo0b57H/+9//vvr27auoqCi9//77ys3NVXl5uYqKiryqkc4fAGA6vuz8Y2Ji3Bfm2e12zZ07t9n3zMvLa/GivMZtz549kr68eO9iLper2f3N+dOf/qTvf//7CgsL89ifnZ2t2267TYMGDdK9996rv/zlL9q4caNKSkq8OHt0/gAAE/LlZ/6VlZUKDw9372+p63/kkUfavLL+qquu0r59+3T06NEmz3322WeKjIxss6533nlHBw8eVEFBQZtjExMTFRoaqsOHDysxMbHN8Y0IfwCApYWHh3uEf0scDoccDkeb41JSUlRbW6u//vWvGjJkiCRp165dqq2tVWpqapuvX7ZsmZKSkjR48OA2x+7fv1/nz59XdHR0m2O/imV/AIDpBPIFfwMGDNDo0aOVnZ2tnTt3aufOncrOztYdd9zhcaV///79tWbNGo/XOp1OvfLKK3rooYeazPvRRx9pzpw52rNnj/7xj3+osLBQ99xzjxISEnTjjTd6VSPhDwAwnUAOf+nLK/JjY2OVlpamtLQ0xcXF6YUXXvAYc/DgQdXW1nrse+mll+RyufSf//mfTebs3Lmz3n77bY0aNUr9+vXTjBkzlJaWpo0bNyo4ONir+lj2BwCYTqDf2CciIkIrV65s4/2bFjBlyhRNmTKl2fExMTHaunWrT+q75M5/+/btGjJkiMLCwuRwODRv3jyfFAQAADrWJYV/YWGhMjIyNH36dO3bt09Tp05VTk6OPvroI1/XBwBAE4G+7B/ovA7/c+fOaerUqVqwYIEmTZqk6667TnPmzFH37t31f//3f7rhhhsUHx+v2NhYPffcc63OVVdXJ6fT6bEBANAWwt8Yrz/z37Rpk86ePavMzEz3vuDgYIWEhMjhcGjr1q3q2rWr/vWvf2nQoEEaN26crrzyymbnmjt3rmbPnn3p1QMAAK953flv3rxZgwcP9riy8MiRIzpx4oQSExPVtWtXSV+uEFy4cKHZCxoa5ebmqra21r1VVlZewiEAAKyGzt8Yrzv/0tJS1dfXe+xbtGiREhMTNWDAAJ08eVIjRozQ4cOH9dRTT7X6gwhGfz8ZAGBNDQ1SO38pt9U5rMrr8C8rK1NDQ4OWL1+ulJQUvfzyy1qyZIneffddSdIVV1yh8vJyHT16VOPGjdP48ePb9XOGAADg8vBq2b+iokI1NTVauXKlnn76acXGxmrdunXasGFDk98UjoyMVFxcnLZt2+bTggEAYNnfGK86/7KyMkVERCg9PV3p6elNnj969Ki6dOmi8PBwOZ1Obdu2TT/4wQ98ViwAABLL/kZ5Ff6lpaWKjY1t8flPPvlEkydPlsvlksvl0iOPPKK4uDjDRQIAAN/xOvxbC/OkpCSVlZUZrQkAgFbR+RvjVfivXbu2g8oAAKD9CH9juLEPAMB0Av3GPoGOW/oCAGAxdP4AANPxxZI9y/4AAJgI4W8My/4AAFgMnT8AwHTo/I0h/AEApkP4G8OyPwAAFkPnDwAwHTp/Ywh/AIDpEP7GsOwPAIDF0PkDAEyHzt8Ywh8AYDqEvzGEPwDAdLixjzF85g8AgMXQ+QMATKehQQoKMjaHlTt/wh8AYDqEvzEs+wMAYDF0/gAA06HzN4bwBwCYDuFvDMv+AABYDJ0/AMB06PyNIfwBAKZD+BvDsj8AABZD+AMATKehwTdbR/mf//kfpaamqmvXrrriiiva9RqXy6W8vDz16tVLXbp00ciRI7V//36PMXV1dfrhD38oh8Ohbt266c4779Qnn3zidX2EPwDAdAI9/Ovr63XPPffoBz/4Qbtf89vf/lbPPPOMFi9erN27dysqKkq33367Tp065R4zc+ZMrVmzRi+99JK2b9+u06dP64477tCFCxe8qo/P/AEApuNyBfZn9rNnz5Yk5efnt2u8y+XS/Pnz9dhjj2ncuHGSpOXLlysyMlJ//vOfNXXqVNXW1mrZsmV64YUXdNttt0mSVq5cqZiYGG3cuFGjRo1qd30BFf6uf/8/6Tx71s+VBD6Xy+nvEkzBefq0v0swh/Pn/V0BviYa/8vk6vBk9sV/A7+cw+n0nMtms8lms/lg/vY7cuSIqqurlZaW5lHHiBEjtGPHDk2dOlV79+7V+fPnPcb06tVLgwYN0o4dO8wb/o1LGzEzZ/q3EFOY6u8CTMGe7O8KAGs6deqU7Ha7z+ft3LmzoqKiVF0d45P5vvGNbygmxnOuxx9/XHl5eT6Zv72qq6slSZGRkR77IyMj9fHHH7vHdO7cWT169GgypvH17RVQ4d+rVy9VVlaqe/fuCjL6HQ4fcTqdiomJUWVlpcLDw/1dTsDiPLUP56l9OE/tE4jnyeVy6dSpU+rVq1eHzB8WFqYjR46ovr7eJ/O5XK4medNS15+Xl+dezm/J7t27lZx86V3HxbU0V9/F2jPmYgEV/p06dVLv3r39XUazwsPDA+ZfrkDGeWofzlP7cJ7aJ9DOU0d0/F8VFhamsLCwDn2P5jzyyCO69957Wx1z1VVXXdLcUVFRkr7s7qOjo937jx075l4NiIqKUn19vU6cOOHR/R87dkypqalevV9AhT8AAIHK4XDI4XB0yNx9+/ZVVFSUioqKlJCQIOnLbwxs3bpVTz75pCQpKSlJoaGhKioq0oQJEyRJVVVVev/99/Xb3/7Wq/cj/AEA8LGKigodP35cFRUVunDhgsrKyiRJ11xzjb7xjW9Ikvr376+5c+cqIyNDQUFBmjlzpp544glde+21uvbaa/XEE0+oa9eumjhxoqQvV1QmT56sn/zkJ7ryyisVERGhWbNmKTY21n31f3sR/m2w2Wx6/PHHL/uVn2bDeWofzlP7cJ7ah/MUuH75y19q+fLl7seN3fzmzZs1cuRISdLBgwdVW1vrHvPf//3fOnv2rKZPn64TJ05o6NCheuutt9S9e3f3mHnz5ikkJEQTJkzQ2bNndeuttyo/P1/BwcFe1Rfk6vjvYwAAgADCL/wBAGAxhD8AABZD+AMAYDGEPwAAFkP4t8NvfvMbpaSk+LuMgMd5AgBzIPzboby8XIMHD/Z3GQGP89S87du3a8iQIQoLC5PD4dC8efP8XVJA4jy1jXMEXyH826G8vFzx8fH+LiPgcZ6aKiwsVEZGhqZPn659+/Zp6tSpysnJ0UcffeTv0gIK56ltnCP4lAsePvjgA9eIESNcYWFhrvj4eNfu3btdwcHBrh07dvi7tIDCeWrb2bNnXb1793atWrXKve+LL75wde/e3bVs2TLXXXfd5briiitcd999tx+r9L/WzlNeXp5rxIgRrgEDBrhiY2NdL7/8sh8r9Z/WztH8+fNdycnJrsGDB7sGDRrkWrp0qR8rhVnQ+X/Fhx9+qKFDhyo5OVnvv/++fvnLX2rs2LFyuVyKi4vzd3kBg/PUPps2bdLZs2eVmZnp3hccHKyQkBDZbDbNmDFDK1as8GOFgaG183TNNddo/vz5+uCDD7Rx40b9+Mc/1pkzZ/xYrX+0do4cDoe2bt2qsrIy7dq1S3PnzlVNTY0fq4UZEP5f8fDDD2vs2LF6+umndfXVVysjI0PDhw/XNddco27duikjI0M9evTQ+PHj/V2qX7V2no4fP66RI0dq4MCBiouL0yuvvOLvcv1m8+bNGjx4sMfPbh45ckQnTpxQYmKibr75Zo+f7bSqts5T40dJPXv2VEREhI4fP+6nSv2nrXPUtWtXSdK5c+d04cIFufjhVrSB3/b/t48//libNm1SSUmJx/7Q0FD3RWwzZszQgw8+6PF7zVbT1nkKCQnR/PnzFR8fr2PHjikxMVHp6enq1q2bnyr2n9LS0ib3HF+0aJESExM1YMAAP1UVeNp7nvbs2aOGhgbFxMRc7hL9rq1zdPLkSY0YMUKHDx/WU0891WF3nsPXB+H/b2VlZQoJCVFsbKzH/pKSEvcdlW6++WZt2bLFD9UFjrbOU3R0tPte1F/t1KwY/mVlZWpoaNDy5cuVkpKil19+WUuWLNG7777r79ICSnvOU01Nje6//349//zzfqzUf9o6R1dccYXKy8t19OhRjRs3TuPHj3ffAx5oDuH/b506dVJDQ4Pq6+sVEvLlaSksLNT+/fu5gv0rvDlPVu7UKioqVFNTo9dff12PPvqoDh06pLi4OG3YsEGJiYn+Li9gtOc81dXVKSMjQ7m5uUpNTfVzxZefN/8sRUZGKi4uTtu2bdM999zjp4phBoT/vyUlJSk0NFSzZs3SrFmz9P777+sHP/iBJPHd9a9o73miUytTRESE0tPTlZ6e7u9yAlZb58nlcmnSpEm65ZZblJWV5YcK/a+tc3T06FF16dJF4eHhcjqd2rZtm/vfSaAlXPD3b7169dLzzz+v1157TcnJyVqwYIEeeOABORwOffOb3/R3eQGjPefJ6p2a9OVntBd/NHKxUaNG6Z577lFhYaF69+6t3bt3X6bqAkdb5+ndd99VQUGB1q5dq/j4eMXHx+u99967jBX6X1vn6JNPPtHw4cM1ePBgfec739EjjzzCt27QpiAXl4V6ZcuWLVq8eLH+8pe/+LuUgORyuTRx4kT169dPeXl5/i7Hb+666y5961vf0sKFC/1dSkDjPLWNc4SOQPh7YdSoUSopKdGZM2cUERGhNWvW6IYbbvB3WQFl+/btGj58uEfn8cILL7TZBQMALh/CHwAAi+EzfwAALIbwBwDAYgh/AAAshvAHAMBiCH8AACyG8AcAwGIIfwAALIbwBwDAYgh/AAAshvAHAMBiCH8AACyG8AcAwGL+HzJTa4CTzqiSAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "from pygimli.frameworks.resolution import modelCovariance\n", "var, MCM = modelCovariance(ves.inv)\n", @@ -226,7 +260,12 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, "source": [ "The model covariance matrix delivers variances and a scaled (dimensionless)\n", "correlation matrix. The latter show the interdependency of the parameters\n", @@ -238,14 +277,35 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 58, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG5CAYAAACHhJ4rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAq+ElEQVR4nO3de3hU9aHu8XdIJjcggRCIhHulWCIJKKQK3ggKIVoQQba7ApIinE2JKKRIvdCitFvqU6VWHaliQdgbPIgWL8BuYBckqFsJAmJNDxtsIAhB5CIRAskkWeePmJExF2bITFbym+/neXicWWvNWm/Cz5mX31oz47AsyxIAAICBWtkdAAAAIFgoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY4XbHcBuVVVVOnLkiNq2bSuHw2F3HAAA4APLsvTNN98oKSlJrVrVP28T8kXnyJEj6tatm90xAADAJTh06JC6du1a7/qQLzpt27aVVP2Lio2NtTkN7OJ2u7Vx40aNGDFCTqfT7jiATxi3Uqm7VMPWDJMkbR6/WTHOGJsToSGBHLMlJSXq1q2b53W8PiFfdGpOV8XGxlJ0Qpjb7VZMTIxiY2ND9gUDLQ/jVgp3hyssOkxS9fM4Rad5C8aYvdhlJ1yMDAAAjBWyRcflcik5OVlpaWl2RwEAAEESskUnOztbBQUFys/PtzsKAAAIkpAtOgAAwHwUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAICvDu+UXvlJ9X8BtAgUHQDw1Sf/VzqwTdqz2u4kAHwU8t91BQAN+rpIKj0hy31e+vvrckiyPn1djv4/lWRJMR2kdt3tTgmgHhQdAGjIMymSpAu/NtBRelx66abvFjx2umkzAfBZyJ664ruuAPhk7BKpVT3/JmwVXr0eQLMVskWH77oC4JPUf5Gm/q3udVP/Vr0eQLMVskUHAPxVZVWfwLK8TmQBaM4oOgBwMa07qqp1J31q9dIj7ntV1XmA1KaT1Lqj3ckAXAQXIwPAxcR10fns3bp9wRZJDs3LelIxYVVSeKTdyQBcBDM6AOCL8Eh53nvlcFBygBaCogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFghW3T4CggAAMwXskWHr4AAAMB8IVt0AACA+Sg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGCtmiwwcGAgBgvpAtOnxgIAAA5gvZogMAAMxnRNFZt26drrjiCv3whz/Uyy+/bHccAADQTITbHaCxKioqlJOToy1btig2NlZXX321xo4dq/j4eLujAQAAm7X4GZ3t27fryiuvVJcuXdS2bVvdeuutys3NtTsWAABoBmwvOnl5eRo1apSSkpLkcDj05ptv1trmhRdeUK9evRQVFaWBAwdq27ZtnnVHjhxRly5dPPe7du2qw4cPN0V0AADQzNl+6urs2bPq37+/fvazn2ncuHG11q9evVqzZs3SCy+8oOuuu04vvviiMjMzVVBQoO7du8uyrFqPcTgc9R6vrKxMZWVlnvslJSWSJLfbLbfbHYCfCC1Rzd89YwD1cbsrLrjtlttR+7mnqTFuqy9fuPC2W6H7u2gJAjlmfd2H7UUnMzNTmZmZ9a5ftGiR7r33Xk2dOlWS9Mwzzyg3N1eLFy/WwoUL1aVLF68ZnC+++ELXXHNNvftbuHChHn/88VrLN27cqJiYmEb8JDDBpk2b7I6AZqqsUqp5yszN3ajIMFvjeAnlcVtulXtu5+bmKsIRYWMa+CoQY7a0tNSn7RxWXVMiNnE4HFq7dq3GjBkjSSovL1dMTIzWrFmjO+64w7PdAw88oN27d2vr1q2qqKhQ37599e6773ouRv7www/VoUOHOo9R14xOt27ddPz4ccXGxgb150Pz5Xa7tWnTJg0fPlxOp9PuOGiGSssr1P83myVJn/xqmGIibP93IuNW0rmKc7ruteskSe//y/uKDo+2OREaEsgxW1JSooSEBJ0+fbrB12/7/09twPHjx1VZWanExESv5YmJiTp69KgkKTw8XE8//bTS09NVVVWluXPn1ltyJCkyMlKRkZG1ljudzpB9osB3GAeoj9P67pR49ThpPk+foTxuLzxVFR4eHrK/h5YmEGPW18c3n/9TG/D9a24sy/JaNnr0aI0ePbqpYwEAgGbO9nddNSQhIUFhYWGe2Zsax44dqzXLAwAA8H3NuuhERERo4MCBtS5a2rRpk4YMGWJTKgAA0FLYfurqzJkz2r9/v+d+YWGhdu/erfj4eHXv3l05OTmaNGmSBg0apMGDB+ull15SUVGRpk+fbmNqAADQEthedHbs2KH09HTP/ZycHEnS5MmT9corr+iuu+7SiRMntGDBAhUXF6tfv37asGGDevToYVdkAADQQthedIYOHVrnh/5daMaMGZoxY0YTJQIAAKZo1tfoAAAANAZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAACNcs+oafXz0Y7tjoJmh6AAAjLGhcIPdEdDMhNsdAABamtLySrsjSJLc7gqVVUql5RVyWg674zS54rNH9GXpl17LNh7cpNt6jZElS+0i26lz6ySb0qFGtDNMDod945OiAwB+GvTb/7Y7wgXCNXf7ZrtD2KJt34dqLTt1/pQm597tuf/NP37XlJFQh0E92mvN9MG2lR1OXQGAD6KdYRrUo73dMXCBc4fvkmV5v4zVvJZaViudO3yXDanwfTsOntI5t32zoMzoAIAPHA6H1kwfbOsT9ve53W7l5m5URsYIOZ1Ou+PYIEOffDVS/+e/f1ZrzYqR/6kfxfe1IRNqlJZXNovZT4oOAPjI4XAoJqL5PG26HZYiw6SYiHA5nc0nV1OKDP/ejI4csmQpyhnWrP6uYB9OXQEAWqz2Ud+dTvxl2i+V3CFZHaI6KD4q3sZUaE6ouwCAFisxJtFze+wPx2pC3wlyV7kVERZhYyo0J8zoAACM4XA4KDnwQtEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAwVsgWHZfLpeTkZKWlpdkdBQAABEnIFp3s7GwVFBQoPz/f7igAACBIQrboAAAA81F0AACAsSg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAwFkUHAAAYK2SLjsvlUnJystLS0uyOAgAAgiRki052drYKCgqUn59vdxQAABAkIVt0AACA+Sg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBY4f5svHfvXr366qvatm2bDhw4oNLSUnXs2FFXXXWVMjIyNG7cOEVGRgYrKwAAgF98mtHZtWuXhg8frv79+ysvL09paWmaNWuWfvOb32jixImyLEuPPvqokpKS9OSTT6qsrCzYuQEAAC7KpxmdMWPG6MEHH9Tq1asVHx9f73b/8z//oz/84Q96+umn9cgjjwQsJAAAwKXwqejs27dPERERF91u8ODBGjx4sMrLyxsdDAAAoLF8OnXlS8lpzPYAAADB4NfFyDW2b9+ud999V8eOHVNVVZXXukWLFgUkGAAAQGP5XXSeeOIJzZs3T1dccYUSExPlcDg86y68DQAAYDe/i84f//hHLV26VFlZWUGIAwAAEDh+f2Bgq1atdN111wUjCwAAQED5XXRmz54tl8sVjCwAAAAB5fepqzlz5ui2227T5ZdfruTkZDmdTq/1f/nLXwIWLphcLpdcLpcqKyvtjgIAAILE7xmdmTNnasuWLerTp486dOiguLg4rz8tRXZ2tgoKCpSfn293FAAAECR+z+isWLFCb7zxhm677bZg5AEAAAgYv2d04uPjdfnllwcjCwAAQED5XXQee+wxzZ8/X6WlpcHIAwAAEDB+n7p69tln9fnnnysxMVE9e/asdTHyzp07AxYOAACgMfwuOmPGjAlCDAAAgMDzu+jMnz8/GDkAAAACzu9rdHxhWVYwdgsAAOAXn4pO3759tWrVKpWXlze43b59+/Tzn/9cTz75ZEDCAQAANIZPp65cLpd++ctfKjs7WyNGjNCgQYOUlJSkqKgonTp1SgUFBXrvvfdUUFCg++67TzNmzAh2bgAAgIvyqegMGzZM+fn5+uCDD7R69WqtWrVKBw4c0Llz55SQkKCrrrpK99xzjyZOnKh27doFOTIAAIBv/LoYeciQIRoyZEiwsgAAAARUUC5GBgAAaA4oOgAAwFgUHQAAYCyKDgAAMBZFBwAAGMvvr4CQpKqqKu3fv1/Hjh1TVVWV17obb7wxIMGCzeVyyeVyqbKy0u4oAAAgSPwuOh9++KHuvvtuHTx4sNZXPTgcjhZTHLKzs5Wdna2SkhLFxcXZHQcAAASB30Vn+vTpGjRokNavX6/OnTvL4XAEIxcAAECj+V109u3bp9dff129e/cORh4AAICA8fti5GuuuUb79+8PRhYAAICA8ntGZ+bMmfrFL36ho0ePKiUlRU6n02t9ampqwMIBAAA0ht9FZ9y4cZKkKVOmeJY5HA5ZltWiLkYGAADm87voFBYWBiMHAABAwPlddHr06BGMHAAAAAHnU9F5++23lZmZKafTqbfffrvBbUePHh2QYAAAAI3lU9EZM2aMjh49qk6dOmnMmDH1bsc1OgAAoDnxqehc+DUP3//KBwAAgOaKL/UEAADGougAAABjhWzRcblcSk5OVlpamt1RAABAkIRs0cnOzlZBQYHy8/PtjgIAAIIkZIsOAAAwn98fGChVv/Nq//79OnbsWK13Yd14440BCQYAANBYfhedDz/8UHfffbcOHjwoy7K81vE5OgAAoDnxu+hMnz5dgwYN0vr169W5c2c5HI5g5AIAAGg0v4vOvn379Prrr6t3797ByAMAABAwfl+MfM0112j//v3ByAIAABBQPs3o7Nmzx3N75syZ+sUvfqGjR48qJSVFTqfTa9vU1NTAJgQAALhEPhWdAQMGyOFweF18PGXKFM/tmnVcjAwAAJoTn4pOYWFhsHMAAAAEnE9Fp0ePHp7beXl5GjJkiMLDvR9aUVGhDz74wGtbAAAAO/l9MXJ6erpOnjxZa/np06eVnp4ekFBNge+6AgDAfH4XnZprcb7vxIkTat26dUBCNQW+6woAAPP5/Dk6Y8eOlVR94XFWVpYiIyM96yorK7Vnzx4NGTIk8AkBAAAukc9FJy4uTlL1jE7btm0VHR3tWRcREaFrr71W06ZNC3xCAACAS+Rz0Vm2bJkkqWfPnpozZ06LOk0FAABCk99fATF//nxJ0rFjx7R37145HA716dNHnTp1Cng4AACAxvD7YuSSkhJNmjRJXbp00U033aQbb7xRXbp00cSJE3X69OlgZAQAALgkfhedqVOn6qOPPtK6dev09ddf6/Tp01q3bp127NjBNToAAKBZ8fvU1fr165Wbm6vrr7/esywjI0NLlizRyJEjAxoOAACgMfye0enQoYPnHVgXiouLU/v27QMSCgAAIBD8Ljrz5s1TTk6OiouLPcuOHj2qBx98UL/61a8CGg4AAKAx/D51tXjxYu3fv189evRQ9+7dJUlFRUWKjIzUV199pRdffNGz7c6dOwOXNMBcLpdcLhfftg4AgMH8LjpjxowJQoyml52drezsbJWUlNR5Kg4AALR8l/w5OgAAAM2d39foSNLXX3+tl19+WQ8//LDnm8x37typw4cPBzQcAABAY/g9o7Nnzx7dcsstiouL04EDBzRt2jTFx8dr7dq1OnjwoFasWBGMnAAAAH7ze0YnJydHWVlZ2rdvn6KiojzLMzMzlZeXF9BwAAAAjeF30cnPz9e//du/1VrepUsXHT16NCChAAAAAsHvohMVFaWSkpJay/fu3auOHTsGJBQAAEAg+F10br/9di1YsEBut1uS5HA4VFRUpIceekjjxo0LeEAAAIBL5XfReeqpp/TVV1+pU6dOOnfunG666Sb17t1bbdu21b//+78HI2NQuFwuJScnKy0tze4oAAAgSPx+11VsbKzee+89bd68WTt37lRVVZWuvvpq3XLLLcHIFzR8YCAAAObzu+jUGDZsmIYNGxbILAAAAAHl16mrqqoqLV26VD/5yU/Ur18/paSkaPTo0VqxYoUsywpWxou644471L59e9155522ZQAAAM2Pz0XHsiyNHj1aU6dO1eHDh5WSkqIrr7xSBw8eVFZWlu64445g5mzQ/fffzwcVAgCAWnw+dfXKK68oLy9Pf/vb35Senu61bvPmzRozZoxWrFihe+65J+AhLyY9PV3vvvtukx8XAAA0bz7P6Lz66qt65JFHapUcqfp6nYceekgrV670O0BeXp5GjRqlpKQkORwOvfnmm7W2eeGFF9SrVy9FRUVp4MCB2rZtm9/HAQAAocfnorNnzx6NHDmy3vWZmZn65JNP/A5w9uxZ9e/fX88//3yd61evXq1Zs2bp0Ucf1a5du3TDDTcoMzNTRUVFfh8LAACEFp9PXZ08eVKJiYn1rk9MTNSpU6f8DpCZmanMzMx61y9atEj33nuvpk6dKkl65plnlJubq8WLF2vhwoV+H6+srExlZWWe+zWf8ux2uz0fgojQU/N3zxhAS8K4lSoqKrxuuxW6v4vmxu2uuOC2W26HFdAx6+s+fC46lZWVCg+vf/OwsDCvARcI5eXl+vjjj/XQQw95LR8xYoQ++OCDS9rnwoUL9fjjj9davnHjRsXExFzSPmGOTZs22R0B8Fsoj9tyq9xzOzc3VxGOCBvT4EJllVJNzcjN3ajIsO/WBWLMlpaW+rSdz0XHsixlZWUpMjKyzvUXzpIEyvHjx1VZWVlrJikxMdHrC0QzMjK0c+dOnT17Vl27dtXatWvr/cTjhx9+WDk5OZ77JSUl6tatm0aMGKHY2NiA/wxoGdxutzZt2qThw4fL6XTaHQfwCeNWOldxTgteWyCp+rUgOjza5kSoUVpeobnbN0uSMjJGKCYiPKBjtq7v3ayLz0Vn8uTJF90mWO+4cjgcXvcty/Jalpub6/O+IiMj6yxrTqczZJ8o8B3GAVqiUB63F56qCg8PD9nfQ3PktL57na4eo+Hfu9+4vytfH+9z0Vm2bNklh7lUCQkJCgsL85q9kaRjx441eL0QAACAdAlf6tmUIiIiNHDgwFrn8jZt2qQhQ4bYlAoAALQUl/xdV4Fy5swZ7d+/33O/sLBQu3fvVnx8vLp3766cnBxNmjRJgwYN0uDBg/XSSy+pqKhI06dPtzE1AABoCWwvOjt27PD6EMKaC4UnT56sV155RXfddZdOnDihBQsWqLi4WP369dOGDRvUo0cPuyIDAIAWwvaiM3To0It+IeiMGTM0Y8aMJkoEAABM0ayv0QEAAGgMig4AADAWRQcAABiLogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QEAAMai6AAAAGNRdAAAQFCkOP6pVc7fqlXxLtsyUHQAAEBQjA3bpiFhBQr79DXbMoTbdmQAAGCW8rPS6S+k0pNyVFRqdNgHkqTwgjek1LGSu0Jtzh1p0kgUHQAAEBhPJHluRkuKdlTfdpSekJaOlFPSzZLcmtpkkTh1BQAAjMWMDgAACIxHvj0tdXSPtHRkrdXue9Ypd3exMpowEkUHAAAERkTr6v+GR0uSqiyHWjksWXLIIUsKj1ZlWGSTRuLUFQAACKzWHVXVupM+tXrpEfe9quo8QGrTSWrdscmjMKMDAAACK66Lzmfv1u0LtkhyaF7Wk4oJq5KsVpJ2N2kUZnQAAEDghUdKqnnblePb+02PogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFgUHQAAYKyQLToul0vJyclKS0uzOwoAAAiSkC062dnZKigoUH5+vt1RAABAkIRs0QEAAOaj6AAAAGNRdAAAgLEoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QEAAMai6AAAAGNRdAAAgLEoOgAAwFgUHQAAYCyKDgAAMFbIFh2Xy6Xk5GSlpaXZHQUAAARJyBad7OxsFRQUKD8/3+4oAAAgSEK26AAAAPNRdAAARvjHiX/YHQHNEEUHAGCEDYUb7I6AZijc7gAA0NxZlqVz7kq7Y9TidleorFIqLa+Q03LYHadJFZ89oq/LvlZ5Zbln2caDm3RbrzGyZKldZDt1bp1kY0KUljeP/2coOgDQAMuydOef/kcfHzxld5R6hGvu9s12h2hybfs+VGvZqfOnNDn3bs/9b/7xu6aMhGaKU1cA0IBz7spmXHJC17nDd8myvF/CHN9OallWK507fJcNqVCXQT3aK9oZZtvxmdEBAB/tmHeLYiLse8L+PrfbrdzcjcrIGCGn02l3nCaWof938navGZwaK0b+p34U39eGTKhLtDNMDod9p1YpOgDgo5iIMMVENJ+nTbfDUmSYFBMRLqez+eRqKlHfmyVwyCFLlqKczevvCfbi1BUAoEWKj4pXh6gOnvs/iv+ROkR1UHxUvI2p0NxQdADgIlIc/9Qq52/VqniX3VFwgctaX6Y3b3/Tc39ZxjJtvHOjLmt9mX2h0OxQdADgIsaGbdOQsAKFffqa3VHwPRFhEZ7bDofD6z4gcY3Od8rPSuXN5yJDNDG3W2GVZdXjwAq1izpRp9NfSKUn5aio1OiwDyRJ4QVvSKljJUtSTLwU19XejIxbqeLcd7fdpZJlSRGt7cuDZoeiU+PpK6TI0PrALXzHKeknkrTH5iBodqIlRX/71OAoPSEtHWlrngsxblX9nvKe3apv/753ddF57LS9mdCscOoKAAAYixmdGr/YK8XG2p0CNgntzyNBvY7uqXsGZ8pfpctSmz7P9zBuVX3qak169e0H90vh0fbmQbND0akR0ZrzuqHM4VZlWGT1GAjVFwzU9u2LZpXlUCuHJevbT2pReHTzeL5g3H73cciS5Iyp/gNcgFNXAFCf1h1V1bqTPrV66RH3varqPEBq00lq3dHuZAB8xIwOANQnrovOZ+/W7Qu2SHJoXtaTigmrksIj7U4GwEfM6ABAQ8IjJdW87cpByQFaGIoOAAAwFkUHAAAYK2SLjsvlUnJystLS0uyOAgAAgiRki052drYKCgqUn59vdxQAABAkIVt0AACA+Sg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QEAAMYK2aLjcrmUnJystLQ0u6MAAIAgCbc7gF2ys7OVnZ2tkpISxcXF2R0HQDMV7QxTwYIMz20ALUvIFh0A8IXD4VBMBE+VQEsVsqeuAACA+Sg6AADAWBQdAABgLIoOAAAwFkUHAAAYi6IDAACMRdEBAADGougAAABjUXQAAICxKDoAAMBYFB0AAGAsig4AADBWyBYdl8ul5ORkpaWl2R0FAAAESch+JW92drays7N1+vRptWvXTiUlJXZHgo3cbrdKS0tVUlIip9NpdxzAJ4xbqdRdqspzlZKkkpISVTgrbE6EhgRyzNa8bluW1eB2DutiWxjuiy++ULdu3eyOAQAALsGhQ4fUtWvXeteHfNGpqqrSkSNH1LZtWzkcjnq3S0tLU35+fhMmC+6xA7HPS92HP48L5LYNrS8pKVG3bt106NAhxcbG+nS8lsCucRus49o1bv19jK/b+7JdqI1bnmsDt49gPNc2pzFrWZa++eYbJSUlqVWr+q/ECdlTVzVatWrVYBOsERYWZtsTSTCOHYh9Xuo+/HlcILf1ZV+xsbHGvGBI9o3bYB3XrnHr72N83d6X7UJt3PJcG7h9BOO5trmN2bi4uItuE7IXI/srOzvbqGMHYp+Xug9/HhfIbe38O7SLXT9zsI5r17j19zG+bu/LdqE2bnmuDdw+gvFc2xLHbMifugKk6unUuLg4nT592ph/GcN8jFu0NHaMWWZ0AEmRkZGaP3++IiMj7Y4C+Ixxi5bGjjHLjA4AADAWMzoAAMBYFB0AAGAsig4AADAWRQcAABiLogMAAIxF0QF8cMcdd6h9+/a688477Y4CXNShQ4c0dOhQJScnKzU1VWvWrLE7EnBR33zzjdLS0jRgwAClpKRoyZIlAdkvby8HfLBlyxadOXNGy5cv1+uvv253HKBBxcXF+vLLLzVgwAAdO3ZMV199tfbu3avWrVvbHQ2oV2VlpcrKyhQTE6PS0lL169dP+fn56tChQ6P2y4wO4IP09HS1bdvW7hiATzp37qwBAwZIkjp16qT4+HidPHnS3lDARYSFhSkmJkaSdP78eVVWVioQczEUHRgvLy9Po0aNUlJSkhwOh958881a27zwwgvq1auXoqKiNHDgQG3btq3pgwLfCuSY3bFjh6qqqtStW7cgp0aoC8S4/frrr9W/f3917dpVc+fOVUJCQqNzUXRgvLNnz6p///56/vnn61y/evVqzZo1S48++qh27dqlG264QZmZmSoqKmripEC1QI3ZEydO6J577tFLL73UFLER4gIxbtu1a6dPPvlEhYWFWrVqlb788svGB7OAECLJWrt2rdeyH//4x9b06dO9lv3oRz+yHnroIa9lW7ZsscaNGxfsiICXSx2z58+ft2644QZrxYoVTRET8NKY59oa06dPt1577bVGZ2FGByGtvLxcH3/8sUaMGOG1fMSIEfrggw9sSgXUz5cxa1mWsrKyNGzYME2aNMmOmIAXX8btl19+qZKSEknV33Kel5enK664otHHDm/0HoAW7Pjx46qsrFRiYqLX8sTERB09etRzPyMjQzt37tTZs2fVtWtXrV27VmlpaU0dF/BpzL7//vtavXq1UlNTPddJ/Md//IdSUlKaOi4gybdx+8UXX+jee++VZVmyLEv33XefUlNTG31sig4gyeFweN23LMtrWW5ublNHAhrU0Ji9/vrrVVVVZUcsoEENjduBAwdq9+7dAT8mp64Q0hISEhQWFuY1eyNJx44dq/UvD6A5YMyiJbJz3FJ0ENIiIiI0cOBAbdq0yWv5pk2bNGTIEJtSAfVjzKIlsnPccuoKxjtz5oz279/vuV9YWKjdu3crPj5e3bt3V05OjiZNmqRBgwZp8ODBeumll1RUVKTp06fbmBqhjDGLlqjZjttGv28LaOa2bNliSar1Z/LkyZ5tXC6X1aNHDysiIsK6+uqrra1bt9oXGCGPMYuWqLmOW77rCgAAGItrdAAAgLEoOgAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAE1m6NChmjVrVpM+vrHHrM+JEyfUqVMnHThwIOD7rs+dd96pRYsWNdnxABPwXVcAJElZWVlavny5JCksLExJSUm67bbb9MQTT6h9+/YBOcbJkyfldDrVtm3bi247dOhQDRgwQM8888wlPb6ux9S1z0s1Z84cnTp1Sn/+85+9lu/Zs0ePPfaYtm7dqrNnz+ryyy9XVlaWZs+erfDw8EYdc8+ePUpPT1dhYaFiY2MbtS8gVDCjA8Bj5MiRKi4u1oEDB/Tyyy/rnXfe0YwZMwK2//j4eL9KSiAe39hj1uXcuXP685//rKlTp3ot37p1q6699lpFR0frrbfe0ieffKK5c+fqqaee0tixY1VVVdWo46ampqpnz55auXJlo/YDhJSgfz86gBZh8uTJ1u233+61LCcnx4qPj/fcr6qqsp588kmrV69eVlRUlJWammqtWbPG6zFr1qyx+vXrZ0VFRVnx8fHWzTffbJ05c8ayLMu66aabrAceeOCi206ePNmS5PWnsLDQ6/F/+tOfrKSkJKuystLr+KNGjbLuuecez/2ax9S3z+XLl1vx8fHW+fPnvfYzduxYa9KkSXX+rt544w0rISHBa1lFRYXVq1cva8KECbW2/+yzzyyn02m9/PLLXrnuu+8+64EHHrDatWtnderUyXrxxRetM2fOWFlZWVabNm2sH/zgB9aGDRu89vXYY49ZN9xwQ525ANTGjA6AOv3zn//UX//6VzmdTs+yefPmadmyZVq8eLE+++wzzZ49WxMnTtTWrVslScXFxfrpT3+qKVOm6B//+IfeffddjR07VlYdZ8gb2vaPf/yjBg8erGnTpqm4uFjFxcXq1q2b1+PHjx+v48ePa8uWLZ5lp06dUm5uriZMmFDrePXtc/z48aqsrNTbb7/t2fb48eNat26dfvazn9X5u8nLy9OgQYO8lm3fvl2FhYV68MEHa22fnJysW2+9VatXr/Zavnz5ciUkJGj79u2aOXOmfv7zn2v8+PEaMmSIdu7cqYyMDE2aNEmlpaWex/z4xz/W9u3bVVZWVmc2AN4ad8IYgFHWrVunNm3aqLKyUufPn5ckz8WvZ8+e1aJFi7R582YNHjxYkvSDH/xA7733nl588UXddNNNKi4uVkVFhcaOHasePXpIklJSUuo81sW2jYiIUExMjC677LI6Hx8fH6+RI0dq1apVuvnmmyVJa9asUXx8vOf+heLi4urcZ3R0tO6++24tW7ZM48ePlyStXLlSXbt21dChQ+s89oEDB5SUlOS1rLCwUJL0wx/+sM7H9OnTR2+99ZbXsv79+2vevHmSpIcffli/+93vlJCQoGnTpkmSfv3rX2vx4sXas2ePrr32WklSly5dVFZWpqNHj3p+bwDqx4wOAI/09HTt3r1bH330kWbOnKmMjAzNnDlTklRQUKDz589r+PDhatOmjefPihUr9Pnnn0uqfuG++eablZKSovHjx2vJkiU6depUncfyZ9v6TJgwQW+88YZndmPlypX613/9V4WFhfm1n2nTpmnjxo06fPiwJGnZsmXKysqSw+Goc/tz584pKirKa1nNxcEnT56s8zGnTp2qdQFxamqq53ZYWJg6dOjgVfYSExMlSceOHfMsi46OliSvWR4A9aPoAPBo3bq1evfurdTUVD377LMqKyvT448/LkmeC2nXr1+v3bt3e/4UFBTo9ddfl1T9Yr1p0yb913/9l5KTk/Xcc8/piiuu8Mx2XMifbeszatQoVVVVaf369Tp06JC2bdumiRMn+v1zX3XVVerfv79WrFihnTt36tNPP1VWVla92yckJNQqZYMHD5bT6dQ777xTa/vKykpt3LhR119/vdfyC08LSpLD4fBaVlO0LryIuaZIdezY0bcfDghxFB0A9Zo/f76eeuopHTlyRMnJyYqMjFRRUZF69+7t9efC62ccDoeuu+46Pf7449q1a5ciIiK0du3aOvff0LYRERGqrKxsMF90dLTGjh2rlStX6tVXX1WfPn00cODAerdvaJ9Tp07VsmXLtHTpUt1yyy21rgm60FVXXaWCggKvZR06dND999+v3/72tzpy5IjXuj/84Q86ceKEZs+e3eDP44u///3v6tq1qxISEhq9LyAUUHQA1Gvo0KG68sor9cQTT6ht27aaM2eOZs+ereXLl+vzzz/Xrl275HK5PJ+/89FHH+mJJ57Qjh07VFRUpL/85S/66quv1Ldv31r7vti2PXv21EcffaQDBw7o+PHj9b41e8KECVq/fr2WLl160dmchvY5YcIEHT58WEuWLNGUKVMa3E9GRoY+++wzr1mdM2fO6P7771evXr2Unp6unTt3SpJ+//vf65FHHtFzzz3nU3m7mG3btmnEiBGN2gcQSig6ABqUk5OjJUuW6NChQ/rNb36jX//611q4cKH69u2rjIwMvfPOO+rVq5ek6utU8vLydOutt6pPnz6aN2+enn76aWVmZtba78W2nTNnjsLCwpScnKyOHTuqqKioznzDhg1TfHy89u7dq7vvvrvBn6WhfcbGxmrcuHFq06aNxowZ0+B+UlJSNGjQIL322mueZU899ZR69Oih999/X//7v/+rZ599VpI0d+5cud1uTZkyRZ07d9ahQ4ca3HdDzp8/r7Vr13ouVgZwcXwyMgB8a/jw4erbt6+npDRkw4YNmjNnjv7+97+rVaum+Tejy+XSW2+9pY0bNzbJ8QAT8PZyACHv5MmT2rhxozZv3qznn3/ep8fceuut2rdvnw4fPtzg9TyB5HQ69dxzzzXJsQBTMKMDIOT17NlTp06d0q9+9SvNmTPH7jgAAoiiAwAAjMXFyAAAwFgUHQAAYCyKDgAAMBZFBwAAGIuiAwAAjEXRAQAAxqLoAAAAY1F0AACAsSg6AADAWBQdAABgLIoOAAAw1v8H8uy2yAuu+coAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "thk = ves.model[:nlay-1]\n", "res = ves.model[nlay-1:]\n", @@ -255,25 +315,43 @@ "fig, ax = plt.subplots()\n", "#ves.showModel(synthModel, ax=ax, label=\"synth\", plot=\"semilogy\", zmax=20)\n", "ves.showModel(ves.model, ax=ax, label=\"model\", zmax=50)\n", + "ax.invert_yaxis()\n", "ax.errorbar(res, mid, marker=\"*\", ls=\"None\", xerr=res*var[nlay-1:])\n", "ax.errorbar(resmean, z, marker=\"*\", ls=\"None\", yerr=thk*var[:nlay-1])" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 63, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Layer Resistivity (Ωm) Res. Uncertainty (Ωm) Thickness (m) Thk. Uncertainty (m)\n", + " 1 62 2 1.4 0.06\n", + " 2 945 60 31.49 39.82\n", + " 3 61 188 ∞ -\n" + ] + } + ], "source": [ - "print(res)\n", - "print(res*var[nlay-1:])\n", - "print(thk)\n", - "print(thk*var[:nlay-1])" + "# Create uncertainty table\n", + "uncertainty_table = pd.DataFrame({\n", + " 'Layer': np.arange(1, nlay+1),\n", + " 'Resistivity (Ωm)': np.round(res, 0).astype(int),\n", + " 'Res. Uncertainty (Ωm)': np.round(res*var[nlay-1:], 0).astype(int),\n", + " 'Thickness (m)': list(np.round(thk, 2)) + [np.inf],\n", + " 'Thk. Uncertainty (m)': list(np.round(thk*var[:nlay-1], 2)) + [np.nan]\n", + "})\n", + "uncertainty_table = uncertainty_table.replace({np.inf: '∞', np.nan: '-'})\n", + "print(uncertainty_table.to_string(index=False))" ] } ], diff --git a/environment.yml b/environment.yml index 270027d..5bcf7e5 100644 --- a/environment.yml +++ b/environment.yml @@ -8,3 +8,4 @@ dependencies: - pygimli>=1.5.0 - jupyterlab - suitesparse=5.10.1 + - pandas