{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with MS/MS data using HDF5\n",
    "\n",
    "The DreaMS project introduces a new data format for MS/MS data. This format is based on HDF5 (Hierarchical Data Format 5) and is particularly suitable for machine learning applications. HDF5 is a binary, non-human-readable format that stores MS/MS data as a collection of tensors. One of these tensors contains MS/MS spectra, while the others contain metadata such as precursor m/z values. There are three main advantages of this data format compared to existing ones such as `.mzML` or `.mgf`. \n",
    "1. It allows for handling huge files containing millions of mass spectra by leveraging the on-disk nature of HDF5. This means that data can be read or processed efficiently in chunks or spectrum-by-spectrum without loading the entire dataset into memory. \n",
    "2. The tensor-based nature of HDF5 allows for seamless integration of MS/MS spectra into machine learning frameworks such as scikit-learn or PyTorch. \n",
    "3. The HDF5 format supports advanced data compression techniques, significantly reducing file size and storage requirements.\n",
    "\n",
    "Notice that the .hdf5 data format can also be used to store complete LC-MS/MS runs, including MS1 data. We omit it in this tutorial for simplicity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Opening files and converting them to .hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DreaMS Python package implements the conversion of common mass spectrometry data formats, such as `.mzML` or `.mgf`, to `.hdf5`. The class responsible for this is named `MSData`. Let's look at an example of opening an `.mzML` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset G73954_1x_BC8_01_17287 into memory (1930 spectra)...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/builder/jenkins/ws/enms_ntly_pyoms_whl_Release3.0.0/OpenMS/src/openms/source/FORMAT/HANDLERS/XMLHandler.cpp(130): While loading '../data/MSV00008490/G73954_1x_BC8_01_17287.mzML': Required attribute 'softwareRef' not present!\n",
      "Warning: Parsing error, \"processingMethod\" is missing the required attribute \"softwareRef\".\n",
      "The software tool which generated this mzML should be fixed. Please notify the maintainers.\n"
     ]
    }
   ],
   "source": [
    "from dreams.utils.data import MSData\n",
    "msdata = MSData.from_mzml('../data/MSV00008490/G73954_1x_BC8_01_17287.mzML')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the `.mzML` file is parsed with `MSData.from_mzml`, it is auto converted to an `.hdf5` format and stored on disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001bc-rw-r--r--@ 1 roman  staff    12M Jul 16 23:29 G73954_1x_BC8_01_17287.mzML\n",
      "-rw-r--r--@ 1 roman  staff   4.0M Jul 20 10:35 G73954_1x_BC8_01_17287.hdf5\n"
     ]
    }
   ],
   "source": [
    "!ls -lhtr ../data/MSV00008490/ | grep G73954_1x_BC8_01_17287"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the data folder now contains a new `.hdf5` file with the same base name. The new file is approximately three times smaller in size, which is a result of two factors: `.hdf5` retains only MS/MS spectra and stores them in a memory-efficient tensor-based format.\n",
    "\n",
    "Similarly, one can open and convert an `.mgf` file. If a file does not follow a standard specification for metadata names, such as naming a precursor m/z as `PEPMASS`, one can specify the correct names using constructor arguments such as `prec_mz_col`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset MassSpecGym into memory (231104 spectra)...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "MSData(pth=../data/MassSpecGym.hdf5, in_mem=True) with 231,104 spectra."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata = MSData.from_mgf('../data/MassSpecGym.mgf', prec_mz_col='PRECURSOR_MZ')\n",
    "print(msdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following sections demonstrate how the content of the `.hdf5` file can be manipulated using the methods of the `MSData` class. Alternatively, one can handle the `.hdf5` file manually using the [h5py library](https://docs.h5py.org/en/stable/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing individual spectra and their metadata\n",
    "\n",
    "Method named `at` retrieves the spectrum and all its metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAADVCAYAAADgkwLjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAylUlEQVR4nO3deVhUZfsH8O/MIKCQCy7DoiJImChKmvmitihqPwFByXDfStMMTSst0FBTSl9Ty8oVX7c3MzC3MM2tzHhxKxdSDJHEZHfBxBFZ5vz+OM0Ass3gGWeG+X6uqyvnzOGc+9w8c7jnOc9zjkwQBAFEREREj0hu7ACIiIiobmBRQURERJJgUUFERESSYFFBREREkmBRQURERJJgUUFERESSYFFBREREkmBRQURERJKwMnYAj0txcTHu3LkDGxsbyOWspYiIiHSlVqvx4MEDNGrUCFZWVZcOFlNU3LlzB1evXjV2GERERGarTZs2aNq0aZXvW0xRYWNjAwBo3bo17OzsjByN8ZSUlODy5ct48sknoVAojB2OUTEXpUpKSpCcnAxPT0+LzgXzUIq5EDEPovv37+Pq1avav6VVsZiiQnPJw9bWFg0aNDByNMajVqvh5OQEOzs7i78MxFyUKikpAQA0aNDAok+czEMp5kLEPJRX07nSaGfSwsJCjBs3Dj/++CMAIDc3F+PHj4e/vz9CQkKQmpqqXXfbtm0YMGAA+vXrhxUrVhgr5DqjsLDQ2CGYDObC8PLzgf/8B5gzB9iwQXxNRHWTUYqKCxcuYPjw4Thz5ox22bx589C7d298//33eOeddzBjxgwAQFJSEqKjo/HNN99g7969OHXqFA4ePGiMsIlIT/HxQJs2wMSJwObNwIQJ4uv//c/YkRGRIRilqNi6dSumTZuGTp06AQCKiopw7NgxhISEAAB69uyJu3fvIiUlBUeOHEG/fv3QsGFDWFtbIyQkBHFxccYIm4j0kJ8PBAcDHToAf/4JXLsGpKYCXl5AUBB7LIjqIqOMqYiKigIAREdHAwDy8vJgbW0Ne3t77TpKpRJZWVnIzs6Gu7t7ueXZ2dm13rdardZeI6stmUz2SD9vTGq1utz/LZmp5EIQBKPuHyi9bvyon42ytm2T4fZtGbZskaF1a3GZqyuwZQvg5ibgm28EjBtn/GMvyxB5MFfMhYh5EOl6/CYxULOqk7pcLq/0hPsof9RTUlJq/bOafTs5OZn9tfj09HRjh2AyjJkLa2trZGZmmkRhAQCJiYmSbSshwRnOzkq0bl3+8+rqCjg7C0hIyIaPT4Zk+5OSlHkwd8yFiHnQjUkUFU2bNkVhYSFUKpV2ZkZOTg6USiUcHR2Rm5urXTcnJweOjo613peHh0e5HpHaqKqo8fPzQ0ZG6UnSysoKrVq1wtChQzF27NhH2qdU1Go10tPT4eLiYhIzHt5//3306tULgYGBtfr5kydPVpnbw4cPw9nZGdevX8cHH3yAc+fOwdnZGeHh4ejZs6c2FxkZGVi0aBGuXr0KT09PzJ8/H0899VSN+540aRIcHBzw8ccfa5ddvHgR8+fPR3JyMjw8PDBv3jx06NCh2u20aNFCv4M2gJKSEiQmJsLb21uyEe6+vjJs3CjDtWvQ9lQAQFoakJEhg6+vEj4+xj/2sgyRB3PFXIiYB5FKpUJycnKN65lEUWFlZYXnn38esbGxGDt2LBISEmBlZQV3d3f06dMHb731FiZOnAg7Ozvs3LkTL7/8cq33JZfLDdYwZDIZIiIi4O/vD0C8i+fx48cxe/ZsNGnSBIMGDTLIfmtDLpebRFFx4sQJzJw5s9axdOnSBb/88ku5ZdOnT0fjxo3RsmVLCIKAsLAweHp64ttvv8WhQ4cwdepUfP/993B0dERmZiYmT56MiRMnIjAwEOvXr0dYWBj2798Pa2vrKve7d+9e/Pzzzxg8eLA2dpVKhcmTJ2PgwIFYtGgRvv76a0yePBkHDx40m2nMCoVCss/HsGHA++8Do0aJlzxcXcWCYvRowMFBhqFDZTDVc7SUeTB3zIXI0vOg67Eb/6/KPyIjI3Hs2DEEBgZiyZIlWL58OWQyGdq3b4/XXnsNI0eORGBgILy8vDBw4EBjh1ulJ554As2bN0fz5s3h5OSEwYMHw9fXFwcOHDB2aCYnNTVVm6/asra21ua7efPmOHHiBJKTk7FgwQIAwPHjx/HXX3/hww8/RNu2bTFp0iT4+Pjg22+/BQDs3r0bnTp1QlhYGNq0aYOIiAjI5fJyU5oflpeXh3//+9/w9vYut/z777+HjY0NZs2ahbZt22L27Nmws7PD/v37a3185szeHtizB7h4EXB3B5ycADc38fWePeL7RFS3GLWnYsuWLdp/t2jRQjtw82FDhw7F0KFDH1dYkrOyskK9evUAAKNHj4anpyd++uknlJSUIC4uDnfv3sX8+fORkJCApk2bIiQkBG+88Ya2Mvz555+xfPlypKamwtXVFeHh4fD19cXnn3+OkydPlstjnz59EBYWhpCQkEr39e233yI6Ohp5eXl48sknERERgWeeeQYAcOXKFXz00Uc4c+YM7OzsMHToUEyZMgVyuRyff/45kpKScOfOHVy+fBlffPEFnn322XLH2a5dO3z66adYsWIFMjIy4OfnhxkzZmD27Nk4d+4cOnTogOXLl0OpVAIAEhIS8K9//Uv7sw9zcXHBkSNHsGPHDoSHh+OPP/6oNs9FRUX49NNPMXnyZDg4OAAAzp07By8vr3I9BV27dsXZs2cBAOfPny/XturXr49Dhw5Vu5/FixcjODgYOTk55ZafO3cOXbt21V4ek8lk6NKlC86ePaud2WRpevQArl4F2rUDCguB9euBV15hQUFUV5lMT0VdVFRUhAMHDiA+Ph5+fn7a5Tt27MCSJUvwxRdfwM7ODmFhYWjatCl27tyJjz/+GN999x1Wr14NALh8+TLeeOMN9OvXD7t370ZgYCCmTJlSbpxJdcru69q1a1iyZIm2+/+ZZ57B9OnToVarcevWLYwYMQItWrRAbGws5s6di//+97/YvHmzdluHDx9GYGAgNm3apJ0O/LAVK1Zg0aJFWLNmDQ4cOIDhw4dj+PDh2LZtG3Jzc7Fu3TrtusePH4evry8A4JdfftH+t3PnTjRs2BCvvvoqAMDf37/CJY7K7Nu3D3fv3sXIkSO1y3JzcyuMWWjatCmysrIAAJmZmbC1tcW0adPQo0cPjBkzptrBvAkJCTh9+jSmTJlS4b2a9mWp7O2Bxo2BFi2A8eNZUBDVZSYxpqIumTt3rrbrvaCgALa2thg7diyCgoK067z44ovo0qULAPGPVEZGBmJjYyGXy+Hu7o733nsP4eHhePPNN7F9+3Z06dJF+0fs9ddfh0qlwt9//61TPGX3dfDgQchkMrRo0QIuLi6YPn06evfuDbVajbi4ONSvXx8LFiyAlZUV2rZti9zcXHz55ZcYN24cAKBZs2YYPnx4tfsbN24cOnfuDABo37493NzcMGDAAABA//79cenSJQDigNHTp0/jo48+AgDtJZDi4mK89dZb6NmzJ0aNGgVAvLW6ra1tjccaExODIUOGlFv3/v37FcZGWFtba2fv3L9/H0uXLkVYWBgmTZqEzZs3Y9y4cfjhhx8qPCPmwYMHmDt3LiIjIyuNp6Z9ERHVdSwqJDZt2jT0798fgPgQs+bNm1cY4OLi4qL995UrV5CXl4euXbtql6nVahQUFOD27dv4888/K8wemD59us7xlN1Xr1694OnpiUmTJqF9+/bo27cvXnnlFVhZWeHKlSvo0KFDuUfaPv3008jNzdUWMGW3VZVWrVpp/21ra1vuZ2xtbbV/YC9evIhWrVrhiSeeKPfzS5Yswc2bN7F27VqdjxEAbt68idOnT+ODDz4ot9zGxgZ5eXnllhUWFmqLAoVCgd69e2P06NEAgAULFuDFF1/EkSNHKozd+eKLL9CxY0c899xzlcZgY2NToYAouy8iorqORYXEmjZtCldX12rXKfuUt+LiYri7u2PlypUV1nviiSeqfW59ZVNbi4uLq9xX/fr1tbc7v3DhAnbs2IGvv/4aO3bsqPTJc5r7h2huelLT0+mAiiOEq5rVkZCQoL30obFv3z5s27YN33zzjd7Tfo8dO4aWLVtWGJuhVCorXM64ceOG9jKFg4MD3NzctO9ZW1vDxcUFmZmZFfaxd+9e3LhxA08//TSA0ueG/PDDDzhz5gyUSiVu3LhR5b6IiOo6jqkwMjc3N2RkZMDBwQGurq5wdXXF9evXsWLFCshkMri6umovGWgMGzYMe/fuRb169XDv3j3t8nv37uHWrVtV7uvMmTNYu3YtfHx88P7772P//v148OABfv31V7i5ueHChQsoKioqt76DgwMaN24s+XEfP35cO0gTEHtsIiIiMGfOHJ3uEfGw8+fPay/zlNW5c2dcuHABBQUF2mW//vpruUs0ZQeAFhYW4q+//kLLli0rbGvLli347rvvsGvXLuzatQt9+vRBnz59sGvXLu2+zpw5o72RlSAI+O2337T7IiKq61hUGFmvXr3g4uKCmTNn4o8//tB24devXx8KhQLDhw/H6dOnsWHDBqSlpWHNmjW4fPkynnnmGXh7e+PSpUvYt28f/vzzT0RGRlZ7vwdbW1usXLkS+/btQ3p6Ovbu3QuVSoV27dph4MCBKCwsRGRkJK5cuYJDhw7h888/x/DhwyW/LXlhYSF+//13bRFw7949TJ06FX5+fvDz80Nubq72v5KSEhQUFNQ4MPXy5cvw8PCosPzZZ5+Fk5MTwsPDcfnyZaxduxbnz5/HkCFDAACDBw/GgQMHsHXrVly9ehUffvghbGxs8OKLLwIA7t69q7184uLioi38XF1dYWdnBzs7O23P1P/93//h77//RlRUFFJSUhAVFYX79+9rx5QQEdV1LCqMTKFQYNWqVVCr1QgNDcXUqVPxwgsvYM6cOQCA1q1b4/PPP8e3336LwMBA/PDDD1i9ejWUSiV8fX0xbtw4REZGYtiwYXjyySer/Vbcvn17LFy4ELGxsfD398fq1auxZMkStG3bFvb29oiOjsa1a9cwaNAgLFiwAGPHjkVYWJjkx3z27Fk89dRT2sspv//+O65cuYLvvvsOvr6+6NWrl/a/zMxMfP/99+jVq1e127xx4wYaNmxYYblCocDKlSuRm5uLkJAQ7NmzB19++SWcnZ21OVm+fDk2b96MgQMH4sqVK4iOjtZOQY2KisLUqVN1Oi57e3usWbMGv/76K0JCQnDu3DmsXbvWbG58RUT0qGSCqTx0wMBUKhWSkpLg6elZYXCgJVGr1fjrr7/QqlUrk7ijpjExF6VKSkpw9uxZ+Pj4GOSugZqxxhcuSL5pSRk6D+aEuRAxDyLN39D27dtX+0XJss+kREREJBkWFURERCQJFhVEREQkCRYVREREJAkWFURERCQJFhVEREQkCRYVREREJAkWFRbo4SdpWjLmgohIOiwqLIwgCMjMzISF3POsWswFEZG0WFRYIP4RLcVcEBFJh0UFERERScLK2AE8bNeuXYiOjoZMJoObmxs++ugjyGQyhIeHIyUlBTKZDAsWLKj0MddERERkPCZVVOTl5WHhwoXYv38/mjVrhjlz5mD9+vW4d+8eHB0dsWLFCiQnJ+P111/HgQMHOMiOiIjIhJjU5Q+1Wo2SkhLcu3cParUaBQUFsLW1xZEjRzBkyBAAgKenJ1q1aoWEhAQjR0tERERlmVRPhYODA6ZPn46AgAA0atQIjRs3xty5c/HFF1/A0dFRu55SqURWVlat9qEpXCyV5tgtOQcazEUpw+dC/s/21QbavjTYJkoxFyLmQaTr8ZtUUZGUlIStW7fi4MGDcHR0xKeffoqZM2dWOkJfLq9dJ0tKSsqjhlknJCYmGjsEk8FclDJULgoKvAAAZ89eNMj2pcY2UYq5EDEPujGpoiI+Ph7du3eHk5MTAGDUqFHo378/HB0dkZubi4YNGwIAcnJyoFQqa7UPDw8P2NvbSxazuSkpKUFiYiK8vb2hUCiMHY5RMRelDJ0LW1vxS4CPj4/k25YS20Qp5kLEPIhUKhWSk5NrXM+kigovLy/ExMTg9u3baNKkCQ4ePAhvb2+0b98eMTEx2hkgqamptZ79IZfLLbphaCgUCubhH8xFKUPnwlzyzDZRirkQWXoedD12kyoqevTogZEjR2L48OGwtraGUqnEokWL0LBhQ0RGRiIgIAAymQyLFy+26N4GIiIiU2RSRQUAjB07FmPHjq2wfNmyZUaIhoiIiHRlUlNKiYiIyHyxqCAiIiJJsKggIiIiSbCoICIiIknoNFBT31ti+/r61ioYIiIiMl86FRWvvvoqnJ2dK72z5cOys7Nx4cKFRw6MiIiIzItORYWtrS0OHz6s0waffvrpRwqIiIiIzJNOYyo2b96s8wb1WZeIiIjqDp16Kry9vcu9vnXrFtasWYM///wTLVu2xMSJE7XP63h4XSIiIrIMtZr9ERERAScnJ4wdOxaOjo4ICwuTOi4iIiIyMzoVFUuXLkV+fr72dU5ODoKCgtCzZ08MHDgQGRkZBguQiIiIzINOlz/at2+P8ePHY9CgQRg2bBjeeOMNBAUFQaFQoKCgAO+9956h4yQiIiITp1NR4e/vj759++K///0vRowYgYkTJ+Lo0aPaR5Rb8uNgiYiISKTzmApra2u8+uqrWLt2LU6ePIkJEyYgMzOTBQUREREB0OOOmsuWLYNCocCsWbMQERGBtLQ0LF26FPXq1cPbb78NFxcXQ8dKREREJkynnorZs2cjLCwMr7/+OmbPng0AcHV1xYoVKzBixAi8//77Bg2SiIiITJ9OPRUqlQotW7aEIAi4f/9+ufe6du2KLVu2GCQ4IiIiMh86FRWzZs3CyJEjUa9ePYSHhxs6JiIiIjJDOhUVzz33HEJCQnTaYG5uLpo3b/5IQREREZH50WlMxUsvvaTzBvVZtzLx8fF4+eWXERgYiClTpiA/Px+5ubkYP348/P39ERISgtTU1EfaBxEREUlPp56KgoICjBw5ssb1BEHAgwcPah1Meno63nnnHWzduhXu7u6YO3cu1q1bh5SUFPTu3RtjxoxBfHw8ZsyYgd27d9d6P0RERCQ9nYqKqKgonTf4yiuv1DqYgwcPws/PD+7u7gCAt99+G/n5+RgwYAAWL14MAOjZsyfu3r2LlJQUeHh41HpfREREJC2diorBgwcbOg4AwNWrV1GvXj1MnToV165dg5eXF95++21YW1vD3t5eu55SqURWVlatigq1Wo2SkhIpwzYrmmO35BxoMBelDJ8L+T/bVxto+9JgmyjFXIiYB5Gux69TUfG4FBcXIz4+Hlu3boVSqcTChQuxbNmySteVy2v1gFWkpKQ8Soh1RmJiorFDMBnMRSlD5aKgwAsAcPbsRYNsX2psE6WYCxHzoBuTKiqaN2+Obt26wcnJCQAQHByMiIgIFBYWQqVSoUGDBgDEp6Qqlcpa7cPDw6Ncr4elKSkpQWJiIry9vS3+FuvMRSl9cpGfD8TGypCaCrRtCwwZIqCmj5StrfglwMfHR6KIDYNtohRzIWIeRCqVCsnJyTWuZ1JFRZ8+fRAWFobs7GwolUocPnwYHTt2hJubG2JjYzF27FgkJCTAyspKO+5CX3K53KIbhoZCoWAe/sFclKopF/HxQHAwcPs24OICpKcD770H7NkD9Oih2/bNAdtEKeZCZOl50PXY9Soqnn/+efj7+8Pf3x+dOnWqVWDV8fb2xsyZMzFhwgQUFxfD1dUVH330EYqLixEREYHY2FhYW1tj+fLlkMlkku+fiKqWny8WFB06AFu2AK1bA2lpwOjRQFAQcPUqauyxIKK6Ta+iYuXKldi3bx+mT58OhUKhLTDatWsnWUCBgYEIDAyssDw6OlqyfRCR/mJixB4KTUEBAK6u4ms3NyA2Fhg/3rgxEpFx6VVUdOzYER07dsTMmTNx/vx5HDx4EDNmzIBcLkdAQAAGDRqkHQ9BRHVLaqp4yUNTUGi4uorLeU86IqrVFIqcnBycP38eiYmJyM3Nhbu7O+7cuYPQ0FD2KBDVUe7u4hiKa9fKL09LE5fXcpgTEdUhevVUbNy4Efv378fFixfh6+uLwYMH48svv4SdnR0AoH///njttdcwYcIEgwRLRMYTGgrMmgWMGiVe8nB1LR1T4eAAPMJ974iojtCrqDh8+DAGDx6M1atXo3HjxhXeb9OmDRYuXChVbERkQuztxVkeQUFir0SLFkB2tlhQ7NnDQZpEpOfljzZt2mDo0KEVCooZM2YAABwcHBAQECBZcERkWnr0EGd5ODoCxcXA+vXia12mkxJR3VdjT0VmZiZiYmIAALt27UKzZs3KvZ+fn4///e9/homOiEyOvT2g+V7B2R5EVFaNRYWTkxMKCgqQl5cHQRCQlZVV7v169epVeSttIiIishw6jal47733AIg3pxoxYoRBAyIiIiLzpFNRsWrVKrzxxhvIycnBZ599Vuk6b731lqSBERERkXnRqajQXPLIzs42aDBERERkvnQqKubPnw8A+Pjjjw0aDBEREZkvvaaUXr9+HXPmzAEA/Pzzz+jWrRv69u2L33//3SDBERERkfnQq6iYO3cuFAoFBEFAVFQUJk+ejNdffx3z5s0zUHhERERkLvS6o2ZSUhLWrVuH1NRUZGZmYuTIkbC1teVlESIiItKvp6JevXq4desW9u/fj65du8LW1hbJycmV3rKbiIiILItePRVjx46Fv78/CgoKsGrVKpw/fx4TJkzgdFIiIiLSr6h49dVX0bdvX9ja2qJFixa4desW1q9fD29vb0PFR0RERGZCr6ICAGxsbHDt2jWkpaVpl506dQrdunWTNDAiIiIyL3oVFWvWrMFnn32GJk2awMqq9EdlMhl++uknqWMjIiIiM6JXUfHVV19hzZo1eO655wwVDwBg48aN2LNnD3bs2IF79+4hPDwcKSkpkMlkWLBgAbp06WLQ/RMREZH+9Jr9UVRUBF9fX0PFAgA4f/48oqOjta8/++wzODo64vvvv8fy5cvx9ttvo7Cw0KAxEBERkf70KipGjhyJTz75BDdv3oRarS73nxT+/vtvzJs3D++++6522ZEjRzBkyBAAgKenJ1q1aoWEhARJ9kdERETS0evyx7Zt23Djxg1s2rRJu0wQBMhkMiQlJT1yMLNnz8aUKVPwxBNPaJdlZ2fD0dFR+1qpVGofcFYbarUaJSUljxSnOdMcuyXnQIO5KKV/LuT/rK/rFwp91zcOtolSzIWIeRDpevx6FRXffPNNrYLRxaZNm9CiRQv07dsXJ06c0C4XBKHCunK5Xh0s5aSkpNT6Z+uSxMREY4dgMpiLUrrmoqDACwBw9uxFg6xvbGwTpZgLEfOgG72KChcXF6hUKhw9ehRZWVkYNmwYUlJSJLlPxXfffYeCggIEBwdDpVIhJycHo0ePhqOjI3Jzc9GwYUMAQE5ODpRKZa334+HhAXt7+0eO11yVlJQgMTER3t7eUCgUxg7HqJiLUvrmwtZWLOx9fHx02r6+6xsL20Qp5kLEPIhUKhWSk5NrXE+vouLChQuYNGkSWrdujaSkJPTt2xdjxozB3LlzMWjQoNrGCgDYvn279t8nTpzA4sWLsWXLFnz88ceIiYnRzgBJTU19pNkfcrncohuGhkKhYB7+wVyU0jcX+ubNXPLMNlGKuRBZeh50PXa9riN8+OGHmDt3LrZu3QorKyu0atUK0dHRWLlyZa2C1MXUqVORm5uLgIAATJ8+HYsXL7bongYiIiJTpVdPRWpqKvz8/ACIN7wCgC5duuDWrVuSBtW9e3fs2LEDAGBvb49ly5ZJun0iIiKSnl49FR4eHoiLiyu37MiRI/Dw8JA0KCIiIjI/evVUzJkzBxMnTsRXX30FlUqF1157DUlJSVi9erWh4iMiIiIzoVdR0aFDBxw4cABHjx5FZmYmmjVrhqVLl6Jx48YGCo+IiIjMhV5FxWuvvYb169cjICCg3PKhQ4ca9B4WREREZPpqLCrS09Px2WefARCnes6aNavc+/n5+UhPTzdMdERERGQ2aiwqXFxc0KlTJ+Tl5QEAWrduXe59a2trvPXWWwYJjoiIiMyHTpc/Ro0aBQDo2LEjXnzxRUPGQ0QWLj8fiIkBUlOBtm2BV14BeGsaIvOg15iKLl26YO3atUhLS6vwZNKPP/5Y0sCoPJ5oyRLExwPBwcDt24CLC5CeDsycCezZA/ToYezoiKgmehUV7777LjIzM/HCCy/AykqvH6VHwBMtWYL8fLGdd+gAbNkCtG4NpKUBo0cDQUHA1asspIlMnV6VwW+//YbDhw+jUaNGhoqHHsITLVmKmBixcNa0cwBwdRVfu7kBsbHA+PHGjZGIqqfXHTUdHR1RUFBgqFioEtWdaG/dEk+0RHVBaqrYE/fQWHC4uorLU1ONExcR6U6vnoqXXnoJ48ePx6BBg+Dg4FDuvSFDhkgaGIl4oiVL4e4uXtq7dq18e09LE5e7uxsvNiLSjV5FxcmTJ9G0aVMcO3as3HKZTMaiwkB4oiVLERoKzJoFjBol9sS5upZe6nNwEAcnE5Fp06uo2LJli6HioCrwREuWwt5eHHwcFCQWyy1aANnZYjvfs4djh4jMgU5FRUJCQo3r+Pr6PnIwVBFPtGRJevQQBx+3awcUFgLr13P6NJE50amomD17drXvy2QyHD58WJKAqCKeaMmS2NsDmmcUcrYHkXnRqag4cuSIoeOgGvBES0REpk6vKaVEREREVWFRQURERJIwuXttb9y4EbGxsZDL5XB1dcWCBQtQXFyMWbNmITs7G7a2tvjkk0/gzrmUREREJsWkiorTp09j+/btiImJgZ2dHZYuXYqlS5fi9u3b6N27N8aMGYP4+HjMmDEDu3fvNna4REREVIZJXf5o0qQJIiMjYWdnBwDw8vJCWloajh07hpCQEABAz549cffuXaSkpBgzVCIiInqISfVUtG3bFm3btgUA5OfnY+XKlQgMDERSUhLsy8yfVCqVyMrKgoeHh977UKvVKCkpkSzmx0usAUtK1DWsVzXNsZtvDqTDXJTSPxf6tkVDry8NtolSzIWIeRDpevwmVVRoZGdn44033oCPjw8GDRqEdevWVVhHLq9dJ4s593AUFHgBAM6evfjI20pMTHzkbdQVzEUpXXOhb1s09PpSY5soxVyImAfdmFxRcenSJUyaNAmhoaF48803UVxcjMLCQqhUKjRo0AAAkJOTA6VSWavte3h4lOv1MCe2tmIh5ePjU+ttlJSUIDExEd7e3lAoFBJFZp6Yi1L65kLftmjo9aXCNlGKuRAxDyKVSoXk5OQa1zOpoiInJwfjx4/HnDlzEBAQAACwsrLC888/j9jYWIwdOxYJCQmwsrKq9ewPuVxu9g1DivgVCoXZ50EqzEUpfXOhb94Mvb5U2CZKMRciS8+DrsduUkXF2rVroVKpsHbtWqxduxYA4ObmhsjISERERCA2NhbW1tZYvnw5ZDKZkaMlqeXnAzEx4uPc27blrciJiMyNSRUVc+bMwZw5cyp9Lzo6+jFHQ49TfDwQHAzcvg24uIiPdZ85U3xoWo8exo6OiIh0YVJTSkl6+fnAf/4DzJkDbNggvjY1+fliQdGhA/Dnn8C1a2JvhZeX+HRWU4yZiIgqMqmeCirv4csBajWgz6SXyr/9y/Hvf9vhMY9/q1ZMjBjjli1A69biMldX8bWbGxAby4eoERGZA/ZUmKj4eKBNG2DiRGDzZmDCBODyZUCl0u3nq/v2/847Hib17T81VSx6NAWFhquruDw11ThxERGRflhUmKCqCgIbG+Cvv3S7HFD1t38Z7txRYPt20xno6u4u9qJcu1Z+eVqauJyPeSEiMg8sKkxQVQWBiwtQUiJeDgCqHy9R3bd/Z2fBpL79h4YCTZoAo0aJhQQg/n/0aMDBQZwFQkREpo9jKkxQVQVBvXqAlZX4fk2zJcp++y+7nbQ0ICNDBnd3QbvM2FM57e3FuIOCxLhbtACys8WCYs8eTislMjXGPmeQ6WJPhQmq6nJAURFQXAw4O9c8W6Lqb/8CGjUqwZAhYlFR2diNNm2A//3vsR4yevQArl4FHB3FY1y/XnzN6aREpsVUzhlkmlhUmKCqCoL0dEChAASh6tkSt26Jl0c03/4vXhSLFCcncSbFxYvA0qUpsLc3vamc9vZA48ZiT8X48fzmQ2RqTO2cQaaHRYUJqqogePAAaNUKyMjQbbZEZd/+r1xRo3PnewCqn8qpKU6IiDR4zqCasKgwUZUVBE8+CTRooN9sieq+/XMqJxHpg+cMqgmLChP2cEGgufGVVLMl9J3KaQ535yQiw+H0b6oJiwozVPbyiJsb0KiROFDqt9+AbdsqH4ugVosFwQcfyLBnT9MaBnNWLE44OIuIOP2basIppSYsP1+8fllUJPYMlL1Nd48ewNdfAwEBwN27YtdjZiYwbFjFh3CpVOI4jIkTARcXGdLTXfHll+J6ZadyNmsG5OQA9esDM2aUj0MzOEtzLVVzIgkKEi/TaAZ+1jTNjFPRiMwXp3/XXVKdm9lTYaI0PQPZ2eL9KR6+TXd+PjB8ONC9u/hH/dIlYNEiccxF//5AVlbpehkZwLPPakZry5CaKtOO1nZ3Bz78ULxb561bYtHSrJm4LU1PxObNNQ/O0qUng70dROaP07/rHinPzSwqTNDD07YyMireplszCvurr8RlbdoA778v/rxKBXh6ig0iJkb84H/1VcXbdd+8Ka4XFgYUFooFysPTxAYMAN59V/xGUtXgrEuXap5mZolT0TgGhcyFvm2V07/rDqnPzSwqTJAut+nWjMJ2cKjYIP78E+jUSWwQly4BSmXFgqBpU7EHpHNnYPFi8d4XW7eW39/q1cC9e0DLluJlkaoGZ+Xm1tyTYWlT0dgrQ+aCbdWySX1uZlFhgnS5TbdmFPaqVZU3iK++EhvEjRviJZSHC4JVq0p7MO7cqbzwOH5cLDZ27658cNaoUWJR06xZzdPMLGkqmiX2ypB5qq6tDh4sh0rFPxF1ndTnZrYYE1TTbbrd3UtHYa9YUX2DaN5cLERGjChfEHz+uaAtJNzdKy88zpwRuzfbt6/8ZlxnzojLn3qq5mlmljQVzdJ6Zch81dRWDx1qYtwAyeCkPjezqDBBNd2mWzMqd88eccbH9etVN4innhKfFXLqlKYgEODmJuDWrdJLGqGhlRceP/xQWmyUHZxVWCj2YCxZIi7XZZqZJU1Fs6ReGTJvNbXV9HQb4wRGj43U52azKSoOHTqEwMBA9O/fH5GRkSgqKjJ2SAbz8H0oHB3L36ZbMyiqRw9x6o9cLs4EqapBaO7CqRmt/cEHafjjD7W2Id28KRYeJ0+W39/Nm+I9MDSNzd4esLMDCgrEMRmjR1eM9+HnjGimmemyTl1hSb0yZN5qaqsuLg+MExg9NlKfm82iqMjNzcW8efOwbt06/PDDD1CpVPjqq6+MHZZBaXoGFApxXETZ23SXZWUlFhWnT1ffIOTy0tHaQUE34ehYviFdvy5eXgHEomH9evHEsndv+cZ2+bJY3Dzc2HSZZmYpU9EsqVeGzFtNbbVv39vGDZAeCynPzWZRVMTHx6NLly5wcnKCTCZDaGgo4uLijB2Wwdnbi0WFlVX523Q/TC4v3xOha4Mo25AEQeytaNdO7PbUTBN7uLE5O4vFTWXb1mWamSVMRbOkXhkyb9W11Z071WjQQG3sEOkxkercbBZ31MzJyYFSqdS+ViqVyNLc3UlParUaJSUlUoX2GIiVREmJ+qF/l39f0xMBAGPGqP9Zr/w6ACAIwj/viW/Wrw80bizHzZulP//wPjTrlFU+hqri1eWYavPz0tDkwBDtoXt34MoVwMtLjsJCYN06AUOGCLC3L/t7MR3650Lf35Oh15eGIduEqaqqrdavX4LExOpyYZzf0eNmWW2i6t+prsdvFkWFWl3xAOVVfW2vQUpKyqOG81h5ez8JADh79jKee87pn39nat9/7jknNGhgj65dS+cpln1fs87DEhMTH2kbD7//8HpVvV/TOrr8vNTK5kJqAweKx+PjkwlzaHq65kLf35Oh15eaIduEqaqqrVaVC2P/jh43S2gTUvxOZYLmq6sJ27VrF3766Sd8+umnAIBTp05hyZIliImJ0XkbKpUKSUlJ8PDwgL0F9z+XlJQgMTER3t7eUCgUxg7HqJiLUsyFiHkoxVyImAeRSqVCcnIy2rdvjwYPD+4rwyx6Knr16oVPPvkE6enpcHZ2xvbt29G7d+9abUsul1t0w9BQKBTMwz+Yi1LMhYh5KMVciCw9D7oeu1kUFc2aNcP8+fMxefJkFBYWonPnznjttdeMHRYRERGVYRZFBQD4+fnBz8/P2GEQERFRFcxiSikRERGZPrPpqXhUmhkkBQUFFn1dTDMtSKVSWXQeAOaiLOZCxDyUYi5EzIPo/v37ACqfjVmWWcz+kMLNmzdx9epVY4dBRERkttq0aYOmTZtW+b7FFBXFxcW4c+cObGxsan2PCyIiIkukVqvx4MEDNGrUCFZWVV/ksJiigoiIiAyLX9mJiIhIEiwqiIiISBIsKoiIiEgSLCqIiIhIEiwqiIiISBIsKoiIiEgSLCqIiIhIEiwqiIiISBIsKoiIiEgSFlFUHDp0CIGBgejfvz8iIyNRVFRk7JAMrrCwEOPGjcOPP/4IAMjNzcX48ePh7++PkJAQpKamatfdtm0bBgwYgH79+mHFihXGCtkgNm7ciICAAAwcOBBhYWG4ffu2ReZiw4YN8Pf3h7+/P8LDw/HgwQPcu3cP06ZNg7+/PwICAvDbb79p17eEz8zGjRsREhICABabi2nTpuGll15CcHAwgoODsXHjRov8fMTHx+Pll19GYGAgpkyZgvz8fIvMgySEOi4nJ0fo2bOnkJGRIajVauGdd94RNmzYYOywDOr3338XQkJChE6dOglHjhwRBEEQpkyZImzatEkQBEH45ZdfhKCgIEEQBOHixYuCn5+fcOfOHeHBgwfCqFGjhAMHDhgtdimdOnVKCAgIEPLz8wVBEIRPPvlEmD17tsXl4ty5c0L//v2Fe/fuCWq1Wpg6daqwYcMGISoqSoiKihIEQRD++OMP4YUXXhAePHhgEZ+Zc+fOCT179hQGDx4sCIJgsbl4/vnnhVu3bpVbZmmfj+vXrwvdu3cXrly5IgiCIERGRgrLli2zuDxIpc73VMTHx6NLly5wcnKCTCZDaGgo4uLijB2WQW3duhXTpk1Dp06dAABFRUU4duyY9ltZz549cffuXaSkpODIkSPo168fGjZsCGtra4SEhNSZ/DRp0gSRkZGws7MDAHh5eSEtLc3ictGpUyfExcWhQYMGuHfvHm7fvo1GjRrhyJEjGDJkCADA09MTrVq1QkJCQp3/zPz999+YN28e3n33Xe0yS8xFRkYGVCoV3nvvPQwcOBBRUVHIz8+3uM/HwYMH4efnB3d3dwDA22+/jdDQUIvLg1TqfFGRk5MDpVKpfa1UKpGVlWXEiAwvKioKL7zwgvZ1Xl4erK2tYW9vr12myUN2dnaF/GRnZz/WeA2lbdu2ePbZZwEA+fn5WLlyJXr16mWRuahXrx527dqF3r17Iy8vD3369EF2djYcHR2162jyUNc/M7Nnz8aUKVPg5OSkXWaJubhx4wZ8fX2xaNEifPvtt8jJycGiRYss7vNx9epV1KtXD1OnTkVwcLA2B5aWB6nU+aJCrVZXWGZpjz6vLAeAmAehkofUymQyQ4f0WGVnZ2PMmDHw8fHBoEGDKl3HEnIxaNAgnDx5Er169UJ4eHilxyuXy+v0Z2bTpk1o0aIF+vbtW265JeaiU6dOWLFiBRwcHGBtbY3XX38dp06dqnTduvz5KC4uxtGjRxEREYGdO3eifv36WLZsWaXr1uU8SKVufDqq4ejoiNzcXO3rnJycct9ILEHTpk1RWFgIlUqlXab5BlbX83Pp0iWEhobCz88PCxYssMhcXLt2DWfOnAEgnvwGDx6MS5cuVXq8dTkPAPDdd9/hxIkTCA4Oxpw5c3DlyhWMHj3aInNx+vRp7UBuQCysrKysLO7z0bx5c3Tr1g1OTk6Qy+UIDg7G+fPnLS4PUqnzRUWvXr1w+vRppKenQxAEbN++Hb179zZ2WI+VlZUVnn/+ecTGxgIAEhISYGVlBXd3d/Tp0wcHDx5EXl4eioqKsHPnzjqTn5ycHIwfPx6zZs3Cm2++CcAyc5GdnY2ZM2ciPz8fABAXF4dnn30Wfn5+iImJAQCkpKQgNTUVXbp0qdOfme3btyMuLg67d+/GwoUL0bZtW2zZssUic1FQUICoqCjcvXsXarUamzZtQr9+/Szu89GnTx+cOHFCewnj8OHD6Nixo8XlQSoyobK+nDrm8OHD+PTTT1FYWIjOnTtj4cKFsLa2NnZYBjd69Gi8+uqr6N27N3JychAREYGsrCxYW1tj4cKF8PLyAgB888032Lx5M4qLi9GnTx/MmjWrTnTnLVy4ELGxsWjTpo12mZubGyIiIiwuF5s2bUJMTAzkcjnatWuHyMhIyOVyREZG4o8//oBMJkN4eDh69uwJwDI+MydOnMDixYuxY8cO5OfnW2QuoqOjsWPHDhQXF6Nbt26YO3cu8vLyLO7zERcXhzVr1qC4uBiurq746KOPUFxcbHF5kIJFFBVERERkeHX+8gcRERE9HiwqiIiISBIsKoiIiEgSLCqIiIhIEiwqiIiISBIsKoiIiEgSLCqIyOS89NJLePDggbHDICI9saggIpOSlJQENzc32NjYGDsUItITb35FRAZ3/fp1jBgxAkOHDsWmTZtgY2ODDz/8ED/++CPi4uLg5OSE5cuXw9PTE6tWrUKjRo1gb2+PuXPnardRVFQEe3t7HD9+3IhHQkTVYU8FET0W2dnZuH//PhISEvDKK6/gzTffRMeOHXH8+HF4eXlhzZo1AICjR4/ixRdfRFBQEM6cOYMzZ87gp59+grOzM959910jHwURVYdFBRE9NuPGjYNCoUC3bt1gY2OD0NBQWFtbo3v37sjKysLt27dx//59ODs7a39GrVbjnXfewb/+9S8MGTLEiNETUU2sjB0AEVmORo0aAQDkcjmeeOIJ7XK5XA61Wo1jx46hV69e5X5m2bJl+Pvvv7Fy5crHGisR6Y9FBRE9NjU9yfHo0aMYNmyY9vW+ffuwc+dO7Nixo848GZSoLuPlDyIyCWq1GufPn8fTTz8NALh06RIiIyOxfPlyKJVKI0dHRLpgTwURmYSzZ89iwIABsLIST0ubNm3C/fv3MXXqVBQWFmrX27t3b7kxF0RkOjillIiIiCTByx9EREQkCRYVREREJAkWFURERCQJFhVEREQkCRYVREREJAkWFURERCQJFhVEREQkCRYVREREJAkWFURERCQJFhVEREQkCRYVREREJAkWFURERCSJ/weF7+VrHRWSsgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAAf5ElEQVR4nO3deVzU5do/8M9sCMgm24CYoCiCqEQsmohL4nIM9NEnLZVJzB6w0knrUTpR4JKGoh7QOkUmiaFGZvmg1TE0999ImiKiR1Z39n0fmJnv748vkscFmAVnhrnef9Ew9zXX6MtP9/3dbg7DMCCEEKIqrrYbIIQQ/UYxSgghaqEYJYQQtVCMEkKIWihGCSFELRSjhOgbhQKHD+OHHyCV4s4dlJUBwMWL2m7LcFGMEqJvPvoIVlZwdcX77yMjAzduAEBysrbbMlx8bTdAiCGRy3H+PCwtMWKE6kWqqhAYCAAyGQB8+y3OnMGff2qmQ6I8mo0S8gytWIHGRkgk+PJLVYZXVqKsrD09AcjlACASISoKPj4aa5IoiWajhDwrzc3g8zF1KgAsXYriYqxcifffh59ft4YfO4awMHh64qOPEBkJgQAhIeDxYGwMAP369WDnpFMUo4Q8KwIBWloAgGGgUGD7dqSmIjUVEydi9WpMnw4O58kDW1rwwQfYvh0MA1dXjBrVvqh/2Lp1Pds8eToO3VNPyLOzYweqq1FdjalTsWwZ+vfH1auorQWAUaNORkUFzJ4tEAgeHnErO9tl/nxkZ0MgwNq1WL0aPJ52midPQTFKyLPV1gaBAElJWLIEAMzM4OeHnJxmMzOzvDx7oTAiImLlypWWlpYMw+zcuXPlypUXfXw8ysuRkkIHQHUTxSghWnL8ODZvxm+/AYCR0RfTp2+7fj0/Px+AtbW1SCTKzMw8deoUh8OJXL78008/hamplhsmT0ExSohWZWVhy5a2M2cs7t2TyuXe3t48Hu/ChQvsL4VC4a5du15++WXt9kg6RzFKiOYUFWHvXpiaYskSSCQYMwYmJvjlF8yY0fm4klu3Nm7blpSU1NjYCIDP5ysUisDAwP379zs6Oj6T1onqKEYJ0ZywMHz+OYqLsXMn+vTBihWwtsby5dixozujKyoqdu3alZCQUFxcDCArK2vkyJE93DHRALr8nhANaWlBv37o2xdDhqCuDgC2bcOGDbh7t5sFbG1tIyMj8/PzPTw8AJSXl/dcs0SDKEYJ0RBj4/b0bG4GlwsA772HqCg895xSZUxNTUeNGgWgtLRU802SHkCX3xOiOaGhiI6GTIZVq3DyJNgrQN3clC0jFApBMao/aDZKiOZMmoSKCpw4AYUCYWEwNweA5cuVLUMxql8oRgnRqKtXcf487t1TpwYboyUlJRrqifQsilFCNEooBAD1JpIODg6g2aj+oBglRKM0EaO0qNcvFKOEaBTFqOGhM/WEaFKZi0v9uHF1DOOtRhGhUMjhcEpLSxUKBZdLcx1dR39DhGiSxNJyyNmzMVlZ6hQxMjKysrKSyWTV1dWaaoz0HIpRQjRJU+tx9iwTnazXCxSjhGiSRq5VKigouHXrlpGR0UXaNlkfUIwSoknnz58HUFpaqvJDf7755htvb+/m5mZTU9N58+ZptDvSIyhGCdGMurq6sLCwBQsWAGhtbf34449r2d1Buq2mpmbBggVvvPFGfX39K6+8kpeXZ2Ji0jPNEk2iB+URogEZGRmhoaH5+fkmJiZmZmbsw5ksLCzCw8NXrFjh5OTUZYXjx48vWrTo/v37FhYWcXFx4eHhPd810RCGEKKO1lb5xx8HDRoEwMPDY9CgQQCMjIy8vb05HA77s0gkys7OflqBlpaWyMhI9sKmMWPG5OfnP8v2ifooRglRQ04O4+fHAHUjRwZNmsRu6jly5MisrCyGYTIzM0UiEZ/PB8DhcIKDg9PT0x8pcO3ateeffx4An8+PiYmRyWTa+BpELRSjhKgqOZkxM2MApn9/xsvr94kTORxOeHh4Y2Pjw++6efOmWCw2fbAhnY+PT3JyskwmUygUiYmJ7OuDBw8+d+6ctr4HURMdGyVEJYmJWLoUAMaOxdWrqK+Hs/P5pKQxL730xLeXl5d//vnnn332WWVlJQBXV1eFQnHz5k0AERER27ZtM6WNP/UWxSghymAY/PIL6uowbRpefhk8Hs6dA4A5c/DVV7Cx6Xy0VCpNTU395JNP8vLyeDyelZVVUlLSzJkzn0XnpMdQjBKijPXrMXYs+vfHtm1YuhRjx6JPHyQkYPHi7teQyWTR0dFtbW1hYWGenp491yx5NujRJIQo4/59TJ4MAKamGDkSSUl48UUMHqxUDT6fv3Hjxh5pj2gDXX5PiDIYBgoFADQ3QyDAwoXKZmiHffv2abIxoj00GyWGor6+XiwWt7W1paSkqF4lIgKrV8PUFEFB4HDU6cfMzEyd4UR30LFRYihkMllKSoqjo+O0adO03QvpVShGCdGOwsJCHo/n7Oys7UaIuujYKDEgubm5P/74o7a7aLdhw4aDBw9quwuiATQbJQZkz549//rXv/bu3ctR77AmIQ+jGCWEELXQop4Q7Vi0aNHAgQMvX76s7UaIumg2SgyLjY3NvXv3dOFxyHl5ecbGxo6OjuwjoIj+ohglhqWlpcXY2FjbXZBehRb1uiUuDuwmZl98gR072l/s+IGoT0cyVCaTrV+//sqVK9puhGgAxahuKS7Grl2Qy5Gbi5yc9hdzc7XaU+/y+++/NzQ0qFXi7l0sW4YNG/DBB6oVKCwsHD9+fHR0tEgkksvlajVDdAAdlNEtXC4WL8aXXwIAw2DDBgCorNRuU72Kk5MTu12H6r78EtHRsLfHjh3IysLduxg8GB4e3Rx94MCB8PDwmpoaBweHzZs383g8tZohOoBmozrH3x/37kEmA4eDqChERcHGBgoFCgu13VmvMGzYMHUfkNzYCAsLALC0RFUV3nwTnp4ICWl/8OjT1dTULFy4cN68eTU1Na+88sq1a9emT5+uVidEN/DWrFmj7R7IX9ra4OaGgQPRpw88PTF0aPvrAgHu30dWFoYN02p/ek4qlS5atCghIcHe3t7NzU3FKqam+PZbmJnhu+9gaoqiIpSX49o1JCXh2LFmoZA/dOjjl/cfP3582rRp586ds7Cw2L59+6ZNm3ThagGiGVrbvoQ8RVAQAzBHjz76emoqc++eNhrqLa5fv+7t7Q2AXdR7e3uzeyKpUis/n/n1V6aqihkyhAEYa2tmwgTG2poBvg0MHDJkSHx8fHNzM/veRzb+zMvL0+S3IjqAYlTnLFzIAMyePf/xYmYmk5zM3L6tpZ50QG1trVgsDg0NLS0tVXasQqHYvn07O/tzcXGZPXu2nZ0dO41wc3P76quvWlpaVOlJLmcOHmRGj2YABmD69m2YMsXjweNHnZycNm/enJGRQRt/9noUozrnvfcYgImL03YfuiQ7O9vLy4vH43E4nD59+ohEohs3bnRzbGlpaXBwMBttIpGovr6eYZiWlpbk5ORhDw6R2Nvbx8TEVFZWqtjfmTNMcDDD4RydOJHL5Xp7ew99cDiGnYS6ublduHBBxeJE51GM6pxNmxiA+d//1XYfmlNeXh4XF/fpp5/K5XJlx8rl8ri4uD59+gB47rnn/P392WDi8Xjz5s27ePFi58N//fVXBwcHAFZWVvv373+8eFpa2ujRo9nIMzMzE4vFd+7cUbbJdllZq8PD2a3qORzO4AfT0o7sJr0VxajO2b2bAZjQUG33oSEXLlxwdXVlL+txdXWNj49vamrq5tji4uK//e1vHWHU0NDAMEx+fr5YLO44PxMQEJCWlqZQKB4Z29TUJBaL2fcEBQXd6/S48pkzZ4KDg9nzQgKBQCQSZWdnq/Z9i4uLY2JirKys2I/28/NTrQ7RIxSjOufYscJx43YuXvxPbTfCMAxTVFS0b9++kydPqjBWJpNt3LiRnZ3Z2dk5OTmxyeLo6BgbG1tTU9P58B9++MHGxoZdcaelpT3y25KSkpiYmH79+rE1R40alZyc3NbWxv72jz/+YE/EGxsbx8bGdnMWfOnSpfnz57N3uHM4nOjoaBW+Nauqqio0NBRASEiIykWIvqAY1TnsI39GjRql7UaYjIwMNvs4HE5wcPD58+e7P/b27dsTJkxgx4aHhzc2NrIraH9/fzb4zM3Nn7aCrq2tDQ8PZ982bdq0oqKip31KXV1dfHz8gAED2De7uLjExcVFRUUZGRkB8PT0zMzMVPZb37x5kz0OMGXKFGXHPkwikQDw9/dXpwjRCxSjOqeoqAiAUChUs05FRcWJEydqa2tVGNvW1hYTE8POy8zNzdlUAjB58uSjj1+K9Zjvv/+enScKhcKff/75kd8+vII2MjJ6ZAUtkUiGDBkCwMTEJD4+/vHV+uOkUunXX3/t7u7ONsnhcLhc7vvvv6/i+XeGCQkJATB79mzVhrMKCwsBODs7q1OE6AWKUZ3T1tbG5XJ5PJ46F8dIJBJra2sA1tbWMTExVVVV3R978+bNcePGsXkkFotbWlo6X0E/rKamhl3MApgzZ05FRcXTPiUzM1MkEnWsoIODg48cORITE8MeRfX19e3+uXiWXC4/dOjQoEGDBg4cePDgQaXGPmLp0qUAZs6cqU6RpqYm9qiCOkWIXqAY1UXsVY0lJSUqjG1tbf373//OhlHH04wsLCxWrVp1//79LocnJyebm5uzp8VPnDjx8K8eX0HHx8ezp31Y586dY09Pm5qaxsfHd6fbwsLCd9555+G7M3k83ocfftja2qrk99YkdpeR4OBgNev84eNzp39/pqujwETfUYzqIl9fX4FAMGrUqCeeg+5EQUHB2LFj2TCKjIyUSqWPr6CvXbv2xLHV1dXz589ns2zu3LlPu4hSKpUmJycPHz6cfaetrW1kZOTt27c7JpL+/v65ublKfd/y8vLo6Ghzc3NbW9tTp04pNbYn5Obm9uvXb8aMGeoWYu9xysnRRFNEd1GM6qKUlBT2SkkAL7zwwv79+7uzwE9OTu7bty97PO706dMP/+rxFfTZs2cffkN6ejp7NsnCwiIxMbHLz5LL5T/99BMb2XhwkTmfz1+zZs0TF/v6RSaTubq6Tp8+Xd1CAQEMwPzn3wXpfShGdYhMJtu6dSv78xNX0I2NjU8cWFZWNnPmzI6JZHV19RPfVlhYKBaLO1bQPj4+ycnJDQ0ND9/xnZ+fr1TPZ86cmTlzJpfLtbGxkUgkSo3VZaNHj546daq6Vf77vxmA+f57TXREdBfFqE7oWLk/EpTsCtrjwbMsbW1tY2JiHjlvc/ToUUdHRwCWlpZ79+7t8rPKysqio6PZSzLZlT4AgUCwceNGlU9qVVdXq3CHki4LDg4OCgpSt8rbbzMAs2OHJjoiuotiVCcEBwdnZWU97bfsFZcdK2j2pvLc3Fz2Rh32uOfkyZPv3r3b/U9saGhITExkZ7t2dnZd3lVpaJYuXTp+/Hh1q/z4I7NiBS3qez2KUW0qKSlhT0nfunWrO+8/ffp0x/kiPp/PLs+NjY23bdum1JmoDk1NTb6+vt38dIOSkpLi7u6ubpXYWGbtWmb5cjrL1LvR0++1adWqVevXrwfg7OzcnfcHBgYePnw4JyeHnYQ2NTU5OztLJJKVK1c+/pzg7jAxMWEY5vr16yqM7d38/Pza2trUKlFUhMZGREcjNhaff66hvoguor2YtCAjI8PKymrYsGFxcXHV1dXKDh86dGhCQgL75M1x48Z1PAVDNRYWFunp6R1PACEsV1dXMzMztUrU14O9YcHUFFKpRroiuolmo1pw6dKlBQsWtLa2CoXCjlsYleXn5xccHKxmhgLo06dPXl6emkV6n+PHj1dUVIjF4ubmZhVLDBmCK1dw5Qq2b8eUKRrtjugYbR9VMEQKhWLXrl1SqVTNOl9++aX6RX738blpZ6dmkd6EfZ4Aj8dj7wGzs7N7/OqIzqSn//VzfT2TlsZcudITfRLdwWEYRttJTlTEblDRcaG+iqZPx4kTqKkB7bAG5ObmhoaGXrhwgcfjBQYGNjQ0XLx4EYCZmdmSJUtWrlzZxVFsmQyjRyM8HBERz6hjogNoUa/HRo8erW6GAqiqApeLq1c10ZF+q9qzx9fH58KFCy4uLidOnDhx4sSFCxfYu2kbGxsTEhIGDx4cEhLCBuvj1qxZc+vePezdC0vLZ9w50S6KUX1WVISCAnWLlJWhtRXHjmmiIb1VXo7Zs60XLdr5wgtz5869dOlSYGAg+5tx48YdPnz4ypUrIpGIx+MdOXLEz8+PfbFjNLuys7e3v3HjBtzd8dprWvoaREu0ekiBqOeDD5iVK9Ut4uDAAMzcuZpoSD/98kv7H0K/fvLU1E7eyG4QYvlgsunt7Z2QkCCVSteuXfu0B74QQ0DHRg2epye4XAwfjtRUbbfyDDEMiothb48DB7BwIRgGkyYhORnPPdfl0Lq6um+++Wbz5s1FRUV8Pj8sLGznzp3PoGWis2hRb/BOnsSyZe1RcvZs+4tnzmi1px7W0oL/+R8cPYqlS+HrCw8PxMTg2LHuZCgACwuLd999Nz8//+2335bJZL/88ktP90t0HMWoPvvHPzBpklqHNRUKvP8+Zs5Enz74xz/+mpD27plpejpmzcLixVi1CmlpyMzEmjXgKvdvwcTEJDo6GoC6NzsR/Ud3MemzESPQ3Az2CvzsbIwYAYbB9evw9OxuhZISDBkCR0c4OuLQIfD52LABACoqeqpnXcAw4PMBgMuFQgGBQLUytra2PB6vsrJSJpOxz3Ilholmo3qrthb79iEkBIcO4fx5JCYCgEyGXbuUKGJl1Z6YjY1g962LikJUFGxte6BjnTF5Mg4cwP/9H7Zuxbx5Kpfh8Xi2trYKhaK8vFyD3RG9Q/8L1VsSCf7rvzByJJycsHUrSkuxYQMUCuWKmJrCzw9r16K+HitWYP/+9tcdHTXerw7p2xdbt2LvXowbh+49FOZphEJhaWlpaWmpY+/+EyOdohjVW5aW+Pe/AaCiAlZWEAoRFYW2NkRGKldHJPrr5w8/bP8hKkpDXeqq5maIxXB0xIN9TFUjFAoBlJaWaqgtopdoUa+3xoxBbi4++QQ7dmDx4vZbOTkcPNgNlHTG3h4cDsrLlZ6//ycHBwcAJSUlGmqL6CWajeotDgexsX/95+bNAMDnY+NGbXWkTwQC9OuHqipUValzIJhmowQ0GyWGSygEAPUSUHdiVCqVrlmzpqqqStuNGCKKUWKoHBwYDqetV8Rodna2v7//2rVrly1bpt1ODBPFKDFQEc7OAi73gHqHNdkYVfPYaHNz88yZM318fL777juFksdqGYb56quvRo8enZWVNWjQoHfeeUedTohqKEaJgTKxtJTL5WpOJNlTTOoUyc7OHjNmzOHDhzMzM+fPn+/m5paQkNDS0tKdsaWlpSEhIREREU1NTSKRKCsrKyAgQOVOiMooRomB0sh63MbGBkBeXl5hYaGyYxmGSUhI8PX1zcrKcnNzCwsLc3FxKSgoWLFihaura1xcXF1dXSfDf/rppxEjRvz88892dnaHDh3as2ePuptHEZVp+QlThGjJunXrAISFhalTpKGhoU+fPhwOh8vlBgcHZ2RkdHPgnTt3Jk2aBIDD4YSHhzc0NDAMI5fL09LS/P392X+b5ubmYrH47t27j4xtamoSi8Xse6ZMmXL//n11vgJRH8UoMTgKhSIxMZHdOGDQoEEq72dVW1vr5eUVHBw8d+5cI/ZWWmDKlCm//fZb5wMPHDhgbW0NwN7ePi0t7fE3sI/cZzfNNjIyEolEHc8zzcjIGDp0KABjY+PY2Fi5XK5a80SDKEaJYbl///7UqVPZyGOfJyIUCmNiYqqqqpSt4+vr6+HhwaZwSUlJTExMx0atXl5eycnJbW1tj4yqra0NDw9n3zNt2rSioqJOPuLixYvz5s3j8XgAuFzuzJkzX331VYFAAGDEiBFXaKc8nUExSgxI7pEj7NFMGxsbPz+/h49uWVpaRkZGdp5rHS5duuTu7m5tbX3s2LGHX6+rq4uPjx8wYABb08XFJT4+vrGxkf2tRCJxdXUFYGJiEh8fr1AouvNZhYWFYrHY1NSUDVMOhyMWi1taWpT97qTnUIwSw9DUxIjFDJ//5ogRL7zwAnt+ydLSMiUlhV1Bs8HHrqCvX7/eSaW0tDQXFxcOhyMSiZ74BqlUmpyc7OHhwda0tbX98MMPly1bxs4rfX19b9y4oWz7paWlc+bMmTBhQmJiorJjSU+jGCUGQCJhhgxhAMbEpGbqVAczMwAvvfTSw2dvLl++LBKJ2GU+e77o3LlzTyz27bffDhgwwNPTkz0v9DTs+aKxY8eyYcrlcnk8XlRUVGtrq4a/HdE2ilHS233yCcPnMwAzbBgzeDAD5EycGBcX98STMw+voAEEBASkpaU9vvo+efLkI8v5Tpw+fdrT09PDw+PQoUPqfheik2hLO9LbbdyIjz5CYCDOn0drKzw8sHcvvL07GVFWVrZjx45//vOf7C3q7u7uBw4cGDFixLPqmOgZilHSS8XHo7YWVVVYtAhbt2LfPnA4WL4csbHtDxXsSkNDw7p167Zs2cLhcM6ePfviiy/2dMtET9FdTKQ3qqvD3buIicGWLfj6a2zaBHd3pKUhIaGbGQrAzMxs3bp1zz33nEKh6PyGImLg6HmjpDdqaoK5OQAIBGAYDBiAa9eU3fsTgLGxsa2t7Z07d6RSqeabJL0FxSjpjRwccO8eLl5EZibGjAGgQoay7OzsALS2tmqwO9LLUIySXmrHDpw4AR+fzs8mdcne3h4AzUZJJ+jYKOmlTEwwY4aaGQqAfVBIN59cRwwTxSghnQkICDAxMaFFPekExSghnfHw8LCwsKAYJZ2gGCWkM8bGxpaWlhSjpBMUo4R0wdramo6Nkk5QjBLSBRcXF5qNkk5QjBLShbFjx7a1tWm7C6K7KEYJ6UJlZeXJkydv3Lih7UaIjqJHk5Deq6wMSUlQKBAWhv79VSjQ3Nz8wQcfbN++ncfjcbnc1157LTIy0tPTU+OdEr1GMUp6r6VLsWkTeDysWIGvv0ZjI/r27f7oixcvhoaG5uTkGBsbOzk53b59WyaTcTicl19+efXq1YGBgT3XONEvtKgnvZqlJczMwOejrAxOTnj9dfz7310OUigUCQkJAQEBOTk5w4cPl0gk+fn5eXl57BOdjxw5Mn78eB8fnz179sjl8mfwJYiu0+5TownpQW++yUilTFsbs2QJs3kzw+UyAMPlMnPmMBLJ0wbdvn17/PjxeLCDfMeGdKzy8vLY2FhHR0f2n4+rq2t8fHxTU1PPfxmiu2hRT3qvggLs3g0OB7NmISgI5uYYOBB//omWFgCJItGAV1+dMWMGux0868CBAxEREdXV1UKhMCkpacaMGU8sLJVKU1NTN2zYkJubC8De3v6tt95666232J3yiMHRdo4T0vMyM9u3tAMYOztm4sRGb28uhwNg6NCh8fHxzc3NNTU1oaGh7D+KOXPmVFRUdFm1ra1t//793g+efmJiYlJQUPAMvg3RNTQbJYZBocDPP2PjRpw/D+Cr8eNTGKagoKCoqAiAvb09j8crLi42NzffsmVLeHi4UrXT09MXLFhQVVX1+uuvf/PNNz3SP9FhdIqJGAYuFyEhkEhw/HjLvHkrJZIzZ86Ul5dPmDDBwcGhrKysuLj4xRdfvHz5srIZCmDKlCmRkZEKhcLa2roneic6jmKUGJiXXjJOTf1/Fy4sXLiQYZhTp06VlJTweLxZs2adPn3a1dVVtarsUdHS0lKN9kr0A8UoMUReXl4pKSn5+fmRkZECgUAul8+dO5fPV30zCIpRQ0YxSgyXs7NzbGzskiVLoHYCsjFaUlKimc6IXqEYJYZu4MCBUDtGHRwc1C9C9BTFKDF0GklAW1tbHo9XWVkpk8k01BfRGxSjxNC94OSUOn58hKmpOkV4PJ6NjY1CoSgvL9dUY0Rf0AbLxNB52dh4nT6N+no167AXTpWWlnbcKkoMBM1GicFj7+BU++wQnaw3WBSjxOAJheBwUF4OhUK9MhSjBopilBg8gQD9+kEmQ2WlOmUUCgXomieDRDFKCODgAHNz1NSoNloul2/atCk1NRVAenq6RjsjeoC3Zs0abfdAiLbZ2MDDA337wtUVhw7B3R337yMnB05OXQ4tLCwMCQnZvXs3AC6XW1BQkJaWZmJi4unpyeXSNMUg0F8zMXjp6WhowEcfQSLB9es4fhwASktx9WqXQ8/s3//8889LJJL+/ft7eXkpFAoOh3P58uVFixYNHz5899dfgza4NwAUo8Tg/fkngoIAYPJkXLqEq1exYQOSkroYVVuLhQsDli8fZGIyZsyY5ubmzMxMe3v7gwcPJicnu7u75+Xlte7ejYEDsWYNqqqewfcg2kIxSgze0KHIzgaAq1fh5oaRIxEVhTfe6GzIsWPw9MS+fdy2toP+/n9kZFRXV8+aNSs7O3v27Nmvv/56dnb2999//6qNDcrLsXYtnJ3x3nu4excAWlvbi3RcGKDeFQJE6+ixzcTgKRTYtg1SKWxssHQpEhMREYGCAty82T5LfZhcjlWrEB8PhoGXFyorce/ewaCgyrlzn/yg0t9/R1wcjh4Fw2DePPj5oaUF5eWYNw8//ohPP4WREZYtw2efPYMvSnoIxSghSnr1Vfz4IwICcPYs5HL4+SElBW5unQ25cgVxcYiIwA8/ICEBCgXefht9++KVV8Dn44svuj6GQHQY3QxKSPesWwceDxUVWL8eJSU4dQo8Hj7+GNHR6PJBpV5eSElBbS2OHAEALhfsPnrV1eDzQbs06zmKUUK6ITMT1tZYtgyFhfjuO3z8McLDsXs3xo9XooilJRoakJ6O3FyMHYvMTAQFwcgIaWk91jd5FihGCemG6mrY2wOAUIiqKgQFIScHAoHSdRISkJGByZPh7o5x49orvPeehrslzxYdGyWkG+rq8O67iInBd9/B3x8vvaTthogOoRglpHuKi3HsGEaMwION6QlhUYwSQoha6PJ7QghRC8UoIYSohWKUEELUQjFKCCFqoRglhBC1/H/tzo/8qf5dHgAAAux6VFh0cmRraXRQS0wgcmRraXQgMjAyMi4wOS41AAB4nL3SW0gUURgH8DNnx7246311dXXb2VXX0b3ohkZlOkciCg3rJaPE6kGzlvSlIO3BVvMSBZZoBlZasSAU2hpavqQ7Vlo9qEReIwkKIzNDowsh2My3oYTaYwOH/2++83FmPjhfetqmkPD4Ie+DhWUXVpKwnJQUHRGSkngT41XJJIv7FL3cCAUsFJaxbut6Z/7zW2ue+PcOJZF5E8u9Sa0kK6REOBKSXnlf8x/xqkY5OiYmRa278R8+oUSUMBvCEmET0T60eJZUhmRympIrkMKXxgolUqqQyg/5+SP/ABQQiAKDUFAwCg5BIWqkDqWxOgyFaGgcFE5TsggUoWW0kSzWRiFdFE3pNqANeqRnWKw3MDojMkaz2BjDSGNpHGtiTHEsNrEonqVxfAJKMCOzhcVmKxNvQ7ZEFtuSGLV4jzQ5wk2i/lwqZA995JMhswz3ii/F5D254NB7RMe0d5OlnXNQLzxVR6pLbFBPfZJFFu1S8IOiBFIzbAEbQraTe7la6I98M8ZFD41tE/3xTD+Xd6WPg6tcm5G+uU5BRI/OH0jbZTkO3thV2Ztf6oIe69kIT4n0G7h1odBzIq8afONii8c5Mwo++LTf86FiL9i1f8FTM9gNdj9T8d8z1eCROiN/6X4TuNO2iVea3OmiJ35k8VZrKdQfNhzmXeNTaaL1Mdl87rQO/qcIJ/PNlnqYxXHdl880OMDN45/Sb9bngLvc4dxPLhRmT9YZessadoPLb/3iKqrTwJOUlDy/XAsuMuUR/87b4JqxAtJY/hLsl68hE+Y5r8uayWvXNNi55RCnSr0Dvtax2DOgaQR/LXNyDfJJcPbj8+TVbBXMxTsqyb6OtzDXnNtOUlrCYZYy111ib4sDz2e1EsZ1DtxeNUjys1vABusLMhs5BB5vaiZ27h146ug8SeyfAe85eZoUDFwFjyyxZMfWEnBKXzH53DAADvsNCQDeABwoI6MAAAQnelRYdE1PTCByZGtpdCAyMDIyLjA5LjUAAHicfVbbbiM3DH3PV8wPWOBNEvnQh02yu1gUmwBt2n/oe/8fPZQdjbNgansImz4jHR5eNA9Hvv54/v2ff4/9kueHh+Og//lExPG3EtHDzyO/HI9fv/94OZ7evjy+e55e/3p5+/OwOMxxD94fsV/eXn++e/h4Orxp9zH0uHAbMtTsoEbrdd4qAM7GKj76cZFmQT69ACqAvc3J0meu6IP7lAJoAFrrhnfkilMG0yiAHUBt7OpdEkgxnLgADgCljSFhY61IJDIL4FwrkjkZoqY2SGJUKzqAEMWHac+/P1MnFq53gn75RQh0CxwjO3CzBwDgyhbmUQF5AamrUIArCED6CpiJgcykYzp+NnGL3itkZgaqaEDJCWSwyaAKmam5YFMx9g4kWSA7FTJzc7EW3mdoSjWDtFKSMzmX0YS5RwKcDatWyMzOZbY+lcOxppk6lWpmei7eZnhE0gMJGaWcsVSC8oMGWcpEoFwWOi0oN0XWKQsXKSerdBJeUBQaKoly35iqVgkl1zRpIpwybCZ3rvpC9HhZqzryBOLotTEzDwXUjtfFVTr1WA2MippRSSA9sY72UEb6UfXIf59lXANQwrY0dCEXkRI5V5n2GHZtOVQpeqlCem6PVpPQbGM0dIw+qrTKtZuiO3XLRdGbqlX8OeyOHAjIlWYdWu9cZlV5tTz+N4R86Y3cbFZEVRYSRSroJ1Rsn45CrJCKkBDy6OaSFWsm3quINPOELuniKf1E9U00V4XsS1DHhMs52xu2nlTO2XEV1KCnreBHoEzL7Seg6GJC31uG5Ix7qt5TX0QFRWrZpPLZ8NZY05uzMxlrTwhbVjMaJ0WacxgChlpuk6utjdfBoUM9a6ANjihbyWQBBelBr4ID6yy7w3LgDeQc1esp/2eJNFtA9EV2Me6Ywr3S3Potj2sy4DgktlGl0TI5GGG2mn2iLcaIUp7rUYTj1zCfDZNEqew086UjNkZK8g7E3Et5AkBtU9egw8DVEVqt+PXl+cOjwPXh4PH15fl8OMi3nI8A+HHoedAzLjuPc8bVz0ObcWHof/v+G58HdDrneQwzLj8PW8YVH28JOHHdn6C8DC/cZsN882+6nHzTnJiknGaT5mSdZvPmJJ5m3J1cvMzmzUk8zabOyT1N3J04nEbuDxZOIzs2WfLyOhzuYxG5+U+P3u61X5B23RfD/YNq0peWsoPIx5I0soPA6JZldhCSQaRnByEZxG0i3y+/amNN1bv5yml0y485Kmn0LJhVMfBs+dVu62z5Mf8kjW7mmszTs5nrtW4wqz5ooUkfRjd9DClJY/fTiJdnE7ckjpFiv9STyc1/epI+jJ0lj4GWxjZ9S93Ts+nbuK2z6VsKD2NbeEvm6dnMszfvOzF/vz/i4/vDfwPPRYIcRQIeAAACPHpUWHRTTUlMRVMgcmRraXQgMjAyMi4wOS41AAB4nH1SO45UMRC8CuGO5LH6/2GEtNIjIGIjIsQ5NuHwVHtgQyYZv3J3u6q6ruu6fl6vr99+vVy3v4e36+XL2+0DPl/4fv4BHfCj5zq/f52n6xz+1t2+3/437tPvl9rqFbHuvENCbT1ys0oBkW1NVevhO5Pl1FSw53rYdnM7NSnBvB66ubR8EOpAk+wI6QMkkeiUkBXZutMOkg6eOZJatmjLDJZBWIt08U5345mTGSZLNgqnwreGlgIIbjwAQKRIFlpZs9cjtrnWAJ2lMoByES9cpHDMK+1hC1KJIQOa3bQYgGpED1eYgudtpyn1kWMEyrjI8piKVAqMFI2GOvhXYSOFxkftgdwpRosLMTwhmNRzks3W1oOQq1ALBoepYhBcJo0stEnZCISF2oWJvJsNi3jcUS3Gg5B1ThcklecgkU2DgBqzNxgVG6oAQWfCEN4GufiunT0FaJSRDSBCMRc7wlOz6juDdUAp6FA9EcYuYqSqGRZ6EIHUSQiEzWTG3qxQAnY6Ft9Z5zwQI1Usz6JCamCHRT4NMjylOe/jgo9D3vD1ZE1h1ShDgKTVZcEH7ph1wBovOnFDssbHCaJOjg1xcM5ZGQ6G4b6pDOyAeLPUglf+jAoeCccVrDITBPpkRaoHEcqQoVRRFD6DYFbyoWRgZDrPRY8pWA1p5cwuVvd1W+8/PoPM4gbrJb2ll2EX9P71c2w4DPdoCbxEzkAEoTBY/PsPhJ7aqw1dGhgAAAAASUVORK5CYII=",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x3203d52a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "{'COLLISION_ENERGY': 30.0,\n",
       " 'FOLD': 'val',\n",
       " 'FORMULA': 'C34H59NO14',\n",
       " 'IDENTIFIER': 'MassSpecGymID0000378',\n",
       " 'INCHIKEY': 'UXDPXZQHTDAXOZ',\n",
       " 'INSTRUMENT_TYPE': 'Orbitrap',\n",
       " 'PARENT_MASS': 705.393024,\n",
       " 'PRECURSOR_FORMULA': 'C34H60NO14',\n",
       " 'SIMULATION_CHALLENGE': 'True',\n",
       " 'adduct': '[M+H]+',\n",
       " 'precursor_mz': 706.4003,\n",
       " 'smiles': 'CCCC[C@@H](C)[C@H]([C@H](C[C@@H](C)CCCCCC[C@H](C[C@@H]([C@H](C)N)O)O)OC(=O)C[C@@H](CC(=O)O)C(=O)O)OC(=O)C[C@@H](CC(=O)O)C(=O)O'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.at(333)  # Random entry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling all MS/MS spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\"spectra\" object is a <class 'numpy.ndarray'> tensor of shape (213548, 2, 128) containing float64 elements.\n"
     ]
    }
   ],
   "source": [
    "spectra = msdata.get_spectra()\n",
    "print(f'\"spectra\" object is a {type(spectra)} tensor of shape {spectra.shape} containing {spectra.dtype} elements.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dimensions of the tensor have the following interpretation. First dimension (213548 elements) is a number of spectra. Second dimension (2 elements) represents m/z and intensity values, where first element is m/z values and second element are intensity values. Third dimension (128 elements) represents the number of signals in each spectrum. One can notice that since a tensor-shaped structre demans each dimension to have a fixed number of elements, each spectrum has the same number of signals. This was achieved during the conversion of .mzML to .hdf5. Let us look, for example, at the first spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 91.0542, 125.0233, 154.0499, 155.0577, 185.0961, 200.107 ,\n",
       "        229.0859, 246.1125,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ],\n",
       "       [  0.2452,   1.    ,   0.0801,   0.3554,   0.3493,   0.045 ,\n",
       "          0.1421,   0.7347,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ,   0.    ,   0.    ,   0.    ,   0.    ,\n",
       "          0.    ,   0.    ]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print spectrum values without scientific notation\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=4, suppress=True)\n",
    "\n",
    "# Show the first spectrum\n",
    "spectra[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectrum has seven signals, with the rest of the elements being zeros. These zeros, often referred to as padding zeros, are used to fill the tensor up to a desired shape (let's denote this as n, where n = 128 in our case). Conversely, when a spectrum has more than n signals, only the n most intense signals are selected and stored in the tensor. The remaining signals are discarded and not stored in the file. The value of n can be controlled using the `n_highest_peaks` argument of the `MSData` constructor. To ensure that no signals are lost, one can set `n_highest_peaks` to the maximum number of peaks in the entire dataset. Since the padding zeros are effectively compressed in HDF5, this will not result in excessively large file sizes.\n",
    "\n",
    "To remove the padding zeros, one can use the `unpad_peak_list` function from the DreaMS package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 91.0542, 125.0233, 154.0499, 155.0577, 185.0961, 200.107 ,\n",
       "        229.0859, 246.1125],\n",
       "       [  0.2452,   1.    ,   0.0801,   0.3554,   0.3493,   0.045 ,\n",
       "          0.1421,   0.7347]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dreams.utils.spectra import unpad_peak_list\n",
    "unpad_peak_list(spectra[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling all metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To retrieve all columns parsed from the input file, one can use the `columns` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ADDUCT',\n",
       " 'COLLISION_ENERGY',\n",
       " 'FOLD',\n",
       " 'FORMULA',\n",
       " 'INCHIKEY',\n",
       " 'INSTRUMENT_TYPE',\n",
       " 'PARENT_MASS',\n",
       " 'PRECURSOR_FORMULA',\n",
       " 'SIMULATION_CHALLENGE',\n",
       " 'SMILES',\n",
       " 'precursor_mz',\n",
       " 'spectrum']"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.columns()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To retrieve a column, one can either use the `get_values` method or square bracket notation. The second option is preferred, as explained in the \"Working with huge files\" chapter of this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['C16H17NO4', 'C16H17NO4', 'C16H17NO4', ..., 'C21H43N5O7',\n",
       "       'C21H43N5O7', 'C21H43N5O7'], dtype='<U17')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_values('FORMULA')  # Or msdata['FORMULA']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To retrieve specific elements, one can pass indices as a second argument to `get_values` or index the column obtained with square brackets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'C16H17NO4'"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_values('FORMULA', 0)   # Or msdata['FORMULA'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['C16H17NO4', 'C16H17NO4', 'C16H17NO4'], dtype='<U9')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_values('FORMULA', [0, 1, 2])  # Or msdata['FORMULA'][0, 1, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MSData` class provides methods for more advanced functionality, such as converting spectra and their metadata to [`matchms`](https://github.com/matchms/matchms) or merging multiple `.hdf5` files. Please refer to the class reference documentation for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with huge files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MSData` class enables working with huge files without loading all their content into memory. Let's look at an example of working with the large-scale [GeMS-A10 dataset](https://huggingface.co/datasets/roman-bushuiev/MassSpecGym/blob/main/data/spectra/GeMS_A10.hdf5) introduced in the DreaMS paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001bc-rw-r--r--  1 roman  staff    14G Jul 20 00:58 ../data/GeMS_A10.hdf5\n"
     ]
    }
   ],
   "source": [
    "pth = '../data/GeMS_A10.hdf5'\n",
    "!ls -lhtr ../data/GeMS_A10.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GeMS-A10 dataset takes 14 GB on disk (the uncompressed HDF5 version would take 47 GB), which is too large to fit into the RAM of an average laptop. Nevertheless, `MSData` allows for efficient handling of the file by loading only the necessary chunks of data at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../data/GeMS_A10.hdf5 contains 23517534 spectra.\n"
     ]
    }
   ],
   "source": [
    "msdata = MSData.from_hdf5(pth, in_mem=False)\n",
    "print(f'{pth} contains {len(msdata)} spectra.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When setting `in_mem=False` in the `MSData` constructor, none of the spectra or their metadata are loaded into memory initially. The data is loaded only when it is accessed. Let's retrieve the first 50 spectra, for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50, 2, 128)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_spectra(range(50)).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or 50 first precursor m/z values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50,)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_prec_mzs(range(50)).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get all precursor m/z values (load all of them into memory), one can use either of the two equivalent options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((23517534,), (23517534,))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msdata.get_prec_mzs().shape, msdata['precursor_mz'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, note that using `msdata['precursor_mz'][:50]` is deprecated for large datasets because `msdata['precursor_mz']` first loads all the values into memory, and only then `[:50]` retrieves a chunk."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dreams",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}