reda_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# REDA Example Notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In case you want to try the newest version from git.\n",
    "# Please be aware that this may break this notebook\n",
    "# !pip install git+git://github.com/geophysics-ubonn/reda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.set_option('display.width', 1000)\n",
    "\n",
    "%matplotlib inline\n",
    "import reda.utils.mpl\n",
    "plt, mpl = reda.utils.mpl.setup()\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(2017)\n",
    "\n",
    "from glob import glob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Usage of reda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get list of data files\n",
    "data_files = sorted(glob('data/pygimli_*.ohm'))\n",
    "# take only each tenth data file\n",
    "data_files = data_files[0:-1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:13,129 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n",
      "2020-12-01 10:47:13,253 - reda.containers.ERT            - %INFO       - Data sized changed from 0 to 741\n",
      "2020-12-01 10:47:13,254 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n",
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:13,382 - reda.containers.ERT            - %INFO       - Data sized changed from 741 to 1482\n",
      "2020-12-01 10:47:13,382 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n",
      "2020-12-01 10:47:13,509 - reda.containers.ERT            - %INFO       - Data sized changed from 1482 to 2223\n",
      "2020-12-01 10:47:13,510 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n",
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:13,646 - reda.containers.ERT            - %INFO       - Data sized changed from 2223 to 2964\n",
      "2020-12-01 10:47:13,646 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n",
      "2020-12-01 10:47:13,786 - reda.containers.ERT            - %INFO       - Data sized changed from 2964 to 3705\n",
      "2020-12-01 10:47:13,787 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n",
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:13,944 - reda.containers.ERT            - %INFO       - Data sized changed from 3705 to 4446\n",
      "2020-12-01 10:47:13,945 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n",
      "2020-12-01 10:47:14,105 - reda.containers.ERT            - %INFO       - Data sized changed from 4446 to 5187\n",
      "2020-12-01 10:47:14,105 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n",
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:14,272 - reda.containers.ERT            - %INFO       - Data sized changed from 5187 to 5928\n",
      "2020-12-01 10:47:14,273 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n",
      "2020-12-01 10:47:14,445 - reda.containers.ERT            - %INFO       - Data sized changed from 5928 to 6669\n",
      "2020-12-01 10:47:14,445 - reda.containers.ERT            - %INFO       - Unified data format (BERT/pyGIMLi) file import\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n",
      "Calculating resistance from apparent resistivity and geometric factors. (r = rhoa_ / k)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:14,629 - reda.containers.ERT            - %INFO       - Data sized changed from 6669 to 7410\n"
     ]
    }
   ],
   "source": [
    "from reda import ERT\n",
    "\n",
    "obj = ERT()\n",
    "\n",
    "for nr, filename in enumerate(data_files):\n",
    "    obj.import_bert(filename, timestep=nr)\n",
    "  \n",
    "obj.compute_K_analytical(spacing=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   a  b  m  n  timestep           k          r       rho_a  valid  id norrec  rdiff  rho_adiff   sigma_a\n",
      "0  4  3  1  2         0   18.849556  16.660311  314.039461    1.0   0    nor    NaN        NaN  0.060023\n",
      "1  5  4  1  2         0   75.398224   4.241190  319.778194    1.0   1    nor    NaN        NaN  0.235783\n",
      "2  6  5  1  2         0  188.495559   1.662108  313.299886    1.0   2    nor    NaN        NaN  0.601646\n",
      "3  7  6  1  2         0  376.991118   0.813218  306.575992    1.0   3    nor    NaN        NaN  1.229682\n",
      "4  8  7  1  2         0  659.734457   0.449058  296.259105    1.0   4    nor    NaN        NaN  2.226883\n"
     ]
    }
   ],
   "source": [
    "# hack together reciprocals\n",
    "import pandas as pd\n",
    "data_rec = obj.data.values\n",
    "dfr = pd.DataFrame(\n",
    "    np.hstack((data_rec[:, 2:4], data_rec[:, 0:2], data_rec[:, 4:])),\n",
    "    columns=obj.data.columns\n",
    ")\n",
    "\n",
    "# need to set the data types accordingly\n",
    "for key, dtype in obj.data.dtypes.items():\n",
    "    dfr[key] = dfr[key].astype(dtype)\n",
    "\n",
    "print(dfr.head(5))\n",
    "\n",
    "# add some noise\n",
    "dfr['rho_a'] = dfr['rho_a'] + 0.05 * dfr['rho_a'].values * np.random.randn(dfr.shape[0])\n",
    "obj.data = pd.concat((obj.data, dfr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   a  b   m   n  timestep            k          r       rho_a  valid  id norrec  rdiff  rho_adiff    sigma_a\n",
      "0  1  2   4   3         0    18.849556  16.660311  314.039461    1.0   0    nor    NaN        NaN   0.060023\n",
      "1  1  2   5   4         0    75.398224   4.241190  319.778194    1.0   1    nor    NaN        NaN   0.235783\n",
      "2  1  2   6   5         0   188.495559   1.662108  313.299886    1.0   2    nor    NaN        NaN   0.601646\n",
      "3  1  2   7   6         0   376.991118   0.813218  306.575992    1.0   3    nor    NaN        NaN   1.229682\n",
      "4  1  2   8   7         0   659.734457   0.449058  296.259105    1.0   4    nor    NaN        NaN   2.226883\n",
      "5  1  2   9   8         0  1055.575132   0.274267  289.509467    1.0   5    nor    NaN        NaN   3.646082\n",
      "6  1  2  10   9         0  1583.362697   0.178931  283.312170    1.0   6    nor    NaN        NaN   5.588756\n",
      "7  1  2  11  10         0  2261.946711   0.121960  275.867204    1.0   7    nor    NaN        NaN   8.199404\n",
      "8  1  2  12  11         0  3110.176727   0.086760  269.838138    1.0   8    nor    NaN        NaN  11.526083\n",
      "9  1  2  13  12         0  4146.902303   0.063809  264.608515    1.0   9    nor    NaN        NaN  15.671840\n"
     ]
    }
   ],
   "source": [
    "print(obj.data.head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])\n"
     ]
    }
   ],
   "source": [
    "# these timesteps where imported\n",
    "print(obj.data.groupby('timestep').groups.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply to all timesteps\n",
    "# obj.filter('K > 1000')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-12-01 10:47:41,530 - reda.containers.ERT            - %INFO       - Data sized changed from 14820 to 14352\n"
     ]
    }
   ],
   "source": [
    "# apply only to timestep 10\n",
    "obj.sub_filter(\n",
    "    'timestep in [3, 7]',\n",
    "    'a > 10 and a < 20 and k > 1000',    \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       a   b   m   n  timestep            k         r       rho_a  valid   id norrec  rdiff  rho_adiff     sigma_a\n",
      "2227   1   2   8   7         3   659.734457  0.289292  190.855950    1.0    4    nor    NaN        NaN    3.456714\n",
      "2228   1   2   9   8         3  1055.575132  0.124387  131.299422    1.0    5    nor    NaN        NaN    8.039450\n",
      "2229   1   2  10   9         3  1583.362697  0.053796   85.178704    1.0    6    nor    NaN        NaN   18.588716\n",
      "2230   1   2  11  10         3  2261.946711  0.022858   51.703278    1.0    7    nor    NaN        NaN   43.748613\n",
      "2231   1   2  12  11         3  3110.176727  0.009497   29.537502    1.0    8    nor    NaN        NaN  105.295861\n",
      "...   ..  ..  ..  ..       ...          ...       ...         ...    ...  ...    ...    ...        ...         ...\n",
      "2941  40  39  32  33         3  1055.575132  0.096384  108.581504    1.0  718    nor    NaN        NaN   10.375209\n",
      "2942  41  40  32  33         3  1583.362697  0.047888   75.130456    1.0  719    nor    NaN        NaN   20.882144\n",
      "2947  40  39  33  34         3   659.734457  0.236961  168.410837    1.0  724    nor    NaN        NaN    4.220112\n",
      "2948  41  40  33  34         3  1055.575132  0.114584  122.195106    1.0  725    nor    NaN        NaN    8.727236\n",
      "2953  41  40  34  35         3   659.734457  0.266089  153.880834    1.0  730    nor    NaN        NaN    3.758146\n",
      "\n",
      "[1008 rows x 14 columns]\n"
     ]
    }
   ],
   "source": [
    "print(obj.data.query('timestep == 3 and rho_a < 200'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histogram plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "a              int64\n",
       "b              int64\n",
       "m              int64\n",
       "n              int64\n",
       "timestep       int64\n",
       "k            float64\n",
       "r            float64\n",
       "rho_a        float64\n",
       "valid        float64\n",
       "id             int64\n",
       "norrec        object\n",
       "rdiff        float64\n",
       "rho_adiff    float64\n",
       "sigma_a      float64\n",
       "dtype: object"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obj.data.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Figure(425.197x141.732)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAACGCAYAAAB0djX4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAANBUlEQVR4nO3de6yk9V3H8fcHCuVSWCi7itx6aJVaqrDgxrS1JkpEqgQlCE3ABvC2tEqVSykXrWxQbDGuBbepuNGyQrC2BiWUGmtC7UVDYnbpUmgjGNuFYmlZYBGolAX26x/Pc+B03cucc2bOM3Oe9ys52Xku88wXPnPme57L/J5UFZIkdWGPrguQJPWXTUiS1BmbkCSpMzYhSVJnbEKSpM7YhCRJnXlV1wVMmiTnA49X1Z1JrgE+WFXPzXOblwJHA3sB7672uvkkAW4EXgS+VlWrk/w0cD5NdpcBj+5unar65nzqU2NE2V8HfB+wL3BuVW1t55v9mBlR/r8A/CzwAnBlVX23nd+b/G1Cs/d2YL/mPcJRwJ5J1gEPA28CvtDOf1VVXZzkPcAxwEHAB4CTgEeq6rMASfYGTqyqX05yYbv9L854rfurak2SW9p13w2cDRwL/BrwuQHW+YMR/v/ok6FmD1BVlwMkWQ0cQvPBMv1aZj9ehv27vyfwW8C9wNPTDWjGa/Uifw/Hzd6/An9TVXduN/9G4A+Bo6vqMuDIJK8BzgX+B3iKptncPPNDiOaDZ3P7+CHgiBnLjgC+0T5+rF03VbVtxrqDrKPhGHb2JDk0yV8AhwNPzFhk9uNn2PkvAw6oqvcDW5KcNGNZb/J3T2j2tu1k/tPAge2/0wL8d1Wt2sX2ngCWto+PAr48Y9kjwPHt42XtutuS7NGu+8iA62g4hp09VfUt4IIk7wfeCny+XWT242fY+T8JTB8u2wIcMGNZb/K3Cc3evcDvJtnt/7uqeibJvydZQ/Om/BhwHPCNqrqrXWdrknuS3AC8GvhokrOA54FPAee0y+5t110L/CXN+aPLaQ7f7G4dDcdQs0/yamA1UMD+wEfMfqyN4nf/C21+BwDv6WP+cew4SVJXPCckSeqMTUiS1BmbkCSpMzYhSVJnbEKSpM5M1CXaS5curampqa7L0Dxs2LDh8apaNpfnmv/km2v+Zj/5dpb9RDWhqakp1q9f33UZmockD831ueY/+eaav9lPvp1l7+E4SVJnJmpPSJKmTV3x6V0u3/ShUxeoEs3HyJpQkmOAq4Dbga3ACcASmqEkLqYZh6mANcAf04ydtL6qPjOqmiRJ42Vkh+Oq6kFgXTt5clVdC9xPM+DekVV1PTDVTn+5Xf6OUdUjSRo/XZ8Tqp38+7IkK5OsT7J+8+bN2y/WImf+/WX2/TCyJpTkUOBM4DRgY5KrgDfTjET7cJKLgE3t9HFJrgT+efvtVNXaqlpRVSuWLZvTlb2aYObfX2bfDyM7J9TeJ+XCnSz+8HbTl42qDknS+Or6cJwkqcdsQpKkztiEJEmdsQlJkjpjE5IkdcYmJEnqjE1IktQZm5AkqTM2IUlSZ2xCkqTOLMj9hJKcBrwVOAj4NvAssH9VXZNkFfAM8FhV3bIQ9UiSxsNC7Qk9B7wW2A84qKpWAyQ5GNjWTp+4oyc6km6/mX9/mX0/LFQTeiPNYKY3AT/VztvtbRzAkXT7zvz7y+z7YaFu7/048Hs0d1a9PsmlAFW1JcmeSS4B7lmgWrSIeItnabItSBOqqk/sYtmqhahBkjR+vDpOktQZm5AkqTM2IUlSZ2xCkqTODNSEktw1yDxJkmZjl1fHJdmH5gumS9svlqZddCBw+IhrkyQtcru7RPsC4CLgMGADrzShp4GPjK4sSVIf7LIJVdUNwA1J3ltVaxaoJklSTwz0ZdWqWpPkbcDUzOdU1c0jqkuS1AMDNaEktwBvADYCL7WzC7AJSZLmbNBhe1YAx1bVDgcZ3Z0kRwPn0dzCYTOwFG/lIEm9N+j3hO4HDp3H61wAPNk+Xu6tHDQo8+8vs++HQZvQUuCrST6T5I7pn1m8zj7AXcB9wOntPG/loN0y//4y+34Y9HDcqnm+zjrgXGAvYJW3cpAkweBXx31+Pi9SVRtpLmrY0bJV89m2JGlyDXp13DO8crhsb5o9mu9U1YGjKkyStPgNuid0wPTjJAF+EXjLqIqSJPXDrEfRrsbtwCnDL0eS1CeDHo47Y8bkHjTfG/ruSCqSJPXGoFfHnTbj8YvAJppDcpIkzdmg54R+ZdSFSJL6Z9Cb2h2R5B+SPNb+3JbkiFEXJ0la3Aa9MOEm4A6a+wodBnyqnSdJ0pwN2oSWVdVNVfVi+7MOcBwNSdK8DHphwhNJ3gV8vJ0+G3hiNi+U5FSagUxvw1G0JUkMvif0q8A7gW8BjwJnAucP+iJJTqAZxPRrzHIUbUnS4jXontA1wHlVtQUgyWuBP6FpToP4OeB/gROAbe28gUbRTrISWAlw1FFHDfhyWizMv7/Mvh8G3RM6broBAVTVkzQNZSBV9UdVdT3wJWDdzFG0gV2Oou1w7v1m/v1l9v0w6J7QHkkO3m5PaNDnvqyqLtrBvFWz3Y4kaXEYtJGsBu5O8nft9FnAtaMpSZLUF4OOmHBzkvXASe2sM6rqq6MrS5LUBwMfUmubjo1HkjQ0s76VgyRJw2ITkiR1xiYkSeqMTUiS1BmbkCSpMzYhSVJnbEKSpM7MeuiduUjyk8DbgB8CPgscDSwBLgcuphnUtKrqhoWoR5I0HhZkT6iqvlhV1wGbgLOq6lrgfuB44Mh2cNOpHT03ycok65Os37x580KUqzFi/v1l9v2wIHtCAEnOAb4OHLKTVXZ4K4eqWgusBVixYsUO19FwTV3x6f83b9OHTu2gEvPvM7Pvh4U6HPdO4Fzgn4CNSa6iORx3C/Bwkoto9pI0BDtqIjvTVXORJFigJlRVnwQ+uZPFH16IGhar2TQcSRo3C3Y4TuPJJiapSzahCWLDkLTY2IQW2KAn/W04kvrAJjQGbDiS+soREyRJnXFPSAMZp+8OabQG2TM3+8Wly8xtQiPkYbbu+YHaX8P4/VuI98bu6lzs70+bkHpvGB8Cw2h2C7UNDc+k/JEzzu8LzwlJkjozFntCSU4EzgD2Az5QVd/puCTpZeP8V6RGa1Kyn5Q6d2QsmhBwNnAlze0eTgZu77QaDcSLFWZnGB8U4/Jh0/fzGLMx38wWe+bj0oRgJ6NoJ1kJrGwnn03ywHarLAUeH2VhQ9KLOnPdbld53ay2tzjyn4QaYYh17uJ9MHD+iyR7mIw6513jXH/3U9X9COlJfgw4neZw3NVV9ewsnru+qlaMqrZhsc7RmIR6J6FGmJw6p01KvZNQZ5c1jsWeUFVtADZ0XYckaWF5dZwkqTOLoQmt7bqAAVnnaExCvZNQI0xOndMmpd5JqLOzGsfinJAkqZ8Ww56QJGlCjcWFCXM1rl9yTXIMcBXN952mgCXAI1X1V0nOo7kccv+quqazIoEkbwZOAY4G1s+sK8kq4Bngsaq6pbsqd8zs58fsR2MS8h+37Cd9T+hsYBVN4Cd3WskMVfUgsK6d3AJsBfZpp5dX1WqAJActeHEzVNVXgG8DR2xX18HAtnb6xA5L3BWznwezH41JyH/csp/oPaHWWJ/Uqqq/BkhySZLXz1zUUUnfo6puTfIU8OvTs3by7zga59rMfrTGuTZgvPMfp+wnvQn9Lc1fRPsBV3dbyiuSHAqcCeybZAlwGHA48AiwMcmlAFX1VGdFAkneARwHvAG4fUZdW5LsmeQS4J4ua9wFs58Hsx+NSch/3LL36jhJUmcm/ZyQJGmC2YQkSZ2xCUmSOmMTkiR1xiYkSeqMTUiS1BmbUIeSTCV5LsnG7eb9Y5IHkjyY5MpZbnPfJBuTbE2ydOhFayjMvr/M/nvZhLr3X1W1HCDJHsBtwI1V9UbgR4EV7W2OB1JVz7Xb++YIatVwmX1/mX3LJjRESX4zyf1JHkry3jls4hRgU1XdAVBVzwMXAu9r/1L6jyTr2r+Ubk3yM0n+Lcl/JvnxYf63aHbMvr/Mfn5sQkOS5JdoBlM8AXgLcHWS2Q6L9Cbg3pkzqupR4EBgb+AHgdXAD7c/5wBvB95HM3KvOmD2/WX28zfpY8eNk98GfqOqXgAeTfICs2/yLwGvmTkjSWjGyHoR+HpV3dfO/wpwV1VVkvtoho1XN8y+v8x+ntwTGoIkewHHt8O4k+QHgCeqamuSVUluSPL7A2zqc8DPt2/AaSfTDCa4DXh+xvyZ09vwD4pOmH1/mf1w2ISG41jgwCSvb08yfhD4sySH07xJngJ+Yncbqap7gS8B1wAk+X7gT1kEu9yLmNn3l9kPwcR30TFxAnAr8HFgf+Dvq2ptko8BvwMsA47c3UaSXAGsAN6V5F9oTk6+DvgocPpoStc8mX1/mf0QeCuHIUhyPXB3VX1iu/mX0pxcPAR4oKrWbLd8Crizqn5kBDVtAlZU1ePD3rZeYfb9ZfbD4Z7QcCwH/nz7mdO3zd2Fl4AlSTZOf2dgvpLsC9wN7EVzzFijtRyz76vlmP28uSckSeqMFyZIkjpjE5IkdcYmJEnqjE1IktQZm5AkqTM2IUlSZ2xCkqTO2IQkSZ35P9lXU0QBZM1LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 425.197x141.732 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import reda.plotters.histograms as RH\n",
    "return_dict = RH.plot_histograms_extra_dims(\n",
    "    obj,\n",
    "    keys=['rho_a', ],\n",
    "    extra_dims=['timestep', ],\n",
    "    Nx=5,\n",
    "    #subquery='timestep in {0}'.format(list(range(0, 100, 10))),\n",
    "    subquery='timestep in [0, 3, 6]'\n",
    ")\n",
    "fig = return_dict[1]['all']\n",
    "print(fig)\n",
    "# print(fig.get_figwidth(), fig.get_figheight())\n",
    "fig.savefig('histogram.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting Pseudosections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAETCAYAAAAoITxqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABe3klEQVR4nO2dd5hkZZX/v6erK3QO0xN6EpNnmAFmyCAiAwoOGFAEVl0Dyoqu4OoaFgy76qq7uirsmtAxof5UREyISJCskgaYwCQm9eTcYbqnp1PV+f1xb9f7nttV1VXVVbf7dp3P89TT773vvW+9t+rc26dOeA8xMxRFURRFUYJE2WhPQFEURVEUJVdUgVEURVEUJXCoAqMoiqIoSuBQBUZRFEVRlMChCoyiKIqiKIGjfLQnoCiKoihK7rz24io+2hrP6Zzn1/Y+wMwrijQlX1EFRlEURVECyNHWOJ59YGZO54SatzQVaTq+owqMoiiKogQQBpBAYrSnMWqoAqMoiqIogYQR59JVYDSIV1EURVECiGOB4Zxew0FEMSJ6lojWENF6Ivq8u/8OItpBRKvd1zJ3PxHRN4hoKxGtJaIzinrRFmqBURRFUZSAUgQXUi+AS5i5i4jCAP5KRH92+z7BzHd7jr8cwHz3dS6A292/RUcVGEVRFEUJIAxGvMD1DNkpkNjlbobdV6Y3uRLAT93zniaieiJqZub9BZ1YCtSFpCiKoigBJQ8XUhMRrbJeN3jHJKIQEa0GcAjAQ8z8jNv1JddNdBsRRd190wDstk7f4+4rOmqBURRFUZQAwgDiWcS1eDjCzGdlHJc5DmAZEdUD+B0RnQLgkwAOAIgAWAngZgD/meubFxK1wCiKoihKQCl0EK8NM7cDeBTACmbezw69AH4M4Bz3sL0AZlinTXf3FR1VYBRFURQlgDCAOHNOr+Egoomu5QVEVAHgUgCbiKjZ3UcA3gTgJfeUewC8y81GOg9Ahx/xL4C6kBRFURQlsBRhFZhmAD8hohAcI8ddzHwvET1CRBMBEIDVAD7gHn8fgCsAbAXQDeA9hZ9SalSBURRFUZQAwuB8YmAyj8m8FsDpKfZfkuZ4BnBjQSeRJarAKIqiKEoQYSBeWP0lUKgCoyiKoigBxFmJt3RRBUZRFEVRAgkhDhrtSYwaqsAoiqIoSgBhAAl1ISmKoiiKEjTUAqMoiqIoSqBwVuJVBUZRFEVRlICRYFVgFEVRFEUJEGqBURRFURQlcDAI8RKuCKQKjKIoiqIEFHUhKYqiKIoSKNSFpCiKoihKACHEWV1IiqIoiqIECKeUgCowiqIoiqIEDHUhKYqiKIoSKJhL24VUuleuKIqiKEpgUQuMoiiKogSUhLqQFEVRFEUJEk4adek6UlSBURRFUZRAUtoxMKrAKIqiKEoA0TRqRVEURVECSVxLCSiKoiiKEiS0mKOiKIqiKIEkoTEwiqIoiqIECc1CUhRFURQlcDBIY2AURVEURQkemoWkKIqiKEqgYEZJrwNTuleuKIqiKIGGkMjxNeyIRDEiepaI1hDReiL6vLt/NhE9Q0RbiehXRBRx90fd7a1u/6ziXrNBFRhFURRFCSAMxwKTyysLegFcwsxLASwDsIKIzgPwFQC3MfM8AG0ArnePvx5Am7v/Nvc4X1AFRlEURVECShxlOb2Ggx263M2w+2IAlwC4293/EwBvcttXuttw+19NRL5EFqsCoyiKoigBhEFIcG6vbCCiEBGtBnAIwEMAtgFoZ+YB95A9AKa57WkAdgOA298BYELhrjI9GsSrKIqiKAElj3VgmoholbW9kplX2gcwcxzAMiKqB/A7AItGNMkioQqMoiiKogQQRl4r8R5h5rOyGp+5nYgeBXA+gHoiKnetLNMB7HUP2wtgBoA9RFQOoA7A0VwnlQ/qQioARDSTiLqIKDTac1HGLypnSrFRGQsahHiOr2FHJJroWl5ARBUALgWwEcCjAK52D3s3gD+47Xvcbbj9jzAzF+4a06MKTB4QUQsRvWZwm5l3MXO1a3bzaw7LiWiPj+9HRPQVIjrqvr7iV6BWqVKicvavRLSdiI4R0T4ius39VacUgVKUMfc9zyCiJ1xl7SARfdjP9y8UgxaYXF5Z0AzgUSJaC+A5AA8x870AbgbwUSLaCifG5Yfu8T8EMMHd/1EAtxT6OtOhDwYlW26AE3W+FM598xCAHQC+O4pzUsYf9wD4sWu6boST1fAvAG4d3Wkp4wUiagJwP4B/hSNfETgukUCSjVUlF5h5LYDTU+zfDuCcFPt7AFxT0ElkiVpgcoSIfgZgJoA/utr7vxHRLCLiwV+KRPQYEX2RiP7uHvNHIppARD93f1k+Zy/2Q0SLiOghImolos1EdK3VdwURbSCiTiLaS0QfJ6IqAH8GMNUdv4uIphJRGRHdQkTbXCvJXe4/AVhzvMH9ZbufiD6ew6W/G8DXmXkPM+8F8HUA143081RSU6pyxszbmLl9cFoAEgDmjfDjVFJQqjIGx0rwADP/nJl7mbmTmTcW4CP1HWYqhgUmMIyvq/EBZn4ngF0A3uCaWv8nzaFvBfBOOClmcwE8BeDHABrh+BM/CwDuDfwQgF8AmOSe9x0iWuyO80MA72fmGgCnwPEvHgdwOYB97hyqmXkfgA/BsZJcBGAqnMWGvu2Z18UA5gO4DMDN5JqPieiVRNSe4dKXAFhjba9x9ylFoITlDET0diI6BuAIHIvf9zIdr+RHCcvYeQBaXaXskKuUzcz8aY1dirCQXWAYX1cztvix+2uyA84vjG3M/Bc3gvvXMCa61wNoYeYfM/MAM78I4DcwJrl+AIuJqJaZ25j5hQzv+QEAn3atJL0APgfgapIxBJ9n5uPMvA7OQ+htAMDMf2Xm+gxjV8PJ7x+kA0A1kcbBjDLjTc7AzL9g5loAC+C4KA9m+2EoRWG8ydh0OBblD8OxQO0A8MssP4sxBQMFLyUQJFSBKR72Q/dEiu1qt30SgHOJqH3wBeAfAUxx+98C4AoAO4nocSI6P8N7ngTgd9Y4GwHEAUy2jtlttXfC+XWTDV0Aaq3tWgBdfkWbK2kZb3KWhJm3AFgP4Du5nqsUlPEmYycA/I6Zn3PjNz4P4BVEVJfl+WMIKmkLjAbx5kch/2nvBvA4M1+a8o2YnwNwJRGFAdwE4C44Ofep5rAbwHuZ+W/eDstPPQPAJrc9E8C+LOe5Ho45/1l3e6m7TykepShnXsrhuC2U4lCKMrbW856B/RHmZCGNL6tKLowvdcw/DgKYU6Cx7gWwgIjeSURh93U2EZ1MRBEi+kciqmPmfgDH4AQ1Ds5hgudXw3cBfImITgKS+fxXet7v34mokoiWAHgPgF9lOc+fwkmhm0ZEUwF8DMAdeV2xki0lJ2dE9E9ENMltLwbwSQAP53vRyrCUnIzBcTe9mYiWucrUvwP4q+siCxyFroUUJMbX1fjHfwP4jGvezCX6fQjM3AknCO2tcH5BHIBTzTPqHvJOAC1uUOMH4Jhkwcyb4Phtt7vzmArg/+CkoT5IRJ0AngZwructHwewFc4/ha8x84MAQEQXElEX0vM9AH8EsA7ASwD+BA2uLDalKGcXAFhHRMcB3Oe+PjWSa1cyUnIyxsyPwJGpP8Gp9TMPwNtHcu2jRbFqIQUF0hCG0sA1u+4AEGZTkEtRCorKmVJsVMYMU5Y08jt+kdJjl5avL7vrec6ylMBYR2NgFEVRFCWAMAPxcWZVyQVVYBRFURQloIw3t1AuqAJTIjBzCzDOFgFQxhwqZ0qxURkzODEwpRvKqgqMoiiKogSUQtdCChKqwCiKoihKACn1dWACocBEKMYxqgIAZFy53ttHaTfksZnOG9JntofIjRhTdrHdV5ZLX+q299ihc7H6vBbGLMcEeTLUKH0flZntMk9fyNoOlSVEX7gsnmyXk+wLwWx7v4bJsflie9euXeju7sbx48ePMPNE5EEkVMEV5e5iw97svIS17ekTmXwFyupbcGahlubInZef3178NxG3iufLLbOEsMzTFwolm4lwSHTFo2S1PUPGjJzVRnpEX2WozxznWdOsIXpyst3X14fNmzdj8uTJ2L17d15yFqEox1CV62lKMRny+Lefx56HZ7mROfbI30DMHBuvlHJUU2lkrqH8uOiri8qScocOHcLBgwfR19eXhYypC2nME6MqnBdeAQCgsGfK1gONQl5hK095HACQ3ecZky0hHdJnjcNROWYiYo7l8jJPnyXcUdlnb9sPYQCIR6yHcgSyL5b+gW0fm/D2xdg6zvPP2NrmsFQoKGq2yyJx0ReJmGzGWKRf9NXGepPt+ugJ0Te10qwd1RiWN3ZNyNz0YZLv9/HFDwAAEokEbrrpJkQiEdx///2oq6vbiTypKK/FK6a4y0H0y2vgXvNPjvv6ZF+fOZbjcp5g6zPMQbl5aNWvsz620Fxads3wB2WDUOilzJOlmIh7EQDFjMBSRYXo4wZTzaK3uUb0dcw2Qt/hqV9duag92b54xhbRd3b1jmQ7RvJ7v2beKgBAS0sLli9fjltvvRU33ngjiCgvOYuhCufSq/M5VSkklmxShv8NXvmjBrPWXt+0BtHXenLMtM+Qz4Hly0yx66ubnhN9b5izKtn+9re/ja9+9avYtGkT5syZk5WMjbf6RrkQCAWGYCknXmGzlZbyDMqNt8/WpMvlmLbS4tWy7e0hfZbSYissgFdJ8fQJJcWrwNhjePqiqY8DpNISj8p/nLbSwp4+W2mhiFRgbKXFVlgAoCJq/qnXROU/+MaYUUwmxuT6UlOj7WZ8z69fu27HJ5fcBy+Dysvq1atx//33o7a2dsgxOcGcVFx4wLO8hK2YZFRSEp6+ElpnKYN1lDJYUhAOy76wEWaukv9ABurM9okmeV73FPMeiWlSUV46eW+yPb9C1oa0LX3XzHt+yNwHlZdPfOITuPHGG4f0KwHAI5uU4X8DVVWadp18pvTNMEpL24KY6LOVltecLqusvGmCkas3zFmXcoqDysujjz6K2bNnpzzGi6ZRK0oAKbjyoigpUOVF8YN8lJdB1IWkKAFClRfFD1R5UfxgJMrLYCmBUkUVGCVQqPKi+IEqL4ofjER5GURjYMY6REk/+ZBYFisGZmif5Wcv98ayZAjizRDnkilQ145tyRgD441zyZA9kbkvdRsAEtH0gbqJmBWn4YlzseNeysOeQN2oiQupGBKoawJuJ8RkMO7UChOo2xyRBV/tiHzvL4kPLnpMzjuN8nLLmrckj5myuP5M5AszeDB4t1/GwNgxMd76YZwooTiXXLADdzMFSkZlABdVmxiEeEO16DsxxcQddM6Q99iJOSb26sxZu0XfRfUvJ9sTy4+JvjfNXSO20ykvj7WYzLcFp0bzlzPFH7IN1LViXgAZ9zIkUHehFah7unw+2nEv75z4N9G3fJYMHAdSKy9/3H5qsn/eKbFhZUzTqBUlAKjlRfEDtbwoflAIy8sgGgOjKGMYVV4UP1DlRfGDQiovYI2BGfsQWWnUnvUkbLdREVKlbZeRM076dOiMa73YqdKZ0qFzSZUWa7143ETh9KnSsPq8qdKhcrMd9qRKR8PGbVQVkanSNWGz1ktdWC4U1lDenWz3s/ysO+ImLfbmxffDSzrl5VV/+UTymA/OMssl3FXWN2SMrEkwuMe9Dq+bKOs06hJzJ2Va6yVTqmqFMcVTtVzYLd5o3EYnJss06s7pZsyuOVI+l8wxqdLnNuwQfTUhk1btdRkB6ZWX53aelGxv6Z2ZbPcmDgwZQxllckmVttZ3GZIqPbU+2W6fJ1Ol204z9/rpS6SMrWgw6dGpXEZAauXlvu2nJPv/a+vrku39vT9JOYYNQ2NgFGVMopYXxQ/U8qL4QUEtLxZqgSkCRBQD8ASAqPs+dzPzZ4noDgAXARiM5ryOmVcXax5KMFHlRfEDVV4UPyiW8qJBvMWjF8AlzNxFRGEAfyWiP7t9n2Dmu4v43kqAUeVF8QNVXhQ/KJbyMogqMEWAnVzTwXXjw+4rvwABQjK+ZWhJADsGxru0v7Udyr4kgEiV9qZDhwuQKu2taZQpVboQNY28qdJRq4CiJ84lbKVOV0ZlqnR11MS5NES7Rd+Uis5k2y4PAABN5aYv7qlW+aFFj4jtdMrL9c9dJ467xaof0h6XaZD5wswmXdqbGm3FuQypd1RKDCluai1j4E1VjZil/u2YF0DGvcQnyJpGJ5pNX8dseb8fW2g++wUL94m+N042sS3Twq2iz7t8ezrlZdvuZnHcL1pflWy39pt5dcU3QhkD5JkqDVHTqF502SUCWpfJZ+fFZ5pU6Ws8NY2umPPSkOmlUl6ebJkrjrmn7RXJ9oGNk5Lt/h5PiY0UlPpCdkXNvyKiEBGtBnAIwEPM/Izb9SUiWktEtxFRNM25NxDRKiJa1ZfoSXWIMs4YDcuLLWf9rHJWCvhteREyht7hT1DGBcW2vAySAOX0Gk8UVYFh5jgzLwMwHcA5RHQKgE8CWATgbACNAG5Oc+5KZj6Lmc+KlMVSHaKMI0bLbWTLWZhUzsY7o+E2EjKGlL/XlHGGX8pLMSCiGUT0KBFtIKL1RPRhd//niGgvEa12X1dY53ySiLYS0WYieq1fc/UlC4mZ24noUQArmPlr7u5eIvoxgI8PPwIZF1AuqdIh203kSYcOpXchiarSntV2E8KF5K0cnaGqdDj1cc526rbzHnZfhlRprwvJriodzpAq7VltNxq2VtsNSxdSTcT8epzgcSFNjBg3UZjkmPb2Bxc+Bi/plJc7t56VPKaqfL445zdHTd+yml3J9oAnTTsnmMGeFXhNXyL1fve8kiFTqnREmrxFqnSNXFFXpEo3yzTqY7PMvXpsnvzcp88/lGxfMmmz6JsVPpxsr5izYcjU0ykvv912erK9u2+JOKe23FjlXmidkWz3JjSBc1TIs6q07TICgL7pVlXphVKptFfYveC0l0Xfuyf9NdnOJVXafpbde3S5OP7JNYuS7aZ15voOyYLqqeGixMAMAPgYM79ARDUAnieih9y+26z/4QAAIloM4K0AlgCYCuAvRLSAmYvuay+aBYaIJhJRvduuAHApgE1E1OzuIwBvAjDUcaiUDBqwq/iBBuwqfuC35WUwCymX17BjMu9n5hfcdieAjQCmZTjlSgB3MnMvM+8AsBXAOSO/uuEppgupGcCjRLQWwHNwYmDuBfBzIloHYB2AJgBfLOIclDGMKi+KH6jyovjBaLmN8lBgmgZjstzXDenGJqJZAE4HMBi/epMbv/ojIho0Y00DYBcf24PMCk/BKGYW0lo4F+7df0mx3lMJDqq8KH6gyoviB6OlvOSZhXSEmc8a7iAiqgbwGwAfYeZjRHQ7gC/AMfx8AcDXAbw31zcvJMFw5BIZ/2amVOlMVaW9adRRq1yAJwYmU1XphIhzyVAuIEOq9JB0aGs7MeQ8tvo8MTDRDKnS1nYo4q0qbWJbvKnSdomA+qh0wk6uMFV8p0RlRV+7qrSX9y18Us47jfJyxRP/Io6bVWV8w0f6ZJzESZUmTfapdpOWeDwu02VzJl2sSynFuXixU1XLPDEI1j03JFXainuJN0kFtXuqWcr92EmeVOkFRl5nLjgo+i5vNmmsC2P7RZ837iWd8vL0zlniuJ6ESV29c7csAHzgsImdiG00cx44NnyKq1IgCpEqPcNTVXqReei2nS7j3i5bZqIa3tz4gujLtqr031vmiGP+3Lo82f7rZhnP1/iCuab6LSbmKtSTIe7OgouQRu2u3fYbAD9n5t8678MHrf7vA7jX3dwLYIZ1+nR3X9Ep3TKWyqiglhfFD9TyovjBWMg2KnQatRuf+kMAG5n5Vmu/vUjSm2HiV+8B8FYiihLRbADzATxbsAvMQDAsMMq4QJUXxQ9UeVH8YCwoL1ycLKQLALwTwDp3HTcA+BSAtxHRMjgupBYA73fmwOuJ6C4AG+BkMN3oRwYSEBQFhpB0HeVbVXrIarui4nT2VaUTYSMsAxkrR8u+RDT1cc772X0eN5HlNvJWnLbdRkOqSltuo4jHTWS7jezVdQGgLmJMmFMrO0Tf5IhxG8XK5Jg2Hz35oSH70ikv3918UfKYT87cI855sPPUZPulNrlC6oUNxpT77BFTMbgvMYI0aqC0XUWD5FTVN32q9MBEo6B2T5fm/WMnmTGPLZQm/IULjPX51ZM2ib45EZMqfdXcF4dMPZ3yYsvZpZVS5t+/ySxb0dklq1+X7zTXV3HYquKeJtteKQD5pko31ou+vulm23YZAUDbmeb5dflSmQj7lkazwu6ls6X8DZJKeVm3a3qy/3uHLhPHP7b6ZDPNF+QzqmGLkcfI3rZkm/qz0wEK7UJi5r8CKU0192U450sAvlTQiWRBMBQYJdCo5UXxA7W8KH4wFiwvhtIuJaAKjFJUVHlR/ECVF8UPxpby4lCMIN6goAqMUjRUeVH8QJUXxQ/GpPICrUY99iEypQAypUpniIEZUhIgQ5xLxqrSmVKlrT5vqvSAlWE6JFXaqiptlwdwxsmuqnTIUy4gU6p0bczEuUyIyfTnpqjZnh5tE32N5V3Jdj/L7+HDJ/9FzjuL8gAA8IYqE+/w956pou+vh016dG1UFlr8wZYLku1lk80Y68s0OGHEeMoF2CU5KCYF264qbce8AMDxmSY+oX2uvDePzzMyefrCFtH3uokmFX5KuF30ZVtV+vubLxTHtQ2Yitfv3Pgu0de5pT7Zjh6V116z29xXVfvNnEO9GitVUAqRKm3FvABA24L0qdJ23Isd8wKkjntJpbw8tGOROObOI69Jtp/YMk/0TXzaXFPdDvksi+w2z1lus+IOB7KIgeHSDtsrZimBGBE9S0Rr3IJQn3f3zyaiZ9zCT78ioshwYynBQi0vih+o5UXxg7FoebHRatTFoRfAJcy8FMAyACuI6DwAX4FTEGoegDYA1xdxDorPqPKi+IEqL4ofjHXlheHEwOTyGk8Us5QAAxj0OYTdFwO4BMDb3f0/AfA5ALdnHIxMNeq8U6XD3lVz86wqLVKec0iVtqtRe1Olw9mlSsObKm25jbyp0rFIplRps8LupGiX6LOrSteVy4rT9SGz/fb5Q9cpSqe8/HG7SYeeEJLzvO+4MbVWlcl5bt8+2bz35E7RF7duRDtFsbP7wSHzUnLDa8IvixphpmpvVen0qdK22+j4Ivndnj7XVBB/46Q1om9OxFScTrXyaTrl5Y6Xz0+213XPEuesazPuyb27Joi+hq1GliqPSLN9xQEz7/Bhc6+UZZniqqQhp1R9k9pOdfIHUd+0+mS71VNVum2Z+Y7s1XUB4B8n/D3ZvnDWtpRTTKW8tOwxyzn8pvU14vi/vGSeQ01/kys1N2wwshNql89V7jDLU/AJa/XzRDYr8ZZ2FlJRV+IlopC7EM4hAA8B2AagnZkHHZK+FX1SiotaXhQ/UMuL4gdj3fJiw5zbazxR1CBedzW+ZURUD+B3ABZlPsPgVsi8AQBi5frPcCwTZOVFyBkqhzlaGU2CqryojAWLICkvQGmnUQ9rgSGiC4ioym2/g4huJaKThjvPhpnbATwK4HwA9UQ0qDilLfrEzCuZ+SxmPitSrjf9WCXIygsg5SyM6PAnKKNCUJUXQGUsSARPedEYmOG4HcBSIloK4GMAfgDgpwAuynQSEU0E0M/M7URUAeBSOAG8jwK4GsCdAN4N4A/DzoAIHHWmOiTOJWKVEsiUKu0tF2DFrwx406jtGBhv5eh8U6UzVZW20qg57E2VNttlnqrS0ZiJJ6mIyNgSO+7Fmyo9KWb8sVOj7aIvROb9wiTfzxv3kk55+dqG14rjLq82laOfPCHTCx9rXZhsrz0o06ivPNMsFb+xfYroa+sxfvHOkPV5jq/7c1SgiKfSslUugOtkDExPs/lx0THLkyq9wFQ2P3fBDtG3vGFzsj2lXJas8Ma9pFNentspf0edZt1j39x6ieg7uq0x2a7bJu/36v0mxTZ2WMbqlB8y8Ql0wuqLZ1cpWLHINlW6QpZzIDtVepqnqvRCI5veVOkLTzMydu0E+exKFfeSSnnZtluWMLlx+7XJ9sZN00Vfw5rUVaUBoPywFedyTMbz8XETE8MD5hoY2fl7NAYmMwNuQO6VAL7FzN8GUDPMOQDQDOBRIloL4DkADzHzvQBuBvBRItoKYAKcqpdKwAi65UUJBkG2vCjBIWiWFxuNgclMJxF9Ek51yguJqAxORlFGmHktgNNT7N8O4JxcJ6qMHVR5UfxAlRfFD4KsvAClHQOTjQLzD3DSnt/LzAeIaCaArxZ3WhIm4zpKRL0r8VpuorDXTZTehSTdRNlXlY5nSpUWfTlUlbbcRORxE5WLqtLSRFodMybtGk+qdEPUmCWnVkgT/aSwMWE2haU5syNuzLcfWPg4vKRTXr740uuTx1xSLV0APzj6ymR7x3GZwtqXMN/naZP3ib493fXJ9oFOafQ70WPp0OWWOb907+WCYVeYBgDUGrdR72TpQuqYZb6HzgVSPs9dtD3ZflOTrBxtVzO/Yo5McQXSKy+2Sf+R40vFOU8enZ9sH90uXQ0N641gxNrkPVa100qPbpfLCvBxk9bK/cYlll2Ka4mTZ1XpoanS5rtsWyRls/UM812+7oy1ou/yBpOen0tV6SdbzArgPzq0Qhy/cbNxG9kuIwBoWmdkJby/XfRxe5pUaXjcRnFLNrOwljDGX1xLLgzrQmLmAwB+AySjz47AyShSSgy1vCh+oJYXxQ+CbnkZhHN8jSeyyUJ6H4C7AXzP3TUNwO+LOCdlDKLKi+IHqrwofjBelBeUeBZSNkG8NwK4AMAxAGDmLQAmFXNSythClRfFD1R5Ufxg3Cgvg5SwCSabGJheZu4j15/pruHi78dAlEyXTkS8adRWSYCwN17Fjo/x9mVIlc7UZ6VHe9OoM6ZKZ6oqbcW5lHtiYOy4F2+qtB33Yse8ALJEwCRPnIuYl0cjv3nx/bI/jfLSv3+uOG5hrD7ZjpGMhdjfY9IgO/s9PuxuE3Ozj6ViNJAw39+Jbvlhx/tNX+ioFQ8zML5+YYwG3oq//U0m/qhrmlwfoHOOkeUF82UMk50q3RiSsSUr5mwQ2+mUl/u3LxbHPdxp4qmeOCBl8NCmicm2N1W64qi5r6p3yXulrM3cH3ZKKwBwrxX30m/dfxoDk5o8q0rbcS/eVOk2K1W6dZn83C847eVk+2pPVelLZm+Gl1TKy8+2nCeO+emeNybbW/bI3+oNq801NW5Kn3JvlwcAgES3JVdx+YyXcS+5/2sdb1aVXMhGgXmciD4FoIKILgXwQQB/LO60lLGAWl4UP1DLi+IH487y4jLeUqNzIRsX0i0ADgNYB+D9AO4D8JliTkoZfVR5UfxAlRfFD8at8oLSjoEZ1gLDzAkA33dfowMB7LqKEkPSofOsKi1SnjOlUcupZEqVjgs3UQ6p0pabKOLpq4waE7Y3VXqitaLupJh0E9luI29V6baBqmT7g4seg5d0yssn1lyTPKa9/4Ih5w1ySpV0JSyt3ZNs/3D3+aIvYbmJQiFpHh7Ya8zMZf3yO6pe2J5s9+y1qiWX8K+RQpFo8KSsNxsT/rHZ8nuomd+WbL9m8kbRN7HcyKDXZQSkV14eazHp0Bt6Z4lz1rSZ2q/tXdINUbPDyFLNbunGrNxnUldDRz0roXaZlaq5R95jdoqrbfrPdpXUcU+mVOmIfHiK9HxrdV1AVpVuW+BxM1tuo/OWvSz63tz0QrKdymUEpFZePvj8O5L9EyOTxfEvbzYrgk98VrrBaluMfET2yeUp0Ga2va5IITseF9KITCgMZ52REiWtAkNE65Ah1oWZTyvKjJRRRS0vih+o5UXxg/FqebEpZRdSJgvM6zP0DQsRzYBTM2kyHEVoJTP/HxF9DsD74LilAOBTzHzfSN5LKQyqvCh+oMqL4geloLwAGHeZRbmQVoFh5p2DbSKaDOBsd/NZZj6UxdgDAD7GzC8QUQ2A54noIbfvNmb+Wr6TVgqPKi+KH6jyovhBySgvJb4S77AxMER0LZzSAY/BWaj9m0T0CWa+O9N5zLwfwH633UlEG+EsgpczTJQsEzCkqnSmOBfRh/R9BUqVtuNeMqVKhzwVp+24F2+qdFXExMDUReQS1E1WqnRTWKap2lWlexOydNV/nCKTyNIpL+945npx3D9NXp1sf3nnFaKvN25EaUpUphD+fofxNk5pkPEHoTIzz/ZuWYW2M25idQYmyM+lt9dcU9Ue812W9UEZIT1TqsS2XWW6b76UwVdPNVWm50QOi75r5j0vttMpL7/eeqY4bpkl199tl/98Nm818QmVLVKuqw4YWarcJSuwh1qtFFdvqrQV98J9Us7ACatp3e8l/Ks3Y6q0FfdClfJ+hqgqXS+62haYB23r6TJG5NzTtibb106SqdJXzZUlKoDUystvt8myfFXl85LtO56V8XwTnjfXVLdNVpUOH7DiXtpkDAx3W2UnBmQM1khTpTNSwrKYTRr1pwGcPWh1IaKJAP4CZ3XerCCiWXAKOz4DZ1G8m4joXQBWwbHStGU4XSkianlR/EAtL4oflI7lxYVLex2YbNKoyzwuo6NZngcAIKJqOLWUPsLMxwDcDmAugGVwLDRfT3PeDUS0iohW9fcfT3WIMkJUefHIGXqHP0HJmVJXXlTG/KHklJdBdCXejNxPRA8A+KW7/Q9w1oIZFiIKw1Fefs7MvwUAZj5o9X8fwL2pzmXmlQBWAkBN/XQerCydU6p0JPVxQP6p0olwBheS7RryuIlst1E4Is2L0bAxW9suI0CusDvBs9ruhLBR7Mo8kllmuZA+vvgBeEmnvPzri/+QPGZBlXy/AwPGBLy7vV70zW08mmy/2D5D9MWtVOnW4zL11b724yfkF5GImWsItUtR7e81Zt72U83nGf9jbneoLWe11DjObu/sWX7Zl5Pt9oXSb9q5yHxHF8zdIfrOrzXmfa/LCEivvNhytrRKfu+f32fck09tkqvtVuw0bqO67fIeq24x90NoSFVpI8tidV1kqAbshRPp+zIQeBnLJVXadhs11ou+3ulmu22BPK9tmfncL1wqK0df3bQq2X7T3DVIRSrl5cBe42786rZ3iuP3tTQl2xNWSTdY/XYrVXpPq+jjTiNXdvo9MIwcFTVVqLAWmAwJOI0AfgVgFoAWANcycxs5y/T/H4ArAHQDuI6ZX0gzdmMWU0gwc3s2c81mHZhPENFVAAbX8F7JzMNWo3Yv6ocANjLzrdb+Zjc+BgDeDOClbCaqFA61vCh+UOqWF8UfStbyMkjhdaN0CTjXAXiYmb9MRLfAWeT2ZgCXA5jvvs6F42U5N83Y+9xXJq0rBGBmNhPNxgIDAH8HEAeQAPDcMMcOcgGAdwJYR0Sr3X2fAvA2IloG52NvgbO6r+ITqrwofqDKi+IHJa+8AAVXYDIk4FwJYLl72E/gJPbc7O7/KTMzgKeJqN5jqLDZyMynp9ifhIiGRmanIZsspH8C8B8AHoHJQvpPZv5RpvOY+a9IrWXpmi+jhCovih+o8qL4gSovyHcl3iYiWmVtr3TdnEPwJOBMtpSSA3BcTICj3Oy2Ttvj7kulwJyfYl8+xwDIzgLzCQCnM/NRACCiCXAsMhkVmELCRMn06ZxSpa0SAYkhfcjQZ8W5hD0xMJmqSlvlAsqGVJU2cQQxT6p0bcz4XOsiMm1vSoVJO54YkSnINSFzbD9LP262VaWvfUoawG6b8USy/epnPyD6Hj7nu8n2wia5FFCkzFzv0R6ZhtvbY+IWBrqlyJ2ImQ/fW0qArc8+fFheX0+d+QzLjlt9JRyRnwt2zAsAtC0yN0T7qVI+LzzVLNH+2sb1ou+d858W2+mUF3vpdgCwCsXjZ3vl82pri1navW61vDmrDlhVpXfKGK1Qqx2f4EmVPmHuFW+qtExxzRDnUkpLnuabKm3FvfTMrBddrScbGes4TX4Hly01kQRvbpThE1fMGRplkEp5uXPrWeKYh9tfk2zv3zpR9DW8ZASwfpsMrI7uNHEvfMxTdsKOpfIzVToDebzVEWY+a7iDvAk4ZMkEMzNR7oVbmLmnEMcMko0CcxSA/S12uvuUgKCWF8UP1PKi+IFaXjwUQVdKlYAD4OCga4iImgEM/ordC8DO3Jju7ss0/llwlmg5CY4eQnD0opxKFGWjwGwF8AwR/QHOR3UlgLVE9FE473hrppOV0UWVF8UPVHlR/ECVlxQU2OqcLgEHwD0A3g3gy+7fP1j7byKiO+EE73akiX+x+Tkc7846OLG1eZGNArPNfQ0yOOmaFMcWhzLjDhrwVI5OZEqVjtl9csiMqdJWenTC02e7jcjjQgpZbiPbZQQAldZ2tbeqdIVVVToq0z9nxIw5M0zSLVUXMubM9y18El7SKS/f3HRJ8pjTamRa8/8efWWyvfykraLvP/a/1oztuWnsdNeJU+QKlfF+Y64l7yrEW4wJur9aftZhuxBwpeyz06q5vIRM+yPAdhvZ5nwA6FhmUosvXCwr/l5lVfxNtfJpOuXljpeNa6i2fKo457nWk8z5z08XfY2W2NXslfdR7JCxLocOyPUvudsy73urSltuo6EprpZMlpKbyKZQqdIzzHbrYo+MnW6+k8uWyOrkb5vwTLKdS1Xp329bmuy/Y+9ycfymDcYo0LhaLl1mu41slxEAsF1V+oRceXr0UqXTk7sjZ1jSJeB8GcBdRHQ9gJ0ArnX77oOTQr0VThr1e7J4j8PMfM9IJ5pNGvXnR/omiv+o5UXxA7W8KH6glpc0FGFxugwJOADw6hTHM4Bcb/7PEtEPADwMmNUdLXdVVqRVYIjof5n5I0T0R6T4iJj5jbm8keIfqrwofqDKi+IHqrxkgoKauPAeAIsAhGFcSAygMAoMgJ+5f7VqdIBQ5UXxA1VeFD9Q5SULgun1PJuZF450kLQKDDM/7/59fKRvMlKYTHzL0JTnTBWnTdtTkBlxO87FmyptV5X2xGzYJQLKPH3lYeMTjYZlip2sKi2zxBoixndfH5bpn9Ey47tPsPTjeuNesq0q/Y6JB5Ltbx2WFsHmCuP/fWrfLNE3q8HEHHhLF1DIbLe2V4u+ck9KuZizXZphkiedcbMJYiqTHyfCx9jqMzJwULqsS5qLX+NJlbbTWJfK5fQvWWKWb7+0QaZKe+Ne0ikv39m0XBz3m/0mBmbbXpnGWrHRfLcTPSUBKg+ZuUUOyeXaqcNKlfZWle61q0p7ygWIStKee1rjXkZQVbpOdLVaZSjsuCoAeN0pJh368gZZEiBV3Esq5eVrG14rjtl6wiStbNwkY6ka1pnnZcMW+cyN7DbPMu44JvrsuJexkiqdkTEyjRz5OxEtZuYNwx+ankwupHXI8NHkmu6kFBe1vCh+oJYXxQ/U8pIDwVRgzgOwhoi2w4mBKXga9evdv4NPqUGX0juQxUeWa0GoXCatSFR5UfxAlRfFD1R5KQlWAIgA6McIVLBMLqSdAEBEl3pqF9xMRC/AKeSUiVwLQqXFcSE57UwVp73upYToyy9VGt5U6agxIYbD0jVSETUm0+qoNJ/WR41ZcnKFNFlOjRq3jb26LgCcUdGSbC+ftQVe0ikvv9hyTvKYTzbLNYU+vevKZPt4v/zQNveaVVBtlxEAbHzaPEzik+X1VdaaeR9vkybn8G4rVbpGfp7xycZEGzos5zJQYbmJ4p40T2u7zM60lV624ECFCcS7+FIrVXqRTGNtP918Z5edJt1E1080KzCfd1LLkHHTKS/3b1+cbL/YJRf33LbJpE5X7ZQuippdRg5qt8rVTss6jSzRMbmsQMJ2G/V7VtTNuqp0MH+yjphcUqWrrNW0G+QPIruq9BAZO9PI2FVL5Yq6H2wyMjZ/RuplQlIpL4+1zE/2t/ScLY5/8DHzr6lut7y+xk1GjiK7PKnS7eYZHIRU6bTkV0pgVCGi9wH4OJylWCoBvAAnTfs7AM7IZaxs1oEhIrqAmf/mbrwCWfybyKMglJIjanlR/EAtL4ofqOUlP4qwDkzRIKKbAZwF4CJmPuDuuwzA95FlBWqbbBSY6wH8iIgGo7XaAbw3lzfJsiCU95wbANwAAOHqhlzeriRQ5aUw2HIWQ+UwR5ceqryMHJWx4VHlZQQESIGBkz59KjMnzafM/CARvQaOBSYnslnI7nkASwcVGGbuGOYUQb4FodzqmCsBoHLijGB9RUVGlZfCYctZLTWqnFmo8lIYVMYyo8pLaWErL9a+g0T0zVzHysYCM/gGOSkuQM4FoTIMZNKl7fIAQOYYmHhs5FWlQxmqSntTpe2q0o0xmf451Yp7mRKVH2VdyPhgy0jOxRv3kk55Of9B6YWbW3dqsn3Yo+DEQta8I9L/W1lurm/vcZkiGZ1vrqF7v0yV7jlqKkucfa6c83P9pszAhOdkLET3FJPf7imoLb8zz2O/Zpf5XnoaLI9m3lU1RoEMaazZkrGq9FL5nHjNqRuT7dd50li9cS/plBd76XYAuPuISZV+ctN80TfhRfO91O2QKfLhVqskQKsnjbXL3DsJLQkwMvJOlTbPDLs8AAC0LTDntS+TMnb5aSZV+soGGQOTKu4llfLy9M5Z4pjfthkZu3f9qaKv3ipy07jZkyq9zyoJ4K0qHbRU6QwEyYUEYBsRvY6Z/2TvJKL/hFPcMacHYdYKTK7kURBKGQa1vCh+oJYXxQ/U8lIgghXE+0EAvyGi6wCsAVANp47SagCpi2BloGgKDHIvCKVkQJUXxQ9UeVH8QJWXAlGEWkjFhJl3EtHZAF4L4GQAxwC8g5nXugG+OTGsAkNElQA+BmAmM7+PiOYDWMjM9w4z0ZwKQmUcq8y4juIya09sZ0yV9vSx5SbKpap0LGK2azyp0rbbaOqQVOl2M77HTRQmY8L84KLH4CWd8mKvSnn5VOk/WxAz5to/t8q1gSpCQ1yQSZ5pmZVsN0+Qrq7uLvNhc4U038ctC/Rzm+QDqXqrEbPOy2RaLG8z6ZpeF1KdVZW4v0r2dU+0UkDtVVbH2o8RO3WVZPIelVnmfU8aayZst1HbQk8a62nmuz13yTbRt6JhXbL9prnShQSkV14SB4xr6Ef75EqoazaaqtINq+UXWL/VuH8ih+T3LlfUle5W4Sbqz2Te1xV1AeQtY1lXlR4iY+Y7ueAU6S5+Y4NZtTnVsg9AauXlvu2nJPt/deQScfzjm4381T8t5zJho3Ebhfe2yzeyq0p7V20OUqr0cARtuk7xx/vdl73/K7mOlc2qGT+Gs1LeoCNyL4Av5vpGSn6o5UXxA7W8KH6glpfCQ5zba1Tn6qwhN+JjBsnGhTSXmf+BiN4GAMzcTVSgVbeUjKjyoviBKi+KH6jyUiSCZYE5mYjWZugnAHUZ+gXZKDB9RFQB92MiorlwLDJKEVHlRfEDVV4UP1DlpYgES4FZlMUxGZbRlmSjwHwWjq9qBhH9HE5w7nXZvkFBsEsJeGJg7HToIanSVho1l3tTpc1n5K2WHI4Y/2ilJwamOmp0N29V6SkVJlVvRkwuXW1XlW4bkAEdn1xyn5x3GuXlkR2y+nh3YkGyfX6V9Dffvu/iZPtQd43o272zKdlumirjXOY3m6z2jS9PE321k03cwrGDMo2arJiYck+Jhb46kyrdf0Sma1Z0WesCeR2aVmyLt89OmbfLCox6DIzXOGnFJAxJYw2b229IGqvF8hXSNdy6OH1V6QsXv5xsX1T/sui7Zt7zYjud8nLnVlkS4M+tFybba7bLir+1G801NGyRv2sih01sC7V5UqWtmIQhlaO1JEBmCiVjVtxLz8x60dVqVy73VJW2K5e/tnGd6Ltizkvwkkp5ebJlrjjm/g5T+uTxl2U6fvUas3ZGvUfGwges9Oh2b1Vp83wOeqp0OsaCWygXBksUFYpsFrJ7yPVJnQfn38OHmflIISehGNTyoviBWl4UP1DLiw8EK426oGSThTRYXGkwrWWmuyrvTmYeSHOakgeqvCh+oMqL4geqvPhEgCwwhSYbF9Jghci1cCwwpwBYD6COiP6ZmR8s4vwAOO6DhGvRHJIqbbmNhqRK224jT6p0Wdhse10e9gq7VRFpPrXdRhOiMv1zcsSYMJvK5cqPPWzcKF849ffwkk55sU2tLf1TxDlrjxkXT5jSm9rLPDbG8mrjzmo/JmuzHDlglKbqSfL6Oo9ZJuiQHLNygzHzDngs1Vxmji3vlL6g2BHT1zVDntdfY35ZVB2Q319/hemrPGquvaxvFO7mTGmstgnfm8Yas5aV9lT8tVOlbZcRABw708jgG5ZIE/5FtWYtKK/LCEivvNhprE90LBPnPPbiycl2/Tr5yGjcZEz60f3ShI/Dxo2aOCHdrbqi7jBkypMokIzZbqOjS6SMdVoydvUpL4q+FXUmBvPS2ZuQilTKi+2a/HvnOeL4P7y4LNlufDYs+upazDM41iJd8yJVeoiMmfMCnyqdgSC5kApNNmnU+wCczsxnMfOZcIoybgdwKYD/KebkSgW1vCh+oJYXxQ/U8uIznONrjEFE/05EfyaiHxBRTg+mbBSYBcy8fnCDmTcAWMTM23OdqDIUVV4UP1DlRfEDVV58Jsc1YMaotWYCgKcBfAnAwmGOFWSjwKwnotuJ6CL39R0AG4goCiDtkq5E9CMiOkREL1n7PkdEe4lotfu6IpfJjjdUeVH8QJUXxQ9UeRklAm6BAdAGp4jjIQCtwxwryCYG5jo4BZg+4m7/DcDH4SgvF6c+BQBwB4BvAfipZ/9tzPy1XCbppFE7n/yQGBi7qnS57LOrStsxLwAQsVKl7fIAAFBtlQiojUq/6sSYSSVujMgYETsOZU9fo+j7/KmyZmU65eXiRz4mjntD8xwzT/JUF55gqgvfe0iWC9jfaZSh2pi8huoqs91+SKZYY8D43o/vk33hNqPvetPZj8+yPs/9UqzsNOeQ/MgQsmJWanfIvsrDZszuJpkeWnXIKvfQZr6vsoEi3aHZLtde7rn2DBV/2a742yyVV7tEgDeN9XWLkwZRnFUtP7RsU6XX7ZLp0F9t+Ydke/sWGWtVu8lckx3zAgDRPe3WpD1prN1WxV9vqrRd/qFUSwLkEudSCBmb6pExq3J55+nye32LFfdyYY2ssZcq7mW4mBcA+P3hM5Ltp9bIVOmGteb+rt8mZSW6q81sZEqV9srYOEmVHpaAXhoRlQEoY+bPE9FUAN8AMDQPPwPZpFGfAPB19+WlK8W+wfOeIKJZuUymVFDLi+IHanlR/EAtL6PLGHULZYSIboKzxlwfER0F8A1mvj7XcYZ1IRHRfCK6m4g2ENH2wVcecx7kJiJa67qYGjK87w1EtIqIVsU9xd6CjCovYwtbzvrH0QLTqryMHcarjAGqvCh58zEApzLzNDiVqS8gos/mOkg2LqQfw9GUboPjMnoPsoudScXtAL4Ax+j1BThWnfemOpCZVwJYCQDRGTN40HVku4wAgMMZUqWtFXbLy2Uane02qvZUla6PGdP35JhMh54aa0+2oySXwakrN6uLfmjRI/CSTnlZv8ukQ/9kgbyGz+5bkWw3x+SquSdX7Eu2tx+dIPq6O415+GiX1BPDTeb6IrXygdrXZtIuOeZZUXeymVv4sBSdsOU26p0tx5z0kEmLPN4sRae/2pjGa3fKz7Or2ZiVbZcRAJDtKhpBaS5bzmqpkdOOle1qp1HpW6Mqk6aeqJcuud6pZjXj9rkyddSuKv2aJRtF3+UNppL0G+bINGogvfKya09zsv29I9L7u71lUrI9YZW8vrodGVKlW9uTTdtlBHjSWBOen4kllCqdVsYyuIkypkrnIGM902wZkynW7ZZr8spTZXXyt9SvSrZfMSv179VUyssft5+a7F/bLd1ET28x7vDGNVLGGl62KpfbbklApkpnkrFxnCqdkWBeZhecmBcw834iuh7AagCfz2WQbBSRCmZ+GAAx805m/hyA1+U2VwdmPsjMcWZOAPg+gHOGO2e8oJYXxQ/U8qL4gVpexghFyELKNQGHiD5JRFuJaDMRvTbLmd8O4NdENM/dngmgO8PxKclGgel1g222ENFNRPRmANXDnZQKImq2Nt+MHAN2gooqL4ofqPKi+IEqL2OMwmch3QFgRYr9tzHzMvd1HwAQ0WIAbwWwxD3nO0QUSnGunDLzdwD8HMAPiKgNwFYAm4noGiKan/lsQzYKzIcBVAL4FwBnAngHgHcNdxIR/RLAUwAWEtEe10T0P0S0zi2nfTGAf812okFFlRfFD1R5UfxAlZcxSIEVGGZ+AtmnM18J4E5m7mXmHXAUkaw8K8z8W2ZeDmAinNX+HwHwCgDfy/K9s4qBmcXMz8HxWb0HAIjoGgDPDDO5t6XY/cNsJyYoM5WlOSpjRMiKc/GmSketStJ2eQAAqI0Zn2t9VPpVJ1cYP/+UqPT5N5SbgGLv8v3euJd0yssXX3q9OO6iahN38k+brhZ9F08xVaZ3nZCxLPt76pLtUJm89soac319+6TPPNFpqmHHq+V5sUNGeY4dFV3oWBy3+mSsSG+DuTMqNsv367HCc2p3yffrbjI6dCIsx4x2mDH7K6Wu3fDSQXNevVXd2xtrkStu7IGIRwAAK+5lSBprhfn+qFpWGo9PMDEJPVNk2YaOWSbupX2JlM9zT9mWbK9okHEu3riXdMrLz7acJ4773aE3JdvPb5T/fBpeNNdUu9OTxmrHvRyWz7WMqdKlksaaC2TFUHljqez4q7CMibJLAlC1lKN4ozGIn2iW8tcxx3yvHYuljC0/xaRHv6n+BdGXKu4lm1Tp9SdM6ZOfrzpX9NW/YGJwGl6WSztE9lkyZsVVAQD3mGeZypiE4GsW0k1E9C4AqwB8jJnbAEyDswjdIHvcfVnj1lRc675+ksu52VhgPpnlPsVCLS+KH6jlRfEDtbyMYXK3wDQNZsW5rxuyeJfbAcwFsAxOYedUy6r4TloLDBFdDuAKANOI6BtWVy0ArUKdAVVeFD9Q5UXxA1VexjD5lQc4wsxnDX+Y9TbMSZM3EX0fwL3u5l4Adhne6e4+X8jkQtoH4HkAb3T/DtIJv2NXiMFuijRFpdvGriTtrSpdabmQKsJyFdvGmHEF2avrAsDUqEnbqwzJlOCaMmP6fN/CJ4dMNZ3y8p1Ny5PH9HK9OGdbn0lhPW9ii+ibHjEm+5ZumSrd1m9cNQNxaUzrbjMrcoZnyeDugW5jnrZTzQGgp9z0xSuliTu232yHO+VdU3nQbPd5Qrz76o1pvG2+nGe03bQT5dJtE+m0zcNyzP7m+mQ71GPr0/mnVIMoad6301YB6TYastppnXET9TfJi++abkz/x2bJaz++wJjDz1zQIvqummhuuVyqStsp+Y+1XyrOebHFPGfqXpIuCnv109je9FWlNY11ZBCVJdOgven4iJjvRFSRBsA1xm3UN0G6iY5bMtYxW8rYiUXmeXXhgq2i721NxvJ/yWy52u4gqZSXR3aYcjU7eueK43/6gnFbNqySMpYxVdpOxy/RqtJ548NHQETNzLzf3bQTcO4B8AsiuhXAVADzATxb/Bk5pFVgmHkNgDVE9P9cH5UyDGp5UfxALS+KH6jlJSAUWIFxE3CWw3E17YGzDtxyIlrmvlsLgPcDADOvJ6K7AGyA45m5kZnjKYYtCplcSOvgfjSUYnEvZj5tyM4SRpUXxQ9UeVH8QJWX4FDoIN5cE3CY+UtwKkn7TiYX0usz9CkWqrwofqDKi+IHqrwEjBL2omVyIe0cbBPRZABnu5vPMvOhYk9MQCb2xRvnEoka71bUE+diV2GuCctYlikVpkTA5Ij0+dslAToGZMrizYvvF9vplJfvbr5IHPfgkSXJ9tI6GePUGjdxE7/fuFT0nTu7Jdne3iFjYFqPmblFPLEskRrjNx7YL68h3GNVhz4hlxcnK8s57inbYN8oA1XSKmdvR1vlefXbzVx6GqTIVRwx39lAhYwJsNOqK/fJ2IvQnsPJdt88a33EfItcwEmdLhtMiY540ljt5drrZAxCT7P5/jpnyPOOzTPtyAJZCuINM19OtpdW7RZ9b523SmynU15e2DlTHHfLjmuS7e0Hm0Rf1Ysmdqduu/QKx3a2m41WOU9NlS4goTKU1boxUx4Z45iJaYvXyjirnokmzqVrmryHOiwZq14k09zffJKpXH565U7Rt2LOhiHTS6W8/Pf6K8Qxa46ZZ9tzO08SfbVrzTXYMS8AEN1ryVWbR8a0qnR+ZL843bgkm2KO18IJyrkGwLUAniGiqzOfVRqo5UXxA7W8KH6glpdgUuhSAkEim4XsPg3g7EGrCxFNBPAXAHcXc2JjHVVeFD9Q5UXxA1VeAsw4U0pyIRsFpszjMjqK7Cw3P4ITR3OImU9x9zUC+BWAWXAima91V/MbZjBGyF1l13YZAUBFJH2q9AQrVbouLFPzplr5u9Uh2TcrciTZvurkF4dMJ53y8mSLSSk80n+yOKfMUn13nWgUfWfUGneW7TICgD1d9cn20XZPfvIeY2bumyFTpfs7LddQlacS90Fjuj7R7KnyHDdum6YXpJvomJUxmfBITrTNXN+AtH7jxARzcF+tJ1W6y4hSb70Uq/pNVnq7x3Q8MGOiNWe7sjHypywEcs37XCnTWAfqjQvpRLO8wI7Z5vo658vPc+4iUzH88snrRd9ZlWa10+WztsBLOuXFlrNfHn2VOGfrAfO51Pxdug4nrDdy/sjDci3KFZM+kGxrqnQRKS9HYpKzona8Urpv+2vNfdnTKG+wrunm3uiaI5+BJy/ak2y/ZYpcUfe8ih3J9qkz9yAVqZQX2zX5TJsMh1y3x6TqV/9Vylj9dmv1873STWS7JlXGCsd4s6rkQjYRA/cT0QNEdB0RXQfgTwDuy+K8OzC0INQtAB5m5vkAHna3A4VaXhQ/UMuL4gdqeRkHFL6YY2AY1gLDzJ8goqsAvNLdtZKZf5fFeU8Q0SzP7ivh5JcDTs2DxwDcnO1kRxtVXhQ/UOVF8QNVXsYB41ApyYVhFRgi+iiAXzHzbwvwfpOt1fwOAJic4X1vAHADAJQ31aU7zDdUeRmf2HIWC9UMc3TxUeVl/CFkLDz6zzJAlZfxAmFEa48HnmxiYGoAPEhErXDiV35t10XIF2ZmovTeO2ZeCWAlAFTMm8oRtyyAHfMCANVRk6pXE5Fpe5OiJoaitlz6XMsstbWfZfruVXNl3Es65eXXW88Ux7XHpybbz7bNEn1zq03a77auiaLvjwfMmoAHO+U/0aqolQ7d61l6vM74wssOyrgMCpnrix72VL092/iiaad8v9ghOyZFvl3lfjOmJ2wIoT7TZ6diA7IkQN0mT1mDOpN2GfUsmJiIGvEMdcnvtr/OxKhEd5swKhrIbRFIW87qKqdyYoLz3Q7UyxiY7kkmXqFzpvw8OxcZmVwyX8YZXD7ppWR7VuSw6PPGvaRTXn708gXiuAePvjrZfuaZhaKvfrP5DOu2yc/MG/dio6nSxcOWsZr66dw907nn+mqkB7+nwWyfmCS60Dfb3HBLT5LLMLxpsnlezYnIFS5Sxb2kUl686dDfOWhk7MUNUsGpX2fuS7sEBQBUtLSbjaMyvJG7TEwiD8g4HpWxEVDCH9ewMTDM/HlmXgLgRgDNAB4nor/k+X4HiagZcGorAPB3PZk8UMuL4gdqeVH8QC0v449STqPOZdmvQ3DcPkcBTBrm2HTcA+DdbvvdAP6Q5zi+oMqL4geqvCh+oMrLOEWDeNNDRB+Es4DdRAC/BvA+Zh66hOPQ81IVhPoygLuI6HoAO91xh6WMGFUxx1RZF5O+i4aocUk0RY+LvnmVxtOVYKmr2W4j7+q6QHrl5f7ti5PHvNi9SJyzqtWkHp47oUX0hcmYSA+fkCu5Hmw1ilFlpTT7N1WYa9rfJvVGnmClHvbL6wsdN9t9DdKnU/uI8cMn6kUXEsajg7inhCdb1aIrD0lXDWdQhTlkzuudLNMuK3YYM/PA3AbRF+ow33WiUq5aGjlimaMrrHTUsvyX4k1EynBiupOqbqd+A0DXNHMN3XOkG/PcxduS7eUNsqrvlLBx171hzroh75lOeblzq6l2/3i7LDv21Hqz9GrNHnm9dTuMTDz2YPZJfprG6g8DMULbAke2+urk59o3wXzulc1dou+SGUbGzquRVaWnhc09lCodH0itvKzbNT3Zv6lvujj+0c0Lku3aTfJesNPxI/s9lcszpUpbbiOVsQJSwh9dNjEwMwB8hJlX5zJwmoJQAPDqNPvHDGp5UfxALS+KH6jlZRwzDt1CuZBNGnX6yL9xiCovih+o8qL4gSovJYAqMAqgyoviD6q8KH6gyktpoBaYMU5ZGaMq4sbARKRf1Y578VaVnlhulujvjMu02A8uekxsp1Nern/uOnHcyVUzkm1v+vWHZj6cbP/kgEx9bes1sR9dPVHRF++0lhD3SOOGHrNUTtmATDMOtZhrGqiScS7lJ6y4E09f10zzHpV75ZidC4xvOtImry/cZadKy3lW7jXfy7G5MsYn1GulX/dK33eixlxD5Y520Qcrrbr8SKfo4pAV+3HMihfo9wTu5EA8QuiY5XwX3VPk9fXPMLFJp86WaaxXTFibbE8ql/O8Ys5LYjud8nLHy+eL46ZYFdGf2SNTXKu3Gnlp3CjjcR67P791Idn+PjUeoWgkKhidS5xnWVW9fJbNaTCxLMvqZfrzq2pMbFVlmYyTSxX3kkp5Wb9rmjjmK/vNQukvHpAxMLEt5r5s2CJlLNxqZJM6pLwnejJUlVYZKw4l/FEGQoEpNmp5UfxALS+KH6jlpbQoZQtM/ikb4wRVXhQ/UOVF8QNVXpRSIhAWmBAlUBdxTJMTPKnSjWGzbacqA8C0cmOSvWS+TG8F0isvtjm/MSJNq1utJTJPih0Vff+19Ypku29AfrSNFcbsGk9IvZF6jKtmoEK6baZPNtewb69MQbYXkbarSAOAbWUmj+upxhRBxoD09qD+JTO3sgGp2pf3WKsXV8prGJhjBqreJ03c5e1W6rtntV3qN98ZhzwrBvcZ0zUflyv4wjJHk+1OGsG62vEY0LHImU90iny/s5pNVenz6neIvkXR/abvpJYh46ZTXh5rmZ9sz5KFifGZLW9Otge2yNWSa/cZl+ATf/q3dJeTG5wY/phSx5bdPH/1Vsd6ccEiJw16VqV8fsyOmpWaJ5ZLd/j8sDl2/oz9SEUq5eW3205P9leSLGOwq7Mx2T6+V8rYpB1GHmIHPauYH25PthOe+1Kk43vczCpjRWAcru2SC4FQYIqBWl4UP1DLi+IHankpYVSBKS1UeVH8QJUXxQ9UeSldCKUdA1NyCowqL4ofqPKi+IEqL4paYHyGiFoAdAKIAxhg5rMyHR8ui2NqpbNE9dRou+iLlpk4CW+5gEtmy7iXdMrLLWveIo57S72Jd3i6c57oq7fSWx89vED0dfeaQIYBT5zL9CYz751H5ZL5XGHiQMqOymCIPbCOjUlJje4y7a5Zsm/ACpep2u2pemsVw450iC5E2804fTUyoITLzHa0w5MOHbFSnjtkDAyXW+/vXeo/YfziZa2eVOkKk27OUyaIvrIOGQuVhPKPS49W9GHBEid9dWmDTJVeEDuQbFd50li9cS/plJdbN14qjptSbr7buw/JyuZ7Npn0+brd8nt49qcfHe5SckfTWh0oQxDVCGRrkIby47h64ioAQGNIlgsIwdwLxxNyqYVUcS+plJdfb5VytLPPxOw9ckRWLt+5qynZblwrr622xcS9hI7IeBzuseS/X6ZYs72MgTfmRWWsKFAJf66jaYG5mJmP+PVmanlR/EAtL4ofqOVFAaBBvKM9AT9Q5UXxA1VeFD9Q5UWx0RgY/2EADxIRA/geM6/0HkBENwC4AQBqmyuS6dLeVOn+hLmEz5xy75A3Sqe8fHfzRcljqkMyhfB3HcYMu+XYRNFXUW5MprtapSuop8uYfaNVchXKJ7YbV9RAt/zYaycZU3LvGjlm4oQZM9Qjzdu2m6hhg+w7blVPrjwgTbn9VaavylNVuqfepDKHuz0uq3YzTsUBT1qzmJh8v1CbcffEJ8pUTjrYmmxzo6ev10qjjsu5xPebSuN2SjV7TNrDYctZQ3MMr5vsrJw7OdwujrPdk2+f/+yQcdIpL8/tNKvo1odkSv4Tx4xJ/4W1c0Vf3Rbzfmu+9a/ZXo6SDTm4ichym4q+RPb5+raMTZ5ajomhYymP60xUJNveFZwHSaW8JA6YdPwwyaUW1nWZ1XfXb5Ur8dZsMO7q2l3yvik/Yp5JfKJH9IlU6QFvyXrr3i9h14avlPDHPFoL2b2Smc8AcDmAG4noVd4DmHklM5/FzGdVNkSHjpAFanlRhsOWs+rG8PAnpEAtL0ombBmrn5D/I1ctL0oqiHN7jSdGRYFh5r3u30MAfgfgnEK/hyovih+o8qL4gSovSlo4x9c4wncFhoiqiKhmsA3gMgCp7aV5osqL4geqvCh+oMqLkpYcrS/ZWGCI6EdEdIiIXrL2NRLRQ0S0xf3b4O4nIvoGEW0lorVEdEbxLnYooxEDMxnA78jxQ5cD+AUz35/phHJKoLHciaNo8iyxfd2Cp8R2OuXlaxteK47r5yqrLZewb+k2Kbs7DjSJvikTTN5xPO7xl4eM/zcclr7hnk7jBgsflq6KgW0m7iVRIyUscsz42hNymrCn3VcrffLlVpZxqM+bYm2ODXfKGJjybnMN5C0l0G35yT0xKRw1k4nXy8rf5SGr5EGPx2fu9aHbY7aZz7qsV8YUJcqN6HKvndac/0+MKPVjTtSJrQlDfi4r5mwQ2+mUl4+tvlYc15swVcmP9FaLvqdeMnFR1Tvll7vmGxr3UlAyxL2QXcLCGwMTth6R1hje8hz5EPLIaqq4l1TKyyM7ZDp0zIoL7GH5bHlym5Gxih1yiYbqveZejx6WMW3UbVWVHhIDY8WmecsFKP5T+K/gDgDfAvBTa98tAB5m5i8T0S3u9s1wwkDmu69zAdzu/vUF3xUYZt4OYGkxxlbLi+IHanlR/EAtL8pwFGMlXmZ+gohmeXZfCWC52/4JgMfgKDBXAvgpMzOAp4monoiamTl1wa4CM27SqFV5UfxAlRfFD1R5UbIm92yvJiJaZW2vTJUJ7GGypZQcgONJAYBpAHZbx+1x96kCM0icCR1xJ8Xw44sfGNKfTnmxV6X8aKN0Pb1rp0l8qgjJFMLF1eaz3ztJpvbOqDbVocs8qu+hDuMi6Dwo3QWhauMq4ZOk22Rgu0mfJM/ilX115j3K+jyVnK1jI8fkXHobzLHxaHpzdzwmTeZl/WacUK90o/TV2yndss9eDTKyU1bZhbVyZ//sKaIr3GV9Tp2eyrYRYw7nOvl5llnmfLtSNXV4/Gw5MMAhHB1w3sfrmgTSKy/f2bQ82a4OSXlp7TOuyqc2yFWdq7eY61v/P+oy8gvyVD2HtZ2xz5JH9OTnQhrgEA7HnefTVXNfTHlMKuXlj9tPTfafHZUrVv+9x6za/K0dl4g+PmTu2Vird1kE8xwq6/S4iU6YlXjttGlnRyJ1G9DU6VEgDwvMkeFWv88EM7O7BMqoM1pp1AVDLS+KH6jlRfEDtbwoOZFrBlL+asdBImoGAPfvIXf/XgAzrOOmu/t8IdAKjCovih+o8qL4gSovSj5QIrdXntwD4N1u+90A/mDtf5ebjXQegA6/4l+AgLiQUqHKi+IHqrwofqDKi5I3BXbmENEv4QTsNhHRHgCfBfBlAHcR0fUAdgIYTLe8D8AVALYC6AbwnsLOJjOBUGCmVszDF079fXI7nfLy3+uvEOe9p94ogr/pahZ91eXGr7vm6FTR98KAsYgd3SGX9j88zcRidO+VcRmCShkjEu8yH3X0oPzYy3qNP71MnobeBiu2xLMCebTN9JWf8KQ1W7a1WJuncnSn5b/3CH9/lfH5l3fK2KDYfrO8eKJCpmuWH/SUtbaITzOp6OFdh+U8rTRq8lSqpiprWfTjJ0Sf7aNHyDpvBNmtE2KLROxLOuXlyRa57P/a469ItmdXyOt7bKu1zPtR+b1HO8aEG3l84k2bttOjvanSdpxLueeRGLVWAbfLCuQpaPXRxbhq7iqxL5Xy8s1NMpZllpUBva5P/lj77RGz9EZfXMbwRFvNtdrxbQBQfty69zzLFNj3JXviWjR1emxRhCykt6XpenWKYxnAqP2yC4QCY6OWF8UP1PKi+IFaXpQRwSjpwOlAKTCqvCh+oMqL4geqvCiFYGzkA40Oo6LAENEKAP8HIATgB8z85eHOSae8XPW3f04ec3a9PGdlm0mjjrM0G2/umJRs9/TLj6Gr21pJtkamPIvVdz0h0OF2s4O7ZGfYWlH3RLN06YQ7zLGhDmmartxvtss91aEjnSYiq8yzoG1vvZVm7DGnh3rNODSQPqqrr0Gu3FnRbczMA5XShVRWZaWCe8zRoX2etGqLxFFTjTo0U1Zrju/ak2x7q96WVRr3UsJK0+a4xweXI+mUly++9PpkuyMuV8t+S+NzyfZXd64QfbTfyFLFfs93ezj9KsSKj9j3R1g+C6jccsmIFXtHthJvKuVly27j5l4SlT/OjsaNu/rZE9KF2RM39+KhI/K8GuP1RdUBeW+Ejlsr6vb0ij7ut2TT6zLSitNjixL+CkajFlIIwLfhLEG8GMDbiGjxcOep5UUpNn19fWp5UYqOWl6UQjG4Em+pVqMeDQvMOQC2uiUFQER3wlmOeEO6E3bt2oVIJKLKi1JUNm/ejFtvvVWVF6VoHDp0SJUXpXAwl7QVbDQUmFRLD2cs/tTd3a3Ki1J0Jk+erMqLUlQOHjyITZs2qfKiFIzxZlXJhTEbxEtENwC4wd3sraurG1qq1cPvhh+2CcCREU2scOOMpbmMjXG25TDO8bQ9J+Xyll45s0vIZ8vXxNbmvK6f6N+8u0b/+wjyOJke6iYUxBnDDtnqGnJ0OrKWM6+MzZmTouz0EA5k6HsBCNr3Eby5jIVxspMxVWB8Jaulh93iUisBgIhWjaR2wyBjaZyxNJfxPM5wFFrOxtr16zjFn8twjPdnWaHGGUtzGYvjpB1fFRhfeQ7AfCKaDUdxeSuAt4/CPBRFURQluDCGZomVEL4rMMw8QEQ3AXgAThr1j5h5vd/zUBRFUZTAU7r6y+jEwDDzfXBqKGTLygK99VgaZyzNZTyP4/d7jrXr13GKO8Zoved4HGcszWUsjpOSUnYhkbfOhaIoiqIoY5+auul85vkfyumcxx+45Xk/Yr/8YMxmISmKoiiKkplStsD4vhJvLhDRCiLaTERbieiWEYzTQkTriGg1Ea0a/ozkeT8iokN2ai0RNRLRQ0S0xf3bkGmMDON8joj2unNaTURXZBrDPWcGET1KRBuIaD0RfTjXOWUYI6f5EFGMiJ4lojXuOJ93988momfc7+xXRBTJc5w7iGiHNZ9lw30++aJyJsYYsYwNM47KWYDlrBAy5p4z7uRsVGSM83iNI8asAkN5lhzIwMXMvCxH09kdAFZ49t0C4GFmng/gYXc7n3EA4DZ3TsvcuKDhGADwMWZeDOA8ADe6n0kuc0o3Rq7z6QVwCTMvBbAMwAoiOg/AV9xx5gFoA3B9nuMAwCes+aweZpy8UDkbQiFkLNM4uc5H5Sw1oyVnqcYARudZlmmcXOdUCDnzXcacUgKc02s8MWYVGFglB5i5D8BgyQHfYOYnALR6dl8J4Cdu+ycA3pTnOPnMZz8zv+C2OwFshLOycdZzyjBGrnNhZh5c9ivsvhjAJQDuzmYuw4zjFypncowRy9gw4+Q6H5WzAlEIORtLz7Jhxsl1PiOWs1GTsUSOr3HEWFZgUpUcyFkwXRjAg0T0PDmrYo6Eycy8320fADB5BGPdRERrXbPssC4CGyKaBeB0AM/kOyfPGDnPh4hCRLQawCEAD8FZS7edmQdL2Wb1nXnHYebB+XzJnc9tRBTN5pryQOUsDYWQsRTj5DwflbMhjEU5G9VnWYpxcp5TIeRsNGRMLTDjn1cy8xlwzLc3EtGrCjEoM4/Eq3g7gLlwTI37AXw92xOJqBrAbwB8hJmP5TOnFGPkPB9mjjPzMjirKZ8DYFG215BpHCI6BcAn3fHOBtAI4OZ8xvaZcSNnhZCxNOOonI2csSZno/osSzPOqMiZ7zKmMTBjlqxKDmQDM+91/x6CUzLpnBHM6yARNQOA+/dQnnM66Ap7AsD3s50TEYXh3Kg/Z+bf5jOnVGPkOx/33HYAjwI4H0A9EQ1mt+X0nVnjrHBNw8zMvQB+nMt8ckTlzEMhZCzdOCpnScaNnI3msyzdOKMtZ/7JGJuK1Nm+xhFjWYFJlhxwI7/fCuCeXAchoioiqhlsA7gMQM4F+yzuAfBut/1uAH/IZ5DBm9TlzdnMiYgIwA8BbGTmW/OZU7oxcp0PEU0konq3XQHgUjj+50cBXJ3NXDKMs8l6iBEcv/NIvrNMqJzJ40csY5nGUTkbf3I2Ws+yTOOMhpyNlowR5/YaT4zZdWC4cCUHJgP4nSM7KAfwC2a+P5sTieiXAJYDaCKiPQA+C+DLAO4iousB7ARwbZ7jLCcnnY4BtAB4fxZTugDAOwGsc/2sAPCpHOeUboy35TifZgA/ISe7ogzAXcx8LxFtAHAnEX0RwItwHi75jPMIEU2EE2i/GsAHhhknL1TOhlAIGcs0jspZgOVsjD3LMo0zGnI2OjI2zqwquaAr8SqKoihKAKmtnsbnLvtgTuf85W+f0ZV4FUVRFEUZZUrYCKEKjKIoiqIEldLVX1SBURRFUZSgMt7WdskFVWAURVEUJaioAqMoiqIoSqBgjLvyALkwlteBUSyI6CNEVJnFccuJ6F63/UbKs+qtPY5SOqicKcVGZaxwEHIrI5Ctu4lSVDynHCuF+4EqMEWGzCqOKbdz4CMAhr3pbZj5Hmb+cp7vpwQIlTOl2KiMlRzeiue5VgovOqrAZAkRvYucYlxriOhn7r47iOhq65gu9+9yInqSiO4BsCHFdoiIvkpEz7ljvt867zEiupuINhHRz8nhXwBMBfAoET2aYm4r3ONfAHCVtf86IvqWNdfvEtEqInqZiF7v7o8R0Y9dbftFIro4xfhV5BREe9Y95kp3/xJ332r3OuYX7AMvUVTOVM6KjcrYOJMx/0oJ5FQp3A80BiYLiGgJgM8AeAUzHyGixixOOwPAKcy8g4iWe7ZvANDBzGeTU5n0b0T0oHve6QCWANgH4G8ALmDmbxDRR+FoxEc8c4vBqfVxCYCtAH6VYU6z4NThmAvnATIPwI1w6qadSkSL4FS5XeA579MAHmHm95KzVPazRPQXOCtK/h8z/5yc5dFDWXwuShpUzlTOio3K2DiUsdyVkiZy3UIuK5l5pXdUOJ8fA/ie21+oyuUFQxWY7LgEwK8Hbzhmbs3inGeZeUea7csAnGb94qkDMB9An3vcHgAgZ2nsWQD+muF9FgHYwcxb3HP+H4Ab0hx7l1vcbAsRbXfPfSWAb7rXtYmIdgLw3vSXAXgjEX3c3Y4BmAngKQCfJqLpAH47OAclb1TOVM6KjcrYeJKx/IJ4j2SxEu8rmXkvEU0C8BARbRJvy8yucjOqqAIzMgbguuGIqAxAxOo77jnW3iYAH2LmB+wD3F83vdauOAr7HXkFLlsBJABvYebNnv0biegZAK8DcB8RvZ+ZHxnpJJUhqJypnBUblbGAylgx1oGxK54T0WDF84NE1MzM+ynPyuWFRmNgsuMRANcQ0QTAicZ297cAONNtvxFAOMvxHgDwz+SUgQcRLSCnsmwmOgHUpNi/CcAsIprrbr8twxjXEFGZe+wcAJsBPAngHwfnAefXiPfmfgDAh4icCnJEdLr7dw6A7cz8DThVWk8b5hqUzKicqZwVG5Wx8SZjBY6BofQVz0dcubzQqAUmC5h5PRF9CcDjRBSHU5X0Ojj+2j8Q0RoA92PoL5V0/ACOOfUF90Y6jOEDolYCuJ+I9jFzMjiNmXtcP/SfiKgbzk2c6uEAALsAPAugFsAH3HO/A+B2IloH51fYdczc697fg3wBwP8CWOv+OtsB4PVwKsW+k4j64fhE/yvL61dSoHKmclZsVMbGm4yNODA3FSkrnhPRc8itUnjR0WrUJQIR3QHgXma+e7TnooxfVM6UYqMyZqiraObz5703p3MeeOm/tBq1oiiKoiijTAmvxKsKTInAzNeN9hyU8Y/KmVJsVMYkWsxRURRFUZTgoQqMoiiKoiiBggEkVIFRFEVRFCVQFCULKTCoAqMoiqIoQUUVGEVRFEVRAocqMIqiKIqiBAqNgVEURVEUJXgwwKW7EIwqMIqiKIoSVNSFpCiKoihKoFAXkqIoiqIogUQtMIqiKIqiBA5VYBRFURRFCRa6kJ2iKIqiKEGDASQ0C0lRFEVRlKChFhhFURRFUQKHKjCKoiiKogQL1jRqRVEURVECBgOsK/EqiqIoihI41AKjKIqiKErg0BgYRVEURVECBbOmUSuKoiiKEkDUAqMoiqIoStBgtcAoiqIoihIstJSAoiiKoihBg6FZSIqiKIqiBBBdB0ZRFEVRlCDBAFgtMIqiKIqiBApmtcAoiqIoihI81AKjKIqiKErwKGELDHEJp2ApiqIoSlAhovsBNOV42hFmXlGM+fiNKjCKoiiKogSOstGegKIoiqIoSq6oAqMoiqIoSuBQBUZRFEVRlMChCoyiKIqiKIFDFRhFURRFUQLH/wfOc2w6sBcTewAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 566.64x282.96 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import reda.plotters.pseudoplots as PS\n",
    "fig = PS.plot_ps_extra(\n",
    "    obj, 'rho_a',\n",
    "    subquery='timestep in [0, 3, 6]'\n",
    ")\n",
    "fig.set_size_inches(7.87, 3.93)\n",
    "fig.tight_layout()\n",
    "fig.savefig('pseudoplots.png', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAADDCAYAAACS9yhvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXNUlEQVR4nO3dfZRcdX3H8ffHZSNRLCnNNiduEpfaIKQgCR3TrLR1S6oC9jS2xZY+8FQ0JzVa4sEq4kGKD8VWD9ZUeYhihTYtJZBqSmkVY9ZWuwluQkjYrMhW0SREWE14sCEJ2Xz7x/0lTMImzCQzd+bOfl7n7Jk7v/v03Tm785l772/uTxGBmZmZFc9LGl2AmZmZHR2HuJmZWUE5xM3MzArKIW5mZlZQDnEzM7OCcoibmZkV1HGNLqCWJk6cGF1dXY0uw8zMrKbWrl3744joOLS9pUK8q6uL/v7+RpdhZmZWU5J+MFq7T6ebmZkVVEsdiZsZPP7442zfvr3RZdREe3s7EydOZMKECY0uxawpOcTNWsz27ds55ZRTaGtra3QpxyQi2LVrF48++qhD3OwwcjudLul4SfdLelDSgKTrUvtSSQ9LekjSFyS1p3ZJWixpSNIGSWflVatZ0R0pwPv64Prrs8dmJonx48c3ugyzquX5P5bnNfHdwDkRcSYwEzhX0hxgKXAqcAYwHnh7Wv48YHr6mQ/clGOtZi2prw/mzoVrrskej+VNpre3ly996UusWLGCZ5999kD7okWLRl1+tPbDLWtWVLX8H6tEbiEemZ+mp+3pJyLi3jQvgPuBKWmZecDtadZqYIKkyXnVa9aKenthzx4YGckee3uPvPxVV13F7t27ueWWWxgYGOBzn/scH/3oR1m+fPmBZdatW8fu3bu55ppruPnmm9m4ceOBeTt37mThwoXcfPPN7Nixgx07drB48WIWLVrEI488wqZNm/jiF7/I8PDwgfahoaH6/PJmOaj2f+xY5do7XVKbpPXAE8B9EbGmbF47cBHwn6mpE9hctvqW1GZmR6mnB8aNg7a27LGn58jLX3DBBSxbtozBwUFOPfVUdu3axZQpU/jGN77xgmWHh4dZsGAB06dPP9C2ceNGSqUSCxYsYPz48QeO2CdOnMjq1auZMWMGl156Kc8999yB9r5mP89vdgTV/o8dq1xDPCJGImIm2dH2bEmnl82+EfiviPjvarYpab6kfkn9w8PDNazWrPV0d8PKlfCRj2SP3d1HXr5UKnHXXXcxZ84cnn32WTZv3kxbWxt79+59wbKTJ09m2bJlBx1Jn3HGGQwMDHDHHXewc+dOHn30UUZGRti3bx8jIyNMmzaNz3zmMwwNDR3UblZU1f6PHStlZ7HzJ+lDwM6I+KSka4FZwO9GxL40/xagNyL+OT1/GOiJiG2H22apVArf7MXGusHBQU477bRGl1Ezrfb7mB0NSWsjonRoe5690zskTUjT44E3At+R9HbgzcAf7g/wZAVwceqlPgd46kgBbmZmNtbk+T3xycBtktrIPjzcGRH3SNoL/ADokwSwPCI+DNwLnA8MATuBy3Ks1azQRkZGWuZ74ja29fVlncN6eup/arqIcgvxiNhAdsr80PZRa0i91RfWu67D8R+OFdVJJ53Ed7/73UaXURPt7e1MnuwvpYxV+7+utWdP1kksj2vMReM7to3CfzhWZJMmTWLSpEmNLsPsmI32dS2/Fx/MA6CMIu/v+ZmZ2Qvl/XWtIvKR+Cj2/+HsPxL3H46ZWf72f13LlzYPzyE+Cv/hmJk1h+5uvwcfiUP8MPyHY2Zmzc7XxM3MzArKIW5mZlZQDnEzM7OCcoibmZkVlEPczMysoBziZmZmBeUQNzMzK6g8hyI9XtL9kh6UNCDputR+sqQ1koYk/Yukcan9pen5UJrflVetZmZmRZDnkfhu4JyIOBOYCZybxgn/a+BTEfGLwA7g8rT85cCO1P6ptJyZmZkluYV4ZH6anrannwDOAe5K7bcBb03T89Jz0vy5SgOOm5mZWc7XxCW1SVoPPAHcB/wv8GRE7E2LbAE603QnsBkgzX8K+Lk86zUzM2tmuYZ4RIxExExgCjAbOPVYtylpvqR+Sf3Dw8PHujkzM7PCaEjv9Ih4ElgFdAMTJO0fiGUKsDVNbwWmAqT5JwI/GWVbSyKiFBGljo6OepduZmbWNPLsnd4haUKaHg+8ERgkC/ML0mKXAF9O0yvSc9L8r0dE5FWvmZlZs8tzKNLJwG2S2sg+PNwZEfdI2gTcIemjwAPArWn5W4F/kDQEbAcuzLFWMzOzppdbiEfEBmDWKO3fI7s+fmj7LuBtOZRmZmZWSL5jm5mZWUE5xM3MzArKIW5mZlZQDnEzM7OCcoibmZkVlEPczMysoBziZmZmBeUQNzMzKyiHuJmZWUE5xM3MzArKIW5mZlZQDnEzM7OCcoibmZkVVJ7jiU+VtErSJkkDkq5I7TMlrZa0XlK/pNmpXZIWSxqStEHSWXnVamZmVgR5jie+F7gyItZJegWwVtJ9wN8A10XEf0g6Pz3vAc4DpqefXwFuSo9mZmZGjkfiEbEtItal6WeAQaATCOBn0mInAo+l6XnA7ZFZDUyQNDmves3MzJpdnkfiB0jqAmYBa4BFwFckfZLsQ8Xr02KdwOay1baktm25FWpmZtbEcu/YJukE4G5gUUQ8DfwZ8J6ImAq8B7i1yu3NT9fS+4eHh2tfsJmZWZPKNcQltZMF+NKIWJ6aLwH2Ty8DZqfprcDUstWnpLaDRMSSiChFRKmjo6M+hZuZmTWhPHuni+woezAibiib9RjwhjR9DvBIml4BXJx6qc8BnooIn0o3MzNL8rwmfjZwEbBR0vrUdjXwDuDTko4DdgHz07x7gfOBIWAncFmOtZqZmTW9Fw1xSSdVsJ19EfHkkRaIiG8COszsXx5l+QAWVrBvMzOzMamSI/HH0s/hAhigDZhWk4rMzMysIpWE+GBEzDrSApIeqFE9ZmZmVqFKOrZ112gZMzMzq6EXDfGI2FWLZczMzKy2Ku6dLqkEfBB4VVpPZP3PXlun2szMzOwIqvmK2VLgL4CNwL76lGNmZmaVqibEhyNiRd0qMTMzs6pUE+LXSvo8sBLYvb+x7PapZmZmlqNqQvwy4FSgnedPpwfP3/fczMzMclRNiL8uIl5Tt0rMzMysKtUMgPI/kmbUrRIzMzOrSjVH4nOAByV9j+yauL9iZmZm1kDVhPi5wDjgObJr4VWRNBW4HZiU1l8SEZ9O895NNtjJCPDvEfG+1P4B4PLU/ucR8ZVq92tmZtaqKgpxSe8A3gu8AngZsI5sGNEbgbMq3Nde4MqIWCfpFcBaSfeRhfo84MyI2C3p59M+ZwAXAr8EvBL4mqRTImKk4t/OzMyshVUyFOn7gRLwhoj4UWp7E/A5qhi5LCK2AdvS9DOSBoFOsvHEPx4Ru9O8J9Iq84A7Uvv3JQ0Bs4G+SvdpZmbWyirp2HYZ8Ef7AxwgIr4K/CbwtaPZqaQuYBawBjgF+DVJayR9Q9Lr0mKdwOay1bakNjMzM6PC0+kR8dwobY9L+rtqdyjpBOBuYFFEPC3pOOAkso5zrwPulPQLVWxvPjAfYNo0D2luZmZjRyVH4v8r6S2HNkr6MNnd2yomqZ0swJeW3eltC7A8MveT3UhmIrAVmFq2+pTUdpCIWBIRpYgodXR0VFOOmZlZoVVyJP5O4G5JlwIPAicA5wPrgYcr3ZEkAbcCgxFxQ9msLwG/AaySdApZD/gfAyuAf5J0A1nHtunA/ZXuz8zMrNW9aIhHxA/Sdeo3A6cBTwN/EhEbUqe3Sp0NXARslLQ+tV0NfAH4gqSHgD3AJRERwICkO4FNZD3bF7pnupmZ2fOU5WVrKJVK0d/f3+gyzMzMakrS2ogoHdr+otfEJa2rxTJmZmZWW5VcEz9N0oYjzBdwYo3qMTMzswpVEuKnVrCMr1WbmZnlrKKObXkUYmZmZtWpZihSMzMzayIOcTMzs4JyiJuZmRVUNeOJHyDpGuD1ZLdBfSAiPlvTqszMzOxFHe2R+M8Bq4GPAa+pXTlmZmZWqaMN8R1AG/AEsL125ZiZmVmlqjqdLuklwEsi4jpJrwQWAw/VpTIzMzM7oopDXNK7gGuBPZJ+AiyOiMvrVpmZmZkdUTWn068EzoiITrIRzc6WdG19yjIzM7MXU02I/5TsGjgRsQ24HPi9SleWNFXSKkmbJA1IuuKQ+VdKCkkT03NJWixpSNIGSWdVUauZmVnLqybEbwKWSfrF9HwasLOK9fcCV0bEDGAOsFDSDMgCHngT8MOy5c8Dpqef+Wn/ZmZmllQc4hFxI7AU+LykHcAQ8LCkt0maXsH62yJiXZp+BhgEOtPsTwHvA8oHN58H3B6Z1cAESZMrrdfMzKzVVfUVs4hYHhE9QAdwFvB1spu+3FLNdiR1AbOANZLmAVsj4sFDFusENpc938LzoW9mZjbmHdUd2yJiL7Ah/dxWzbqSTgDuBhaRnWK/muxU+lGRNJ/sdDvTpk072s20hL4+6O2Fnh7o7m50NWZmVm9HFeJHS1I7WYAvjYjlks4ATgYelAQwBVgnaTbZLV2nlq0+JbUdJCKWAEsASqVSHDp/rOjrg7lzYc8eGDcOVq50kJuZtbrcBkBRltK3AoMRcQNARGyMiJ+PiK6I6CI7ZX5WRPwIWAFcnHqpzwGeSr3ibRS9vVmAj4xkj729ja7IzMzqLc8j8bOBi4CNktantqsj4t7DLH8vcD5ZB7qdwGV1r7DAenqyI/D9R+I9PY2uyMzM6i23EI+IbwJ6kWW6yqYDWFjnslpGd3d2Ct3XxM3Mxo5cr4lbfXV3O7zNzMaS3K6Jm5mZWW05xM3MzArKIW5mZlZQDnEzM7OCcoibmZkVlEPczMysoBziZmZmBeUQNzMzKyiHuJmZWUE5xM3MzArKIW5mZlZQeQ5FOlXSKkmbJA1IuiK1f0LSdyRtkPSvkiaUrfMBSUOSHpb05rxqNTMzK4I8j8T3AldGxAxgDrBQ0gzgPuD0iHgt8F3gAwBp3oXALwHnAjdKasuxXjMzs6aWW4hHxLaIWJemnwEGgc6I+GpE7E2LrQampOl5wB0RsTsivk82rvjsvOo1MzNrdg25Ji6pC5gFrDlk1p8C/5GmO4HNZfO2pDYzMzOjASEu6QTgbmBRRDxd1v5BslPuS6vc3nxJ/ZL6h4eHa1usmZlZE8s1xCW1kwX40ohYXtZ+KfBbwB9HRKTmrcDUstWnpLaDRMSSiChFRKmjo6NutZuZmTWbPHunC7gVGIyIG8razwXeB/x2ROwsW2UFcKGkl0o6GZgO3J9XvWZmZs3uuBz3dTZwEbBR0vrUdjWwGHgpcF+W86yOiAURMSDpTmAT2Wn2hRExkmO9ZmZmTS23EI+IbwIaZda9R1jnY8DH6laUmZlZgfmObWZmZgXlEDczG0P6+uD667NHK748r4mbmVkD9fXB3LmwZw+MGwcrV0J3d6OrsmPhI3EzszGitzcL8JGR7LG3t9EV2bFyiJuZjRE9PdkReFtb9tjT0+iK7Fj5dLqZ2RjR3Z2dQu/tzQLcp9KLzyFuZjaGdHc7vFuJT6ebmZkVlEPczMysoBziZmZmBeUQNzMzKyiHuJmZWUE5xM3MzAoqz/HEp0paJWmTpAFJV6T2kyTdJ+mR9PizqV2SFksakrRB0ll51WpmZlYEeR6J7wWujIgZwBxgoaQZwFXAyoiYDqxMzwHOA6ann/nATTnWamZm1vRyC/GI2BYR69L0M8Ag0AnMA25Li90GvDVNzwNuj8xqYIKkyXnVa2Zm1uwack1cUhcwC1gDTIqIbWnWj4BJaboT2Fy22pbUdui25kvql9Q/PDxcv6LNzMyaTO4hLukE4G5gUUQ8XT4vIgKIarYXEUsiohQRpY6OjhpWamZm1txyDXFJ7WQBvjQilqfmx/efJk+PT6T2rcDUstWnpDYzMzMj397pAm4FBiPihrJZK4BL0vQlwJfL2i9OvdTnAE+VnXa3FtDXB9dfnz0WQdHqNbPWl+coZmcDFwEbJa1PbVcDHwfulHQ58APg99O8e4HzgSFgJ3BZjrVanfX1wdy5sGdPNq7xypXNPbJS0erdr6+vWMNOul6z6uQW4hHxTUCHmT13lOUDWFjXoqxhenuzQBwZyR57e5v7TbBo9ULxPni4XrPq+Y5t1hA9PdkbX1tb9tjT0+iKjqxo9cLoHzyames1q16ep9PNDujuzo5cinIqsmj1wvMfPPYfKTb7Bw/Xa1Y9ZWetW0OpVIr+/v5Gl2HWNIp2zdb1mo1O0tqIKL2g3SFuZmbW3A4X4r4mbmZmVlAOcTMzs4JyiJuZmRVUS10TlzRMdsOYWpkI/LiG27MX8mtcf36N68uvb/35NYZXRcQLBghpqRCvNUn9o3UksNrxa1x/fo3ry69v/fk1PjyfTjczMysoh7iZmVlBOcSPbEmjCxgD/BrXn1/j+vLrW39+jQ/D18TNzMwKykfiZmZmBeUQPwxJ50p6WNKQpKsaXU8rkTRV0ipJmyQNSLqi0TW1Kkltkh6QdE+ja2lFkiZIukvSdyQNSvId1GtM0nvS+8RDkv5Z0vGNrqmZOMRHIakN+CxwHjAD+ENJMxpbVUvZC1wZETOAOcBCv751cwUw2OgiWtingf+MiFOBM/FrXVOSOoE/B0oRcTrQBlzY2Kqai0N8dLOBoYj4XkTsAe4A5jW4ppYREdsiYl2afobsja+zsVW1HklTgLcAn290La1I0onArwO3AkTEnoh4sqFFtabjgPGSjgNeBjzW4HqaikN8dJ3A5rLnW3DI1IWkLmAWsKbBpbSivwXeB+xrcB2t6mRgGPj7dMni85Je3uiiWklEbAU+CfwQ2AY8FRFfbWxVzcUhbg0j6QTgbmBRRDzd6HpaiaTfAp6IiLWNrqWFHQecBdwUEbOA/wPcf6aGJP0s2VnQk4FXAi+X9CeNraq5OMRHtxWYWvZ8SmqzGpHUThbgSyNieaPraUFnA78t6VGyy0HnSPrHxpbUcrYAWyJi/1mku8hC3WrnN4HvR8RwRDwHLAde3+CamopDfHTfBqZLOlnSOLKOFCsaXFPLkCSy64iDEXFDo+tpRRHxgYiYEhFdZH+/X48IH8HUUET8CNgs6TWpaS6wqYEltaIfAnMkvSy9b8zFnQcPclyjC2hGEbFX0ruAr5D1hvxCRAw0uKxWcjZwEbBR0vrUdnVE3Nu4ksyOyruBpenD/veAyxpcT0uJiDWS7gLWkX2r5QF897aD+I5tZmZmBeXT6WZmZgXlEDczMysoh7iZmVlBOcTNzMwKyiFuZmZWUA5xsxaXRtp6Z5p+ZfrKTr32NVPS+fXavpkdzCFu1vomAO8EiIjHIuKCOu5rJuAQN8uJvydu1uIk7R+F72HgEeC0iDhd0qXAW4GXA9PJBpoYR3Yjnt3A+RGxXdKryYbm7QB2Au+IiO9IehtwLTACPEV2i8whYDzZbYqvB+4B/g44HWgH/jIivpz2/TvAiWSDC/1jRFxX31fCrPX4jm1mre8q4PSImJlGjbunbN7pZKPIHU8WwO+PiFmSPgVcTDYS2hJgQUQ8IulXgBuBc4APAW+OiK2SJkTEHkkfIhv7+V0Akv6K7JavfyppAnC/pK+lfc9O+98JfFvSv0dEfx1fB7OW4xA3G9tWpTHdn5H0FPBvqX0j8No00tzrgWXZrasBeGl6/BbwRUl3kg1MMZo3kQ3E8t70/HhgWpq+LyJ+AiBpOfCrgEPcrAoOcbOxbXfZ9L6y5/vI3h9eAjwZETMPXTEiFqQj87cAayX98ijbF/B7EfHwQY3Zeodey/O1PbMquWObWet7BnjF0ayYxnn/frr+jTJnpulXR8SaiPgQMEw2fO+h+/oK8O40AhWSZpXNe6OkkySNJ7s2/62jqdFsLHOIm7W4dMr6W5IeAj5xFJv4Y+BySQ8CA2Sd5AA+IWlj2u7/AA8Cq4AZktZL+gPgI2Qd2jZIGkjP97ufbEz5DcDdvh5uVj33Tjez3KXe6Qc6wJnZ0fGRuJmZWUH5SNzMzKygfCRuZmZWUA5xMzOzgnKIm5mZFZRD3MzMrKAc4mZmZgXlEDczMyuo/wfWUqucHpvWUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 566.929x198.425 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import reda.plotters.time_series as TS\n",
    "\n",
    "fig, ax = TS.plot_quadpole_evolution(\n",
    "    obj,\n",
    "    [10, 11, 15, 14],\n",
    "    'rho_a',\n",
    "    threshold=0.05,\n",
    "    rolling=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add some noise to the data\n",
    "import numpy as np\n",
    "np.random.seed(2017)\n",
    "obj.data['rho_a'] = obj.data['rho_a'] + (np.abs(obj.data['rho_a'] * 0.05 * np.random.randn(obj.data.shape[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAADECAYAAACP8hjXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXQUlEQVR4nO3df5RcZX3H8ffHJYEo1DTNnjQmWUMlCBEksVPMiuKeRCuCp9gWLbZSUWyaGlppEUU8imCV9mihKpU2igVtCvIjxZRiK8YsFUmCmxASkgVcIjbEELYk/GowPzbf/nGfhUnYTWaSmTtzZz+vc+bMvc/9Md+9J7vfPD/meRQRmJmZWfG8rNEBmJmZ2cFxEjczMysoJ3EzM7OCchI3MzMrKCdxMzOzgnISNzMzK6jckrikIyTdK+l+SeskXZbKr5P0M0mr02tGKpekr0jqk7RG0hvyitXMzKwIDsvxs3YAsyPiOUmjgLslfS8duygibtnn/HcC09LrjcA16X1Y48ePj6lTp9Y2ajMzswZbuXLl/0ZE+77luSXxyGaVeS7tjkqv/c00cybwrXTdckljJU2MiM3DXTB16lR6enpqFrOZmVkzkPTzocpz7ROX1CZpNfAEcGdErEiHPp+azK+SdHgqmwRsLLv8sVRmZmZm5NucTkQMADMkjQX+TdIJwCeBx4HRwALgE8Dlld5T0lxgLkBHR0etQzYrnC1btrB169ZGh1Ez48aNY8KECY0Ow6wp5ZrEB0XEU5KWAqdFxJdS8Q5J/wx8LO1vAqaUXTY5le17rwVkyZ9SqeSJ4G3E27p1K8ceeyxtbW2NDuWQDQwM8PDDDzuJmw0jz9Hp7akGjqQxwNuBByVNTGUC3g08kC5ZDPxxGqU+C3h6f/3hZvai/SXwZcvgiiuy92bXCv8RsZEnz9+xPPvEJwJLJa0BfkLWJ347sFDSWmAtMB7463T+HcAGoA/4OvCRHGM1a0nLlsGcOfDpT2fvh/JHpru7m9tuu43Fixfz/PPPv1B+wQUXDHn+UOXDnWtWVLX8HatEnqPT1wAzhyifPcz5Acyvd1xmI0l3N+zcCQMD2Xt3N3R2Dn/+xRdfzGWXXcZ1113Hm9/8Zu655x62bNnC9OnTGTduHACrVq3i1FNP5Qtf+AKTJk1i7dq1L1y/fft2LrroIk488US2bdvGtm3b+Pa3v82GDRuYP38+69ev57rrruOMM87ghhtuYMOGDZx//vkcc8wx9X0QZnVS7e/YofKMbWYjSFcXjB4NbW3Ze1fX/s8/66yzuPnmm+nt7eW4447jl7/8JZMnT+auu+56ybn9/f3MmzePadOmvVC2du1aSqUS8+bNY8yYMS/U2MePH8/y5cuZPn065557Lrt27XqhfFkR2vnNhlHt79ihchI3G0E6O2HJEvjc57L3A9UQSqUSt9xyC7NmzeL5559n48aNtLW1sXv37pecO3HiRG6++Wb6+vpeKDvxxBNZt24dN954I9u3b+fRRx9lYGCAPXv2MDAwQEdHB1dffTV9fX17lZsVVbW/Y4dKWat1ayiVSuHJXmyk6+3t5fjjj290GDXTaj+P2cGQtDIiSvuWuyZuZmZWUE7iZi2oVZqkW+XnMKuXhkz2Ymb1M27cOB5++OFGh1Ezg6PgzeylnMTNWsyECRM8w5nZCOHmdDMzs4JyEjczMysoJ/FhFGl+aTMzG5ncJz6Ewblvd+7MZtzJ4wv7ZmZm1XJNfAhDzX1rZmbWbJzEh5D33LdmZmYHw83pQxic+7a7O0vgbko3M7Nm5CQ+jM5OJ28zM2tuuTWnSzpC0r2S7pe0TtJlqXyhpIckPSDpm5JGpfIuSU9LWp1en8krVjMzsyLIsya+A5gdEc+lRH23pO8BC4H3p3P+FfgwcE3a/1FEvCvHGM3MzAojtyQe2Zqnz6XdUekVEXHH4DmS7gUm5xWTmZlZkeU6Ol1Sm6TVwBPAnRGxouzYKOAc4D/LLulMze/fk/S6Ye45V1KPpJ7+/v56hm9mZtZUck3iETEQETPIatsnSzqh7PDXgP+OiB+l/VXAqyPiJOCrwG3D3HNBRJQiotTe3l6/4M3MzJpMQ74nHhFPAUuB0wAkXQq0A39Vds4zEfFc2r4DGCVpfP7RmplZo3gK7P3LrU9cUjuwKyKekjQGeDvwt5I+DLwDmBMRe8rO/3VgS0SEpJPJ/sPxZF7xmplZY3kK7APLc3T6ROB6SW1kCfmmiLhd0m7g58AySQCLIuJy4Czgz9Lx54Gz0+A4MzMbAYaaAttJfG95jk5fA8wconzIGCLiauDqesdlZmbNaXAK7MGauKfAfinP2GZmZk3JU2AfmJO4mZk1LU+BvX9exczMzKygnMTNzMwKyknczMysoJzEzczMCspJ3MzMrKCcxM3MzArKSdzMzKygnMTNzMwKyknczMysoJzEzczMCspJ3MzMrKCcxM3MzAoqtyQu6QhJ90q6X9I6SZel8qMlrZDUJ+k7kkan8sPTfl86PjWvWM3MzIogz5r4DmB2RJwEzABOkzQL+Fvgqog4BtgGnJfOPw/YlsqvSueZmZlZklsSj8xzaXdUegUwG7gllV8PvDttn5n2ScfnSFI+0ZqZmTW/XPvEJbVJWg08AdwJPAI8FRG70ymPAZPS9iRgI0A6/jTwa3nGa2Zm1sxyTeIRMRARM4DJwMnAcYd6T0lzJfVI6unv7z/U25mZmRVGQ0anR8RTwFKgExgr6bB0aDKwKW1vAqYApOOvBJ4c4l4LIqIUEaX29vZ6h25mZtY08hyd3i5pbNoeA7wd6CVL5mel0z4AfDdtL077pOM/jIjIK14zM7Nmd9iBT6mZicD1ktrI/vNwU0TcLmk9cKOkvwbuA65N518LfFtSH7AVODvHWM3MzJpebkk8ItYAM4co30DWP75v+S+B9+QQmpmZWSF5xjYzM7OCchI3MzMrqAM2p0saV8F99qQR52ZmZpaTSvrEf5Fe+5strQ3oqElEZmZmVpFKknhvRLxkQFo5SffVKB4zMzOrUCV94p01OsfMzMxq6IBJPH3V65DPMTMzs9qq+HvikkrAp4BXp+tEtjjZ6+sUm5mZme1HNZO9LAQuAtYCe+oTjpmZmVWqmiTeHxGL6xaJmZmZVaWaJH6ppG8AS4Adg4URsajmUZmZmdkBVZPEP0i2/vcoXmxOD8BJ3MzMrAGqSeK/FRGvrVskZmZmVpVq5k6/R9L0ukViZmZmVammJj4LuF/SBrI+cX/FzMzMrIGqSeKnAaOBXWR94VWRNAX4FjAhXb8gIr4s6TvAYDP9WOCpiJghaSrQCzyUji2PiHnVfq6ZmVmrqiiJS/oT4GPAUcDLgVXAJcDXgDdU+Fm7gQsjYpWko4CVku6MiD8o+5y/A54uu+aRiJhR4f3NzMxGlEqWIv0EUALeGhGPp7LfBr5OFSuXRcRmYHPaflZSLzAJWJ/uKeC9wOwqfwYzM7MRqZKBbR8E/nAwgQNExPeBtwE/OJgPTU3lM4EVZcVvAbZExE/Lyo6WdJ+kuyS9ZZh7zZXUI6mnv7//YMIxMzMrpIpGp0fEriHKtgBfrfYDJR0J3ApcEBHPlB16H3BD2f5moCMtg/pXwL9K+pUh4lgQEaWIKLW3t1cbjpmZWWFVksQfkXTGvoWSLiebva1ikkaRJfCF5TO9SToM+D3gO4NlEbEjIp5M2yuBR4Bjq/k8MzOzVlbJwLaPALdKOhe4HzgSOB1YzYsjxw8o9XlfC/RGxJX7HH4b8GBEPFZ2fjuwNSIGJP0GMA3YUOnnmZmZtboDJvGI+Lmk3wLeARwPPAO8PyLWpEFvlToFOAdYK2l1KrskIu4AzmbvpnSAU4HLJe0im+Z1XkRsreLzzMzMWpoiqv7Kd9MqlUrR09PT6DDMzMxqStLKiCjtW37APnFJq2pxjpmZmdVWJX3ix0tas5/jAl5Zo3jMzMysQpUk8eMqOGfgUAMxMzOz6lQ0sC2PQMzMzKw61SxFamZmZk3ESdzMzKygnMTNzMwKqpr1xF8g6dPAm4BNwH0R8Q81jcrMzMwO6GBr4r8GLAc+D7y2duGYmZlZpQ42iW8D2oAnAE+FamZm1gBVNadLehnwsoi4TNKrgK8AD9QlMjMzM9uvipO4pPOBS4Gdkp4EvhIR59UtMjMzM9uvaprTLwROjIhJZCuanSLp0vqEZWZmZgdSTRJ/jqwPnIjYDJwH/H49gjIzM7MDqyaJXwPcLOmYtN8BbK/0YklTJC2VtF7SOkkfTeWflbRJ0ur0Or3smk9K6pP0kKR3VBGrmZlZy6u4TzwivibpceAbkk4CjgIWSnoPsDoifnqAW+wGLoyIVZKOAlZKujMduyoivlR+sqTpwNnA64BXAT+QdGxEeLEVMzMzqvyKWUQsioguoB14A/BDsklf/qmCazdHxKq0/SzQC0zazyVnAjdGxI6I+BnQB5xcTbxmZmat7KC+Jx4RuyNiTURcHxF/GRGzq7le0lRgJrAiFZ0vaY2kb0r61VQ2CdhYdtlj7D/pm5mZjSi5z50u6UjgVuCCiHiGrK/9NcAMYDPwd1Xeb66kHkk9/f39tQ7XzMysaeWaxCWNIkvgCyNiEUBEbImIgYjYA3ydF5vMNwFTyi6fnMr2EhELIqIUEaX29vb6/gBmZmZNJLckLknAtUBvRFxZVj6x7LTf5cUZ4BYDZ0s6XNLRwDTg3rziNTMza3YHtYrZQToFOAdYK2l1KrsEeJ+kGUAAjwJ/ChAR6yTdBKwnG9k+3yPT92/ZMujuhq4u6OxsdDRmZlZvuSXxiLgb0BCH7tjPNZ8nWynNDmDZMpgzB3buhNGjYckSJ3Izs1aX+8A2q4/u7iyBDwxk793djY7IzMzqzUm8RXR1ZTXwtrbsvaur0RGZmVm95dknbnXU2Zk1obtP3Mxs5HASbyGdnU7eZmYjiZvTzczMCspJ3MzMrKCcxM3MzArKSdzMzKygnMTNzMwKyknczMysoJzEzczMCspJ3MzMrKCcxM3MzArKSdzMzKygckvikqZIWippvaR1kj6ayr8o6UFJayT9m6SxqXyqpOclrU6vf8wrVjMzsyLIsya+G7gwIqYDs4D5kqYDdwInRMTrgYeBT5Zd80hEzEiveTnGambWkpYtgyuuyN6t+HJbACUiNgOb0/azknqBSRHx/bLTlgNn5RWTmdlIsmwZzJkDO3dmSxYvWeJFk4quIX3ikqYCM4EV+xz6EPC9sv2jJd0n6S5Jb8krPjOzVtTdnSXwgYHsvbu70RHZocp9KVJJRwK3AhdExDNl5Z8ia3JfmIo2Ax0R8aSk3wRuk/S68mvSdXOBuQAdHR15/AhmZoXU1ZXVwAdr4l1djY7IDlWuSVzSKLIEvjAiFpWVnwu8C5gTEQEQETuAHWl7paRHgGOBnvJ7RsQCYAFAqVSKHH4MM7NC6uzMmtC7u7ME7qb04sstiUsScC3QGxFXlpWfBnwceGtEbC8rbwe2RsSApN8ApgEb8orXzKwVdXY6ebeSPGvipwDnAGslrU5llwBfAQ4H7szyPMvTSPRTgcsl7QL2APMiYmuO8ZqZmTW1PEen3w1oiEN3DHP+rWRN72ZmZjYEz9hmZmZWUE7iZmZmBeUkbmZmVlBO4mZmZgXlJG5mTcPzeptVJ/cZ28zMhuJ5vc2q55q4mTUFz+ttVj0ncTNrCoPzere1eV5vs0q5Od3MmoLn9TarnpO4mTUNz+ttVh03p5tVyCOnzazZuCZuVgGPnDazZuSauFkFPHLazJqRk7hZBTxy2syakZvTzSpQ1JHTy5YVL+ai8TO2RsotiUuaAnwLmAAEsCAivixpHPAdYCrwKPDeiNgmScCXgdOB7cC5EbEqr3jN9lW0kdPux68/P2NrtDyb03cDF0bEdGAWMF/SdOBiYElETAOWpH2AdwLT0msucE2OsVoOPNq7vtyPX39+xtZoudXEI2IzsDltPyupF5gEnAl0pdOuB7qBT6Tyb0VEAMsljZU0Md3HCs41mPob7McffMbux689P2NrtIb0iUuaCswEVgATyhLz42TN7ZAl+I1llz2WyvZK4pLmktXU6ejoqF/QVlND1WCcxGurqP34ReJnbI2WexKXdCRwK3BBRDyTdX1nIiIkRTX3i4gFwAKAUqlU1bXWOK7B5KNo/fhF5GdsjZRrEpc0iiyBL4yIRal4y2AzuaSJwBOpfBMwpezyyanMWoBrMGZmhy7P0ekCrgV6I+LKskOLgQ8Af5Pev1tWfr6kG4E3Ak+7P7y1uAZjZnZo8qyJnwKcA6yVtDqVXUKWvG+SdB7wc+C96dgdZF8v6yP7itkHc4zVzMys6eU5Ov1uQMMcnjPE+QHMr2tQZmZmBeZpV83MzApKWYW3NUjqJ2uSr5XxwP/W8H72Un7G9ednXH9+xvXl5wuvjoj2fQtbKonXmqSeiCg1Oo5W5mdcf37G9ednXF9+vsNzc7qZmVlBOYmbmZkVlJP4/i1odAAjgJ9x/fkZ15+fcX35+Q7DfeJmZmYF5Zq4mZlZQTmJD0PSaZIektQn6eIDX2HVkDRF0lJJ6yWtk/TRRsfUiiS1SbpP0u2NjqUVpSWSb5H0oKReSZ5IuMYk/WX6G/GApBskHdHomJqJk/gQJLUB/wC8E5gOvE/S9MZG1XJ2AxdGxHRgFjDfz7guPgr0NjqIFvZl4D8j4jjgJPysa0rSJOAvgFJEnAC0AWc3Nqrm4iQ+tJOBvojYEBE7gRuBMxscU0uJiM0RsSptP0v2x29SY6NqLZImA2cA32h0LK1I0iuBU8kWdiIidkbEUw0NqjUdBoyRdBjwcuAXDY6nqTiJD20SsLFs/zGcYOpG0lRgJrCiwaG0mr8HPg7saXAcrepooB/459Rl8Q1Jr2h0UK0kIjYBXwL+B9hMtprl9xsbVXNxEreGknQk2RrzF0TEM42Op1VIehfwRESsbHQsLeww4A3ANRExE/g/wONnakjSr5K1gh4NvAp4haT3Nzaq5uIkPrRNwJSy/cmpzGpI0iiyBL4wIhY1Op4WcwrwO5IeJesOmi3pXxobUst5DHgsIgZbkG4hS+pWO28DfhYR/RGxC1gEvKnBMTUVJ/Gh/QSYJuloSaPJBlIsbnBMLUWSyPoSeyPiykbH02oi4pMRMTkippL9+/1hRLgGU0MR8TiwUdJrU9EcYH0DQ2pF/wPMkvTy9DdjDh48uJfc1hMvkojYLel84L/IRkN+MyLWNTisVnMKcA6wVtLqVHZJRNzRuJDMqvbnwML0n/0NwAcbHE9LiYgVkm4BVpF9o+U+PHvbXjxjm5mZWUG5Od3MzKygnMTNzMwKyknczMysoJzEzczMCspJ3MzMrKCcxM1aXFpp6yNp+1XpKzv1+qwZkk6v1/3NbG9O4matbyzwEYCI+EVEnFXHz5oBOImb5cTfEzdrcZIGV+F7CPgpcHxEnCDpXODdwCuAaWQLTYwmm4RnB3B6RGyV9BqypXnbge3An0TEg5LeA1wKDABPk02R2QeMIZum+ArgduCrwAnAKOCzEfHd9Nm/C7ySbHGhf4mIy+r7JMxaj2dsM2t9FwMnRMSMtGLc7WXHTiBbQe4IsgT8iYiYKekq4I/JVkJbAMyLiJ9KeiPwNWA28BngHRGxSdLYiNgp6TNkaz+fDyDpC2RTvn5I0ljgXkk/SJ99cvr87cBPJP1HRPTU8TmYtRwncbORbWlaz/1ZSU8D/57K1wKvT6vMvQm4OZu6GoDD0/uPgesk3US2MMVQfptsIZaPpf0jgI60fWdEPAkgaRHwZsBJ3KwKTuJmI9uOsu09Zft7yP4+vAx4KiJm7HthRMxLNfMzgJWSfnOI+wv4/Yh4aK/C7Lp9+/Lct2dWJQ9sM2t9zwJHHcyFaY33n6X+b5Q5KW2/JiJWRMRngH6y5Xv3/az/Av48rUCFpJllx94uaZykMWR98z8+mBjNRjIncbMWl5qsfyzpAeCLB3GLPwLOk3Q/sI5skBzAFyWtTfe9B7gfWApMl7Ra0h8AnyMb0LZG0rq0P+hesvXk1wC3uj/crHoenW5muUuj018YAGdmB8c1cTMzs4JyTdzMzKygXBM3MzMrKCdxMzOzgnISNzMzKygncTMzs4JyEjczMysoJ3EzM7OC+n8qArCwRGcxVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 566.929x198.425 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = TS.plot_quadpole_evolution(\n",
    "    obj,\n",
    "    [10, 11, 15, 14],\n",
    "    'rho_a',\n",
    "    threshold=0.05,\n",
    "    rolling=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Data Journal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "--- Data Journal Start ---\n",
      "2020-12-01 10:48:34.535947\n",
      "Data was imported from file pygimli_paper_data000.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data010.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data020.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data030.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data040.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data050.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data060.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data070.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data080.ohm (741 data points)\n",
      "Data was imported from file pygimli_paper_data090.ohm (741 data points)\n",
      "A filter was applied with query \"a > 10 and a < 20 and k > 1000\". In total 468 records were removed\n",
      "--- Data Journal End ---\n",
      "\n"
     ]
    }
   ],
   "source": [
    "obj.print_data_journal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}