diff --git a/doc/user-guide/filtering-and-resampling.ipynb b/doc/user-guide/filtering-and-resampling.ipynb new file mode 100644 index 0000000..9fbe14d --- /dev/null +++ b/doc/user-guide/filtering-and-resampling.ipynb @@ -0,0 +1,685 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b314e777-7ffb-4e62-b4c5-ce8a785c5181", + "metadata": {}, + "source": [ + "# Filtering and resampling Xarray datasets with xbatcher\n", + "\n", + "There are many cases in machine learning where we want to discard invalid observations or modify the distribution of a target variable. This notebook demonstrates how `BatchGenerators` can be used to make filtered or resampled datasets by passing functions that identify usable data or assign a sample weight to each patch." + ] + }, + { + "cell_type": "markdown", + "id": "7158f5f3-42f5-4dcd-87ee-045d9d0e85f5", + "metadata": {}, + "source": [ + "### Libraries and toy data" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5d912ff0-d808-4704-8dea-b9e1b5a53bf1", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "import xbatcher as xb" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7fb892c1-50fd-48c8-8567-b150946b53c9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 31MB\n",
+       "Dimensions:  (lat: 25, time: 2920, lon: 53)\n",
+       "Coordinates:\n",
+       "  * lat      (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0\n",
+       "  * lon      (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0\n",
+       "  * time     (time) datetime64[ns] 23kB 2013-01-01 ... 2014-12-31T18:00:00\n",
+       "Data variables:\n",
+       "    air      (time, lat, lon) float64 31MB ...\n",
+       "Attributes:\n",
+       "    Conventions:  COARDS\n",
+       "    title:        4x daily NMC reanalysis (1948)\n",
+       "    description:  Data is from NMC initialized reanalysis\\n(4x/day).  These a...\n",
+       "    platform:     Model\n",
+       "    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...
" + ], + "text/plain": [ + " Size: 31MB\n", + "Dimensions: (lat: 25, time: 2920, lon: 53)\n", + "Coordinates:\n", + " * lat (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0\n", + " * lon (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0\n", + " * time (time) datetime64[ns] 23kB 2013-01-01 ... 2014-12-31T18:00:00\n", + "Data variables:\n", + " air (time, lat, lon) float64 31MB ...\n", + "Attributes:\n", + " Conventions: COARDS\n", + " title: 4x daily NMC reanalysis (1948)\n", + " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n", + " platform: Model\n", + " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly..." + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds = xr.tutorial.open_dataset('air_temperature')\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "85fbfe7b-c006-4052-b276-41e2354263a4", + "metadata": {}, + "source": [ + "### Filtering\n", + "\n", + "Here we add a QA variable to the air temperature dataset. Suppose that 1% of the time there is an instrument failure, and we do not want any cells with the QA flag set to go into a model." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "3e448fe0-af4c-49e9-be8d-529b29cb12fc", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAHHCAYAAACY6dMIAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWBFJREFUeJzt3Xl4TNf/B/D3ZF8nC7KRxb4LjYqxlJJK8VVbW0UrTVVLY021pUXoIlq19huUEl0oVXSzlViqRBB7S2JPShZLs5LFzPn94Zf7NZKQ2TKZmffree7zJOeee+/nJO3k42xXJoQQICIiIjIBVsYOgIiIiKiqmLgQERGRyWDiQkRERCaDiQsRERGZDCYuREREZDKYuBAREZHJYOJCREREJoOJCxEREZkMJi5ERERkMpi4kFkLCgrCq6++auwwiIhIT5i4kMk7ePAgZs6ciZycHGOHUq2OHDmCsWPHomXLlnB2dkZAQABefPFFpKamVlj/7NmzePbZZ+Hi4gJPT0+88soruHHjRrl6n3zyCZ577jl4e3tDJpNh5syZFd5v8+bNCA8Ph5+fH+zt7VGvXj08//zzOHPmjEbtWLlyJZo3bw4HBwc0btwYX3zxRbk6KSkpmDRpEjp16gQHBwfIZDJcuXJFo+dU97OIyDCYuJDJO3jwIGbNmlVh4pKSkoIVK1ZUf1DV4NNPP8XGjRvRs2dPLFq0CG+88Qb++OMPPPHEE+WSh3/++QdPPfUULly4gNmzZ2Py5MnYsmULnnnmGZSUlKjVnTZtGo4cOYJ27do98vmnT5+Gh4cHJkyYgCVLlmDMmDE4fvw4OnTogJMnT1apDV9++SVef/11tGzZEl988QUUCgXGjx+PTz/9VK1eYmIiFi9ejPz8fDRv3rxK9zbms4jIgASRiZs7d64AIC5fvmzsUKrVgQMHRHFxsVpZamqqsLe3F8OHD1crHzNmjHB0dBRXr16Vynbu3CkAiC+//FKtbtnP8caNGwKAiImJqXJMmZmZwsbGRrz55puPrXvnzh1Rq1Yt0bdvX7Xy4cOHC2dnZ3H79m2p7NatWyIvL08Iod3vuzqfRUSGxR4XMmkzZ87EO++8AwCoX78+ZDKZWtf+w3NcVq9eDZlMhj///BPjx49HnTp14O7ujjfffBMlJSXIycnBiBEj4OHhAQ8PD7z77rsQD71AXaVSYeHChWjZsiUcHBzg7e2NN998E//++291NRsA0KlTJ9jZ2amVNW7cGC1btsTZs2fVyjdu3Ij//Oc/CAgIkMrCwsLQpEkT/PDDD2p1g4KCtI7Jy8sLTk5OVRq227NnD27duoW33npLrTwqKgqFhYXYsmWLVObp6QlXV1et46rOZxGRYdkYOwAiXQwaNAipqan4/vvvsWDBAtSuXRsAUKdOnUdeN27cOPj4+GDWrFk4dOgQli9fDnd3dxw8eBABAQGYPXs2tm7dirlz56JVq1YYMWKEdO2bb76J1atXIzIyEuPHj8fly5fx3//+F8ePH8eBAwdga2tb6XOLi4uRn59fpbaVtUUTQghkZWWhZcuWUtm1a9eQnZ2N9u3bl6vfoUMHbN26VePnPCgnJwelpaXIzMzEwoULkZeXh549ez72uuPHjwNAubhCQkJgZWWF48eP4+WXX9YpNmM8i4gMi4kLmbQ2bdrgiSeewPfff48BAwZUubfA29sbW7duhUwmw1tvvYULFy5g7ty5ePPNN7F06VIAwBtvvIGgoCCsWrVKSlz+/PNPfPXVV1izZg2GDRsm3e/pp5/Gs88+iw0bNqiVP+z7779HZGRklWJ8uKenKtasWYNr167hww8/lMoyMjIAAL6+vuXq+/r64vbt2yguLoa9vb3GzwOAjh07IiUlBQDg4uKCadOmYeTIkY+9LiMjA9bW1vDy8lIrt7OzQ61atXD9+nWt4jH2s4jIsJi4kEUaOXIkZDKZ9H1oaCgSExPV/uBaW1ujffv2SE5Olso2bNgANzc3PPPMM7h586ZUHhISAhcXF+zZs+eRiUt4eDh27typ59bcd+7cOURFRUGhUCAiIkIqv3v3LgBUmJg4ODhIdbRNXOLj45GXl4dLly4hPj4ed+/ehVKphJXVo0ei7969W26o68G4yuLWh+p8FhEZFhMXskgPzvUAADc3NwCAv79/ufIH566cP38eubm55f7lXiY7O/uRz/X19a2w50NXmZmZ6Nu3L9zc3PDjjz/C2tpaOufo6Ajg/jDVw4qKitTqaEOhUEhfv/TSS9JKnM8//xwAcOPGDSiVSqmOi4sLXFxc4OjoWG5F04NxaRNTdT6LiIyDiQtZpAf/sD+u/MEhG5VKBS8vL6xZs6bC6x83t+bu3bvIzc2tUow+Pj5Vqpebm4vevXsjJycH+/fvh5+fn9r5skSpbMjoQRkZGfD09NS6t+VhHh4e6NGjB9asWSMlLk8++SSuXr0q1YmJicHMmTPh6+sLpVKJ7OxstUSwpKQEt27dKteOqqjOZxGRcTBxIZP34JCPoTVs2BC7du1C586dtfpX+vr16/U6x6WoqAj9+vVDamoqdu3ahRYtWpSrU7duXdSpUwdHjx4td+7w4cNo27ZtleKpqoeTszVr1qgNxTRo0AAApOcePXoUffr0kc4fPXoUKpVKq7iq81lEZBxMXMjkOTs7A0C17Jz74osvYsmSJfjoo48we/ZstXP37t1DQUEB3N3dK71en3NclEolhgwZgsTERPz8889qQzYPGzx4ML7++mukp6dLw2EJCQlITU3FpEmTtHr+w70XAHDlyhUkJCSord7p3Llzhdf36NEDnp6eWLp0qVoysXTpUjg5OaFv374ax1SdzyIi42DiQiYvJCQEAPDBBx/gpZdegq2tLfr16yclNPrUrVs3vPnmm4iNjcWJEyfQq1cv2Nra4vz589iwYQMWLVqE559/vtLr9TnH5e2338Yvv/yCfv364fbt2/juu+/Uzj+4vPf999/Hhg0b8PTTT2PChAkoKCjA3Llz0bp163I9QN9++y2uXr2KO3fuAAD++OMPfPzxxwCAV155BYGBgQCA1q1bo2fPnmjbti08PDxw/vx5rFy5EqWlpZgzZ85j43d0dMRHH32EqKgovPDCCwgPD8f+/fvx3Xff4ZNPPoGnp6dUNzc3V9qe/8CBAwCA//73v3B3d4e7uzvGjh1bY55FRAZm1O3viPTko48+EnXr1hVWVlZqO50GBgaKiIgIqV58fLwAII4cOaJ2fUxMjAAgbty4oVYeEREhnJ2dyz1v+fLlIiQkRDg6OgpXV1fRunVr8e6774rr16/rvW2V6datmwBQ6fGwM2fOiF69egknJyfh7u4uhg8fLjIzMzW67549e6R6MTExon379sLDw0PY2NgIPz8/8dJLL4lTp05p1I7ly5eLpk2bCjs7O9GwYUOxYMECoVKp1Opcvny50pgCAwNr5LOIyDBkQmixWQQRERGREXDLfyIiIjIZTFyIiIjIZDBxISIiIpPBxIWIiMgC/PHHH+jXrx/8/Pwgk8nw008/PfaavXv34oknnoC9vT0aNWqE1atXGzzOx2HiQkREZAEKCwsRHByMuLi4KtW/fPky+vbti6effhonTpzAxIkT8frrr2PHjh0GjvTRuKqIiIjIwshkMmzevBkDBgyotM57772HLVu24MyZM1LZSy+9hJycHGzfvr0aoqyY2W9Ap1KpcP36dbi6ulbr1vBERGR6hBDIz8+Hn5/fY99wrouioqJKX/ypCSFEub9t9vb2enn/WGJiIsLCwtTKwsPDMXHiRJ3vrQuzT1yuX79e7o2/REREj5Keno569eoZ5N5FRUWoH+iCzGzl4ys/houLCwoKCtTKyl4uqqvMzEx4e3urlXl7eyMvLw9379412lvVzT5xcXV1BQB0QR/YwNbI0VB125x6+pHnBzZpXU2REJEpuIdS/Imt0t8OQygpKUFmthKXkwMhd9W+VycvX4X6IVeRnp4OuVwulevrbe81ldknLmVdaDawhY2MiYulkbtaP/I8/5sgIjX/P+uzOqYWyF2tdEpcpPvI5WqJi774+PggKytLrSwrKwtyudxovS2ABSQuRERENZFSqKDUYXmMUqj0F0wFFAoFtm7dqla2c+fOR76JvjpwOTQREZERqCB0PjRRUFCAEydO4MSJEwDuL3c+ceIE0tLSAABTp07FiBEjpPqjR4/GpUuX8O677+LcuXNYsmQJfvjhB0yaNElvPwNtMHEhIiKyAEePHkW7du3Qrl07AEB0dDTatWuHGTNmAAAyMjKkJAYA6tevjy1btmDnzp0IDg7GvHnz8NVXXyE8PNwo8ZfhUBEREZERqKCCLoM9ml7dvXt3PGrrtop2xe3evTuOHz+uaWgGxcSFiIjICJRCQKnDHrC6XGvKOFREREREJoM9LkREREagzQTbh6+3RExciIiIjEAFASUTF41xqIiIiIhMBntciIiIjIBDRdph4kJERGQEXFWkHSYuRERERqD6/0OX6y0RExcya+F+wcYOgYiI9IiJCxERkREodVxVpMu1poyJCxERkREoBXR8O7T+YjElXA5NREREJoM9LkREREbAybnaYeJCRERkBCrIoIRMp+stEYeKiIiIyGSwx4WIiMgIVOL+ocv1loiJCxERkREodRwq0uVaU8ahIiIiIjIZ7HEhIiIyAva4aIeJCxERkRGohAwqocOqIh2uNWVMXIiIiIyAPS7a4RwXIiIiMhnscSEiIjICJayg1KH/QKnHWEwJExciAgDsuH7ykefD/YKrKRIiyyB0nOMiLHSOC4eKiIiIyGSwx4WIiMgIODlXO0btcQkKCoJMJit3REVFAQCKiooQFRWFWrVqwcXFBYMHD0ZWVpYxQyYiItILpbDS+bBERm31kSNHkJGRIR07d+4EALzwwgsAgEmTJuHXX3/Fhg0bsG/fPly/fh2DBg0yZshERERkREYdKqpTp47a93PmzEHDhg3RrVs35ObmYuXKlVi7di169OgBAIiPj0fz5s1x6NAhdOzY0RghExER6YUKMqh06D9QwTLfslhj+plKSkrw3Xff4bXXXoNMJkNycjJKS0sRFhYm1WnWrBkCAgKQmJhY6X2Ki4uRl5endhAREdU0ZXNcdDksUY1JXH766Sfk5OTg1VdfBQBkZmbCzs4O7u7uavW8vb2RmZlZ6X1iY2Ph5uYmHf7+/gaMmoiIiKpTjUlcVq5cid69e8PPz0+n+0ydOhW5ubnSkZ6erqcIiYiI9IeTc7VTI5ZDX716Fbt27cKmTZukMh8fH5SUlCAnJ0et1yUrKws+Pj6V3sve3h729vaGDJeIiEhn9+e46PCSRQ4VGU98fDy8vLzQt29fqSwkJAS2trZISEiQylJSUpCWlgaFQmGMMImIiPRG9f9b/mt76DKx15QZvcdFpVIhPj4eERERsLH5Xzhubm4YOXIkoqOj4enpCblcjnHjxkGhUHBFERERkYUyeuKya9cupKWl4bXXXit3bsGCBbCyssLgwYNRXFyM8PBwLFmyxAhREhER6Zeu81SUwjKXQxs9cenVqxdEJT98BwcHxMXFIS4urpqjIiIiMiyVjsM9lrqPi9ETFyKqGYz19me+lZqINMHEhYiIyAiUQgal0OElizpca8qYuBARERlB2eog7a+3zKEiy1xLRURERCaJPS5ERERGoBJWUOmwqkjFVUVERERUXThUpB0OFREREZHJYI8LERGREaig28oglf5CMSlMXIiIiIxA9w3oLHPQhIkLERGREei+5b9lJi6W2WoiIiIySexxISIiMgIVZFBBlzku3DmXiIiIqgmHirRjma0mIiIik8QeFyIiIiPQfQM6y+x7YOJCREYV7hds7BCIjEIlZFDpso+Lhb4d2jLTNSIiIjJJ7HEhIiIyApWOQ0XcgI6IiIiqje5vh7bMxMUyW01EREQmiT0uRERERqCEDEodNpHT5VpTxsSFiIjICDhUpB0mLkREREaghG69Jkr9hWJSLDNdIyIiIpPEHhciIiIj4FCRdpi4EBERGQFfsqgdy2w1ERGRhYqLi0NQUBAcHBwQGhqKw4cPP7L+woUL0bRpUzg6OsLf3x+TJk1CUVFRNUVbHhMXIiIiIxCQQaXDIbSY2Lt+/XpER0cjJiYGx44dQ3BwMMLDw5GdnV1h/bVr12LKlCmIiYnB2bNnsXLlSqxfvx7vv/++rs3XGhMXIiIiIygbKtLl0NT8+fMxatQoREZGokWLFli2bBmcnJywatWqCusfPHgQnTt3xrBhwxAUFIRevXph6NChj+2lMSQmLkRERCYsLy9P7SguLq6wXklJCZKTkxEWFiaVWVlZISwsDImJiRVe06lTJyQnJ0uJyqVLl7B161b06dNH/w2pIk7OJSIiMgKVkEEltN/Hpexaf39/tfKYmBjMnDmzXP2bN29CqVTC29tbrdzb2xvnzp2r8BnDhg3DzZs30aVLFwghcO/ePYwePdqoQ0VMXIiIiIxAqePbocuuTU9Ph1wul8rt7e11jq3M3r17MXv2bCxZsgShoaG4cOECJkyYgI8++gjTp0/X23M0wcSFiIjIhMnlcrXEpTK1a9eGtbU1srKy1MqzsrLg4+NT4TXTp0/HK6+8gtdffx0A0Lp1axQWFuKNN97ABx98ACur6p9xwjkuRERERlA2VKTLoQk7OzuEhIQgISHhfzGoVEhISIBCoajwmjt37pRLTqytrQEAQggNW6wf7HEhIiIyAhWsoNKh/0Cba6OjoxEREYH27dujQ4cOWLhwIQoLCxEZGQkAGDFiBOrWrYvY2FgAQL9+/TB//ny0a9dOGiqaPn06+vXrJyUw1Y2JCxERkREohQxKHSbnanPtkCFDcOPGDcyYMQOZmZlo27Yttm/fLk3YTUtLU+thmTZtGmQyGaZNm4Zr166hTp066NevHz755BOt49aVTBirr6ea5OXlwc3NDd3RHzYyW2OHQ0RENdg9UYq9+Bm5ublVmjeijbK/S2P2D4K9i/Z/l4oLSrG06yaDxloTsceFiIjICPS1HNrSMHEhIiIyAqHj26EFX7JIREREVLOxx4WIiMgIlJBBqcWLEh+83hIxcSEiIjICldBtnorKrJfWVI5DRURERGQy2ONCRERkBCodJ+fqcq0pY+JCRERkBCrIoNJhnoou15oyJi5ERERGYIydc82BZfYzERERkUlijwsREZERcI6Ldoze6mvXruHll19GrVq14OjoiNatW+Po0aPSeSEEZsyYAV9fXzg6OiIsLAznz583YsRERES6U0Embfuv1WGhc1yMmrj8+++/6Ny5M2xtbbFt2zb8/fffmDdvHjw8PKQ6n332GRYvXoxly5YhKSkJzs7OCA8PR1FRkREjJyIiImMw6lDRp59+Cn9/f8THx0tl9evXl74WQmDhwoWYNm0a+vfvDwD45ptv4O3tjZ9++gkvvfRStcdMRESkD0LHVUWCPS7V75dffkH79u3xwgsvwMvLC+3atcOKFSuk85cvX0ZmZibCwsKkMjc3N4SGhiIxMbHCexYXFyMvL0/tICIiqml0GibS8c3SpsyoiculS5ewdOlSNG7cGDt27MCYMWMwfvx4fP311wCAzMxMAIC3t7fadd7e3tK5h8XGxsLNzU06/P39DdsIIiIiqjZGHSpSqVRo3749Zs+eDQBo164dzpw5g2XLliEiIkKre06dOhXR0dHS93l5eUxeiIioxuGqIu0YtdW+vr5o0aKFWlnz5s2RlpYGAPDx8QEAZGVlqdXJysqSzj3M3t4ecrlc7SAiIqppOFSkHaMmLp07d0ZKSopaWWpqKgIDAwHcn6jr4+ODhIQE6XxeXh6SkpKgUCiqNVYiIiIyPqMOFU2aNAmdOnXC7Nmz8eKLL+Lw4cNYvnw5li9fDgCQyWSYOHEiPv74YzRu3Bj169fH9OnT4efnhwEDBhgzdCIiIp3wXUXaqVLiMmjQII1vvGzZMnh5eT2yzpNPPonNmzdj6tSp+PDDD1G/fn0sXLgQw4cPl+q8++67KCwsxBtvvIGcnBx06dIF27dvh4ODg8YxERER1RS6DvdY6lCRTAghHlfJysoKL774IhwdHat007Vr1+Ls2bNo0KCBzgHqKi8vD25ubuiO/rCR2Ro7HCIiqsHuiVLsxc/Izc012BzJsr9LvbePgq2zndb3KS0swbZnVxg01pqoykNFixcvfmwPSpkff/xR64CIqGbacf1kpefC/YKrMRKqyKN+PwB/R2Q+qpS47NmzB56enlW+6bZt21C3bl2tgyIiIjJ3HCrSTpUSl27duml00y5dumgVDBERkaVg4qIdrVYVqVQqXLhwAdnZ2VCpVGrnnnrqKb0ERkRERPQwjROXQ4cOYdiwYbh69Soentcrk8mgVCr1FhwREZG5EtBtSfNjV9aYKY0Tl9GjR6N9+/bYsmULfH19IZNZZlcVERGRLjhUpB2NE5fz58/jxx9/RKNGjQwRDxEREVGlNN7yPzQ0FBcuXDBELERERBaD7yrSTpV6XE6dOiV9PW7cOLz99tvIzMxE69atYWurvqlbmzZt9BshERGRGeJQkXaqlLi0bdsWMplMbTLua6+9Jn1ddo6Tc4mIiMiQqpS4XL582dBxEBERWRT2uGinSolLYGCg9PUff/yBTp06wcZG/dJ79+7h4MGDanWJiIioYkLIIHRIPnS51pRpPDn36aefxu3bt8uV5+bm4umnn9ZLUEREROZOBZnOhyXSOHEpm8vysFu3bsHZ2VkvQRERERFVpMr7uAwaNAjA/Ym4r776Kuzt7aVzSqUSp06dQqdOnfQfIRERkRniHBftVDlxcXNzA3C/x8XV1RWOjo7SOTs7O3Ts2BGjRo3Sf4REpGbH9ZOVngv3CzbYcw15b9JdTf39GOu/V1PAOS7aqXLiEh8fLy2H/uKLL+Di4mKwoIiIiIgqotEcFyEE1qxZg4yMDEPFQ0REZBG4c652NEpcrKys0LhxY9y6dctQ8RAREVmEsqEiXQ5LpPGqojlz5uCdd97BmTNnDBEPERERUaU0fjv0iBEjcOfOHQQHB8POzk5tki6ACvd4ISIiInVCx+EeS+1x0ThxWbhwoQHCICIisiwCwAOvANTqekukceISERFhiDiIiIiIHkvjxAW4v+HcTz/9hLNnzwIAWrZsieeeew7W1tZ6DY6IiMhcqSCDTIdt+y11y3+NE5cLFy6gT58+uHbtGpo2bQoAiI2Nhb+/P7Zs2YKGDRvqPUgiIiJzww3otKPxqqLx48ejYcOGSE9Px7Fjx3Ds2DGkpaWhfv36GD9+vCFiJCIiMjvcx0U7Gve47Nu3D4cOHYKnp6dUVqtWLcyZMwedO3fWa3BERERED9I4cbG3t0d+fn658oKCAtjZ2eklKCIiInMnhI6riix0WZHGQ0X/+c9/8MYbbyApKQlCCAghcOjQIYwePRrPPfecIWIkIiIyO9w5VzsaJy6LFy9Gw4YNoVAo4ODgAAcHB3Tu3BmNGjXCokWLDBEjEREREQAthorc3d3x888/4/z58zh37hwAoHnz5mjUqJHegyOi8sL9go0dAlGV8b/XynFVkXa02scFABo3bozGjRvrMxYiIiKLoRIyyHRIPriqqIqUSiVWr16NhIQEZGdnQ6VSqZ3fvXu33oIjIiIiepDGicuECROwevVq9O3bF61atYJMZpkZHxERkS64qkg7Gicu69atww8//IA+ffoYIh4iIiKLcD9x0WWOix6DMSEaryqys7PjRFwiIiIyCo0Tl7fffhuLFi2CsNRUj4iISA+4j4t2NB4q+vPPP7Fnzx5s27YNLVu2hK2trdr5TZs26S04IiIicyX+/9Dlekuk1T4uAwcONEQsREREFoP7uGhH48QlPj6+SvUOHDiA9u3bw97eXuOgiIiIiCqi8RyXqurduzeuXbtmqNsTERGZNqGHwwIZLHHh5F0iIqJH0HVirpZDRXFxcQgKCoKDgwNCQ0Nx+PDhR9bPyclBVFQUfH19YW9vjyZNmmDr1q1aPVsftN7yn4iIiEzL+vXrER0djWXLliE0NBQLFy5EeHg4UlJS4OXlVa5+SUkJnnnmGXh5eeHHH39E3bp1cfXqVbi7u1d/8P+PiQsREZERGGPn3Pnz52PUqFGIjIwEACxbtgxbtmzBqlWrMGXKlHL1V61ahdu3b+PgwYPSKuKgoCDtg9YDgw0VERERUeX0tY9LXl6e2lFcXFzh80pKSpCcnIywsDCpzMrKCmFhYUhMTKzwml9++QUKhQJRUVHw9vZGq1atMHv2bCiVSv3/QKrIYIkL32FERERkeP7+/nBzc5OO2NjYCuvdvHkTSqUS3t7eauXe3t7IzMys8JpLly7hxx9/hFKpxNatWzF9+nTMmzcPH3/8sd7bUVUGGyri5FwiIqJH0GGCrXQ9gPT0dMjlcqlYn9uQqFQqeHl5Yfny5bC2tkZISAiuXbuGuXPnIiYmRm/P0YTGicvdu3chhICTkxMA4OrVq9i8eTNatGiBXr16SfXy8/P1FyUREZGZ0dccF7lcrpa4VKZ27dqwtrZGVlaWWnlWVhZ8fHwqvMbX1xe2trawtraWypo3b47MzEyUlJTAzs5O+wZoSeOhov79++Obb74BcH+JVGhoKObNm4f+/ftj6dKlGt1r5syZkMlkakezZs2k80VFRYiKikKtWrXg4uKCwYMHl/uBExER0ePZ2dkhJCQECQkJUplKpUJCQgIUCkWF13Tu3BkXLlyASqWSylJTU+Hr62uUpAXQInE5duwYunbtCgD48ccf4e3tjatXr+Kbb77B4sWLNQ6gZcuWyMjIkI4///xTOjdp0iT8+uuv2LBhA/bt24fr169j0KBBGj+DiIioxjHCBnTR0dFYsWIFvv76a5w9exZjxoxBYWGhtMpoxIgRmDp1qlR/zJgxuH37NiZMmIDU1FRs2bIFs2fPRlRUlLat1pnGQ0V37tyBq6srAOD333/HoEGDYGVlhY4dO+Lq1auaB2BjU2EXVW5uLlauXIm1a9eiR48eAO6/bqB58+Y4dOgQOnbsqPGziIiIagpjvKtoyJAhuHHjBmbMmIHMzEy0bdsW27dvlybspqWlwcrqf30a/v7+2LFjByZNmoQ2bdqgbt26mDBhAt577z2t49aVxolLo0aN8NNPP2HgwIFSYwAgOzu7SmNsDzt//jz8/Pzg4OAAhUKB2NhYBAQEIDk5GaWlpWrLtpo1a4aAgAAkJiZWmrgUFxerLQXLy8vTOCYiIqJqYYR1LGPHjsXYsWMrPLd3795yZQqFAocOHTJwVFWn8VDRjBkzMHnyZAQFBaFDhw7SuNjvv/+Odu3aaXSv0NBQrF69Gtu3b8fSpUtx+fJldO3aFfn5+cjMzISdnV253fketWwLAGJjY9WWhfn7+2vaRCIiIqqhNO5xef7559GlSxdkZGQgODhYKu/ZsycGDhyo0b169+4tfd2mTRuEhoYiMDAQP/zwAxwdHTUNDQAwdepUREdHS9/n5eUxeSEiohrHGENF5kCrDeh8fHzg6uqKnTt34u7duwCAJ598Um1FkDbc3d3RpEkTXLhwAT4+PigpKUFOTo5anUct2wLur18vWxpW1SViRERE1Y5vh9aKxonLrVu30LNnTzRp0gR9+vRBRkYGAGDkyJF4++23dQqmoKAAFy9ehK+vL0JCQmBra6u2bCslJQVpaWmVLtsiIiIi86Zx4jJp0iTY2toiLS1N2oQOuD9Tefv27Rrda/Lkydi3bx+uXLmCgwcPYuDAgbC2tsbQoUPh5uaGkSNHIjo6Gnv27EFycjIiIyOhUCi4ooiIiMyATA+H5dF4jsvvv/+OHTt2oF69emrljRs31ng59D///IOhQ4fi1q1bqFOnDrp06YJDhw6hTp06AIAFCxbAysoKgwcPRnFxMcLDw7FkyRJNQyYiIqp5dB3usdChIo0Tl8LCQrWeljK3b9/W+P0I69ate+R5BwcHxMXFIS4uTqP7EhERkXnSeKioa9eu0pb/wP23QKtUKnz22Wd4+umn9RocERGR2eLkXK1o3OPy2WefoWfPnjh69ChKSkrw7rvv4q+//sLt27dx4MABQ8RIRESkkR3XT1Z6LtwvuNJz1UpPb4e2NBr3uLRq1Qqpqano3Lkz+vfvj8LCQgwaNAjHjx9Hw4YNDREjEREREQAtelwAwM3NDdOmTdN3LERERBZDiPuHLtdbIq02oNu/fz9efvlldOrUCdeuXQMAfPvtt2pvdiYiIqJH4BwXrWicuGzcuBHh4eFwdHTEsWPHpBca5ubmYvbs2XoPkIiIyCyVzXHR5bBAGicuH3/8MZYtW4YVK1bA1tZWKu/cuTOOHTum1+CIiIiIHqTxHJeUlBQ89dRT5crd3NzKvVeIiIiIKiYT9w9drrdEGve4+Pj44MKFC+XK//zzTzRo0EAvQREREZk9znHRisaJy6hRozBhwgQkJSVBJpPh+vXrWLNmDSZPnowxY8YYIkYiIiIiAFoMFU2ZMgUqlQo9e/bEnTt38NRTT8He3h6TJ0/GuHHjDBEjERGR+eEGdFrRKHFRKpU4cOAAoqKi8M477+DChQsoKChAixYt4OLiYqgYiYiIzA9fsqgVjRIXa2tr9OrVC2fPnoW7uztatGhhqLiIiIiIytFqy/9Lly4ZIhYiIiLLwcm5WtFqH5fJkyfjt99+Q0ZGBvLy8tQOIiIiqgImLlrReHJunz59AADPPfccZLL/TQwSQkAmk0GpVOovOiIiIi3UmDdAk95pnLjs2bPHEHEQERFZFq4q0orGiUu3bt0MEQcREZFFsaSdc+/cuYO0tDSUlJSolbdp00bje2mcuJw6darCcplMBgcHBwQEBMDe3l7jQIiIiCyKBSyHvnHjBiIjI7Ft27YKz2szvUTjxKVt27Zqc1seZmtriyFDhuDLL7+Eg4ODxgERERGReZg4cSJycnKQlJSE7t27Y/PmzcjKysLHH3+MefPmaXVPjVcVbd68GY0bN8by5ctx4sQJnDhxAsuXL0fTpk2xdu1arFy5Ert378a0adO0CoiIiIjMw+7duzF//ny0b98eVlZWCAwMxMsvv4zPPvsMsbGxWt1T4x6XTz75BIsWLUJ4eLhU1rp1a9SrVw/Tp0/H4cOH4ezsjLfffhuff/65VkERERGZOxl0nOOit0gMp7CwEF5eXgAADw8P3LhxA02aNEHr1q1x7Ngxre6pcY/L6dOnERgYWK48MDAQp0+fBnB/OCkjI0OrgIiIiMg8NG3aFCkpKQCA4OBgfPnll7h27RqWLVsGX19fre6pceLSrFkzzJkzR21mcGlpKebMmYNmzZoBAK5duwZvb2+tAiIiIrIIZcuhdTlquAkTJkgdGTExMdi2bRv8/f2xaNEizJ49W6t7ajxUFBcXh+eeew716tWTljGdPn0aSqUSv/32GwDg0qVLeOutt7QKiIiIyCJYwKqil19+Wfo6JCQEV69exblz5xAQEIDatWtrdU+NE5dOnTrh8uXLWLNmDVJTUwEAL7zwAoYNGwZXV1cAwCuvvKJVMERERGQ+oqOjq1x3/vz5VaqnceICAK6urhg9erQ2lxIRERFgET0ux48fx7Fjx3Dv3j00bdoUAJCamgpra2s88cQTUr1HbbPyMK0Sl2+//RZffvklLl26hMTERAQGBmLBggVo0KAB+vfvr80tiYiILIol7Jzbr18/uLq64uuvv4aHhwcA4N9//0VkZCS6du2Kt99+W+N7ajw5d+nSpYiOjkbv3r3x77//SrveeXh4YOHChRoHQEREROZp3rx5iI2NlZIW4H6+UK0b0H3xxRdYsWIFPvjgA9jY/K/Dpn379tJyaCIiInoMoYejhsvLy8ONGzfKld+4cQP5+fla3VPjoaLLly+jXbt25crt7e1RWFioVRBUfXZcP/nI83wVPGmD/10RacEC5rgMHDgQkZGRmDdvHjp06AAASEpKwjvvvINBgwZpdU+NE5f69evjxIkT5Tah2759O5o3b65VEERERJbGEua4LFu2DJMnT8awYcNQWloKALCxscHIkSMxd+5cre6pceISHR2NqKgoFBUVQQiBw4cP4/vvv0dsbCy++uorrYIgIiIi8+Pk5IQlS5Zg7ty5uHjxIgCgYcOGcHZ21vqeGicur7/+OhwdHTFt2jTcuXMHw4YNg5+fHxYtWoSXXnpJ60CIiIgsiq6735rAzrllnJ2dpU1rdaXVcujhw4dj+PDhuHPnDgoKCqQXKBEREVEVWcAcF0PQKnEp4+TkBCcnJ33FQkRERPRIVUpc2rVrV+Vd7bR9TTUREZElsYTJuYZQpcRlwIAB0tdFRUVYsmQJWrRoAYVCAQA4dOgQ/vrrL75YkYiIqKo4VKSVKiUuMTEx0tevv/46xo8fj48++qhcnfT0dP1GR0RERPQAjXfO3bBhA0aMGFGu/OWXX8bGjRv1EhQREZHZE/8bLtLmsNQeF40TF0dHRxw4cKBc+YEDB+Dg4KCXoIiIiMyeBWz5bwgaryqaOHEixowZg2PHjqlt37tq1SpMnz5d7wESERERldE4cZkyZQoaNGiARYsW4bvvvgMANG/eHPHx8XjxxRf1HiAREZFZ4uRcrWi1j8uLL77IJIWIiEgHXA6tHY3nuBAREREZS5V6XDw9PZGamoratWtX6aYBAQHYv39/uTdIk/GF+wUbOwQyQzX1v6sd109Wes6QMRvruUSWoEqJS05ODrZt2wY3N7cq3fTWrVtQKpU6BUZERGTWOMdFK1We4xIREWHIOIiIiCwK57hop0pzXFQqlcZHgwYNNA5mzpw5kMlkmDhxolRWVFSEqKgo1KpVCy4uLhg8eDCysrI0vjcRERGZvhozOffIkSP48ssv0aZNG7XySZMm4ddff8WGDRuwb98+XL9+HYMGDTJSlERERHrEzec0ViMSl4KCAgwfPhwrVqyAh4eHVJ6bm4uVK1di/vz56NGjB0JCQhAfH4+DBw/i0KFDRoyYiIhIR9w5Vys1InGJiopC3759ERYWplaenJyM0tJStfJmzZohICAAiYmJFd6ruLgYeXl5agcRERGZB602oNOndevW4dixYzhy5Ei5c5mZmbCzs4O7u7taube3NzIzMyu8X2xsLGbNmmWIUImIiPSGk3O1Y9Qel/T0dEyYMAFr1qzR2wsap06ditzcXOlIT0/Xy32JiIj0ikNFWtEqcbl48SKmTZuGoUOHIjs7GwCwbds2/PXXXxrdJzk5GdnZ2XjiiSdgY2MDGxsb7Nu3D4sXL4aNjQ28vb1RUlKCnJwcteuysrLg4+NT4T3t7e0hl8vVDiIiIjIPGicu+/btQ+vWrZGUlIRNmzahoKAAAHDy5EnExMRodK+ePXvi9OnTOHHihHS0b98ew4cPl762tbVFQkKCdE1KSgrS0tKgUCg0DZ2IiKjGKBsq0uWwRBonLlOmTMHHH3+MnTt3ws7OTirv0aOHxit9XF1d0apVK7XD2dkZtWrVQqtWreDm5oaRI0ciOjoae/bsQXJyMiIjI6FQKNCxY0dNQyciIqo5jDRUFBcXh6CgIDg4OCA0NBSHDx+u0nXr1q2DTCbDgAEDtHuwnmicuJw+fRoDBw4sV+7l5YWbN2/qJagHLViwAP/5z38wePBgPPXUU/Dx8cGmTZv0/hwiIqJqZYTEZf369YiOjkZMTAyOHTuG4OBghIeHS9M+KnPlyhVMnjwZXbt21fyheqZx4uLu7o6MjIxy5cePH0fdunV1Dmjv3r1YuHCh9L2DgwPi4uJw+/ZtFBYWYtOmTZXObyEiIqLKzZ8/H6NGjUJkZCRatGiBZcuWwcnJCatWrar0GqVSieHDh2PWrFla7YqvbxonLi+99BLee+89ZGZmQiaTQaVS4cCBA5g8eTJGjBhhiBiJiIjMjr7muDy8d1lxcXGFzyspKUFycrLa3mhWVlYICwurdG80APjwww/h5eWFkSNH6rX92tJ4H5fZs2cjKioK/v7+UCqVaNGiBZRKJYYNG4Zp06YZIkYiIq2E+wVb1HPJxOi6pPn/r/X391crjomJwcyZM8tVv3nzJpRKJby9vdXKvb29ce7cuQof8eeff2LlypU4ceKEDoHql8aJi52dHVasWIHp06fjzJkzKCgoQLt27dC4cWNDxEdERESPkJ6errb1h729vV7um5+fj1deeQUrVqxA7dq19XJPfdB659yAgAAEBAToMxYiIiLLoacel6ruWVa7dm1YW1sjKytLrbyyvdEuXryIK1euoF+/flKZSqUCANjY2CAlJQUNGzbUoQHaqVLiEh0dXeUbzp8/X+tgiIiILEV1b/lvZ2eHkJAQJCQkSEuaVSoVEhISMHbs2HL1mzVrhtOnT6uVTZs2Dfn5+Vi0aFG5IarqUqXE5fjx42rfHzt2DPfu3UPTpk0BAKmpqbC2tkZISIj+IyQiIiK9iI6ORkREBNq3b48OHTpg4cKFKCwsRGRkJABgxIgRqFu3LmJjY+Hg4IBWrVqpXV/27sCHy6tTlRKXPXv2SF/Pnz8frq6u+Prrr+Hh4QEA+PfffxEZGVkj1ncTERGZBD0NFWliyJAhuHHjBmbMmIHMzEy0bdsW27dvlybspqWlwcrKqK8xfCyZEEKjptetWxe///47WrZsqVZ+5swZ9OrVC9evX9drgLrKy8uDm5sbuqM/bGS2xg6HiIhqsHuiFHvxM3Jzcw32rruyv0vNx86Gtb32LxhWFhfh7H/fN2isNZHGaVVeXh5u3LhRrvzGjRvIz8/XS1BEREREFdE4cRk4cCAiIyOxadMm/PPPP/jnn3+wceNGjBw5EoMGDTJEjERERObHSO8qMnUaL4detmwZJk+ejGHDhqG0tPT+TWxsMHLkSMydO1fvARIREZklI8xxMQcaJy5OTk5YsmQJ5s6di4sXLwIAGjZsCGdnZ70HR0REZK5k/3/ocr0l0noDOmdnZ7Rp00afsRARERE9ksaJy9NPPw2ZrPI8b/fu3ToFREREZBE4VKQVjROXtm3bqn1fWlqKEydO4MyZM4iIiNBXXERERGatunfONRcaJy4LFiyosHzmzJkoKCjQOSAiIiKiyuhte7yXX34Zq1at0tftiIiIzBuXQ2tF68m5D0tMTISDg/Y7ABIREVkcC00+dKFx4vLwJnNCCGRkZODo0aOYPn263gIjIiIiepjGiYtcLldbVWRlZYWmTZviww8/RK9evfQaHBERkbni5FztaJy4rF692gBhEBERWRguh9aKxpNzGzRogFu3bpUrz8nJQYMGDfQSFBEREVFFNO5xuXLlCpRKZbny4uJiXLt2TS9BERERmTsOFWmnyonLL7/8In29Y8cOuLm5Sd8rlUokJCQgKChIr8ERERGZLQ4VaaXKicuAAQMAADKZrNwOuba2tggKCsK8efP0GhwREZG5Yo+LdqqcuKhUKgBA/fr1ceTIEdSuXdtgQRERERFVROM5LpcvXzZEHERERJaFQ0VaqVLisnjxYrzxxhtwcHDA4sWLH1l3/PjxegmMiIjIrDFx0UqVEpcFCxZg+PDhcHBwqPQli8D9+S9MXIiIiMhQqpS4PDg8xKEiIiIi3XFyrnY0nuPy4YcfYvLkyXByclIrv3v3LubOnYsZM2boLTgiIqKK7Lh+8pHnw/2CqykSHXCoSCsa75w7a9YsFBQUlCu/c+cOZs2apZegiIiIiCqicY+LEELtJYtlTp48CU9PT70ERUREZO5kQkAmtO820eVaU1blxMXDwwMymQwymQxNmjRRS16USiUKCgowevRogwRJRERkdjhUpJUqJy4LFy6EEAKvvfYaZs2apbblv52dHYKCgqBQKAwSJBERERGgQeJSts1//fr10alTJ9ja2hosKCIiInPHVUXa0XiOS7du3aSvi4qKUFJSonZeLpfrHhUREZG541CRVjReVXTnzh2MHTsWXl5ecHZ2hoeHh9pBREREj1fW46LLYYk0Tlzeeecd7N69G0uXLoW9vT2++uorzJo1C35+fvjmm28MESMRERERAC2Gin799Vd888036N69OyIjI9G1a1c0atQIgYGBWLNmDYYPH26IOImIiMwLh4q0onGPy+3bt9GgQQMA9+ez3L59GwDQpUsX/PHHH/qNjoiIyExxqEg7GicuDRo0kN5X1KxZM/zwww8A7vfEuLu76zU4IiIiogdpnLhERkbi5Mn774iYMmUK4uLi4ODggEmTJuGdd97Re4BERERmSejhsEAaz3GZNGmS9HVYWBjOnTuH5ORkNGrUCG3atNFrcERERObMUod7dKFxj8vDAgMDMWjQIHh6euKNN97QR0xEREREFdK4x6Uyt27dwsqVK7F8+XJ93ZKIiKhC4X7Bxg5Bd0LcP3S53gLpLXEhIiKiquOW/9rReaiIiIiIqLoYNXFZunQp2rRpA7lcDrlcDoVCgW3btknni4qKEBUVhVq1asHFxQWDBw9GVlaWESMmIiLSE64q0kqVh4oGDRr0yPM5OTkaP7xevXqYM2cOGjduDCEEvv76a/Tv3x/Hjx9Hy5YtMWnSJGzZsgUbNmyAm5sbxo4di0GDBuHAgQMaP4uIiKgmkanuH7pcb4mqnLi4ubk99vyIESM0eni/fv3Uvv/kk0+wdOlSHDp0CPXq1cPKlSuxdu1a9OjRAwAQHx+P5s2b49ChQ+jYsaNGzyIiIqpRuOW/VqqcuMTHxxsyDiiVSmzYsAGFhYVQKBRITk5GaWkpwsLCpDrNmjVDQEAAEhMTK01ciouLUVxcLH2fl5dn0LiJiIio+hh9cu7p06fh4uICe3t7jB49Gps3b0aLFi2QmZkJOzu7cq8R8Pb2RmZmZqX3i42NhZubm3T4+/sbuAVERESa47uKtGP0xKVp06Y4ceIEkpKSMGbMGERERODvv//W+n5Tp05Fbm6udKSnp+sxWiIiIj0p28dFl8MCGX0fFzs7OzRq1AgAEBISgiNHjmDRokUYMmQISkpKkJOTo9brkpWVBR8fn0rvZ29vD3t7e0OHTUREREZg9B6Xh6lUKhQXFyMkJAS2trZISEiQzqWkpCAtLQ0KhcKIERIREemOQ0XaMWqPy9SpU9G7d28EBAQgPz8fa9euxd69e7Fjxw64ublh5MiRiI6OhqenJ+RyOcaNGweFQsEVRUREZPq4qkgrRk1csrOzMWLECGRkZMDNzQ1t2rTBjh078MwzzwAAFixYACsrKwwePBjFxcUIDw/HkiVLjBkyERERGZFRE5eVK1c+8ryDgwPi4uIQFxdXTRERERFVD76rSDtGn5xLRERkkfh2aK0wcSEiIoPZcf1kpefC/YKrMRIyF0xciIiIjIBDRdph4kJERGQMXFWkFSYuRERERsAeF+3UuA3oiIiIiCrDHhciIiJjUIn7hy7XWyAmLkRERMbAOS5a4VARERERmQz2uBARERmBDDpOztVbJKaFiQsREZExcOdcrXCoiIiIiEwGExciIiIjKNvHRZdDG3FxcQgKCoKDgwNCQ0Nx+PDhSuuuWLECXbt2hYeHBzw8PBAWFvbI+tWBiQsREZExCD0cGlq/fj2io6MRExODY8eOITg4GOHh4cjOzq6w/t69ezF06FDs2bMHiYmJ8Pf3R69evXDt2jXNH64nTFyIiIgsxPz58zFq1ChERkaiRYsWWLZsGZycnLBq1aoK669ZswZvvfUW2rZti2bNmuGrr76CSqVCQkJCNUf+P0xciIiIjEAmhM4HAOTl5akdxcXFFT6vpKQEycnJCAsLk8qsrKwQFhaGxMTEKsV8584dlJaWwtPTU/cfgJa4qoiISM92XD9Z6blwv+BqjMT4LK29GlH9/6HL9QD8/f3VimNiYjBz5sxy1W/evAmlUglvb2+1cm9vb5w7d65Kj3zvvffg5+enlvxUNyYuRERERvBgr4m21wNAeno65HK5VG5vb69zbBWZM2cO1q1bh71798LBwcEgz6gKJi5EREQmTC6XqyUulalduzasra2RlZWlVp6VlQUfH59HXvv5559jzpw52LVrF9q0aaNTvLriHBciIiJjqOZVRXZ2dggJCVGbWFs20VahUFR63WeffYaPPvoI27dvR/v27TV7qAGwx4WIiMgYjLBzbnR0NCIiItC+fXt06NABCxcuRGFhISIjIwEAI0aMQN26dREbGwsA+PTTTzFjxgysXbsWQUFByMzMBAC4uLjAxcVF+9h1wMSFiIjIQgwZMgQ3btzAjBkzkJmZibZt22L79u3ShN20tDRYWf1vMGbp0qUoKSnB888/r3afyiYAVwcmLkREREagy+63ZddrY+zYsRg7dmyF5/bu3av2/ZUrV7R7iAExcSEiIjIGvmRRK5ycS0RERCaDPS5ERERGIFPdP3S53hIxcSEiIjIGDhVphUNFREREZDLY40JERGQMWmwiV+56C8TEhYiIyAj09a4iS8PEhYhIz/hG5JrtUW/vzstXwqNJNQXCOS5a4RwXIiIiMhnscSEiIjIGAUCXJc2W2eHCxIWIiMgYOMdFOxwqIiIiIpPBHhciIiJjENBxcq7eIjEpTFyIiIiMgauKtMKhIiIiIjIZ7HEhIiIyBhUAmY7XWyAmLkREREbAVUXaYeJCRERkDJzjohXOcSEiIiKTwR4XIiIiY2CPi1aYuBARERkDExetcKiIiIiITAZ7XPDoV5zz9fREROblUZ/r90QpgEvVEwiXQ2uFiQsREZERcDm0djhURERERCbDqIlLbGwsnnzySbi6usLLywsDBgxASkqKWp2ioiJERUWhVq1acHFxweDBg5GVlWWkiImIiPSkbHKuLocFMmrism/fPkRFReHQoUPYuXMnSktL0atXLxQWFkp1Jk2ahF9//RUbNmzAvn37cP36dQwaNMiIURMREemBSuh+WCCjznHZvn272verV6+Gl5cXkpOT8dRTTyE3NxcrV67E2rVr0aNHDwBAfHw8mjdvjkOHDqFjx47GCJuIiIiMpEbNccnNzQUAeHp6AgCSk5NRWlqKsLAwqU6zZs0QEBCAxMREo8RIRESkFxwq0kqNWVWkUqkwceJEdO7cGa1atQIAZGZmws7ODu7u7mp1vb29kZmZWeF9iouLUVxcLH2fl5dnsJiJiIi0p2vyYZmJS43pcYmKisKZM2ewbt06ne4TGxsLNzc36fD399dThERERHrEHhet1IjEZezYsfjtt9+wZ88e1KtXTyr38fFBSUkJcnJy1OpnZWXBx8enwntNnToVubm50pGenm7I0ImIiKgaGTVxEUJg7Nix2Lx5M3bv3o369eurnQ8JCYGtrS0SEhKkspSUFKSlpUGhUFR4T3t7e8jlcrWDiIioxuGqIq0YdY5LVFQU1q5di59//hmurq7SvBU3Nzc4OjrCzc0NI0eORHR0NDw9PSGXyzFu3DgoFAquKCIiItMmVPcPXa63QEZNXJYuXQoA6N69u1p5fHw8Xn31VQDAggULYGVlhcGDB6O4uBjh4eFYsmRJNUdKRERENYFRExdRhYlFDg4OiIuLQ1xcXDVEREREVE10nWBroZNza8xyaCIiIouiEtBpSTPnuFiuR73inIiIiGoOJi5ERETGwKEirTBxISIiMgYBHRMXvUViUmrEBnREREREVcEeFyIiImPgUJFWmLgQEREZg0oFQIdN5FTcgI6IiIiqC3tctMI5LkRERGQy2ONCRERkDOxx0QoTFyIiImPgzrla4VARERERmQz2uBARERmBECoIof3KIF2uNWVMXIiIiIxBCN2Geyx0jguHioiIiMhksMeFiIjIGISOk3MttMeFiQsREZExqFSATId5KhY6x4VDRURERGQy2ONCRERkDBwq0goTFyIiIiMQKhWEDkNFXA5NRERE1Yc9LlrhHBciIiIyGexxISIiMgaVAGTscdEUExciIiJjEAKALsuhLTNx4VARERERmQz2uBARERmBUAkIHYaKBHtciIiIqNoIle6HFuLi4hAUFAQHBweEhobi8OHDj6y/YcMGNGvWDA4ODmjdujW2bt2q1XP1hYkLERGRhVi/fj2io6MRExODY8eOITg4GOHh4cjOzq6w/sGDBzF06FCMHDkSx48fx4ABAzBgwACcOXOmmiP/H5kw876mvLw8uLm5oTv6w0Zma+xwiIioBrsnSrEXPyM3Nxdyudwgz5D+LskG6vR36Z4oxV6xWaNYQ0ND8eSTT+K///0vAEClUsHf3x/jxo3DlClTytUfMmQICgsL8dtvv0llHTt2RNu2bbFs2TKtY9cFe1yIiIiMoZqHikpKSpCcnIywsDCpzMrKCmFhYUhMTKzwmsTERLX6ABAeHl5p/epg9pNzyzqU7qFUpw0KiYjI/N1DKYDqmfiq69+lsljz8vLUyu3t7WFvb1+u/s2bN6FUKuHt7a1W7u3tjXPnzlX4jMzMzArrZ2Zmah+4jsw+ccnPzwcA/AnjTiYiIiLTkZ+fDzc3N4Pc287ODj4+PvgzU/e/Sy4uLvD391cri4mJwcyZM3W+d01l9omLn58f0tPT4erqCplMhry8PPj7+yM9Pd1g45c1CdtrviyprQDba+5qSnuFEMjPz4efn5/BnuHg4IDLly+jpKRE53sJISCTydTKKuptAYDatWvD2toaWVlZauVZWVnw8fGp8BofHx+N6lcHs09crKysUK9evXLlcrncIj4MyrC95suS2gqwveauJrTXUD0tD3JwcICDg4PBn/MgOzs7hISEICEhAQMGDABwf3JuQkICxo4dW+E1CoUCCQkJmDhxolS2c+dOKBSKaoi4YmafuBAREdF90dHRiIiIQPv27dGhQwcsXLgQhYWFiIyMBACMGDECdevWRWxsLABgwoQJ6NatG+bNm4e+ffti3bp1OHr0KJYvX260NjBxISIishBDhgzBjRs3MGPGDGRmZqJt27bYvn27NAE3LS0NVlb/W3DcqVMnrF27FtOmTcP777+Pxo0b46effkKrVq2M1QTLS1zs7e0RExNT6RiguWF7zZcltRVge82dpbXXmMaOHVvp0NDevXvLlb3wwgt44YUXDBxV1Zn9BnRERERkPrgBHREREZkMJi5ERERkMpi4EBERkclg4kJEREQmwywSl9jYWDz55JNwdXWFl5cXBgwYgJSUFLU6RUVFiIqKQq1ateDi4oLBgweX2w0wLS0Nffv2hZOTE7y8vPDOO+/g3r171dmUKnlce2/fvo1x48ahadOmcHR0REBAAMaPH4/c3Fy1+5hCe6vyuy0jhEDv3r0hk8nw008/qZ0zhbYCVW9vYmIievToAWdnZ8jlcjz11FO4e/eudP727dsYPnw45HI53N3dMXLkSBQUFFRnU6qkKu3NzMzEK6+8Ah8fHzg7O+OJJ57Axo0b1eqYSnuXLl2KNm3aSJusKRQKbNu2TTpvTp9TwKPba06fU1TNhBkIDw8X8fHx4syZM+LEiROiT58+IiAgQBQUFEh1Ro8eLfz9/UVCQoI4evSo6Nixo+jUqZN0/t69e6JVq1YiLCxMHD9+XGzdulXUrl1bTJ061RhNeqTHtff06dNi0KBB4pdffhEXLlwQCQkJonHjxmLw4MHSPUylvVX53ZaZP3++6N27twAgNm/eLJWbSluFqFp7Dx48KORyuYiNjRVnzpwR586dE+vXrxdFRUVSnWeffVYEBweLQ4cOif3794tGjRqJoUOHGqNJj1SV9j7zzDPiySefFElJSeLixYvio48+ElZWVuLYsWNSHVNp7y+//CK2bNkiUlNTRUpKinj//feFra2tOHPmjBDCvD6nhHh0e83pc4qql1kkLg/Lzs4WAMS+ffuEEELk5OQIW1tbsWHDBqnO2bNnBQCRmJgohBBi69atwsrKSmRmZkp1li5dKuRyuSguLq7eBmjo4fZW5IcffhB2dnaitLRUCGG67a2srcePHxd169YVGRkZ5RIXU22rEBW3NzQ0VEybNq3Sa/7++28BQBw5ckQq27Ztm5DJZOLatWsGjVdXFbXX2dlZfPPNN2r1PD09xYoVK4QQpt1eIYTw8PAQX331ldl/TpUpa29FzOVzigzLLIaKHlbW1ejp6QkASE5ORmlpKcLCwqQ6zZo1Q0BAABITEwHc73pv3bq12uu7w8PDkZeXh7/++qsao9fcw+2trI5cLoeNzf09B021vRW19c6dOxg2bBji4uIqfPGXqbYVKN/e7OxsJCUlwcvLC506dYK3tze6deuGP//8U7omMTER7u7uaN++vVQWFhYGKysrJCUlVW8DNFTR77dTp05Yv349bt++DZVKhXXr1qGoqAjdu3cHYLrtVSqVWLduHQoLC6FQKMz+c+rh9lbEXD6nyLDMbudclUqFiRMnonPnztKWxJmZmbCzs4O7u7taXW9vb2RmZkp1Hvyfo+x82bmaqqL2PuzmzZv46KOP8MYbb0hlptjeyto6adIkdOrUCf3796/wOlNsK1Bxey9dugQAmDlzJj7//HO0bdsW33zzDXr27IkzZ86gcePGyMzMhJeXl9q9bGxs4OnpaXLtBYAffvgBQ4YMQa1atWBjYwMnJyds3rwZjRo1AgCTa+/p06ehUChQVFQEFxcXbN68GS1atMCJEyfM8nOqsvY+zFw+p8jwzC5xiYqKwpkzZ9T+BWrOHtfevLw89O3bFy1atMDMmTOrNzg9q6itv/zyC3bv3o3jx48bMTLDqKi9KpUKAPDmm29KL0Vr164dEhISsGrVKunFaKaosv+Wp0+fjpycHOzatQu1a9fGTz/9hBdffBH79+9H69atjRSt9po2bYoTJ04gNzcXP/74IyIiIrBv3z5jh2UwlbX3weTFnD6nyPDMaqho7Nix+O2337Bnzx7Uq1dPKvfx8UFJSQlycnLU6mdlZUlDCz4+PuVm75d9X9HwQ01QWXvL5Ofn49lnn4Wrqys2b94MW1tb6Zyptbeytu7evRsXL16Eu7s7bGxspC7mwYMHS0MJptZWoPL2+vr6AkC5f7E2b94caWlpAO63KTs7W+38vXv3cPv2bZNr78WLF/Hf//4Xq1atQs+ePREcHIyYmBi0b98ecXFxAEyvvXZ2dmjUqBFCQkIQGxuL4OBgLFq0yGw/pyprbxlz+pyi6mEWiYsQAmPHjsXmzZuxe/du1K9fX+18SEgIbG1tkZCQIJWlpKQgLS1NGmtVKBQ4ffq02gfgzp07IZfLK+zWNKbHtRe4/y+YXr16wc7ODr/88gscHBzUzptKex/X1ilTpuDUqVM4ceKEdADAggULEB8fD8B02go8vr1BQUHw8/Mrt2Q4NTUVgYGBAO63NycnB8nJydL53bt3Q6VSITQ01PCN0MDj2nvnzh0AUHtbLQBYW1tLvU+m1N6KqFQqFBcXm93nVGXK2guYz+cUVTOjTg3WkzFjxgg3Nzexd+9ekZGRIR137tyR6owePVoEBASI3bt3i6NHjwqFQiEUCoV0vmzZXa9evcSJEyfE9u3bRZ06dWrksrvHtTc3N1eEhoaK1q1biwsXLqjVuXfvnhDCdNpbld/tw1DJcuia3lYhqtbeBQsWCLlcLjZs2CDOnz8vpk2bJhwcHMSFCxekOs8++6xo166dSEpKEn/++ado3LhxjVwe/Lj2lpSUiEaNGomuXbuKpKQkceHCBfH5558LmUwmtmzZIt3HVNo7ZcoUsW/fPnH58mVx6tQpMWXKFCGTycTvv/8uhDCvzykhHt1ec/qcouplFokLgAqP+Ph4qc7du3fFW2+9JTw8PISTk5MYOHCgyMjIULvPlStXRO/evYWjo6OoXbu2ePvtt6VleTXJ49q7Z8+eSutcvnxZuo8ptLcqv9uKrnkwcRHCNNoqRNXbGxsbK+rVqyecnJyEQqEQ+/fvVzt/69YtMXToUOHi4iLkcrmIjIwU+fn51diSqqlKe1NTU8WgQYOEl5eXcHJyEm3atCm3PNpU2vvaa6+JwMBAYWdnJ+rUqSN69uwpJS1CmNfnlBCPbq85fU5R9ZIJIYRh+nKIiIiI9Mss5rgQERGRZWDiQkRERCaDiQsRERGZDCYuREREZDKYuBAREZHJYOJCREREJoOJCxEREZkMJi5Ej3HlyhXIZDLpdQL6JpPJ8NNPP2l9/d69eyGTySCTyTBgwIBH1u3evTsmTpyo9bPo0cp+Dw+/4ZmI9IeJC9Vor7766mP/GBuav78/MjIy0KpVKwD/SxQefhmesaWkpGD16tXGDsMiVPbfZUZGBhYuXFjt8RBZEiYuRI9hbW0NHx8f6c3TNZWXl1eN+Jd+SUmJsUMwGh8fH7i5uRk7DCKzxsSFTNq+ffvQoUMH2Nvbw9fXF1OmTMG9e/ek8927d8f48ePx7rvvwtPTEz4+Ppg5c6baPc6dO4cuXbrAwcEBLVq0wK5du9SGbx4cKrpy5QqefvppAICHhwdkMhleffVVAPff3Pzwv7bbtm2r9rzz58/jqaeekp61c+fOcm1KT0/Hiy++CHd3d3h6eqJ///64cuWKxj+bwsJCjBgxAi4uLvD19cW8efPK1SkuLsbkyZNRt25dODs7IzQ0FHv37lWrs2LFCvj7+8PJyQkDBw7E/Pnz1RKkmTNnom3btvjqq69Qv3596Q2/OTk5eP3111GnTh3I5XL06NEDJ0+eVLv3zz//jCeeeAIODg5o0KABZs2aJf3+hBCYOXMmAgICYG9vDz8/P4wfP75KbX9cu27duoWhQ4eibt26cHJyQuvWrfH999+r3ePHH39E69at4ejoiFq1aiEsLAyFhYWYOXMmvv76a/z888/S0NDDPzMiMpya/U9Ioke4du0a+vTpg1dffRXffPMNzp07h1GjRsHBwUEtWfj6668RHR2NpKQkJCYm4tVXX0Xnzp3xzDPPQKlUYsCAAQgICEBSUhLy8/Px9ttvV/pMf39/bNy4EYMHD0ZKSgrkcjkcHR2rFK9KpcKgQYPg7e2NpKQk5ObmlptvUlpaivDwcCgUCuzfvx82Njb4+OOP8eyzz+LUqVOws7Or8s/nnXfewb59+/Dzzz/Dy8sL77//Po4dO4a2bdtKdcaOHYu///4b69atg5+fHzZv3oxnn30Wp0+fRuPGjXHgwAGMHj0an376KZ577jns2rUL06dPL/esCxcuYOPGjdi0aROsra0BAC+88AIcHR2xbds2uLm54csvv0TPnj2RmpoKT09P7N+/HyNGjMDixYvRtWtXXLx4EW+88QYAICYmBhs3bsSCBQuwbt06tGzZEpmZmeUSn8o8rl1FRUUICQnBe++9B7lcji1btuCVV15Bw4YN0aFDB2RkZGDo0KH47LPPMHDgQOTn52P//v0QQmDy5Mk4e/Ys8vLyEB8fDwDw9PSs8u+FiHRk3Hc8Ej1aRESE6N+/f4Xn3n//fdG0aVOhUqmksri4OOHi4iKUSqUQQohu3bqJLl26qF335JNPivfee08IIcS2bduEjY2N2ht4d+7cqfaG6cuXLwsA4vjx40KI/73V9t9//1W7b2BgoFiwYIFaWXBwsIiJiRFCCLFjxw5hY2Mjrl27Jp3ftm2b2rO+/fbbcm0qLi4Wjo6OYseOHRX+HCqKJz8/X9jZ2YkffvhBKrt165ZwdHQUEyZMEEIIcfXqVWFtba0WjxBC9OzZU0ydOlUIIcSQIUNE37591c4PHz5cuLm5Sd/HxMQIW1tbkZ2dLZXt379fyOVyUVRUpHZtw4YNxZdffik9Z/bs2Wrnv/32W+Hr6yuEEGLevHmiSZMmoqSkpMJ2V6Yq7apI3759xdtvvy2EECI5OVkAEFeuXKmw7qP+u4yPj1f7+RCRfrHHhUzW2bNnoVAoIJPJpLLOnTujoKAA//zzDwICAgAAbdq0UbvO19cX2dnZAO5PaPX394ePj490vkOHDgaL19/fH35+flKZQqFQq3Py5ElcuHABrq6uauVFRUW4ePFilZ918eJFlJSUIDQ0VCrz9PRE06ZNpe9Pnz4NpVKJJk2aqF1bXFyMWrVqAbj/8xk4cKDa+Q4dOuC3335TKwsMDESdOnXU2lFQUCDdp8zdu3eldpw8eRIHDhzAJ598Ip1XKpUoKirCnTt38MILL2DhwoVo0KABnn32WfTp0wf9+vV77FyjqrRLqVRi9uzZ+OGHH3Dt2jWUlJSguLgYTk5OAIDg4GD07NkTrVu3Rnh4OHr16oXnn38eHh4ej3w2ERkeExcye7a2tmrfy2QyqFQqvT/HysoKQgi1stLSUo3uUVBQgJCQEKxZs6bcuQcTA30oKCiAtbU1kpOTpeGdMi4uLhrdy9nZudy9fX19K5z7UTY/pqCgALNmzcKgQYPK1XFwcIC/vz9SUlKwa9cu7Ny5E2+99Rbmzp2Lffv2lfudatquuXPnYtGiRVi4cCFat24NZ2dnTJw4UZpYbG1tjZ07d+LgwYP4/fff8cUXX+CDDz5AUlIS6tevr8mPhoj0jIkLmazmzZtj48aNEEJIvS4HDhyAq6sr6tWrV6V7NG3aFOnp6cjKyoK3tzcA4MiRI4+8pmyeiVKpVCuvU6cOMjIypO/z8vJw+fJltXjT09ORkZEBX19fAMChQ4fU7vHEE09g/fr18PLyglwur1IbKtKwYUPY2toiKSlJ6nn6999/kZqaim7dugEA2rVrB6VSiezsbHTt2rXC+zRt2rTcz+NxP5+ydmRmZsLGxgZBQUGV1klJSUGjRo0qvY+joyP69euHfv36ISoqCs2aNcPp06fxxBNPVHpNVdp14MAB9O/fHy+//DKA+/OPUlNT0aJFC6mOTCZD586d0blzZ8yYMQOBgYHYvHkzoqOjYWdnV+73T0TVg6uKqMbLzc3FiRMn1I709HS89dZbSE9Px7hx43Du3Dn8/PPPiImJQXR0NKysqvaf9jPPPIOGDRsiIiICp06dwoEDBzBt2jQAUBuCelBgYCBkMhl+++033LhxAwUFBQCAHj164Ntvv8X+/ftx+vRpREREqP2LPywsDE2aNEFERAROnjyJ/fv344MPPlC79/Dhw1G7dm30798f+/fvx+XLl7F3716MHz8e//zzT5V/Zi4uLhg5ciTeeecd7N69G2fOnMGrr76q9nNp0qQJhg8fjhEjRmDTpk24fPkyDh8+jNjYWGzZsgUAMG7cOGzduhXz58/H+fPn8eWXX2Lbtm2V/mwebKtCocCAAQPw+++/48qVKzh48CA++OADHD16FAAwY8YMfPPNN5g1axb++usvnD17FuvWrZN+/qtXr8bKlStx5swZXLp0Cd999x0cHR0RGBj4yGdXpV2NGzeWelTOnj2LN998E1lZWdI9kpKSMHv2bBw9ehRpaWnYtGkTbty4gebNmwO4v4Ls1KlTSElJwc2bNzXuWSMiHRh7kg3Ro0RERAgA5Y6RI0cKIYTYu3evePLJJ4WdnZ3w8fER7733nigtLZWu79atmzQZtUz//v1FRESE9P3Zs2dF586dhZ2dnWjWrJn49ddfBQCxfft2IUT5yblCCPHhhx8KHx8fIZPJpHvl5uaKIUOGCLlcLvz9/cXq1avVJucKIURKSoro0qWLsLOzE02aNBHbt29Xm5wrhBAZGRlixIgRonbt2sLe3l40aNBAjBo1SuTm5lb4M6pssnB+fr54+eWXhZOTk/D29hafffZZuZ9HSUmJmDFjhggKChK2trbC19dXDBw4UJw6dUqqs3z5clG3bl3h6OgoBgwYID7++GPh4+MjnY+JiRHBwcHl4srLyxPjxo0Tfn5+wtbWVvj7+4vhw4eLtLQ0qc727dtFp06dhKOjo5DL5aJDhw5i+fLlQgghNm/eLEJDQ4VcLhfOzs6iY8eOYteuXRX+DB72uHbdunVL9O/fX7i4uAgvLy8xbdo0MWLECGnC7d9//y3Cw8NFnTp1hL29vWjSpIn44osvpPtnZ2eLZ555Rri4uAgAYs+ePdI5Ts4lMiyZEA8NyhNZuAMHDqBLly64cOECGjZsaOxwHmvv3r14+umn8e+//1bLBnSjRo3CuXPnsH//foM/yxStXr0aEydOrHE7KxOZC85xIYu3efNmuLi4oHHjxrhw4QImTJiAzp07m0TS8qB69eqhX79+5TZS09Xnn3+OZ555Bs7Ozti2bRu+/vprLFmyRK/PMBcuLi64d++etAkfEekfExeyePn5+XjvvfeQlpaG2rVrIywsrMJdZmuq0NBQnD9/HoDmq4Gq4vDhw/jss8+Qn5+PBg0aYPHixXj99df1/pyq2r9/P3r37l3p+bI5R8ZQ9iLOh1czEZH+cKiIiEzK3bt3ce3atUrPP2qVEhGZPiYuREREZDK4HJqIiIhMBhMXIiIiMhlMXIiIiMhkMHEhIiIik8HEhYiIiEwGExciIiIyGUxciIiIyGQwcSEiIiKT8X9UWe8Vk5FLPwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ds['qa'] = (('time', 'lat', 'lon'), np.random.rand(*ds.air.shape) < 0.01)\n", + "ds['qa'].isel(time=0).plot()" + ] + }, + { + "cell_type": "markdown", + "id": "a9e888d1-4c4d-4ec7-aa47-8ec7c17c3b63", + "metadata": {}, + "source": [ + "Define a small function to determine which patches to keep. The function should take the underlying dataset as its first argument, and a dictionary of slice objects as the second argument. Each dictionary corresponds to one batch from the `BatchGenerator`. Batches for which the function returns True are retained in the `BatchGenerator`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "000a8888-aa69-4bc0-b932-a7879f43fa85", + "metadata": {}, + "outputs": [], + "source": [ + "def myfilter(ds, batch):\n", + " return (ds.isel(**batch).qa == 0).all()" + ] + }, + { + "cell_type": "markdown", + "id": "e615808d-5c5b-46c5-955e-bbe95f01aa35", + "metadata": {}, + "source": [ + "Now we pass the filter function to the batch generator and verify that none of the anomalous pixels make it into resulting batches." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "9afa1d84-8991-4572-a15e-f67b278b169c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Original generator: 29200 batches\n", + "Filtered generator: 8344 batches\n" + ] + } + ], + "source": [ + "bgen_original = xb.BatchGenerator(ds, {'lat': 5, 'lon': 5, 'time': 5})\n", + "\n", + "bgen_filtered = xb.BatchGenerator(\n", + " ds, {'lat': 5, 'lon': 5, 'time': 5}, filter_fn=myfilter\n", + ")\n", + "\n", + "print('Original generator:', len(bgen_original), 'batches')\n", + "print('Filtered generator:', len(bgen_filtered), 'batches')\n", + "\n", + "for batch in bgen_filtered:\n", + " assert (batch.qa == 0).all()" + ] + }, + { + "cell_type": "markdown", + "id": "32b49b74-3e46-435e-84d6-0c5b0e8610cf", + "metadata": {}, + "source": [ + "### Resampling\n", + "\n", + "Now we show how to resample a `BatchGenerator`. Note that this approach only supports *undersampling*. That is, you can only remove batches, not duplicate them. Now our task is to define a function with the same signature as the filter, but this time returning a non-negative float that indicates the relative sample weight of this patch. This functionality uses `np.random.choice` to select batches, so use `np.random.seed` to ensure reproducibility.\n", + "\n", + "Suppose we want to sample the dataset to emphasize batches with low air temperature. One option is to return a higher sample weight for batches with mean air temperature below a certain threshold." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "26ffc331-558b-49fb-a7fe-2da91bf8276b", + "metadata": {}, + "outputs": [], + "source": [ + "threshold = 270\n", + "\n", + "\n", + "def myresample(ds, batch):\n", + " window_mean = ds.isel(**batch).air.mean()\n", + " if window_mean > threshold:\n", + " return 1\n", + " else:\n", + " return 4\n", + "\n", + "\n", + "bgen_resampled = xb.BatchGenerator(\n", + " ds, {'lat': 5, 'lon': 5, 'time': 5}, resample_fn=myresample, resample_n=5_000\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6c6b1ec5-8e63-47b0-b8c3-82e188e0109f", + "metadata": {}, + "source": [ + "Now, we can compare the distribution of batch mean air temperature between the resampled and original generator." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e3ef2917-080a-40f3-8b47-97794b112cfa", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHVCAYAAAB8NLYkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQDJJREFUeJzt3Xl0VFW+9vGnMpOESghkIE0YFFqI4AAoRBQHuAQM0miUi4CAIioEFGhR0hdBaRSlW7ARBAcEWuWq0OIAAsYgOBAGQVoGRRExaSEJNiYlQypD7fcPXs61mqhBU6lw6vtZ66x1ap9dVb/aVMKTfSaHMcYIAAAAZ70gfxcAAACA2kGwAwAAsAmCHQAAgE0Q7AAAAGyCYAcAAGATBDsAAACbINgBAADYBMEOAADAJgh2AAAANkGwA4D/b/HixXI4HDpw4MAZP3f9+vVyOBxav359rdf1Yw6HQw8++KBP3wPA2YtgBwAAYBMh/i4AAOqLW265RQMHDlR4ePgZP7d79+46ceKEwsLCfFAZANQMM3YAAt6xY8ckScHBwYqIiJDD4Tjj1wgKClJERISCgvi1CsB/+A0EwFY++eQT9enTR06nU9HR0erRo4c2bdpkbT91HN2GDRs0evRoJSQkqFmzZl7bfnyMncfj0YMPPqjk5GRFRkbq6quv1p49e9SyZUsNHz7c6lfdMXZXXXWV2rdvrz179ujqq69WZGSkfve732nmzJleNZeXl2vKlCnq1KmTYmJiFBUVpSuuuELvvfeeT8YIgH2xKxaAbezevVtXXHGFnE6n7rvvPoWGhurpp5/WVVddpQ0bNqhLly5W39GjRys+Pl5TpkyxZuyqk52drZkzZ+q6665Tenq6/vnPfyo9PV1lZWU1qun7779X7969dcMNN2jAgAFavny57r//fnXo0EF9+vSRJLlcLj333HO6+eabNXLkSP3www9auHCh0tPTtWXLFl100UW/aVwABA6CHQDbmDx5sioqKvThhx/qnHPOkSQNHTpU5513nu677z5t2LDB6hsXF6fc3FwFBwf/5OsVFRVp1qxZ6t+/v1asWGG1P/TQQzU+M/XgwYP6+9//rltuuUWSNGLECLVo0UILFy60gl2jRo104MABr+PzRo4cqbZt2+rJJ5/UwoULazwGAAIbu2IB2EJVVZXeeecd9e/f3wp1ktS0aVMNGjRIH374oVwul9U+cuTInw11kpSbm6vKykqNHj3aq33s2LE1ris6OlpDhgyxHoeFhenSSy/V/v37rbbg4GAr1Hk8Hh05ckSVlZXq3Lmztm/fXuP3AgCCHQBbOHz4sI4fP67zzjvvtG3t2rWTx+NRQUGB1daqVatffM1vvvlGktS6dWuv9ri4ODVq1KhGdTVr1uy0kzEaNWqk77//3qttyZIluuCCCxQREaHGjRsrPj5eq1atUmlpaY3eBwAkgh2AANWgQYM6eZ+fmhU0xljrL774ooYPH65zzz1XCxcu1Jo1a5STk6NrrrlGHo+nTuoEYA8cYwfAFuLj4xUZGam9e/eetu3zzz9XUFCQUlJStHXr1hq/ZosWLSRJ+/bt85rh+/e//33ajNtvsXz5cp1zzjl67bXXvGb3pk6dWmvvASAwMGMHwBaCg4PVq1cvvfHGG16XKykqKtLSpUt1+eWXy+l0ntFr9ujRQyEhIZo/f75X+9y5c2ujZMupWb0fz+Jt3rxZeXl5tfo+AOyPGTsAtjF9+nTl5OTo8ssv1+jRoxUSEqKnn35abrf7tGvH1URiYqLuuecePf744+rXr5969+6tf/7zn1q9erWaNGnyqy5kXJ2+ffvqtdde0/XXX6+MjAx9/fXXWrBggVJTU3X06NFaeQ8AgYFgB8A2zj//fH3wwQfKzs7WjBkz5PF41KVLF7344ote17A7E4899pgiIyP17LPP6t1331VaWpreeecdXX755YqIiKiVuocPH67CwkI9/fTTWrt2rVJTU/Xiiy9q2bJlXhc8BoBf4jA/nvsHAPyikpISNWrUSNOnT9f//M//+LscALBwjB0A/IwTJ06c1vbEE09IOnnLMACoT9gVCwA/45VXXtHixYt17bXXKjo6Wh9++KH+93//V7169VK3bt38XR4AeCHYAcDPuOCCCxQSEqKZM2fK5XJZJ1RMnz7d36UBwGk4xg4AAMAmOMYOAADAJmy7K9bj8ejgwYNq2LBhrV1rCgAAoK4ZY/TDDz8oOTlZQUE/Pydn22B38OBBpaSk+LsMAACAWlFQUKBmzZr9bB/bBruGDRtKOjkIZ3obIQD4rY4dO6bk5GRJJ//QjIqK8nNFAM5WLpdLKSkpVrb5ObYNdqd2vzqdToIdgDoXGRmpV199VZIUHx+vkBDb/roFUEdqcmgZv2kAwAdCQkJ00003+bsMAAGGs2IBAABsghk7APCByspKrVixQpJ0/fXXsysWQJ3gNw0A+IDb7daAAQMkSUePHiXYAagT7IoFAACwCYIdAACATbBvAACAs0DLSav8XUKtOfBohr9LsC1m7AAAAGyCYAcAAGATBDsAAACb4Bg7APCBsLAwLVq0yFoHgLpAsAMAHwgNDdXw4cP9XQaAAMOuWAAAAJtgxg4AfKCyslJr166VJKWnp3PnCQB1gt80AOADbrdbffv2lcQtxQDUHXbFAgAA2ATBDgAAwCbOONh9++23GjJkiBo3bqwGDRqoQ4cO+vjjj63txhhNmTJFTZs2VYMGDdSzZ099+eWXXq9x5MgRDR48WE6nU7GxsRoxYoSOHj3q1efTTz/VFVdcoYiICKWkpGjmzJm/8iMCAAAEhjMKdt9//726deum0NBQrV69Wnv27NHjjz+uRo0aWX1mzpypOXPmaMGCBdq8ebOioqKUnp6usrIyq8/gwYO1e/du5eTkaOXKlXr//fd1xx13WNtdLpd69eqlFi1aaNu2bfrLX/6iBx98UM8880wtfGQAAAB7chhjTE07T5o0SR999JE++OCDarcbY5ScnKw//vGPuvfeeyVJpaWlSkxM1OLFizVw4EB99tlnSk1N1datW9W5c2dJ0po1a3TttdfqX//6l5KTkzV//nz9z//8jwoLC60Le06aNEmvv/66Pv/88xrV6nK5FBMTo9LSUjmdzpp+RACoFceOHVN0dLSkkydPREVF+bkinO1aTlrl7xJqzYFHM/xdwlnlTDLNGc3Yvfnmm+rcubNuuukmJSQk6OKLL9azzz5rbf/6669VWFionj17Wm0xMTHq0qWL8vLyJEl5eXmKjY21Qp0k9ezZU0FBQdq8ebPVp3v37l5Xa09PT9fevXv1/fffV1ub2+2Wy+XyWgAAAALJGQW7/fv3a/78+WrTpo3Wrl2rUaNG6e6779aSJUskSYWFhZKkxMREr+clJiZa2woLC5WQkOC1PSQkRHFxcV59qnuNH7/Hf5oxY4ZiYmKsJSUl5Uw+GgDUqrCwMM2dO1dz587llmIA6swZXVjJ4/Goc+fOeuSRRyRJF198sXbt2qUFCxZo2LBhPimwprKzszVhwgTrscvlItwB8JvQ0FBlZWX5uwwAAeaMZuyaNm2q1NRUr7Z27dopPz9fkpSUlCRJKioq8upTVFRkbUtKSlJxcbHX9srKSh05csSrT3Wv8eP3+E/h4eFyOp1eCwAAQCA5o2DXrVs37d2716vtiy++UIsWLSRJrVq1UlJSknJzc63tLpdLmzdvVlpamiQpLS1NJSUl2rZtm9Vn3bp18ng86tKli9Xn/fffV0VFhdUnJydH5513ntcZuABQX1VVVWn9+vVav369qqqq/F0OgABxRsFu/Pjx2rRpkx555BHt27dPS5cu1TPPPGPtbnA4HBo3bpymT5+uN998Uzt37tTQoUOVnJys/v37Szo5w9e7d2+NHDlSW7Zs0UcffaQxY8Zo4MCBSk5OliQNGjRIYWFhGjFihHbv3q1XXnlFf/vb37x2tQJAfVZWVqarr75aV199tdflngDAl87oGLtLLrlEK1asUHZ2tqZNm6ZWrVrpiSee0ODBg60+9913n44dO6Y77rhDJSUluvzyy7VmzRpFRERYfV566SWNGTNGPXr0UFBQkDIzMzVnzhxre0xMjN555x1lZWWpU6dOatKkiaZMmeJ1rTsAAAB4O6Pr2J1NuI4dAH/iOnaobVzHLnD57Dp2AAAAqL8IdgAAADZBsAMAALAJgh0AAIBNnNFZsQCAmgkNDdXMmTOtdQCoCwQ7APCBsLAwTZw40d9lAAgw7IoFAACwCWbsAMAHqqqqtH37dklSx44dFRwc7OeKAAQCgh0A+EBZWZkuvfRSSVygGEDdYVcsAACATRDsAAAAbIJgBwAAYBMEOwAAAJsg2AEAANgEwQ4AAMAmuNwJAPhAaGiopk6daq0DQF0g2AGAD4SFhenBBx/0dxkAAgy7YgEAAGyCGTsA8AGPx6PPPvtMktSuXTsFBfF3NADfI9gBgA+cOHFC7du3l8QtxQDUHf6EBAAAsAmCHQAAgE0Q7AAAAGyCYAcAAGATBDsAAACbINgBAADYBJc7AQAfCA0N1b333mutA0BdcBhjjL+L8AWXy6WYmBiVlpbK6XT6uxwAAH6TlpNW+bsEVOPAoxk+f48zyTTsigUAALAJdsUCgA94PB7l5+dLkpo3b84txQDUCYIdAPjAiRMn1KpVK0ncUgxA3eFPSAAAAJtgxg4AYEucbIBAxIwdAACATRDsAAAAbIJgBwAAYBMEOwAAAJvg5AkA8IGQkBCNHj3aWgeAusBvGwDwgfDwcM2bN8/fZQAIMOyKBQAAsAlm7ADAB4wx+u677yRJTZo0kcPh8HNFAALBb5qxe/TRR+VwODRu3DirraysTFlZWWrcuLGio6OVmZmpoqIir+fl5+crIyNDkZGRSkhI0MSJE1VZWenVZ/369erYsaPCw8PVunVrLV68+LeUCgB16vjx40pISFBCQoKOHz/u73IABIhfHey2bt2qp59+WhdccIFX+/jx4/XWW29p2bJl2rBhgw4ePKgbbrjB2l5VVaWMjAyVl5dr48aNWrJkiRYvXqwpU6ZYfb7++mtlZGTo6quv1o4dOzRu3DjdfvvtWrt27a8tFwAAwPZ+VbA7evSoBg8erGeffVaNGjWy2ktLS7Vw4ULNmjVL11xzjTp16qRFixZp48aN2rRpkyTpnXfe0Z49e/Tiiy/qoosuUp8+ffTnP/9Z8+bNU3l5uSRpwYIFatWqlR5//HG1a9dOY8aM0Y033qjZs2fXwkcGAACwp191jF1WVpYyMjLUs2dPTZ8+3Wrftm2bKioq1LNnT6utbdu2at68ufLy8tS1a1fl5eWpQ4cOSkxMtPqkp6dr1KhR2r17ty6++GLl5eV5vcapPj/e5fuf3G633G639djlcv2ajwYAAY37qwJntzMOdi+//LK2b9+urVu3nratsLBQYWFhio2N9WpPTExUYWGh1efHoe7U9lPbfq6Py+XSiRMn1KBBg9Pee8aMGXrooYfO9OMAAADYxhntii0oKNA999yjl156SREREb6q6VfJzs5WaWmptRQUFPi7JAAAgDp1RsFu27ZtKi4uVseOHRUSEqKQkBBt2LBBc+bMUUhIiBITE1VeXq6SkhKv5xUVFSkpKUmSlJSUdNpZsqce/1Ifp9NZ7WyddPJioE6n02sBAAAIJGe0K7ZHjx7auXOnV9utt96qtm3b6v7771dKSopCQ0OVm5urzMxMSdLevXuVn5+vtLQ0SVJaWpoefvhhFRcXKyEhQZKUk5Mjp9Op1NRUq8/bb7/t9T45OTnWawBAfRcSEqJhw4ZZ6wBQF87ot03Dhg3Vvn17r7aoqCg1btzYah8xYoQmTJiguLg4OZ1OjR07VmlpaerataskqVevXkpNTdUtt9yimTNnqrCwUJMnT1ZWVpbCw8MlSXfddZfmzp2r++67T7fddpvWrVunV199VatWcVAvgLNDeHg4198EUOdq/c/I2bNnKygoSJmZmXK73UpPT9dTTz1lbQ8ODtbKlSs1atQopaWlKSoqSsOGDdO0adOsPq1atdKqVas0fvx4/e1vf1OzZs303HPPKT09vbbLBQAAsA2HMcb4uwhfcLlciomJUWlpKcfbAahzxhjrjhORkZFnzS3FuNwJcGYOPJrh8/c4k0zzm24pBgCo3vHjxxUdHa3o6GhuKQagzhDsAAAAbIJgBwAAYBMEOwAAAJsg2AEAANgEwQ4AAMAmCHYAAAA2wX1uAMAHgoODdeONN1rrAFAXCHYA4AMRERFatmyZv8sAEGAIdgDwG3CnBgD1CcfYAQAA2ATBDgB8wFNepm8e66tvHusrT3mZv8sBECAIdgAAADZBsAMAALAJgh0AAIBNEOwAAABsgmAHAABgEwQ7AAAAm+ACxQDgA46gIDU4p7O1DgB1gWAHAD7gCAlTwk0P+rsMAAGGPyMBAABsgmAHAABgEwQ7APABT3mZ8mdlKn9WJrcUA1BnOMYOAHzEVLj9XQKAAMOMHQAAgE0Q7AAAAGyCYAcAAGATBDsAAACbINgBAADYBGfFAoAvOBwKT2lvrQNAXSDYAYAPBIWGK2nQo/4uA0CAYVcsAACATRDsAAAAbIJgBwA+4CkvU8GcQSqYM4hbigGoMxxjBwA+4jnh8ncJAAIMM3YAAAA2QbADAACwCYIdAACATRDsAAAAbIJgBwAAYBOcFQsAvuBwKCypjbUOAHXhjGbsZsyYoUsuuUQNGzZUQkKC+vfvr71793r1KSsrU1ZWlho3bqzo6GhlZmaqqKjIq09+fr4yMjIUGRmphIQETZw4UZWVlV591q9fr44dOyo8PFytW7fW4sWLf90nBAA/CAoNV9Nhs9V02GwFhYb7uxwAAeKMgt2GDRuUlZWlTZs2KScnRxUVFerVq5eOHTtm9Rk/frzeeustLVu2TBs2bNDBgwd1ww03WNurqqqUkZGh8vJybdy4UUuWLNHixYs1ZcoUq8/XX3+tjIwMXX311dqxY4fGjRun22+/XWvXrq2FjwwAAGBPDmOM+bVPPnz4sBISErRhwwZ1795dpaWlio+P19KlS3XjjTdKkj7//HO1a9dOeXl56tq1q1avXq2+ffvq4MGDSkxMlCQtWLBA999/vw4fPqywsDDdf//9WrVqlXbt2mW918CBA1VSUqI1a9bUqDaXy6WYmBiVlpbK6XT+2o8IAD+r5aRV/i4BgB8deDTD5+9xJpnmN508UVpaKkmKi4uTJG3btk0VFRXq2bOn1adt27Zq3ry58vLyJEl5eXnq0KGDFeokKT09XS6XS7t377b6/Pg1TvU59RrVcbvdcrlcXgsA+Iunokz/mn+b/jX/NnkquKUYgLrxq4Odx+PRuHHj1K1bN7Vv316SVFhYqLCwMMXGxnr1TUxMVGFhodXnx6Hu1PZT236uj8vl0okTJ6qtZ8aMGYqJibGWlJSUX/vRAOC3M1KVq1hVrmLpV+8XAYAz86uDXVZWlnbt2qWXX365Nuv51bKzs1VaWmotBQUF/i4JAACgTv2qy52MGTNGK1eu1Pvvv69mzZpZ7UlJSSovL1dJSYnXrF1RUZGSkpKsPlu2bPF6vVNnzf64z3+eSVtUVCSn06kGDRpUW1N4eLjCwznzDAAABK4zCnbGGI0dO1YrVqzQ+vXr1apVK6/tnTp1UmhoqHJzc5WZmSlJ2rt3r/Lz85WWliZJSktL08MPP6zi4mIlJCRIknJycuR0OpWammr1efvtt71eOycnx3oNAGc3TjgAAN84o2CXlZWlpUuX6o033lDDhg2tY+JiYmLUoEEDxcTEaMSIEZowYYLi4uLkdDo1duxYpaWlqWvXrpKkXr16KTU1VbfccotmzpypwsJCTZ48WVlZWdaM21133aW5c+fqvvvu02233aZ169bp1Vdf1apV/GcAAADwU87oGLv58+ertLRUV111lZo2bWotr7zyitVn9uzZ6tu3rzIzM9W9e3clJSXptddes7YHBwdr5cqVCg4OVlpamoYMGaKhQ4dq2rRpVp9WrVpp1apVysnJ0YUXXqjHH39czz33nNLT02vhIwMAANjTb7qOXX3GdeyA+isQdsV6KspUuGSCJClp2CwFhUb4uSIAvlDfrmPHvWIBwAeCQiOUfPtT/i4DQID5TRcoBgAAQP1BsAMAALAJgh0A+ICnokwHnxutg8+N5pZiAOoMx9gBgC8YqeLf+dY6ANQFZuwAAABsgmAHAABgEwQ7AAAAm+AYO+AsEQgX9QUA/DYEO9gWQQgAEGgIdgDgCw4p2JlgrQNAXSDYAYAPBIVGqNmo5/1dBoAAw8kTAAAANkGwAwAAsAmCHQD4gKfCrUNLxuvQkvHyVLj9XQ6AAMExdgDgC8aovPBLax0A6gIzdgAAADbBjB28cO03AADOXszYAQAA2ATBDgAAwCYIdgAAADbBMXYA4CNBDZz+LgFAgCHYAYAPBIVFKOXupf4uA0CAYVcsAACATRDsAAAAbIJgBwA+4Klwq3DpJBUuncQtxQDUGY6xAwBfMEbugl3WOgDUBWbsAAAAbIJgBwAAYBMEOwAAAJsg2AEAANgEwQ4AAMAmOCsWAHzEERru7xIABBiCHQD4QFBYhJpP+Ie/ywAQYNgVCwAAYBMEOwAAAJtgV2wtaDlplb9LAFDPmMpyHV7xiCQp/vo/yRES5ueKAAQCgh0A+IDxeHRi/8fWusPP9QAIDOyKBQAAsAmCHQAAgE0Q7AAAAGyiXge7efPmqWXLloqIiFCXLl20ZcsWf5cEAABQb9XbYPfKK69owoQJmjp1qrZv364LL7xQ6enpKi4u9ndpAAAA9VK9PSt21qxZGjlypG699VZJ0oIFC7Rq1So9//zzmjRp0mn93W633G639bi0tFSS5HK5fF6rx33c5+8B4OziKS/7v3X3ccl4/FgNAF+pi5xx6j2MMb/c2dRDbrfbBAcHmxUrVni1Dx061PTr16/a50ydOtVIYmFhYWFhYWGx5VJQUPCLGapezth99913qqqqUmJiold7YmKiPv/882qfk52drQkTJliPPR6Pjhw5osaNG8vh4ApSLpdLKSkpKigokNPp9Hc59RbjVDOMU80wTjXDONUM41QzdhwnY4x++OEHJScn/2Lfehnsfo3w8HCFh4d7tcXGxvqnmHrM6XTa5ovuS4xTzTBONcM41QzjVDOMU83YbZxiYmJq1K9enjzRpEkTBQcHq6ioyKu9qKhISUlJfqoKAACgfquXwS4sLEydOnVSbm6u1ebxeJSbm6u0tDQ/VgYAAFB/1dtdsRMmTNCwYcPUuXNnXXrppXriiSd07Ngx6yxZnJnw8HBNnTr1tN3V8MY41QzjVDOMU80wTjXDONVMoI+Tw5ianDvrH3PnztVf/vIXFRYW6qKLLtKcOXPUpUsXf5cFAABQL9XrYAcAAICaq5fH2AEAAODMEewAAABsgmAHAABgEwQ7AAAAmyDYAQAA2ATBDgAAwCYIdgAAADZBsAMAALAJgh0AAIBNEOwAAABsgmAHAABgEwQ7AAAAmyDYAQAA2ATBDgAAwCYIdgAAADZBsAMAALAJgh0AAIBNEOwAoB5xOBx68MEHa+311q9fL4fDofXr19faawKovwh2AAAANkGwAwAAsAmCHYCzxrFjx/xdAgDUawQ7APXSgw8+KIfDoT179mjQoEFq1KiRLr/8cknSiy++qE6dOqlBgwaKi4vTwIEDVVBQ4PX8L7/8UpmZmUpKSlJERISaNWumgQMHqrS01OqzaNEiXXPNNUpISFB4eLhSU1M1f/7802pp2bKl+vbtq/Xr16tz585q0KCBOnToYB239tprr6lDhw6KiIhQp06d9Mknn3g9f/jw4YqOjtb+/fuVnp6uqKgoJScna9q0aTLG/OJYfPvtt7rtttuUmJio8PBwnX/++Xr++edP6/evf/1L/fv3V1RUlBISEjR+/Hi53e5ffH0A9hHi7wIA4OfcdNNNatOmjR555BEZY/Twww/rgQce0IABA3T77bfr8OHDevLJJ9W9e3d98sknio2NVXl5udLT0+V2uzV27FglJSXp22+/1cqVK1VSUqKYmBhJ0vz583X++eerX79+CgkJ0VtvvaXRo0fL4/EoKyvLq459+/Zp0KBBuvPOOzVkyBD99a9/1XXXXacFCxboT3/6k0aPHi1JmjFjhgYMGKC9e/cqKOj//nauqqpS79691bVrV82cOVNr1qzR1KlTVVlZqWnTpv3k5y8qKlLXrl3lcDg0ZswYxcfHa/Xq1RoxYoRcLpfGjRsnSTpx4oR69Oih/Px83X333UpOTtYLL7ygdevW1fK/CIB6zQBAPTR16lQjydx8881W24EDB0xwcLB5+OGHvfru3LnThISEWO2ffPKJkWSWLVv2s+9x/Pjx09rS09PNOeec49XWokULI8ls3LjRalu7dq2RZBo0aGC++eYbq/3pp582ksx7771ntQ0bNsxIMmPHjrXaPB6PycjIMGFhYebw4cNWuyQzdepU6/GIESNM06ZNzXfffedV08CBA01MTIz1GZ544gkjybz66qtWn2PHjpnWrVufVg8A+2JXLIB67a677rLWX3vtNXk8Hg0YMEDfffedtSQlJalNmzZ67733JMmakVu7dq2OHz/+k6/doEEDa720tFTfffedrrzySu3fv99rl60kpaamKi0tzXrcpUsXSdI111yj5s2bn9a+f//+095vzJgx1vqpGbjy8nK9++671dZnjNE//vEPXXfddTLGeH3m9PR0lZaWavv27ZKkt99+W02bNtWNN95oPT8yMlJ33HHHT35+APbDrlgA9VqrVq2s9S+//FLGGLVp06bavqGhodZzJkyYoFmzZumll17SFVdcoX79+mnIkCFW6JOkjz76SFOnTlVeXt5pAbC0tNSr74/Dm/R/4TElJaXa9u+//96rPSgoSOecc45X2+9//3tJ0oEDB6r9PIcPH1ZJSYmeeeYZPfPMM9X2KS4uliR98803at26tRwOh9f28847r9rnAbAngh2Aeu3Hs2oej0cOh0OrV69WcHDwaX2jo6Ot9ccff1zDhw/XG2+8oXfeeUd33323ZsyYoU2bNqlZs2b66quv1KNHD7Vt21azZs1SSkqKwsLC9Pbbb2v27NnyeDxer13d+/1cu6nBSRG/5FQNQ4YM0bBhw6rtc8EFF/zm9wFgHwQ7AGeNc889V8YYtWrVyprt+jkdOnRQhw4dNHnyZG3cuFHdunXTggULNH36dL311ltyu9168803vWbjTu3OrW0ej0f79+/3qvuLL76QdPKs2+rEx8erYcOGqqqqUs+ePX/29Vu0aKFdu3bJGOM1a7d3797fXjyAswbH2AE4a9xwww0KDg7WQw89dNqMmDFG//73vyVJLpdLlZWVXts7dOigoKAg6/Ifp2bafvw6paWlWrRokc/qnzt3rle9c+fOVWhoqHr06FFt/+DgYGVmZuof//iHdu3addr2w4cPW+vXXnutDh48qOXLl1ttx48f/8lduADsiRk7AGeNc889V9OnT1d2drYOHDig/v37q2HDhvr666+1YsUK3XHHHbr33nu1bt06jRkzRjfddJN+//vfq7KyUi+88IIVlCSpV69eCgsL03XXXac777xTR48e1bPPPquEhAQdOnSo1muPiIjQmjVrNGzYMHXp0kWrV6/WqlWr9Kc//Unx8fE/+bxHH31U7733nrp06aKRI0cqNTVVR44c0fbt2/Xuu+/qyJEjkqSRI0dq7ty5Gjp0qLZt26amTZvqhRdeUGRkZK1/FgD1F8EOwFll0qRJ+v3vf6/Zs2froYceknTyBIZevXqpX79+kqQLL7xQ6enpeuutt/Ttt98qMjJSF154oVavXq2uXbtKOnlSwfLlyzV58mTde++9SkpK0qhRoxQfH6/bbrut1usODg7WmjVrNGrUKE2cOFENGzbU1KlTNWXKlJ99XmJiorZs2aJp06bptdde01NPPaXGjRvr/PPP12OPPWb1i4yMVG5ursaOHasnn3xSkZGRGjx4sPr06aPevXvX+ucBUD85TG0c4QsA+EnDhw/X8uXLdfToUX+XAsDmOMYOAADAJgh2AAAANkGwAwAAsAmOsQMAALAJZuwAAABswraXO/F4PDp48KAaNmx42r0TAQAAzhbGGP3www9KTk5WUNDPz8nZNtgdPHjwtJtzAwAAnK0KCgrUrFmzn+1j22DXsGFDSScHwel0+rkaAIHm2LFjSk5OlnTyD82oqCg/VwTgbOVyuZSSkmJlm59j22B3aver0+kk2AGoc5GRkXr11VclSfHx8QoJse2vWwB1pCaHlvGbBgB8ICQkRDfddJO/ywAQYDgrFgAAwCaYsQMAH6isrNSKFSskSddffz27YgHUCX7TAIAPuN1uDRgwQJJ09OhRgh2AOsGuWAAAAJsg2AEAANgE+wYAAMBv0nLSKn+XUCsOPJrh7xJ+M2bsAAAAbIJgBwAAYBMEOwAAAJvgGDsA8IGwsDAtWrTIWgeAukCwAwAfCA0N1fDhw/1dBoAAw65YAAAAm2DGDgB8oLKyUmvXrpUkpaenc+cJAHWC3zQA4ANut1t9+/aVxC3FANQddsUCAADYBMEOAADAJgh2AAAANkGwAwAAsAmCHQAAgE1wmhYAn2o5aZW/S6gVBx7N8HcJAPCLCHYA4ANhYWGaO3eutQ5Uxy5/+KD+INgBgA+EhoYqKyvL32UACDAcYwcAAGATzNgBgA9UVVXpgw8+kCRdccUVCg4O9nNFAAIBwQ4AfKCsrExXX321pJO3FIuKivJzRQACAbtiAQAAbIJgBwAAYBM+CXbffvuthgwZosaNG6tBgwbq0KGDPv74Y2u7MUZTpkxR06ZN1aBBA/Xs2VNffvml12scOXJEgwcPltPpVGxsrEaMGKGjR4/6olwAAABbqPVg9/3336tbt24KDQ3V6tWrtWfPHj3++ONq1KiR1WfmzJmaM2eOFixYoM2bNysqKkrp6ekqKyuz+gwePFi7d+9WTk6OVq5cqffff1933HFHbZcLAABgG7V+8sRjjz2mlJQULVq0yGpr1aqVtW6M0RNPPKHJkyfrD3/4gyTp73//uxITE/X6669r4MCB+uyzz7RmzRpt3bpVnTt3liQ9+eSTuvbaa/XXv/5VycnJtV02AADAWa/WZ+zefPNNde7cWTfddJMSEhJ08cUX69lnn7W2f/311yosLFTPnj2ttpiYGHXp0kV5eXmSpLy8PMXGxlqhTpJ69uypoKAgbd68udr3dbvdcrlcXgsAAEAgqfVgt3//fs2fP19t2rTR2rVrNWrUKN19991asmSJJKmwsFCSlJiY6PW8xMREa1thYaESEhK8toeEhCguLs7q859mzJihmJgYa0lJSantjwYANRYaGqqZM2dq5syZCg0N9Xc5AAJEre+K9Xg86ty5sx555BFJ0sUXX6xdu3ZpwYIFGjZsWG2/nSU7O1sTJkywHrtcLsIdAL8JCwvTxIkT/V0GgABT6zN2TZs2VWpqqldbu3btlJ+fL0lKSkqSJBUVFXn1KSoqsrYlJSWpuLjYa3tlZaWOHDli9flP4eHhcjqdXgsAAEAgqfVg161bN+3du9er7YsvvlCLFi0knTyRIikpSbm5udZ2l8ulzZs3Ky0tTZKUlpamkpISbdu2zeqzbt06eTwedenSpbZLBoBaV1VVpa1bt2rr1q2qqqrydzkAAkSt74odP368LrvsMj3yyCMaMGCAtmzZomeeeUbPPPOMJMnhcGjcuHGaPn262rRpo1atWumBBx5QcnKy+vfvL+nkDF/v3r01cuRILViwQBUVFRozZowGDhzIGbEAzgplZWW69NJLJXFLMQB1p9aD3SWXXKIVK1YoOztb06ZNU6tWrfTEE09o8ODBVp/77rtPx44d0x133KGSkhJdfvnlWrNmjSIiIqw+L730ksaMGaMePXooKChImZmZmjNnTm2XCwAAYBsOY4zxdxG+4HK5FBMTo9LSUo63A/yo5aRV/i6hVhx4NOOM+h87dkzR0dGSmLHDT7PLz4ddnOnPeV05k0zDvWIBAABsgmAHAABgEwQ7AAAAm6j1kycAwI7O9FgoT3mZtd7ugTUKCov4md51p74eQwSgdhDsAMAHHMHBiul2s7UOAHWBYAcAPuAIDlXs5YN/uSMA1CKOsQMAALAJZuwAwAeM8ajiuwJJUmiTFDkc/B0NwPcIdgDgA6aiXIeez5IkpYxfLkc9OXkCgL3xJyQAAIBNMGMH1FPcaggAcKaYsQMAALAJgh0AAIBNEOwAAABsgmAHAABgE5w8AQA+4AgOlvPSG6x1AKgLBDvYDmeToj5wBIeq0dW3+bsMAAGGXbEAAAA2wYwdAPiAMR5VuQ5LkoKd8dxSDECdINgBgA+YinJ9u2CEJG4pBqDu8CckAACATRDsAAAAbIJgBwAAYBMcYwcAOOtwWSOgeszYAQAA2ATBDgAAwCbYFQsAPuAIClb0xRnWOgDUBYIdAPiAIyRUjXuN8ncZAAIMu2IBAABsghk7APABY4w8J1ySpKAGTjkcDj9XBCAQMGMHAD5gKtz615OD9a8nB8tUuP1dDoAAQbADAACwCYIdAACATRDsAAAAbIJgBwAAYBMEOwAAAJsg2AEAANiEz4Pdo48+KofDoXHjxlltZWVlysrKUuPGjRUdHa3MzEwVFRV5PS8/P18ZGRmKjIxUQkKCJk6cqMrKSl+XCwC1whEUrKj2PRTVvge3FANQZ3x6geKtW7fq6aef1gUXXODVPn78eK1atUrLli1TTEyMxowZoxtuuEEfffSRJKmqqkoZGRlKSkrSxo0bdejQIQ0dOlShoaF65JFHfFlyQGs5aZW/SwBswxESqiYZ4/1dxmn4OQfszWczdkePHtXgwYP17LPPqlGjRlZ7aWmpFi5cqFmzZumaa65Rp06dtGjRIm3cuFGbNm2SJL3zzjvas2ePXnzxRV100UXq06eP/vznP2vevHkqLy/3VckAAABnNZ8Fu6ysLGVkZKhnz55e7du2bVNFRYVXe9u2bdW8eXPl5eVJkvLy8tShQwclJiZafdLT0+VyubR79+5q38/tdsvlcnktAOAvxhh5ysvkKS+TMcbf5QAIED7ZFfvyyy9r+/bt2rp162nbCgsLFRYWptjYWK/2xMREFRYWWn1+HOpObT+1rTozZszQQw89VAvVA8BvZyrcKph9oyQpZfxyOcIi/FwRgEBQ6zN2BQUFuueee/TSSy8pIqLufpFlZ2ertLTUWgoKCursvQEAAOqDWg9227ZtU3FxsTp27KiQkBCFhIRow4YNmjNnjkJCQpSYmKjy8nKVlJR4Pa+oqEhJSUmSpKSkpNPOkj31+FSf/xQeHi6n0+m1AAAABJJaD3Y9evTQzp07tWPHDmvp3LmzBg8ebK2HhoYqNzfXes7evXuVn5+vtLQ0SVJaWpp27typ4uJiq09OTo6cTqdSU1Nru2QAAABbqPVj7Bo2bKj27dt7tUVFRalx48ZW+4gRIzRhwgTFxcXJ6XRq7NixSktLU9euXSVJvXr1Umpqqm655RbNnDlThYWFmjx5srKyshQeHl7bJQMAANiCT69j91Nmz56toKAgZWZmyu12Kz09XU899ZS1PTg4WCtXrtSoUaOUlpamqKgoDRs2TNOmTfNHuQAAAGeFOgl269ev93ocERGhefPmad68eT/5nBYtWujtt9/2cWUAAAD24ZcZOwCwO0dQkCLP62atA0BdINgBgA84QsIU3z/b32UACDD8GQkAAGATBDsAAACbINgBgA94ysv0zWN99c1jfeUpL/N3OQACBMEOAADAJgh2AAAANkGwAwAAsAmCHQAAgE0Q7AAAAGyCYAcAAGAT3HkCAHzAERSkBud0ttYBoC4Q7ADABxwhYUq46UF/lwEgwPBnJAAAgE0Q7AAAAGyCYAcAPuApL1P+rEzlz8rklmIA6gzH2AGAj5gKt79LABBgmLEDAACwCYIdAACATRDsAAAAbIJgBwAAYBMEOwAAAJvgrFgA8AWHQ+Ep7a11AKgLBDsA8IGg0HAlDXrU32UACDDsigUAALAJgh0AAIBNEOwAwAc85WUqmDNIBXMGcUsxAHWGY+wAwEc8J1z+LgFAgGHGDgAAwCYIdgAAADZBsAMAALAJgh0AAIBNEOwAAABsgrNiAcAXHA6FJbWx1gGgLhDsAMAHgkLD1XTYbH+XASDAsCsWAADAJgh2AAAANkGwAwAf8FSU6V/zb9O/5t8mTwW3FANQN2o92M2YMUOXXHKJGjZsqISEBPXv31979+716lNWVqasrCw1btxY0dHRyszMVFFRkVef/Px8ZWRkKDIyUgkJCZo4caIqKytru1wA8A0jVbmKVeUqloy/iwEQKGo92G3YsEFZWVnatGmTcnJyVFFRoV69eunYsWNWn/Hjx+utt97SsmXLtGHDBh08eFA33HCDtb2qqkoZGRkqLy/Xxo0btWTJEi1evFhTpkyp7XIBAABsw2GM8enfkocPH1ZCQoI2bNig7t27q7S0VPHx8Vq6dKluvPFGSdLnn3+udu3aKS8vT127dtXq1avVt29fHTx4UImJiZKkBQsW6P7779fhw4cVFhZ22vu43W653W7rscvlUkpKikpLS+V0On35EW2j5aRV/i4BsA1PeZkKZp/8HZcyfrmCwiL8XBGAX3Lg0Qx/l1Atl8ulmJiYGmUanx9jV1paKkmKi4uTJG3btk0VFRXq2bOn1adt27Zq3ry58vLyJEl5eXnq0KGDFeokKT09XS6XS7t37672fWbMmKGYmBhrSUlJ8dVHAgAAqJd8Guw8Ho/GjRunbt26qX379pKkwsJChYWFKTY21qtvYmKiCgsLrT4/DnWntp/aVp3s7GyVlpZaS0FBQS1/GgAAgPrNpxcozsrK0q5du/Thhx/68m0kSeHh4QoPD/f5+wAAANRXPpuxGzNmjFauXKn33ntPzZo1s9qTkpJUXl6ukpISr/5FRUVKSkqy+vznWbKnHp/qAwD1mkMKbdxcoY2bS9xRDEAdqfVgZ4zRmDFjtGLFCq1bt06tWrXy2t6pUyeFhoYqNzfXatu7d6/y8/OVlpYmSUpLS9POnTtVXFxs9cnJyZHT6VRqamptlwwAtS4oNELJtz+l5NufUlAoJ04AqBu1vis2KytLS5cu1RtvvKGGDRtax8TFxMSoQYMGiomJ0YgRIzRhwgTFxcXJ6XRq7NixSktLU9euXSVJvXr1Umpqqm655RbNnDlThYWFmjx5srKystjdCgAA8BNqPdjNnz9fknTVVVd5tS9atEjDhw+XJM2ePVtBQUHKzMyU2+1Wenq6nnrqKatvcHCwVq5cqVGjRiktLU1RUVEaNmyYpk2bVtvlAgAA2IbPr2PnL2dyzRecxHXsgNrjqShT4ZIJkqSkYbPYHQucBexwHTufnhULAAHLSBX/zrfWAaAu+PwCxQAAAKgbBDsAAACbINgBAADYBMEOAADAJgh2AAAANsFZsQDgCw4p2JlgrQNAXSDYAYAPBIVGqNmo5/1dBoAAw65YAAAAmyDYAQAA2ATBDgB8wFPh1qEl43VoyXh5Ktz+LgdAgOAYu1rAPVYBnMYYlRd+aa0DQF1gxg4AAMAmCHYAAAA2QbADAACwCYIdAACATRDsAAAAbIKzYgHAR4IaOP1dAoAAQ7ADAB8ICotQyt1L/V0GgADDrlgAAACbINgBAADYBMEOAHzAU+FW4dJJKlw6iVuKAagzHGMHAL5gjNwFu6x1AKgLzNgBAADYBMEOAADAJgh2AAAANkGwAwAAsAmCHQAAgE1wViwA+IgjNNzfJQAIMAQ7APCBoLAINZ/wD3+XASDAsCsWAADAJgh2AAAANsGuWADwAVNZrsMrHpEkxV//JzlCwvxcEYBAQLADAB8wHo9O7P/YWnf4uR4AgYFdsQAAADZBsAMAALAJgh0AAIBN1OtgN2/ePLVs2VIRERHq0qWLtmzZ4u+SAAAA6q16G+xeeeUVTZgwQVOnTtX27dt14YUXKj09XcXFxf4uDQAAoF6qt2fFzpo1SyNHjtStt94qSVqwYIFWrVql559/XpMmTTqtv9vtltvtth6XlpZKklwul89r9biP+/w9AJxdPOVl/7fuPi4Zjx+rAVATdZEZfo1TdRljfrmzqYfcbrcJDg42K1as8GofOnSo6devX7XPmTp1qpHEwsLCwsLCwmLLpaCg4BczVL2csfvuu+9UVVWlxMREr/bExER9/vnn1T4nOztbEyZMsB57PB4dOXJEjRs3lsPBFaRcLpdSUlJUUFAgp9Pp73LqLcapZhinmmGcaoZxqhnGqWbsOE7GGP3www9KTk7+xb71Mtj9GuHh4QoPD/dqi42N9U8x9ZjT6bTNF92XGKeaYZxqhnGqGcapZhinmrHbOMXExNSoX708eaJJkyYKDg5WUVGRV3tRUZGSkpL8VBUAAED9Vi+DXVhYmDp16qTc3FyrzePxKDc3V2lpaX6sDAAAoP6qt7tiJ0yYoGHDhqlz58669NJL9cQTT+jYsWPWWbI4M+Hh4Zo6deppu6vhjXGqGcapZhinmmGcaoZxqplAHyeHMTU5d9Y/5s6dq7/85S8qLCzURRddpDlz5qhLly7+LgsAAKBeqtfBDgAAADVXL4+xAwAAwJkj2AEAANgEwQ4AAMAmCHYAAAA2QbA7i82YMUOXXHKJGjZsqISEBPXv31979+6ttq8xRn369JHD4dDrr7/utS0/P18ZGRmKjIxUQkKCJk6cqMrKyjr4BHWjJuN01VVXyeFweC133XWXVx/G6aS8vDxdc801ioqKktPpVPfu3XXixAlr+5EjRzR48GA5nU7FxsZqxIgROnr0aF1+FJ/6pXE6cODAad+lU8uyZcusfnyfpMLCQt1yyy1KSkpSVFSUOnbsqH/84x9efQL9+yRJX331la6//nrFx8fL6XRqwIABp13A3+7jNH/+fF1wwQXW3STS0tK0evVqa3tZWZmysrLUuHFjRUdHKzMz87QxsvvPnOUX7yaLeis9Pd0sWrTI7Nq1y+zYscNce+21pnnz5ubo0aOn9Z01a5bp06ePkWRWrFhhtVdWVpr27dubnj17mk8++cS8/fbbpkmTJiY7O7sOP4lv1WScrrzySjNy5Ehz6NAhayktLbW2M04nbdy40TidTjNjxgyza9cu8/nnn5tXXnnFlJWVWX169+5tLrzwQrNp0ybzwQcfmNatW5ubb77ZHx/JJ35pnCorK72+R4cOHTIPPfSQiY6ONj/88IPVh++TMf/1X/9lLrnkErN582bz1VdfmT//+c8mKCjIbN++3eoT6N+no0ePmnPOOcdcf/315tNPPzWffvqp+cMf/mAuueQSU1VVZb2O3cfpzTffNKtWrTJffPGF2bt3r/nTn/5kQkNDza5du4wxxtx1110mJSXF5Obmmo8//th07drVXHbZZdbzA+Fn7hSCnY0UFxcbSWbDhg1e7Z988on53e9+Zw4dOnRasHv77bdNUFCQKSwstNrmz59vnE6ncbvddVV6napunK688kpzzz33/ORzGKeTunTpYiZPnvyTz9mzZ4+RZLZu3Wq1rV692jgcDvPtt9/6tF5/+amfux+76KKLzG233WY95vt0UlRUlPn73//u1S8uLs48++yzxhi+T8YYs3btWhMUFOT1h2ZJSYlxOBwmJyfHGBOY42SMMY0aNTLPPfecKSkpMaGhoWbZsmXWts8++8xIMnl5ecaYwPqZY1esjZSWlkqS4uLirLbjx49r0KBBmjdvXrX32c3Ly1OHDh2UmJhotaWnp8vlcmn37t2+L9oPqhsnSXrppZfUpEkTtW/fXtnZ2Tp+/Li1jXGSiouLtXnzZiUkJOiyyy5TYmKirrzySn344YfWc/Ly8hQbG6vOnTtbbT179lRQUJA2b95ctx+gjvzU9+mUbdu2aceOHRoxYoTVxvfppMsuu0yvvPKKjhw5Io/Ho5dfflllZWW66qqrJPF9kiS32y2Hw+F1F4WIiAgFBQVZP3uBNk5VVVV6+eWXdezYMaWlpWnbtm2qqKhQz549rT5t27ZV8+bNlZeXJymwfuYIdjbh8Xg0btw4devWTe3bt7fax48fr8suu0x/+MMfqn1eYWGh1xddkvW4sLDQdwX7yU+N06BBg/Tiiy/qvffeU3Z2tl544QUNGTLE2s44Sfv375ckPfjggxo5cqTWrFmjjh07qkePHvryyy8lnRyLhIQEr9cKCQlRXFxcwIzTf1q4cKHatWunyy67zGrj+3TSq6++qoqKCjVu3Fjh4eG68847tWLFCrVu3VoS3ydJ6tq1q6KionT//ffr+PHjOnbsmO69915VVVXp0KFDkgJnnHbu3Kno6GiFh4frrrvu0ooVK5SamqrCwkKFhYUpNjbWq39iYqL1+QPpZ67e3isWZyYrK0u7du3ymj158803tW7dOn3yySd+rKx+qW6cJOmOO+6w1jt06KCmTZuqR48e+uqrr3TuuefWdZl+V904eTweSdKdd95p3bP54osvVm5urp5//nnNmDHDL7X60099n045ceKEli5dqgceeKCOK6tffmqcHnjgAZWUlOjdd99VkyZN9Prrr2vAgAH64IMP1KFDBz9V6z/VjVN8fLyWLVumUaNGac6cOQoKCtLNN9+sjh07KigosOZmzjvvPO3YsUOlpaVavny5hg0bpg0bNvi7rHqHYGcDY8aM0cqVK/X++++rWbNmVvu6dev01VdfnfZXTGZmpq644gqtX79eSUlJ2rJli9f2U2cSVbfr9mz2U+NUnVP3JN63b5/OPfdcxklS06ZNJUmpqale/du1a6f8/HxJJ8eiuLjYa3tlZaWOHDkSMOP0Y8uXL9fx48c1dOhQr3a+TyfP9Jw7d6527dql888/X5J04YUX6oMPPtC8efO0YMECvk//X69evfTVV1/pu+++U0hIiGJjY5WUlKRzzjlHUuD83IWFhVmzuZ06ddLWrVv1t7/9Tf/93/+t8vJylZSUeP1/V1RUZH3+QPqZC6y4bzPGGI0ZM0YrVqzQunXr1KpVK6/tkyZN0qeffqodO3ZYiyTNnj1bixYtkiSlpaVp586dXr8UcnJy5HQ6T/sP/Gz1S+NUnVNjdSrMME5Sy5YtlZycfNqlGL744gu1aNFC0slxKikp0bZt26zt69atk8fjscLy2e5Mvk8LFy5Uv379FB8f79XO90nWMaz/OesUHBxszQ7zffLWpEkTxcbGat26dSouLla/fv0kBcY4Vcfj8cjtdqtTp04KDQ1Vbm6utW3v3r3Kz89XWlqapMD4mbP488wN/DajRo0yMTExZv369V6XVjh+/PhPPkc/cbmTXr16mR07dpg1a9aY+Ph4W50C/kvjtG/fPjNt2jTz8ccfm6+//tq88cYb5pxzzjHdu3e3XoNxOmn27NnG6XSaZcuWmS+//NJMnjzZREREmH379ll9evfubS6++GKzefNm8+GHH5o2bdrY6rILNf25+/LLL43D4TCrV68+7TX4PhlTXl5uWrduba644gqzefNms2/fPvPXv/7VOBwOs2rVKut1+D4Z8/zzz5u8vDyzb98+88ILL5i4uDgzYcIEr9ex+zhNmjTJbNiwwXz99dfm008/NZMmTTIOh8O88847xpiTlztp3ry5Wbdunfn4449NWlqaSUtLs54fCD9zpxDszmKSql0WLVr0s8/5cbAzxpgDBw6YPn36mAYNGpgmTZqYP/7xj6aiosK3xdehXxqn/Px80717dxMXF2fCw8NN69atzcSJE70uL2AM43TKjBkzTLNmzUxkZKRJS0szH3zwgdf2f//73+bmm2820dHRxul0mltvvdW6fpsd1HScsrOzTUpKite1xn6M75MxX3zxhbnhhhtMQkKCiYyMNBdccMFplz/h+2TM/fffbxITE01oaKhp06aNefzxx43H4/F6HbuP02233WZatGhhwsLCTHx8vOnRo4cV6owx5sSJE2b06NGmUaNGJjIy0lx//fXm0KFDXq9h95+5UxzGGOPLGUEAAADUDY6xAwAAsAmCHQAAgE0Q7AAAAGyCYAcAAGATBDsAAACbINgBAADYBMEOAADAJgh2AAAANkGwAwAAsAmCHQAAgE0Q7AAAAGzi/wHPo5Whpvm0qgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "original_mean = np.array([batch.air.mean() for batch in bgen_original])\n", + "filtered_mean = np.array([batch.air.mean() for batch in bgen_resampled])\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\n", + "\n", + "line = dict(x=threshold, color='black', linestyle='dashed')\n", + "\n", + "ax1.hist(original_mean)\n", + "ax1.axvline(**line)\n", + "ax1.set_title('original')\n", + "ax2.hist(filtered_mean)\n", + "ax2.axvline(**line)\n", + "ax2.set_title('resampled')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10.17" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/xbatcher/generators.py b/xbatcher/generators.py index 2616ab3..e6f5895 100644 --- a/xbatcher/generators.py +++ b/xbatcher/generators.py @@ -10,8 +10,9 @@ import numpy as np import xarray as xr -PatchGenerator = Iterator[dict[Hashable, slice]] -BatchSelector = list[dict[Hashable, slice]] +Selector = dict[Hashable, slice] +PatchGenerator = Iterator[Selector] +BatchSelector = list[Selector] BatchSelectorSet = dict[int, BatchSelector] @@ -414,6 +415,18 @@ class BatchGenerator: cache_preprocess: callable, optional A function to apply to batches prior to caching. Note: The caching API is experimental and subject to change. + filter_fn: callable, optional + Function that determines whether a batch is removed. This function should + take a ``Dataset`` or ``DataArray`` as its first argument, a Selector + object as its second argument, and return ``True`` for batches that should be + kept. + resample_fn: callable, optional + Function that determines the relative importance of this batch for + resampling. This function should have the same signature as ``filter_fn``, + but return a float. + resample_n: int + Number of batches to keep after resampling. Must be larger than zero and + less than the number of batches available after filtering. Yields ------ @@ -431,6 +444,9 @@ def __init__( preload_batch: bool = True, cache: dict[str, Any] | None = None, cache_preprocess: Callable | None = None, + filter_fn: Callable[..., bool] | None = None, + resample_fn: Callable[..., float] | None = None, + resample_n: int | None = None, ): if input_overlap is None: input_overlap = {} @@ -439,6 +455,9 @@ def __init__( self.ds = ds self.cache = cache self.cache_preprocess = cache_preprocess + self.filter_fn = filter_fn + self.resample_fn = resample_fn + self.resample_n = resample_n self._batch_selectors: BatchSchema = BatchSchema( ds, @@ -449,6 +468,52 @@ def __init__( preload_batch=preload_batch, ) + # Extract the list of selectors for filtering/resampling. Both steps + # can only remove batches, so if this list gets shorter we know + # we have to re-enumerate the selectors property. + if self._batch_selectors.concat_input_dims: + batches = [s for s in self._batch_selectors.selectors[0]] + else: + batches = [s[0] for s in self._batch_selectors.selectors.values()] + + n_initial_batches = len(batches) + + if self.filter_fn is not None: + batches = [b for b in batches if self.filter_fn(self.ds, b)] + if len(batches) == 0: + warnings.warn('Filtering resulted in no batches.') + + if self.resample_fn is not None: + assert ( + self.resample_n is not None + ), 'resample_n must be provided to resample batches.' + assert len(batches) >= self.resample_n, ( + f'Cannot sample {self.resample_n} slices from this dataset ' + f'when there are {len(batches)} available.' + ) + + weight = np.array([self.resample_fn(self.ds, s) for s in batches]) + assert np.any( + weight > 0 + ), 'Sample weight vector does not have any positive values.' + weight = weight / np.sum(weight) + + batches_to_keep = np.random.choice( + len(batches), self.resample_n, replace=False, p=weight + ) + + batches = [batches[i] for i in batches_to_keep] + + # Re-enumerate the list of batches only if filtering or resampling + # occurred. + if len(batches) < n_initial_batches: + if self._batch_selectors.concat_input_dims: + self._batch_selectors.selectors = {0: [b for b in batches]} + else: + self._batch_selectors.selectors = { + i: [b] for i, b in enumerate(batches) + } + @property def input_dims(self): return self._batch_selectors.input_dims diff --git a/xbatcher/tests/test_generators.py b/xbatcher/tests/test_generators.py index 166648c..59ec960 100644 --- a/xbatcher/tests/test_generators.py +++ b/xbatcher/tests/test_generators.py @@ -49,6 +49,30 @@ def sample_ds_3d(): return ds +@pytest.fixture(scope='module') +def sample_filter_fn(): + """ + Sample filter function for testing. + """ + + def myfilter(ds, patch): + return ds.isel(patch).bar.mean() > 4.5 + + return myfilter + + +@pytest.fixture(scope='module') +def sample_resample_fn(): + """ + Sample resample function for testing. + """ + + def myresample(ds, patch): + return ds.isel(patch).foo.mean() + + return myresample + + def test_constructor_dataarray(): """ Test that the xarray.DataArray passed to the batch generator is stored @@ -96,7 +120,7 @@ def test_batch_1d(sample_ds_1d, input_size): validate_generator_length(bg) expected_dims = get_batch_dimensions(bg) for n, ds_batch in enumerate(bg): - assert ds_batch.dims['x'] == input_size + assert ds_batch.sizes['x'] == input_size expected_slice = slice(input_size * n, input_size * (n + 1)) ds_batch_expected = sample_ds_1d.isel(x=expected_slice) xr.testing.assert_identical(ds_batch_expected, ds_batch) @@ -146,7 +170,7 @@ def test_batch_1d_no_coordinate(sample_ds_1d, input_size): validate_generator_length(bg) expected_dims = get_batch_dimensions(bg) for n, ds_batch in enumerate(bg): - assert ds_batch.dims['x'] == input_size + assert ds_batch.sizes['x'] == input_size expected_slice = slice(input_size * n, input_size * (n + 1)) ds_batch_expected = ds_dropped.isel(x=expected_slice) xr.testing.assert_identical(ds_batch_expected, ds_batch) @@ -187,7 +211,7 @@ def test_batch_1d_overlap(sample_ds_1d, input_overlap): expected_dims = get_batch_dimensions(bg) stride = input_size - input_overlap for n, ds_batch in enumerate(bg): - assert ds_batch.dims['x'] == input_size + assert ds_batch.sizes['x'] == input_size expected_slice = slice(stride * n, stride * n + input_size) ds_batch_expected = sample_ds_1d.isel(x=expected_slice) xr.testing.assert_identical(ds_batch_expected, ds_batch) @@ -204,11 +228,11 @@ def test_batch_3d_1d_input(sample_ds_3d, input_size): validate_generator_length(bg) expected_dims = get_batch_dimensions(bg) for n, ds_batch in enumerate(bg): - assert ds_batch.dims['x'] == input_size + assert ds_batch.sizes['x'] == input_size # time and y should be collapsed into batch dimension assert ( - ds_batch.dims['sample'] - == sample_ds_3d.dims['y'] * sample_ds_3d.dims['time'] + ds_batch.sizes['sample'] + == sample_ds_3d.sizes['y'] * sample_ds_3d.sizes['time'] ) expected_slice = slice(input_size * n, input_size * (n + 1)) ds_batch_expected = ( @@ -279,8 +303,8 @@ def test_batch_3d_2d_input(sample_ds_3d, input_size): yn, xn = np.unravel_index( n, ( - (sample_ds_3d.dims['y'] // input_size), - (sample_ds_3d.dims['x'] // x_input_size), + (sample_ds_3d.sizes['y'] // input_size), + (sample_ds_3d.sizes['x'] // x_input_size), ), ) expected_xslice = slice(x_input_size * xn, x_input_size * (xn + 1)) @@ -427,3 +451,123 @@ def preproc(ds): ds_cache = bg[1] xr.testing.assert_equal(ds_no_cache, ds_cache) xr.testing.assert_identical(ds_no_cache, ds_cache) + + +def test_filter_1d(sample_ds_1d, sample_filter_fn): + bg = BatchGenerator(sample_ds_1d, input_dims={'x': 5}) + + bg_filter = BatchGenerator( + sample_ds_1d, input_dims={'x': 5}, filter_fn=sample_filter_fn + ) + + assert len(bg_filter) < len(bg) + + for batch in bg_filter: + assert batch.bar.mean() > 4.5 + + +def test_filter_3d(sample_ds_3d, sample_filter_fn): + bg = BatchGenerator(sample_ds_3d, input_dims={'x': 5, 'y': 5, 'time': 5}) + + bg_filter = BatchGenerator( + sample_ds_3d, input_dims={'x': 5, 'y': 5, 'time': 5}, filter_fn=sample_filter_fn + ) + + assert len(bg_filter) < len(bg) + + for batch in bg_filter: + assert batch.bar.mean() > 4.5 + + +def test_filter_3d_concat(sample_ds_3d, sample_filter_fn): + bg = BatchGenerator( + sample_ds_3d, input_dims={'x': 5, 'y': 5, 'time': 5}, concat_input_dims=True + ) + + bg_filter = BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + filter_fn=sample_filter_fn, + concat_input_dims=True, + ) + + assert bg_filter[0].sizes['input_batch'] < bg[0].sizes['input_batch'] + + assert (bg_filter[0].bar.mean(dim=['x_input', 'y_input', 'time_input']) > 4.5).all() + + +@pytest.mark.parametrize('n', [5, 10]) +def test_resample_1d(sample_ds_1d, sample_resample_fn, n): + bg = BatchGenerator( + sample_ds_1d, input_dims={'x': 5}, resample_fn=sample_resample_fn, resample_n=n + ) + assert len(bg) == n + + +@pytest.mark.parametrize('n', [10, 50, 100]) +def test_resample_3d(sample_ds_3d, sample_resample_fn, n): + bg = BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + resample_fn=sample_resample_fn, + resample_n=n, + ) + assert len(bg) == n + + +def test_filter_prevents_resample(sample_ds_3d, sample_resample_fn): + def strict_filter(*args): + return False + + with pytest.raises(AssertionError, match='Cannot sample 1000 slices'): + BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + resample_fn=sample_resample_fn, + resample_n=1000, + ) + + +def test_error_missing_resample_n(sample_ds_3d): + with pytest.raises(AssertionError, match='resample_n must be provided'): + BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + resample_fn=sample_resample_fn, + resample_n=None, + ) + + +def test_error_large_resample_n(sample_ds_3d, sample_resample_fn): + with pytest.raises(AssertionError, match='Cannot sample 999999999 slices'): + BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + resample_fn=sample_resample_fn, + resample_n=999999999, + ) + + +def test_error_all_zero_resample_weight(sample_ds_3d): + def zero(*args): + return 0 + + with pytest.raises(AssertionError, match='Sample weight vector'): + BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + resample_fn=zero, + resample_n=100, + ) + + +def test_warning_empty_filter(sample_ds_3d): + def strict_filter(*args): + return False + + with pytest.warns(UserWarning, match='no batches'): + BatchGenerator( + sample_ds_3d, + input_dims={'x': 5, 'y': 5, 'time': 5}, + filter_fn=strict_filter, + )