{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exploration\n", "The first step in the `CORBASS` algorithm is an exploration of the parameter space. This way the region of interest, i.e. the region where the probability mass is concentrated, is found. Here we give an example how to run the `corbass.exploration` module and how to interpret the output. First we import the function `explore` from the module:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import sys\n", "import os\n", "# relative import\n", "sys.path.append(os.path.abspath('') + '/../')\n", "from corbass.exploration import explore" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After the import we load the example parameter file:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from corbass.utils import load\n", "pars = load('./Example_Parfile.py')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we use the `explore` function to explore the parameter space. This is done via grid evaluation. The model parameter bounds and the number of gridpoints are specified in the parameter file. Then for each point in the grid, i.e. each choice of model parameters, the posterior probability is calculated. `explore` returns an instance of `ExplorationResult`, containing the posterior on the grid and a scale. We use the `write` member function to write the results to the disk, to use them in [the next example](https://gitext.gfz-potsdam.de/arthus/corbass/blob/master/examples/Example_1_Exploration.ipynb)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exploring 1450, this may take a while...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100% (15625 of 15625) |##################| Elapsed Time: 0:43:00 Time: 0:43:00\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Exploring 1550, this may take a while...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100% (15625 of 15625) |##################| Elapsed Time: 0:58:53 Time: 0:58:53\n" ] } ], "source": [ "for bin_name, data in pars.data:\n", " print(f\"Exploring {bin_name}, this may take a while...\")\n", " rslt = explore(pars.bin_fname(bin_name), data, pars.n_ex, pars.bounds,\n", " r_ref=pars.r_ref)\n", " rslt.write()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can load the results for one of the bins and analyse them. The standard output name can be accessed from the `ExplorationResult` class, the output location is specified in the parameter file:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from corbass.exploration import ExplorationResult\n", "import numpy as np\n", "\n", "with np.load(f'{pars.bin_fname(1450)}{ExplorationResult.suffix}') as fh:\n", " posterior = fh['posterior']\n", " scale = fh['scale']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To demonstrate how this result is used, we focus on the marginal posterior for the error level and the residual scaling. First we get the corresponding parameter space arrays and build a meshgrid:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from corbass.utils import bounds2grid\n", "lams, errs, ress = bounds2grid(pars.bounds, pars.n_ex)\n", "ER, RE = np.meshgrid(errs[0], ress[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we integrate out the non-dipole scaling, by means of a simple Riemann sum:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "marginal = posterior.sum(axis=0)*lams[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the marginal distribution:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEGCAYAAACgt3iRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAabklEQVR4nO3df9Rd1V3n8feHECAFUoj8MBJsaY11oCqUDEVbK4gWrFqwIzWzqsRazZRBrc5oC85MrePSqWPtUuwCZVVLaHEwxVYyrNKaRIrTGQqGlpZfRWKDEEFif0FaMZDwnT/uDr2GJ8+9zz25eZ775P1a66y7z757n7PPcyDfe87Ze59UFZIkdXHQbDdAkjT5DCaSpM4MJpKkzgwmkqTODCaSpM4Onu0GjMshObQO4/DZboYkTZTtfPkLVXXsTOvN22ByGIfz8pwz282QpImyoa7/+1HqeZtLktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLU2ViDSZKjklyf5HNJ7kvyXUmWJFmf5IH2eXRf+cuSbE5yf5Jz+/JPT3JX++7yJBlnuyVJMzPuK5PfBz5aVd8GfCdwH3ApsLGqlgMb2zpJTgZWAqcA5wFXJFnQtnMlsBpY3pbzxtxuSdIMjC2YJFkMvAr4Y4CqeqqqvgKcD6xpxdYAF7T0+cB1VbWjqrYAm4EzkiwFFlfVrVVVwDV9dSRJc8A4r0xeBPwT8L4kn07y3iSHA8dX1aMA7fO4Vv4E4OG++ltb3gktvWf+cyRZnWRTkk1Ps2PfHo0kaa/GGUwOBl4GXFlVpwFfo93S2oupnoPUNPnPzay6qqpWVNWKhRw60/ZKkkY0zmCyFdhaVbe19evpBZfH2q0r2ue2vvIn9tVfBjzS8pdNkS9JmiPGFkyq6h+Bh5O8pGWdA9wLrANWtbxVwA0tvQ5YmeTQJCfRe9B+e7sVtj3Jma0X10V9dSRJc8DBY97+zwPXJjkE+DzwRnoBbG2SNwEPARcCVNU9SdbSCzg7gUuqalfbzsXA1cAi4Ka2SJLmiPQ6SM0/i7OkXp5zZrsZkjRRNtT1d1TVipnWcwS8JKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKmzsQaTJA8muSvJnUk2tbwlSdYneaB9Ht1X/rIkm5Pcn+TcvvzT23Y2J7k8ScbZbknSzOyPK5Ozq+rUqlrR1i8FNlbVcmBjWyfJycBK4BTgPOCKJAtanSuB1cDytpy3H9otSRrSbNzmOh9Y09JrgAv68q+rqh1VtQXYDJyRZCmwuKpuraoCrumrI0maA8YdTAr4yyR3JFnd8o6vqkcB2udxLf8E4OG+ultb3gktvWf+cyRZnWRTkk1Ps2MfHoYkaToHj3n7r6iqR5IcB6xP8rlpyk71HKSmyX9uZtVVwFUAi7NkyjKSpH1vrFcmVfVI+9wGfBg4A3is3bqifW5rxbcCJ/ZVXwY80vKXTZEvSZojxhZMkhye5MjdaeDVwN3AOmBVK7YKuKGl1wErkxya5CR6D9pvb7fCtic5s/XiuqivjiRpDpj2NleS1w2xjX+pqo9MkX888OHWi/dg4E+r6qNJ/gZYm+RNwEPAhQBVdU+StcC9wE7gkqra1bZ1MXA1sAi4qS2SpDkivQ5Se/ky+SK9q4DpxnW8qqpevK8b1tXiLKmX55zZboYkTZQNdf0dfUM5hjboAfxNVfXT0xVI8oGZ7lSSNL8MembynkEbqKqf2EdtkSRNqEHB5Ir90gpJ0kRzokdJUmeDnpm8KMm6vX1ZVa/dx+2RJE2gQcHkn4Df3R8NkSRNrkHBZHtV3bJfWiJJmliDgsmD+6MRmp8OOuSQ0SoePOKUcbueGanaM089Ndr+arT9SfPRtP/XVtWzI+CTfDfwwv46VXXN2FomSZoYQ/0ETPJ+4MXAncDuKU52v1tEknSAG/Z+wgrg5Jpu7hVJ0gFr2HEmdwPfOM6GSJIm17BXJscA9ya5Hb7+CkPHmUiSYPhg8o5xNkKSNNmGCiaONZEkTWfaZyZJbhy0gWHKSJLmt0FXJq+cbm4uei/NOnkftkfzyYiDD/OCZSPV23XUopHqHfzY4yPVe+bRx0ar9+STI9WT5rJB/7efP8Q2Rhw+LEmaLwaNgPdZiSRpIN9nIknqzGAiSerMYCJJ6myoYJLkh5N8OsmXkjyRZHuSJ8bdOEnSZBi27+bvAa8D7nKyR0nSnoa9zfUwcLeBRJI0lWGvTN4KfCTJLfzriR7fPZZWaV6oHTsGF5rCQSP+Znn41YePVO/JE0d7I+Q33rx0pHpL/mrLSPV2bfunkerVrl2DC0kdDRtMfhP4KnAYMOK7WCVJ89WwwWRJVb16rC2RJE2sYZ+ZbEhiMJEkTWnYYHIJ8NEkT9o1WJK0p4HBJEmAU6rqoKpaVFWLq+rIqlo8zA6SLGhjVG5s60uSrE/yQPs8uq/sZUk2J7k/ybl9+acnuat9d3lrkyRpjhgYTFp34A932MdbgPv61i8FNlbVcmBjWyfJycBK4BTgPOCKJAtanSuB1cDytpzXoT2SpH1s2Ntcn0zyb2e68STLgB8C3tuXfT6wpqXXABf05V9XVTuqaguwGTgjyVJgcVXd2gLbNX11JElzwLC9uc4G3pzkQeBr9F6KVVX1HQPq/R69MSpH9uUdX1WP0tvAo0mOa/knAJ/sK7e15T3d0nvmP0eS1fSuYDiM5w0+KknSPjFsMPnBmW44yQ8D26rqjiRnDVNliryaJv+5mVVXAVcBLM4SR+tL0n4yVDCpqr9P8kpgeVW9L8mxwBEDqr0CeG2S19Ab7Lg4yQeAx5IsbVclS4FtrfxW4MS++suAR1r+sinyNceNOvL6mS0Pj1Tv2E8fPbjQFF73+g0j1XvLax4cqd63XPfmkeq95KojBxeawq7ND864Tu18eqR96cA17KzBvwa8DbisZS0EPjBdnaq6rKqWVdUL6T1Y/6uq+glgHbCqFVsF3NDS64CVSQ5NchK9B+23t1ti25Oc2XpxXdRXR5I0Bwx7m+tHgdOATwFU1SNJRvuZBO8E1iZ5E/AQcGHb5j1J1gL3AjuBS6pq90/bi4GrgUXATW2RJM0RwwaTp6qqkhRAkhnNqFdVHwc+3tJfBM7ZS7nfpDcP2J75m4CXzmSfkqT9Z9iuwWuT/BFwVJKfBTbwr7v7SpIOYMM+gH9Xkh8AngBeAry9qtaPtWWSpIkxVDBJ8ttV9TZg/RR5kqQD3LC3uX5girwZjz2RJM1P016ZJLkY+I/Ai5J8tu+rI4H/O86GSZImx6DbXH9Krxvu/6BNyNhsr6ovja1VOqA9s+NfRqr3vI99dnChKXxoxFf1bP+Nvx6p3i0X/s5I9c7e8daR6i3/k2dmXGfXlodG2peDHQ9c0waTqnoceBz49/unOZKkSTTsMxNJkvbKYCJJ6sxgIknqbFBvru1MPd377veZDPXqXknS/DboAfyokzlKkg4gw070CEB7K+Jhu9erarT+g5KkeWXY95m8NskDwBbgFuBBnAZektQMe2XyG8CZwIaqOi3J2Tj2RHPMqIMdF924aaR6t37ulJHqffAHzx6p3kHfMFI1/vlblsy4zvO+9JWR9rXzS18eqR4184GVmluG7c31dHsPyUFJDqqqm4FTx9guSdIEGfbK5CtJjgD+Grg2yTZ6b0OUJGnoK5PzgSeBXwI+Cvwd8CPjapQkabIM+3Ksr/WtrhlTWyRJE2rYl2P1D148BFgIfM1Bi5IkGP7K5F8NXkxyAXDGWFokSZo4I83NVVV/AXzfPm6LJGlCDXub63V9qwcBK5h6zi5J0gFo2K7B/T23dtIbAX/+Pm+NNAtq166R6u28f/NI9b7pwa0j1TvomJkPPgTg0ENGqzeCLFgwUr3a6aDFSTfsM5M3jrshkqTJNWgK+j9gmttZVfUL+7xFkqSJM+gB/CbgDnozBb8MeKAtpwKj3RuQJM07g95nsgYgyU8BZ1fV0239D4G/HHvrJEkTYdiuwd8E9I81OaLlSZI0dG+udwKfTnJzW/9e4B1jaZEkaeIM25vrfUluAl7esi6tqn8cX7MkSZNk2ttcSb6tfb6M3m2th9vyTS1vurqHJbk9yWeS3JPk11v+kiTrkzzQPo/uq3NZks1J7k9ybl/+6Unuat9dniSjH7IkaV8bdGXyn4DVwO9O8V0x/ZQqO4Dvq6qvJlkIfKJd3bwO2FhV70xyKXAp8LYkJwMrgVPoBa4NSb61qnYBV7Z2fBL4CHAevjZYE2rUN0I+8w+PjLbDzHzWpBw02u+1esaJMQ5Ug3pzrW6fM37PaFUV8NW2urAtRW/k/Fktfw3wceBtLf+6qtoBbEmyGTgjyYPA4qq6FSDJNcAFGEwkac4Y6idLkguTHNnS/zXJh5KcNkS9BUnuBLYB66vqNuD4qnoUoH0e14qfQO8W2m5bW94JLb1n/lT7W51kU5JNT7NjmEOTJO0Dw17//req2p7klcC59K4o/nBQparaVVWnAsvoXWW8dJriU11X1zT5U+3vqqpaUVUrFnLooOZJkvaRYYPJ7tHuPwRcWVU30HtJ1lCq6iv0bmedBzyWZClA+9zWim0FTuyrtgx4pOUvmyJfkjRHDBtM/iHJHwGvBz6S5NBBdZMcm+Soll4EfD/wOWAdsKoVWwXc0NLrgJVJDk1yErAcuL3dCtue5MzWi+uivjqSpDlg2EGLr6d3VfGuqvpKu6L4lQF1lgJrkiygF3jWVtWNSW4F1iZ5E/AQcCFAVd2TZC1wL71p7i9pPbkALgauBhbRe/Duw3dJmkPS63Q1RMHe85LlbQDjscARVbVlrK3rYHGW1Mtzzmw3Q5p9k9A1uHyfyVyxoa6/o6pWzLTesL25fo1e993LWtZC4AMz3ZkkaX4a9jbXjwKnAZ8CqKpHdncVljTHjfCrv3zBhGZo2Ovfp9ogxAJIcvj4miRJmjTDBpO1rTfXUUl+FtgAvHd8zZIkTZJhZw1+V5IfAJ4AXgK8varWj7VlkqSJMewzE1rwWA/PTpPyhqq6dmwtkyRNjEEDDxe3aeHfk+TV6fk54PP0xp5IkjTwyuT9wJeBW4GfoTdQ8RDg/Kq6c8xtkyRNiEHB5EVV9e0ASd4LfAH45qraPvaWSZImxqDeXE/vTrSpTbYYSCRJexp0ZfKdSZ5o6QCL2nrovf9q8VhbJ0maCIPetLhgfzVEkjS5Zj4DnCRJezCYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjobWzBJcmKSm5Pcl+SeJG9p+UuSrE/yQPs8uq/OZUk2J7k/ybl9+acnuat9d3mSjKvdkqSZG+eVyU7gP1fVvwHOBC5JcjJwKbCxqpYDG9s67buVwCnAecAVSXa/g/5KYDWwvC3njbHdkqQZGlswqapHq+pTLb0duA84ATgfWNOKrQEuaOnzgeuqakdVbQE2A2ckWQosrqpbq6qAa/rqSJLmgP3yzCTJC4HTgNuA46vqUegFHOC4VuwE4OG+altb3gktvWf+VPtZnWRTkk1Ps2NfHoIkaRpjDyZJjgD+HPjFqnpiuqJT5NU0+c/NrLqqqlZU1YqFHDrzxkqSRjLWYJJkIb1Acm1VfahlP9ZuXdE+t7X8rcCJfdWXAY+0/GVT5EuS5ohx9uYK8MfAfVX17r6v1gGrWnoVcENf/sokhyY5id6D9tvbrbDtSc5s27yor44kaQ44eIzbfgXwk8BdSe5seb8KvBNYm+RNwEPAhQBVdU+StcC99HqCXVJVu1q9i4GrgUXATW2RJM0R6XWQmn8WZ0m9POfMdjMkaaJsqOvvqKoVM63nCHhJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmdjCyZJ/iTJtiR39+UtSbI+yQPt8+i+7y5LsjnJ/UnO7cs/Pcld7bvLk2RcbZYkjWacVyZXA+ftkXcpsLGqlgMb2zpJTgZWAqe0OlckWdDqXAmsBpa3Zc9tSpJm2diCSVX9NfClPbLPB9a09Brggr7866pqR1VtATYDZyRZCiyuqlurqoBr+upIkuaIg/fz/o6vqkcBqurRJMe1/BOAT/aV29rynm7pPfOnlGQ1vasYgK9uqOvv31cN34+OAb4w243YDw6E4zwQjhE8zvnmJaNU2t/BZG+meg5S0+RPqaquAq7aV42aDUk2VdWK2W7HuB0Ix3kgHCN4nPNNkk2j1Nvfvbkea7euaJ/bWv5W4MS+csuAR1r+sinyJUlzyP4OJuuAVS29CrihL39lkkOTnETvQfvt7ZbY9iRntl5cF/XVkSTNEWO7zZXkfwFnAcck2Qr8GvBOYG2SNwEPARcCVNU9SdYC9wI7gUuqalfb1MX0eoYtAm5qy3w20bfpZuBAOM4D4RjB45xvRjrO9DpJSZI0OkfAS5I6M5hIkjozmMySJOe1qWM2J7l0iu/PSvJ4kjvb8vbZaGcXU02ps8f3aVPkbE7y2SQv299t7GqIY5z48wiQ5MQkNye5L8k9Sd4yRZn5cD6HOc6JPqdJDktye5LPtGP89SnKzPxcVpXLfl6ABcDfAS8CDgE+A5y8R5mzgBtnu60dj/NVwMuAu/fy/WvodagIcCZw22y3eQzHOPHnsR3HUuBlLX0k8LdT/Dc7H87nMMc50ee0nZ8jWnohcBtwZtdz6ZXJ7DgD2FxVn6+qp4Dr6E0pM6/U1FPq9DsfuKZ6PgkctXsc0qQY4hjnhap6tKo+1dLbgft47mwU8+F8DnOcE62dn6+21YVt2bMn1ozPpcFkdpwAPNy3vrdpYr6rXYrelOSU/dO0/WrYv8Okm1fnMckLgdPo/aLtN6/O5zTHCRN+TpMsSHInvYHj66uq87mcK9OpHGiGmSbmU8ALquqrSV4D/AW9wZzzyYymy5lQ8+o8JjkC+HPgF6vqiT2/nqLKRJ7PAcc58ee0euP4Tk1yFPDhJC+tqv7nfjM+l16ZzI69TR/zrKp6YvelaFV9BFiY5Jj918T9YuDfYdLNp/OYZCG9f2CvraoPTVFkXpzPQcc5n85pVX0F+DjPfbXHjM+lwWR2/A2wPMlJSQ6h9y6Xdf0Fknzj7heBJTmD3rn64n5v6XitAy5qPUfOBB6vNqv0fDFfzmM7hj8G7quqd++l2MSfz2GOc9LPaZJj2xUJSRYB3w98bo9iMz6X3uaaBVW1M8nPAR+j17PrT6o3pcyb2/d/CPwYcHGSncCTwMpq3SwmRaaeUmchPHuMH6HXa2Qz8M/AG2enpaMb4hgn/jw2rwB+Erir3WsH+FXgm2H+nE+GO85JP6dLgTXpvYDwIGBtVd24x78/Mz6XTqciSerM21ySpM4MJpKkzgwmkqTODCaSpM4MJpKkzgwmmnVJdvXNwHpnpphFeR/u66wkN871bfZt+6eSvKel35zkon203e9pM8be2cYazKTuBUlO3hft0PzhOBPNBU9W1anTFUiyoL7+KmeSHFxVOwdteNhyk6D1/99X3gC8q6reN0LdC4Ab6b1meyjz6Txoal6ZaM5K8mCStyf5BHBhko8n+a0ktwBvSfKCJBvb+xY2JvnmVu/qJO9OcjPw29Ns//D03kfyN0k+neT8ln9b/+R9bb+n7638NNs/Jb33RtzZ2ri85V/U1j+T5P0t70fafj+dZEOS46fY3juS/HJfm367bf9vk3xPy39ekrVt+3/Wtrlij+38DPB64O1Jrk1yRPv7fSrJXf3HtWdbk3w38Frgd9pxvTjJqUk+2cp9OMnRfW189nxN97fSPDDbc+u7uAC7gDv7lh9v+Q8Cb+0r93Hgir71/w2saumfBv6ipa+m98t5wRT7Oov2Lgrgt4CfaOmj6L274nDgl4Bfb/lLgb8dUP7Zbe6xrz8A3tDShwCLgFOA+4FjWv6S9nk0Xx9E/DPA77b0TwHvael3AL/c97fYXeY1wIaW/mXgj1r6pcBOYMUUbbsa+LGWPhhY3NLH0Bv1nGna+mzdtv5Z4Htb+r8DvzfV+XKZ34u3uTQXTHeb68+mWf8u4HUt/X7gf/Z998Hquy22F68GXrv71z5wGL1pM9YC6+lNjfJ64IMDyu/NrcB/SbIM+FBVPZDk+4Drq+oLAFW1+10oy4A/S++dEYcAWwa0HWD3JIR3AC9s6VcCv9+2fXeSzw6xnQC/leRVwDP0pho/HthbW79eMXk+cFRV3dKy1vD1vxc89/xpnjKYaK772oD1fv1zA01XbrcA/66q7n/OF8kXk3wH8OPAf5iu/FS3pACq6k+T3Ab8EPCxdnspTD2V9x8A766qdUnOoncVMsiO9rmLr/+/PNXU4YO8ATgWOL2qnk7yIL1Aube2zsQw50HzgM9MNMn+H70Zl6H3D+InZlj/Y8DPJ8/OAHta33fXAW8Fnl9Vdw1R/jmSvAj4fFVdTm8W1u8ANgKvT/INrcySVvz5wD+09KoZHke/T9C7mqL1uPr2Ieo8H9jWAsnZwAta/t7aup3eK22pqseBL+9+ZkNvksRb0AHHYKK5YNEeXYPfOWS9XwDe2G7l/CQzf8j7G/Rm+P1skrvb+m7X0wtUa4csP5UfB+5Ob/bZb6P3GtR7gN8EbknyGWD3NOfvAD6Y5P8AX5jhcfS7Aji2/U3eRu95xuMD6lwLrEiyiV5Q/hzANG29DviV1lngxfSC3++0fZ5K77mJDjDOGizNI+lNK76wqv6l/UO/EfjWqnpqlpumec5nJtL88jzg5vTeFhjgYgOJ9gevTCRJnfnMRJLUmcFEktSZwUSS1JnBRJLUmcFEktTZ/wfWw/4XCJdnIAAAAABJRU5ErkJggg==\n", "text/plain": [ "