From cbeeb8031e71ab99e92929067187b0c6f5c52e8b Mon Sep 17 00:00:00 2001 From: Alexis Jeandet Date: Tue, 23 Apr 2024 22:54:10 +0200 Subject: [PATCH] Some documentation rework and resampling notebook example with CM Signed-off-by: Alexis Jeandet --- docs/examples/Resampling.ipynb | 183 +++++++++++++++++++++++++++ speasy/products/variable.py | 20 +++ speasy/signal/filtering/__init__.py | 55 ++++++-- speasy/signal/resampling/__init__.py | 21 ++- 4 files changed, 264 insertions(+), 15 deletions(-) create mode 100644 docs/examples/Resampling.ipynb diff --git a/docs/examples/Resampling.ipynb b/docs/examples/Resampling.ipynb new file mode 100644 index 00000000..2c5b282f --- /dev/null +++ b/docs/examples/Resampling.ipynb @@ -0,0 +1,183 @@ +{ + "cells": [ + { + "metadata": {}, + "cell_type": "markdown", + "source": "# Resampling and Interpolation example", + "id": "71d9f362d398650b" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "## Only for Google Colab users:", + "id": "5ba4bb6f40f25b2e" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": "%pip install --upgrade ipympl speasy", + "id": "f71c66b3f08e988f" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "try:\n", + " from google.colab import output\n", + "\n", + " output.enable_custom_widget_manager()\n", + "except:\n", + " print(\"Not running inside Google Collab\")" + ], + "id": "4475cf6b42ee3d5e" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "## For all users:", + "id": "6389dd46e0178598" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-23T20:45:02.685401Z", + "start_time": "2024-04-23T20:44:57.549277Z" + } + }, + "cell_type": "code", + "source": [ + "import speasy as spz\n", + "from speasy.signal.resampling import interpolate\n", + "import numpy as np\n", + "%matplotlib widget\n", + "mms1_products = spz.inventories.tree.cda.MMS.MMS1\n", + "# Use this instead if you are not using jupyterlab yet\n", + "#%matplotlib notebook\n", + "import matplotlib.pyplot as plt\n", + "from datetime import datetime\n", + "import scipy.constants as cst\n" + ], + "id": "56b9bd9bc5d9af1b", + "outputs": [], + "execution_count": 1 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-23T20:52:39.321020Z", + "start_time": "2024-04-23T20:52:39.317849Z" + } + }, + "cell_type": "code", + "source": [ + "def mms1_mirror_cm(start: datetime, stop: datetime) -> spz.SpeasyVariable:\n", + " b, Tperp, Tpara, N = spz.get_data(\n", + " [\n", + " mms1_products.FGM.MMS1_FGM_SRVY_L2.mms1_fgm_b_gsm_srvy_l2,\n", + " mms1_products.DIS.MMS1_FPI_FAST_L2_DIS_MOMS.mms1_dis_tempperp_fast,\n", + " mms1_products.DIS.MMS1_FPI_FAST_L2_DIS_MOMS.mms1_dis_temppara_fast,\n", + " mms1_products.DIS.MMS1_FPI_FAST_L2_DIS_MOMS.mms1_dis_numberdensity_fast,\n", + " ],\n", + " start,\n", + " stop\n", + " )\n", + " anisotropy = Tperp[\"eT_perp\"] / Tpara[\"eT_para\"]\n", + " bt = b[\"Bt\"]\n", + " Pperp = np.multiply(N[\"N\"], Tperp[\"eT_perp\"])\n", + " Pperp, anisotropy = interpolate(bt, [Pperp, anisotropy])\n", + " betaperp = Pperp * 1e6 * cst.e * 2 * cst.mu_0 / (bt * 1e-9) ** 2\n", + " cm = np.multiply(betaperp, anisotropy - 1.)\n", + " cm.columns[0] = \"Mirror mode instability criterion\"\n", + " return cm" + ], + "id": "71b3ff9a04f247c3", + "outputs": [], + "execution_count": 12 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-23T20:52:40.293571Z", + "start_time": "2024-04-23T20:52:39.802095Z" + } + }, + "cell_type": "code", + "source": [ + "plt.figure()\n", + "mms1_mirror_cm(datetime(2021, 1, 1), datetime(2021, 1, 2)).plot()" + ], + "id": "73b10678e0df3daf", + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABsaElEQVR4nO3dd3xN9/8H8NfN3oPsCEEkQhIr9gwlNqWlRYjtZ+/SIqhVtaqKrqCt0VqtKkoVRc3YewVBIiQkiMz7+f1Rud9cNycSuSP33tfz8cjDveece+77nRv3vu9nHZkQQoCIiIiIjIaJrgMgIiIiIu1iAUhERERkZFgAEhERERkZFoBERERERoYFIBEREZGRYQFIREREZGRYABIREREZGRaAREREREaGBSARERGRkWEBSERERGRkWAASERERGRkWgERERERGhgUgERERkZFhAUhERERkZFgAEhERERkZFoBERERERoYFIBEREZGRMctvoxACBw4cwMGDB3H79m2kpaXB1dUVNWrUwDvvvAMfHx/tR0pEREREaqHUAvjy5UvMmTMHPj4+aNOmDf744w88ffoUpqamuHHjBqKiolC+fHm0bdsWR48e1V3URERERPTWlFoA/f39UbduXaxcuRLh4eEwNzdXecCdO3ewbt06dO/eHVOmTMHAgQO1GS8RERERFZNMCCFy71y4cAFBQUGFemBmZibu3LmDSpUqaTI+IiIiIlIzpQKQiIiIiAyf5CzgnJwcpfvHjh3DP//8g6ysLG3ERUREREQaolIAxsfHo1GjRrC0tETTpk3x5MkTtG/fHvXr10ezZs0QFBSE+Ph43URLRERERMWmUgB+9NFHEEJg69at8PT0RPv27ZGamoq4uDjcuXMH7u7umD17tm6iJSIiIqJiUxkD6OXlhS1btqBevXpITk6Gi4sL9uzZgxYtWgAA9u3bhwEDBuDmzZu6ipmIiIiIikGlBfDJkyfw9vYGAJQqVQo2NjYoV66cYn/FihXZBUxERESkx1QKQDc3N6UCb/jw4ShVqpTi/pMnT2Bra6u9CImIiIhIrVQKwOrVq+PIkSOK+/PmzVMqAA8dOoSQkBDtRUhEREREalXkdQBPnDgBa2vrQi8YTUREREQlCxeCJiIiIjIykgtB5+Xg4IBbt25pPhoiIiIi0rhCFYBsJCQiIiIyHIUqAImIiIjIcBSqAOzVqxccHBw0Hw0RERERaRwngRAREREZGckWwLS0NO1GQkRERERakW8BmJycrLj2LxEREREZFpUC8MGDB2jSpAlq1Kihm4iIiIiISKOUCsDr16+jQYMGqFWrFpYvX667qIiIiIhIY5QmgXh4eKBx48b4+eefYWLCFWKIiIiIDJFSlff8+XN4e3uz+CMiIiIyYGZ57+zZswft27eHvb09Pv30U91FRUREREQao1QA1q9fHwcOHEDr1q1ha2uLSZMm6S4yIiIiItKIfBeCvn37NsLDw3H16lXdREVEREREGiN5JZBHjx7B1dVV+xERERERkUbxUnBERERERobTfYmIiIiMjJnUjqSkJEybNg379u1DYmIi5HK50v7k5GRtxEdEREREaiZZAPbq1Qs3b95E//794e7uDplMpt3IiIiIiEgjJMcA2tvb49ChQ6hWrZr2oyIiIiIijZEcA1i5cmW8fPlSu9EQERERkcZJtgCeOHECkyZNwrRp0xAUFARzc3Ol/Q4ODtqKkYiIiIjUSHIMoJOTE1JSUtC8eXOl7UIIyGQy5OTkaCM+IiIiIlIzyQKwZ8+esLCwwLp16zgJhIiIiMiASHYB29jY4PTp0wgICNB+VERERESkMZItgKGhoYiLi2MBWAC5XI4HDx7A3t6eLaRERER6QgiBZ8+ewcvLCyYmxnlNDMkWwI0bN2L69OmYMGECgoODVSaBhISEaCvGEuvevXvw8fHRdRhERET0FuLi4lCmTBldh6ETkgVgfhWxTCbjJJA8UlJS4OTkhLi4OM6KJiIi0hOpqanw8fHB06dP4ejoqOtwdEKyCzg2Nla7keih3G5fBwcHFoBERER6xpiHb0kWgOXKldNuJERERESkFZIjH+fOnYvo6GiV7dHR0fjss880HRcRERERaYhkAfj111+jcuXKKturVq2KlStXajouIiIiItIQyQIwISEBnp6eKttdXV0RHx+v6biIiIiISEMkC0AfHx8cPnxYZfvhw4fh5eWl6biIiIiISEMkJ4EMGDAAo0ePRlZWluJ6wHv37sXEiRMxbtw4bcZIRERERGokWQBOnDgRycnJGDp0KDIzMwEAVlZW+OijjzB58mRtxkhEREREaiS5EHSu58+f4/Lly7C2tkalSpVgaWmpvehKuNTUVDg6OiIlJYXrABIREekJfn4X0AKYy87ODrVr19ZONERERESkcUqTQIYMGYK4uLhCPfDnn3/G2rVrNRUXEREREWmIUgugq6srgoKC0KBBA3Ts2BGhoaHw8vKClZUVnjx5gkuXLuHQoUPYsGEDvL298c033+guciIiIiJ6KypjABMTE/H9999jw4YNuHDhgtLB9vb2eOeddzBo0CC0atVK27GWOBxDQEREpH/4+f2GSSBPnz7FnTt38PLlS7i4uKBixYpGfeHk1/EPiIjexvZzD5CYmoF+jcrrOhQio8TP7zdMAnFycoKTk5P2oiEiMgLD150GADSu5IJK7va6DoeIjJDklUCIiEizkl9k6joEIjJSLACJiIiIjAwLQCIiIiIjwwKQiIiIyMhIFoDNmzfH06dPVbanpqaiefPmmo6LiIiIiDREsgDcv38/MjNVByinp6fj4MGDmo6LiIiIiDREZRmYc+fOKW5funQJCQkJivs5OTnYtWsXvL29tRchEREREamVSgFYvXp1yGQyyGSyfLt6ra2t8eWXX2orPiIiIiJSM5UCMDY2FkIIVKhQAcePH4erq6tin4WFBdzc3GBqaqrtOImIiIhITVQKwHLlygEA5HK5LuIhIiIiIg2TnASyZs0a/PHHH4r7EydOhJOTExo0aIA7d+5oKz4iIiIiUjPJAnDOnDmwtrYGABw5cgTLli3D/Pnz4eLigjFjxmgzRiIiIiJSI5Uu4FxxcXHw8/MDAPz666947733MGjQIDRs2BDNmjXTZoxERAZJ6DoAIjJaki2AdnZ2SEpKAgDs3r0b77zzDgDAysoKL1++1F6ERERERKRWki2ALVu2xIABA1CjRg1cu3YN7dq1AwBcvHgRvr6+2oyRiIiIiNRIsgXwq6++Qv369fHo0SNs3rwZpUuXBgDExMTgww8/1GaMRERERKRGki2ATk5OWLZsmcr2GTNmaDomIiIiItIgyQIwV1paGu7evatyXeCQkBBNxkVEREREGiJZAD569AiRkZHYtWtXvvtzcnI0GRcRERERaYjkGMDRo0fj6dOnOHr0KKytrbFr1y6sWbMGlSpVwrZt27QbJRERERGpjWQL4N9//43ffvsNtWvXhomJCcqVK4eWLVvCwcEBc+fOVcwKJiIiIiL9ItkC+OLFC7i5uQEASpUqhUePHgEAgoODcerUKe1FSERERERqJVkABgQE4OrVqwCA6tWr4+uvv8b9+/excuVKeHp6ajNGIiIiIlIjyS7g0aNHIz4+HgAQFRWF8PBwrF27FhYWFli9erU2YyQiIiIiNVIqAFNTU+Hg4AAA6Nmzp2J7jRo1cPv2bVy5cgVly5aFi4uL9iMlIiIiIrVQ6gJ2dnZGYmIiAKB58+Z4+vSpYp+NjQ1q1qzJ4o+ISE2E0HUERGSslApAOzs7JCUlAQD279+PrKwsXcVFRERERBqi1AX8zjvvICwsDIGBgQCAd999FxYWFvk+8O+//9ZOhERERESkVkoF4E8//YQ1a9bg5s2bOHDgAKpWrQobGxvdRUdEREREaqdUAFpbW2PIkCEAgJMnT+Kzzz6Dk5OTrmIjIiIiIg2QXAdw3759hSr+HBwccOvWLbUFNHfuXNSuXRv29vZwc3ND586dFesR5hJCYPr06fDy8oK1tTWaNWuGixcvKh2TkZGBESNGwMXFBba2tujYsSPu3bundMyTJ08QEREBR0dHODo6IiIiQmniCxEREZEhkiwAC0uoeRrbgQMHMGzYMBw9ehR79uxBdnY2WrVqhRcvXiiOmT9/PhYtWoRly5bhxIkT8PDwQMuWLfHs2TPFMaNHj8bWrVuxYcMGHDp0CM+fP0f79u2Rk5OjOKZHjx44c+YMdu3ahV27duHMmTOIiIhQaz5EREREJY4oJjs7O3Hz5s3inkZSYmKiACAOHDgghBBCLpcLDw8PMW/ePMUx6enpwtHRUaxcuVIIIcTTp0+Fubm52LBhg+KY+/fvCxMTE7Fr1y4hhBCXLl0SAMTRo0cVxxw5ckQAEFeuXClUbCkpKQKASElJUVu+RGT4yn20XZT7aLv498ZjXYdCZJT4+S1EsVsANS0lJQV4dT1iAIiNjUVCQgJatWqlOMbS0hJNmzbFv//+CwCIiYlBVlaW0jFeXl4ICgpSHHPkyBE4Ojqibt26imPq1asHR0dHxTFEREREhkjyUnAlgRACY8eORaNGjRAUFAQASEhIAAC4u7srHevu7o47d+4ojrGwsICzs7PKMbmPT0hIgJubm8pzurm5KY55XUZGBjIyMhT3U1NTi50jERERkbYVuwVQJpOpJ5J8DB8+HOfOncP69evf+LxCiDfG8vox+R1f0Hnmzp2rmDDi6OgIHx+fImRDREREVDKUuEkguUaMGIFt27Zh3759KFOmjGK7h4cHkKclMFdiYqKiVdDDwwOZmZl48uRJgcc8fPhQ5XkfPXqk0rqYa/LkyUhJSVH8xMXFqSFTIiIiIu2SLAD3799fqBPs3LkT3t7eagtICIHhw4djy5Yt+Pvvv1G+fHml/eXLl4eHhwf27Nmj2JaZmYkDBw6gQYMGAIBatWrB3Nxc6Zj4+HhcuHBBcUz9+vWRkpKC48ePK445duwYUlJSFMe8ztLSEg4ODko/RERvS4AXAyYi3ZAcA9i6dWt4e3ujb9++6NOnj2R3Z6NGjdQa0LBhw7Bu3Tr89ttvsLe3V7T0OTo6wtraGjKZDKNHj8acOXNQqVIlVKpUCXPmzIGNjQ169OihOLZ///4YN24cSpcujVKlSmH8+PEIDg7GO++8AwAIDAxE69atMXDgQHz99dcAgEGDBqF9+/YICAhQa05EREREJYlkC+CDBw8watQobNmyBeXLl0d4eDh++eUXZGZmajSgFStWICUlBc2aNYOnp6fi5+eff1YcM3HiRIwePRpDhw5FaGgo7t+/j927d8Pe3l5xzOLFi9G5c2d069YNDRs2hI2NDX7//XeYmpoqjlm7di2Cg4PRqlUrtGrVCiEhIfjxxx81mh8RERGRrslEIQbxnTlzBtHR0Vi/fj3kcjl69uyJ/v37o1q1atqJsoRKTU2Fo6MjUlJS2B1MRIXmO+kPAMC6gXXRoKKLrsMhMjr8/C7kJJDq1atj0qRJGDZsGF68eIHo6GjUqlULjRs3VrkEGxERERGVbAUWgFlZWdi0aRPatm2LcuXK4c8//8SyZcvw8OFDxMbGwsfHB++//772oiUiIiKiYpOcBDJixAjF+nu9evXC/PnzFYsxA4CtrS3mzZsHX19f7URKRERERGohWQBeunQJX375Jbp27QoLC4t8j/Hy8sK+ffs0GR8RERERqZlkF3BUVBTef/99leIvOzsb//zzDwDAzMwMTZs21XyURERERKQ2kgVgWFgYkpOTVbanpKQgLCxM03ERERERkYZIFoBS18RNSkqCra2tpuMiIiIiIg1RGQPYpUsXAIBMJkNkZCQsLS0V+3JycnDu3DnJS6URERERUcmnUgA6OjoCr1oA7e3tYW1trdhnYWGBevXqYeDAgdqNkoiIiIjURqUAXLVqFQDA19cX48ePZ3cvEZGmvPE6TEREmiG5DExUVJR2IyEiIiIirVAqAGvWrIm9e/fC2dkZNWrUyHcSSK5Tp05pIz4iIiIiUjOlArBTp06KSR+dO3fWVUxEREREpEFKBWDebl92ARMREREZJsl1AOPi4nDv3j3F/ePHj2P06NH45ptvtBUbEREREWmAZAHYo0cPxXV+ExIS8M477+D48eP4+OOPMXPmTG3GSERERERqJFkAXrhwAXXq1AEA/PLLLwgODsa///6LdevWYfXq1dqMkYiIiIjUSLIAzMrKUkwI+euvv9CxY0cAQOXKlREfH6+9CImIiIhIrSQLwKpVq2LlypU4ePAg9uzZg9atWwMAHjx4gNKlS2szRiIiIiJSI8kC8LPPPsPXX3+NZs2a4cMPP0S1atUAANu2bVN0DRMRERGR/pG8EkizZs3w+PFjpKamwtnZWbF90KBBsLGx0VZ8RERERKRmkgUgAJiamioVf3h1jWAiIiIi0l+SXcAPHz5EREQEvLy8YGZmBlNTU6UfIiIqHqHrAIjIaEm2AEZGRuLu3buYOnUqPD09C7wuMBERERHpD8kC8NChQzh48CCqV6+u3YiIiIiISKMku4B9fHwgBDsoiIg07UVGNhJS0nUdBhEZEckCcMmSJZg0aRJu376t3YiIiIxM7dl/od7cvYhPeanrUIjISEh2AXfv3h1paWmoWLEibGxsYG5urrQ/OTlZG/ERERm8tMwcAMDx2GR0qu6t63CIyAhIFoBLlizRbiREREaOo26ISFskC8A+ffpoNxIiIiMnuDAMEWmJ5BhAALh58yamTJmCDz/8EImJiQCAXbt24eLFi9qKj4jIaLAFkIi0RbIAPHDgAIKDg3Hs2DFs2bIFz58/BwCcO3cOUVFR2oyRiMgosAAkIm2RLAAnTZqEWbNmYc+ePbCwsFBsDwsLw5EjR7QVHxGR0dD3+i87R44fj9zG9YfPdB0KEb2BZAF4/vx5vPvuuyrbXV1dkZSUpOm4iIiMjr6vvbr22F1M/e0iWi7+R9ehENEbSBaATk5OiI+PV9l++vRpeHtzmQIiInXT7/IPOBP3VNchEFEhSRaAPXr0wEcffYSEhATIZDLI5XIcPnwY48ePR+/evbUbJRGRAVJp8NP3CpCI9IZkATh79myULVsW3t7eeP78OapUqYImTZqgQYMGmDJlinajJCIyAvq+DIxM1wEQUaFJrgNobm6OtWvXYubMmTh9+jTkcjlq1KiBSpUqaTdCIiIjoedDAIlIj0gWgLkqVqyIihUraicaIiIjxvqPiLRFqQAcO3ZsoR+4aNEiTcRDRGS09L4FkH3ARHpDqQA8ffq00s6YmBjk5OQgICAAAHDt2jWYmpqiVq1a2o2SiMgI6PsYQCLSH0oF4L59+xS3Fy1aBHt7e6xZswbOzs4AgCdPnqBv375o3Lix9iMlIjJw+t4CKGMTIJHekJwFvHDhQsydO1dR/AGAs7MzZs2ahYULF2orPiIio6Hn9R8R6RHJAjA1NRUPHz5U2Z6YmIhnz3iZHyIitdPzJkAZGwCJ9IZkAfjuu++ib9++2LRpE+7du4d79+5h06ZN6N+/P7p06aLdKImIjIB+l39EpE8kl4FZuXIlxo8fj169eiErK+u/g83M0L9/f3z++efajJGIyCiUlAbA6EOxuP/0Jaa0C4SMzXpEBinfAjAnJwcnTpzArFmz8Pnnn+PmzZsQQsDPzw+2trbaj5KIyAiIElIBztx+CQDQubo3gss4FvpxLBWJ9Ee+BaCpqSnCw8Nx+fJllC9fHiEhIdqPjIjIyJSM8u9/nmVk6ToEItIQyTGAwcHBuHXrlnajISIyIq+v+1dCGgD/p4jxsLeYSH9IFoCzZ8/G+PHjsX37dsTHxyM1NVXph4iI1EvP6z8i0iOSk0Bat24NAOjYsaPSIGAhBGQyGXJycrQTIRER6USJa5EkIrWRLADzXhWEiIg0r6RMAslV1EvT8UogRPpDsgBs2rSpdiMhIiIiIq2QHAMIAAcPHkSvXr3QoEED3L9/HwDw448/4tChQ9qKj4hIL/1z7RFaL/kH5+49LfRjitMAeCfpBQb9cBKn7z55+5MUMx5OAiHSH5IF4ObNmxEeHg5ra2ucOnUKGRkZAIBnz55hzpw52oyRiEjv9I4+jisJzxC56kShH1PULte8Bv8Yg92XHuLd5f++9TlU4yEiQyVZAM6aNQsrV67Et99+C3Nzc8X2Bg0a4NSpU9qKj4hIrz1Pzy70scVpAYxLTnv7B0so6phEtgAS6Q/JAvDq1ato0qSJynYHBwc8fVr4Lg0iImNm8upddteFeAz5MQYpL6UXVy5Oi5uJBqovtgASGS7JSSCenp64ceMGfH19lbYfOnQIFSpU0EZsRER6L3dm7JCf/us58Xa2ljy2OC2AJiYaaH5jBUhksCRbAAcPHoxRo0bh2LFjkMlkePDgAdauXYvx48dj6NCh2o2SiEhPvV6XPXmRKXlsccYAaqb+K2o87AMm0heSLYATJ05ESkoKwsLCkJ6ejiZNmsDS0hLjx4/H8OHDtRslEZGeer1rtqCWuvxaAJ+8yMSx2GS0CHSDuan0wg0a6QJmCyCRwSpwGZjZs2fj8ePHOH78OI4ePYpHjx7h008/1V50RET67rW6zPS1Qu15RsGTRLqu/BdDforBsr9vFPw0Gmh8m7jpXJGO5yQQIv1RYAEIADY2NnB3d4eXlxfs7Oy0ExURkYF4vSbK2wIoBPDr6ft57qs2ud169AIAsON8fMHPU8Tq68HTl+gdfRz7riZKHpNUQHc1Eek3yQIwOzsbU6dOhaOjI3x9fVGuXDk4OjpiypQpyMqSnsVGRET/k/raMjCv9wDnLfqK2uWa97FFHQP48dbz+OfaI/QtwjqFb8IGQCL9ITkGcPjw4di6dSvmz5+P+vXrAwCOHDmC6dOn4/Hjx1i5cqU24ySiEkAuF5qZbWpETAv4/WXlyCX3Zb62Ly0zG22/OIh6FUpjXteQIo8BfPQso9DHZuXICxx/SET6R/J/9Pr167F69WoMHjwYISEhCAkJweDBgxEdHY3169drN0oi0rnrD5+h2ozd+GpfwWPRqGAFFWpHbyVL7ruT9L+FnvddTUTw9N24nZSGDSfi3nhevCre87YYFrZe3HgyDpU+2YldFxIU5yEi/SdZAFpZWamsAQgAvr6+sLCw0HRcRFTCzPj9Ep5lZOPzP6/qOhS99nqhlrecOn47Gb+fffDGc/RddQI5rxViJgU00N17koYKH+9Ai0UH3nh1j9eff8KriSBDforBiPWn0Xj+PqRl5j9xJW9q+68mYt2xu0h6XviWRiLSHsm3jGHDhuHTTz9VXAMYADIyMjB79mwuA0NkhDjDUz3y9qT2jj6Oab9dVNo/Yv1pZGZLdwVLKagFsNFn+4BXE0peZuUAABJS8i/MRqw/LXme388+wP2nL7Hn0sM3xhO56gQ+3nq+SNdCJiLtURoD2KVLF6Wdf/31F8qUKYNq1aoBAM6ePYvMzEy0aNFCu1ESEZVwWTly9Pz2GELKOGJK+ypK+87fS1HcLswYSv8pO3F7XrtCP/fLzBxkFbJozG0AfJynZW7CxrOo5uOEXvXKFeocozacQSM/F5S2s1TaLstnGsj5+ykq24hI95QKQEdHR6WdXbt2Vbrv4+OjnaiIqMS59CBV1yGUaH9fScTx28k4fjsZw5v7Ke3rsOyQ4vbr6wBKSc/KgZW5aaGODZy2q9BxyvPpAt4Ycw8bY+7h/dAyhT5PrVl/4ZfB9VGnfKk3Hvv9oVj0a+hb5KVqiEhzlArAVatW6S4SHVm+fDk+//xzxMfHo2rVqliyZAkaN26s67CISpyMt+iWLIyHqen48cgd9KxXFp6O0tfJLenydtvO+P2S5HHHY6UneuR1I/E5qno5qGzffTHhLSP8T0FzOHp+e6xI5+r29REcnBgGn1I2QAHDBD7dfgkVXW3RLMCtSOcnIs0x6nn9P//8M0aPHo1PPvkEp0+fRuPGjdGmTRvcvXtX16ERlThvmjzwtgb9cBLL9t1A7++Pa+T82iCE8lVzC2otPXnnSaHO2f7LQ5i85bzSts0x9zDox5i3jhMAxv1yFn9JjOErbGx5XUl4prid9Fx64ejfzjxA568O489iFrBEpB4yIfGunpSUhGnTpmHfvn1ITEyEXK787T85uXDfYkuyunXrombNmlixYoViW2BgIDp37oy5c+e+8fGpqalwdHRESkoKHBxUv6kT6drzjGzM2XEZ7UM80aCiS7HOVXXaLrzI/G8CgdT4tKwcOfZceoiNJ+OQlSPQsboXuoUWPHTEd9IfittS583IzoGFqYnWuxCFELj1+AVc7S3hYGWe7zF3kl6g6ef7YWFmomgFrOxhr1QYaVtkA19M71gVsY9fIGrbRfxz7ZHGnsunlDX+HN0ET9Oy0GDe34V6zJp+ddDU31VjMRG9CT+/CygA27Rpg5s3b6J///5wd3dXeePt06ePtmLUiMzMTNjY2GDjxo149913FdtHjRqFM2fO4MCBA288h6b+gGLuJCPm1TdxIfKf3Ze3veFpWhYcrc0hk/1vELamPifz/rXkjSG/wd+vk8k0d3F5TZ47l4CAEP8t2yEDkJyWCRdbSyi3/byKB7I3/n4K8xr993xCkZvI8xrk3Y48LXS5MQLAoj3XFPuntAuEEEX725j1x+VCHzuyuR+WvuF6tZr2Xq0y2BRzT2nbyl41MeSnU5KP8S1tg9t51tgriK4LO0MzsrkfHqSkq7xmeZmayNA+xBMmMhm2nr6P5pXd0KBiabXFwHGJJV+Nsk6oWdZZredkAVhAAWhvb49Dhw4pZgAbmgcPHsDb2xuHDx9GgwYNFNvnzJmDNWvW4OpV1bXOMjIylJbFSU1NhY+Pj9r/gBbvuYYv9l5X2/mIiIj01eh3KmH0O/5qPScLwAIuBVe5cmW8fPlSu9HowOvf/oQQkt8I586dixkzZmg8pkBPB3Sp4f2/mF79K3vttqJRSQCPnmfA1c4yn7YozSrsd2d5nt9rnrAV91+/nd/rIPVcb8o57+PEa8/3Jvk+p+y/Vlcn6/y7BPM+T35xFuX5leKQ5d6WvWrtfXVf9r9tULQsyl61igqsPx6HQE8HBLjbQSaTSY7le/33vfX0/SJGSVQ07YI98cf5+DceZ2VuAmtzUzxJ++869O/meX/MVdAYVak9hWkRL2qrOalfgLu9rkMwSJItgCdOnMCkSZMwbdo0BAUFwdxc+cNO3yvmt+kC1lYLIFFJcCPxOd5Z9N//gwszwvHbmfv4ZOsFxX07y/99f4z4/hgOXn8MABjUpAKGN/fLd8xcVo4cOXKBylOVly3xcLCCiQyIbOiLfg3Lw8zUBLsuJODaw2dK3dgXZoTjzwsJGLfxLBZ3r4Zxv5xVmdVa1csBvw5riNWHb+Pxiwx8feAWAKCK53/bn6ZlotFn+1SurZvX0ckt0HXFv2hQsTSCvB3Ru345yS+G07ddxOp/bytt61zdC7+eefMVPTSpU3Uv/JZPDONb+WN480pIfpGJLssPK7q/R7aoBFc7C0x9bWHqguwb3wxhC/YDAG7NaQsTExl+O3MfozacUTnW09EKjtbmii70m3PaFnhdZCJNYgtgAS2ATk5OSElJQfPmzZW257bM5OTkaCM+jbGwsECtWrWwZ88epQJwz5496NSpU76PsbS0hKWlZb77iAyNn5sdto9oBDd7S9hZmildAzZv8QcAP/avW6hzmpuawNwUuD67DU7EJsPV3hK/nXmAPg184Wqv/H+rdZAHWgd5YFCTCniRka1YdLhrrTLoWuu/9ereCXTHvzeT0NTfFdvPxePWo+cY1yoApiYyDGxSAQAwuU2g0nndHKxwbXYbHL2VhFl/XEKwtxPWH1ee+e/haIXDk5Tf+6TkXiM3ryUf1EBqejZO3k7GoUnN4WBlDiEETsc9RUUXO1hZmCBgSuHX7iuqLUMboGZZZ6UCsLavM6a1r4rgMv+t91rK1gK+LraKAnBsy/+62D6sUxZ+n+yUPPfm/2uAiq62cLL575Kgvw5rCDMT2RsXuHawMkegp4OiAGTxR6RbkgVgz549YWFhgXXr1uU7CcQQjB07FhEREQgNDUX9+vXxzTff4O7duxgyZIiuQyMqEYK8/7c4/OvXni0Oc1MTNPD7b1by+PCAAo+1MjeVXBDZ3soc4VU9gFcTQIqiXoXS2D6iMe4mpakUgEVhaZ7/alrf9wlFtlzA/NW132QymdJA9h/61UHv6P+Wvnkn0B1/XX7z5dVyzehYFVHb/tdSV6d8KaX1Bat4/tei4WRjjqevuk03Dmmgcp78XlKzPNeqC3C3x9WHr096EYriDwCq+zgp75X4MzHAjxAivSZZAF64cAGnT59GQEDBb876rHv37khKSsLMmTMRHx+PoKAg7NixA+XKFe5ySETGJEfbA0y1JFtevAWuG1R0wZ0k1QJSJpPB3FS66mni74rjn7SApakpbCxNUamAVrfX9Wngq1QAbhhYDxU+3qG4n1swL+5WHX1Xn0CExCXe3rS2Y35F25u+B0h9USho/CkRaZ9kARgaGoq4uDiDLgABYOjQoRg6dKiuwyAq8SzMDHPd+OI2bLYP8XzrFkQ3eysAQHae8YjlStvgTgHL0nzbO1TpvqejlWT3a1hlN5yd1goO1vm/1Xs7FXzllfx6ft5Uw+VITTIq+GFEpGWSBeCIESMwatQoTJgwAcHBwSqTQEJCQrQRHxGVEO/XKoNtZ+4b3AK+FV1t0biSi2ISS1E19HPB7jFNcO9JGvqtPqn2+HKFV3XHFx/UUOkOr1+h4DXxHG2kZ6tPalMZ6Vk5ijGVr8uvaHtTK57UfhPD/P5ApLckC8Du3bsDAPr166fYltuEbwiTQIioaKzMTfMdR6bvZDIZfuxfF00/31dgy1tB/N3t4edqhxaV3eDnbqf2GPFq0kR+YyHlrwquuuVL4VhscpEWSXayscCSD2pI7s+vC/hNDaZSk6sLs1g8EWmPZAEYGxur3UiIiHSouJNcTExk+D6ydrHjyFsm7R/fDM1eLbMiVUDlRr2yVy38fu4BOlbzKnYMuUxkMqVL3CFPwSlFaj8n/RKVLJIFICdCEJExkatxlrO62FnleYuWKKByw3a2tUDv+r5qfX6ZDPhjRCMM+OHk/1pH3/Br6lDNC1N+vZD/yYioxJAsAH/44YcCH9i7d29NxENEpBNSkxd0qbStRb63tUUGoJK7PZZ+UAOdvjoMFKIL2NHaHEu6V8fon5UXgzaRAeVdNNM9TkRFJ1kAjho1Sul+VlYW0tLSYGFhARsbGxaARGRQPqxTFkv+uv7GSRXalHcWrlT7mUbb1V49f7U8a/29qQsYr65C8jQtEzXKOisKRxmAwU0rIC0zG+9Ucddg0ERUGJIF4JMnT1S2Xb9+Hf/3f/+HCRMmaDouIiKtGh7mhzrlS6ksbGzM8hu3V5iGUplMhsiG5V8713+TWCa3DZR8HBFpT5Em5leqVAnz5s1TaR0kItJ3ZqYmaFDRBTYWkt+LtULqqkvS2zUYSz7b3rajnEMAiUqWIq/MZGpqigcPdHuRcyIiQ1WS6qT8F4IuWgmYuyxNhJonqBBR8Uh+1d22bZvSfSEE4uPjsWzZMjRs2FAbsRER0SshZRwLcZR6vW0XcF6r+9bB3eQX8HOzV1tcRFR8kgVg586dle7LZDK4urqiefPmWLhwoTZiIyIyervHNMGpO0/Qubp3vvs12WJYxtlGZZsoYiewhZkJiz+iEkiyAJQX8wLpRERUfP7u9vB3124BtW5AXWw4EYep7auo7CuBq+UQ0VvQ7WhnIiJSblMrAYMAG/i5oIGfi67DICINkiwAc3JysHr1auzduxeJiYkqLYJ///23NuIjIjJ4xWlVk5odTERUkAIXgl69ejXatWuHoKAgvskQEWlIUcfVEREVl2QBuGHDBvzyyy9o27atdiMiIjIyxWoBVGcgRGQ0JNcBtLCwgJ+fn3ajISIyQhamJrA2NwUAeDpa6TocIjICkgXguHHj8MUXXxR50U8iIioaExMZTk9ricszW8PMpGjr89etUEpjcRGR4ZLsAj506BD27duHnTt3omrVqjA3N1fav2XLFm3ER0RkFKxetQAW9iv3wYlhOHX3CTqEeGk0LiIyTJIFoJOTE959913tRkNEZOQK2+viU8oGPqVUF2omIioMyQJw1apVhTrB4cOHERoaCktLS3XGRUREREQaUrTBJvlo06YN7t+/r55oiIiIiEjjil0AcpIIEZH68C2ViLSh2AUgERGpDxeFJiJtYAFIRFSCsAWQiLSBBSARUQnCApCItKHYBSCvEUxERESkXzgJhIioBOEYQCLSBsl1AF++fAkhBGxs/lto9M6dO9i6dSuqVKmCVq1aKY579uyZdiIlIjIC/E5NRNog2QLYqVMn/PDDDwCAp0+fom7duli4cCE6deqEFStWaDNGIiKjEeTtqOsQiMgISBaAp06dQuPGjQEAmzZtgru7O+7cuYMffvgBS5cu1WaMRERGY1wrf4xqUQk7RzXWdSj54rBvIsMg2QWclpYGe3t7AMDu3bvRpUsXmJiYoF69erhz5442YyQiMho2FmYY09Jf12EQkYGTbAH08/PDr7/+iri4OPz555+KcX+JiYlwcHDQZoxEREREpEaSBeC0adMwfvx4+Pr6om7duqhfvz7wqjWwRo0a2oyRiIiIiNRIsgv4vffeQ6NGjRAfH49q1aoptrdo0QLvvvuutuIjIiIiIjWTLAABwMPDAx4eHkrb6tSpo+mYiIiIiEiDlArALl26YPXq1XBwcECXLl0KfOCWLVs0HRsRERERaYBSAejo6Ki4tJujI9eiIiIiIjJESgXgqlWr8r1NRERERIaj2NcCJiIiIiL9otQCWKNGDUUX8JucOnVKUzERERERkQYpFYCdO3dW3E5PT8fy5ctRpUoVxRqAR48excWLFzF06FDtR0pEREREaqFUAEZFRSluDxgwACNHjsSnn36K14+Ji4vTXoREREREpFaSYwA3btyI3r17q2zv1asXNm/erOm4iIiIiEhDJAtAa2trHDp0SGX7oUOHYGVlpem4iIiIiEhDJK8EMnr0aPzf//0fYmJiUK9ePeDVGMDo6GhMmzZNmzESERERkRpJFoCTJk1ChQoV8MUXX2DdunUAgMDAQKxevRrdunXTZoxEREREpEYFXgu4W7dubyz21q9fj44dO8LW1lbdsRERERGRBhR7IejBgwfj4cOH6omGiIiIiDSu2AWgEEI9kRARUYknQ+EuFkBEJRsvBUdERERkZFgAEhERERkZFoBERERERoYFIBEREZGRKXYBWK5cOZibm6snGiIiIiLSuALXASyMCxcuqCcSIiIiItKKIrcAnj17FqamppqJhoiIiIg07q26gLn2HxEREZH+UukC7tKlS4EPSElJgUzGhUCJiIiI9JVKAfj777+jZcuWcHd3z/cBOTk52oiLiIiIiDREpQAMDAxE165d0b9//3wfcObMGWzfvl0bsRERERGRBqiMAaxVqxZOnTol+QBLS0uULVtW03ERERERkYaotACuXLmywG7ewMBAxMbGajouIiIqgQQ4CZDIEKgUgJaWlrqJhIiIiIi0QnIZmFOnTuH8+fOK+7/99hs6d+6Mjz/+GJmZmdqKj4iIiIjUTLIAHDx4MK5duwYAuHXrFj744APY2Nhg48aNmDhxojZjJCIiIiI1kiwAr127hurVqwMANm7ciCZNmmDdunVYvXo1Nm/erM0YiYiIiEiNJAtAIQTkcjkA4K+//kLbtm0BAD4+Pnj8+LH2IiQiohJDBl4IgMgQSBaAoaGhmDVrFn788UccOHAA7dq1AwDExsZKLhJNRERERCWfZAG4ZMkSnDp1CsOHD8cnn3wCPz8/AMCmTZvQoEEDjQRz+/Zt9O/fH+XLl4e1tTUqVqyIqKgolUknd+/eRYcOHWBrawsXFxeMHDlS5Zjz58+jadOmsLa2hre3N2bOnKlyDeMDBw6gVq1asLKyQoUKFbBy5UqN5EVERERUkqgsA5MrJCREaRZwrs8//xympqYaCebKlSuQy+X4+uuv4efnhwsXLmDgwIF48eIFFixYALy6FF27du3g6uqKQ4cOISkpCX369IEQAl9++SUAIDU1FS1btkRYWBhOnDiBa9euITIyEra2thg3bhzwqiWzbdu2GDhwIH766SccPnwYQ4cOhaurK7p27aqR/IiIiIhKApl4vVksH8+fP1eMB8zl4OCgybgUPv/8c6xYsQK3bt0CAOzcuRPt27dHXFwcvLy8AAAbNmxAZGQkEhMT4eDggBUrVmDy5Ml4+PChYl3DefPm4csvv8S9e/cgk8nw0UcfYdu2bbh8+bLiuYYMGYKzZ8/iyJEjhYotNTUVjo6OSElJ0drvg4hIF3wn/QEAWBVZG2GV3XQdDlGx8PO7gC7g2NhYtGvXDra2tnB0dISzszOcnZ3h5OQEZ2dnrQWYkpKCUqVKKe4fOXIEQUFBiuIPAMLDw5GRkYGYmBjFMU2bNlVa1Do8PBwPHjzA7du3Fce0atVK6bnCw8Nx8uRJZGVl5RtLRkYGUlNTlX6IiIiI9I1kF3DPnj0BANHR0XB3d4dMpv2ZXzdv3sSXX36JhQsXKrYlJCSoTEJxdnaGhYUFEhISFMf4+voqHZP7mISEBJQvXz7f87i7uyM7OxuPHz+Gp6enSjxz587FjBkz1JojERERkbZJFoDnzp1DTEwMAgICiv0k06dPf2PhdOLECYSGhiruP3jwAK1bt8b777+PAQMGKB2bXzEqhFDa/voxuT3dRT0mr8mTJ2Ps2LGK+6mpqfDx8SkwLyIiIqKSRrIArF27NuLi4tRSAA4fPhwffPBBgcfkbbF78OABwsLCUL9+fXzzzTdKx3l4eODYsWNK2548eYKsrCxFi56Hh4eiNTBXYmIikKclUOoYMzMzlC5dOt8YLS0tea1kIiIi0nuSBeB3332HIUOG4P79+wgKCoK5ubnS/pCQkEI/iYuLC1xcXAp17P379xEWFoZatWph1apVMDFRHqZYv359zJ49G/Hx8Ypu2t27d8PS0hK1atVSHJN7zWILCwvFMV5eXopCs379+vj999+Vzr17926Ehoaq5EpERERkSCQLwEePHuHmzZvo27evYptMJlN0tebk5Kg9mAcPHqBZs2YoW7YsFixYgEePHin2eXh4AABatWqFKlWqICIiAp9//jmSk5Mxfvx4DBw4UDGTp0ePHpgxYwYiIyPx8ccf4/r165gzZw6mTZum6N4dMmQIli1bhrFjx2LgwIE4cuQIvv/+e6xfv17teRERERGVJJIFYL9+/VCjRg2sX79ea5NAdu/ejRs3buDGjRsoU6aM0r7c8Xmmpqb4448/MHToUDRs2BDW1tbo0aOHYp1AAHB0dMSePXswbNgwhIaGwtnZGWPHjlUav1e+fHns2LEDY8aMwVdffQUvLy8sXbqUawASERGRwZNcB9DW1hZnz55VXAGEVHEdISIyFlwHkAwJP78LWAewefPmOHv2rHajISIiIiKNk+wC7tChA8aMGYPz588jODhYZWJEx44dtREfEREREamZZAE4ZMgQAMDMmTNV9mlqEggRERERaZ5kAfj6tX+JiIiIyDBIjgEkIiIiIsMk2QIIAMePH8f+/fuRmJio0iK4aNEiTcdGRERERBogWQDOmTMHU6ZMQUBAgMo6gNpYE5CIiEogvv0TGQTJAvCLL75AdHQ0IiMjtRsREREREWmU5BhAExMTNGzYULvREBEREZHGSRaAuZdIIyIiIiLDItkFPH78eLRr1w4VK1ZElSpVVBaC3rJlizbiIyIiIiI1kywAR4wYgX379iEsLAylS5fmxA8iIiIiAyFZAP7www/YvHkz2rVrp92IiIiIiEijJMcAlipVChUrVtRuNERERESkcZIF4PTp0xEVFYW0tDTtRkREREREGiXZBbx06VLcvHkT7u7u8PX1VZkEcurUKW3ER0RERERqJlkAdu7cWbuREBEREZFWSBaAUVFRhTrB+vXr0bFjR9ja2qozLiIiIiLSEMkxgIU1ePBgPHz4UD3REBEREZHGFbsAFEKoJxIiIiIi0opiF4BEREREpF9YABIREREZGRaAREREREaGBSARERGRkSl2AViuXDmVRaKJiMgwyXQdABGphWQBGBkZiX/++eeNJ7hw4QJ8fHzUHRcRERERaYhkAfjs2TO0atUKlSpVwpw5c3D//n3tRkZEREREGiFZAG7evBn379/H8OHDsXHjRvj6+qJNmzbYtGkTsrKytBslEREREalNgWMAS5cujVGjRuH06dM4fvw4/Pz8EBERAS8vL4wZMwbXr1/XXqREREREpBaFmgQSHx+P3bt3Y/fu3TA1NUXbtm1x8eJFVKlSBYsXL9Z8lERERESkNpIFYFZWFjZv3oz27dujXLly2LhxI8aMGYP4+HisWbMGu3fvxo8//oiZM2dqN2IiIiIiKhYzqR2enp6Qy+X48MMPcfz4cVSvXl3lmPDwcDg5OWk6RiIiIiJSI8kCcPHixXj//fdhZWUl+WBnZ2fExsZqKjYiIiIi0gDJAjAiIkK7kRARUYkndB0AEamFUgHYpUuXQj9wy5YtmoiHiIiIiDRMqQB0dHTUXSREREREpBVKBeCqVat0FwkRERERaYXkGMBciYmJuHr1KmQyGfz9/eHm5qadyIiIiIhIIyTXAUxNTUVERAS8vb3RtGlTNGnSBN7e3ujVqxdSUlK0GyURERERqY1kAThgwAAcO3YM27dvx9OnT5GSkoLt27fj5MmTGDhwoHajJCIiIiK1kewC/uOPP/Dnn3+iUaNGim3h4eH49ttv0bp1a23FR0RERERqJtkCWLp06XxnBTs6OsLZ2VnTcRERERGRhkgWgFOmTMHYsWMRHx+v2JaQkIAJEyZg6tSp2oqPiIiIiNRMsgt4xYoVuHHjBsqVK4eyZcsCAO7evQtLS0s8evQIX3/9teLYU6dOaSdaIiLSKZmuAyAitZAsADt37qzdSIiIiIhIKyQLwKioKO1GQkRERERaITkGkIiIiIgMk1ILYKlSpXDt2jW4uLjA2dkZMpn0aI/k5GRtxEdEREREaqZUAC5evBj29vaK2wUVgERERESkn5QKwD59+ihuR0ZG6iIeIiIiItIwyTGApqamSExMVNmelJQEU1NTTcdFRERERBoiWQAKIfLdnpGRAQsLC03GREREREQapLIMzNKlSwEAMpkM3333Hezs7BT7cnJy8M8//6By5crajZKIiIiI1EalAFy8eDHwqgVw5cqVSt29FhYW8PX1xcqVK7UbJRERERGpjUoBGBsbCwAICwvDli1b4OzsrIu4iIiIiEhDJK8Esm/fPu1GQkRERERaoVQAjh07ttAPXLRokSbiISIiIiINUyoAT58+XagHcYFoIiIiIv2lVACy25eIiIjI8EmuA0hEREREhklyEkhYWFiBXb1///23pmIiIiIiIg2SLACrV6+udD8rKwtnzpzBhQsXlK4ZTERExoNjwIkMg2QBmLsg9OumT5+O58+fazImIiIiItKgIo8B7NWrF6KjozUTDRERERFpXJELwCNHjsDKykoz0RARERGRxkl2AXfp0kXpvhAC8fHxOHnyJKZOnaqN2IiIiIhIAyQLQEdHR6X7JiYmCAgIwMyZM9GqVSttxEZEREREGiBZAK5atUq7kRARERGRVkgWgHk9f/4ccrlcaZuDg4OmYiIiIiIiDZKcBBIbG4t27drB1tYWjo6OcHZ2hrOzM5ycnODs7KzxwDIyMlC9enXIZDKcOXNGad/du3fRoUMH2NrawsXFBSNHjkRmZqbSMefPn0fTpk1hbW0Nb29vzJw5E0IIpWMOHDiAWrVqwcrKChUqVMDKlSs1nhcRERGRrkm2APbs2RMAEB0dDXd3d60v/jlx4kR4eXnh7NmzSttzcnLQrl07uLq64tChQ0hKSkKfPn0ghMCXX34JAEhNTUXLli0RFhaGEydO4Nq1a4iMjIStrS3GjRsHvCpw27Zti4EDB+Knn37C4cOHMXToULi6uqJr165azZWIiIhImyQLwHPnziEmJgYBAQHajQjAzp07sXv3bmzevBk7d+5U2rd7925cunQJcXFx8PLyAgAsXLgQkZGRmD17NhwcHLB27Vqkp6dj9erVsLS0RFBQEK5du4ZFixZh7NixkMlkWLlyJcqWLYslS5YAAAIDA3Hy5EksWLCABSAREREZNMku4Nq1ayMuLk670QB4+PAhBg4ciB9//BE2NjYq+48cOYKgoCBF8QcA4eHhyMjIQExMjOKYpk2bwtLSUumYBw8e4Pbt24pjXp/NHB4ejpMnTyIrKyvf2DIyMpCamqr0Q0RERKRvJFsAv/vuOwwZMgT3799HUFAQzM3NlfaHhISoPRghBCIjIzFkyBCEhoYqirW8EhIS4O7urrTN2dkZFhYWSEhIUBzj6+urdEzuYxISElC+fPl8z+Pu7o7s7Gw8fvwYnp6eKs89d+5czJgxQy25EhEREemKZAH46NEj3Lx5E3379lVsk8lkEEJAJpMhJyen0E8yffr0NxZOJ06cwL///ovU1FRMnjy5wGPzG4+YG5fUMbkTQIp6TF6TJ0/G2LFjFfdTU1Ph4+NTYKxEREREJY1kAdivXz/UqFED69evL/YkkOHDh+ODDz4o8BhfX1/MmjULR48eVeq6BYDQ0FD07NkTa9asgYeHB44dO6a0/8mTJ8jKylK06Hl4eChaA3MlJiYCeVoCpY4xMzND6dKl843R0tJSJTYiIiIifSNZAN65cwfbtm2Dn59fsZ/ExcUFLi4ubzxu6dKlmDVrluL+gwcPEB4ejp9//hl169YFANSvXx+zZ89GfHy8opt29+7dsLS0RK1atRTHfPzxx8jMzISFhYXiGC8vL0XXcP369fH7778rPf/u3bsRGhqq0t1NREREZEgkJ4E0b95cZQkWTStbtiyCgoIUP/7+/gCAihUrokyZMgCAVq1aoUqVKoiIiMDp06exd+9ejB8/HgMHDlQsTt2jRw9YWloiMjISFy5cwNatWzFnzhzFDGAAGDJkCO7cuYOxY8fi8uXLiI6Oxvfff4/x48drNWciIiIibZNsAezQoQPGjBmD8+fPIzg4WKVVrGPHjtqIT4WpqSn++OMPDB06FA0bNoS1tTV69OiBBQsWKI5xdHTEnj17MGzYMISGhsLZ2Rljx45VGr9Xvnx57NixA2PGjMFXX30FLy8vLF26lEvAEBERkcGTidcvj/GKiYlk42CRJ4EYqtTUVDg6OiIlJYWXxiMig+Y76Q8AwA/96qCJv6uuwyEqFn5+F9AC+Pq1f4mIiIjIMEg38xVScHCwThaMJiIiIqK3U+wC8Pbt25JXziAiIiKikqfYBSARERmPfAeNE5HeYQFIREREZGRYABIREREZGRaAREREREaGBSARERGRkSlUAZieni657+uvv4a7u7s6YyIiIiIiDZIsAOVyOT799FN4e3vDzs4Ot27dAgBMnToV33//veK4Hj16wNbWVjvREhEREVGxSRaAs2bNwurVqzF//nxYWFgotgcHB+O7777TVnxEREREpGaSBeAPP/yAb775Bj179oSpqalie0hICK5cuaKt+IiIiIhIzSQLwPv378PPz09lu1wu55U/iIiIiPSYZAFYtWpVHDx4UGX7xo0bUaNGDU3HRUREREQaYia1IyoqChEREbh//z7kcjm2bNmCq1ev4ocffsD27du1GyURERERqY1kC2CHDh3w888/Y8eOHZDJZJg2bRouX76M33//HS1bttRulERERESkNpItgAAQHh6O8PBw7UVDRERERBrHK4EQERERGRmlFkBnZ2fIZLJCPTA5OVlTMRERUQlVuE8IIirplArAJUuWKG4nJSVh1qxZCA8PR/369QEAR44cwZ9//ompU6dqP1IiIiIiUgulArBPnz6K2127dsXMmTMxfPhwxbaRI0di2bJl+OuvvzBmzBjtRkpEREREaiE5BvDPP/9E69atVbaHh4fjr7/+0nRcRERERKQhkgVg6dKlsXXrVpXtv/76K0qXLq3puIiIiIhIQySXgZkxYwb69++P/fv3K8YAHj16FLt27cJ3332nzRiJiIiISI0kC8DIyEgEBgZi6dKl2LJlC4QQqFKlCg4fPoy6detqN0oiIiIiUpsCF4KuW7cu1q5dq71oiIiIiEjjlArA1NRUODg4KG4XJPc4IiIiItIvKgtBx8fHw83NDU5OTvkuCi2EgEwmQ05OjjbjJCIiIiI1USoA//77b5QqVUpxu7BXBSEiIiIi/aFUADZt2lRxu1mzZrqIh4iIiIg0THIdwAoVKmDq1Km4evWqdiMiIiIiIo2SLACHDx+OXbt2ITAwELVq1cKSJUsQHx+v3eiIiIiISO0kC8CxY8fixIkTuHLlCtq3b48VK1agbNmyaNWqFX744QftRklEREREaiNZAOby9/fHjBkzcPXqVRw8eBCPHj1C3759tRMdEREREaldgQtB5zp+/DjWrVuHn3/+GSkpKXjvvfc0HxkRERERaYRkAXjt2jWsXbsW69atw+3btxEWFoZ58+ahS5cusLe3126URERERKQ2kgVg5cqVERoaimHDhuGDDz6Ah4eHdiMjIqISh8vDEhkGyQLwypUr8Pf31240RERERKRxkgUgiz/1ycnJQVZWlq7DICIjYW5uDlNTU12HQUQlmFIBWKpUKVy7dg0uLi5wdnYu8FJwycnJ2ohPrwkhkJCQgKdPn+o6FCIyMk5OTvDw8OAlPYkoX0oF4OLFixUTPJYsWaKrmAxGbvHn5uYGGxsbvhETkcYJIZCWlobExEQAgKenp65DIqISSKkA7NOnT763qehycnIUxV/p0qV1HQ4RGRFra2sAQGJiItzc3NgdTEQqlArA1NTUQj/QwcFBE/EYjNwxfzY2NroOhYiMUO57T1ZWFgtAIlKhVAA6OTkVupsyJydHUzEZFHb7EpEu8L2HiAqiVADu27dPcfv27duYNGkSIiMjUb9+fQDAkSNHsGbNGsydO1f7kRIRERGRWihdC7hp06aKnx9++AGLFi3C3Llz0bFjR3Ts2BFz587FggULsGrVKt1FTBrVrFkzyGQyyGQynDlzRi3nXL16NZycnNRyLn3h6+tbrIlUt2/fVutr8Lb2798PmUxW4Ez211/f6dOno3r16or7kZGR6Ny5s8ZjLYzC5FMc2sh1+vTpiv+jnKxHRG/LRGrHkSNHEBoaqrI9NDQUx48f13RcpEMDBw5EfHw8goKCgDzFSK7cD1FnZ2ekp6crPfb48eOKD6dc3bt3x7Vr17SYgf7z8fFReg2KS5NF+Jte3y+++AKrV69W3G/WrBlGjx6tkVjepEGDBoiPj4ejoyOggd/L67mqw+rVq9GsWTPF/fHjxyM+Ph5lypRR6/MQkXGRLAB9fHywcuVKle1ff/01fHx8NB0X6ZCNjQ08PDxgZia5TjgAwN7eHlu3blXaFh0djbJlyypts7a2hpubm+R5pBbJftvFsw1h0W1TU9NCvQYlwZteX0dHxxLRApyVlQULCwuNrI2Xk5MDuVyulVzt7Ozg4eHBiR1EVCySBeDixYuxfPlyBAUFYcCAARgwYACCgoKwfPlyLF68WLtRUonUp08fREdHK+6/fPkSGzZsUFlCSKqLMDo6GhUqVIClpSWEEJDJZFi5ciU6deoEW1tbzJo1CwCwYsUKVKxYERYWFggICMCPP/6odH6px73O19cXs2bNQu/evWFnZ4dy5crht99+w6NHj9CpUyfY2dkhODgYJ0+eVHrc5s2bUbVqVVhaWsLX1xcLFy5U2p+YmIgOHTrA2toa5cuXx9q1a1WeOyUlBYMGDYKbmxscHBzQvHlznD17VvJ3+3oXcG6r6969exEaGgobGxs0aNAAV69eVTzm7NmzCAsLg729PRwcHFCrVi2cPHkS+/fvR9++fZGSkqJonZ0+fToA4KeffkJoaCjs7e3h4eGBHj16KNaPy+vw4cOoVq0arKysULduXZw/f17y9X1d3m7RyMhIHDhwAF988YUiltjYWPj5+WHBggVKj7tw4QJMTExw8+ZNyXNHR0crXhtPT08MHz5csS+/v4u8XcAF/V4yMzMxceJEeHt7w9bWFnXr1sX+/ftVct6+fTuqVKkCS0tL3LlzR6ULOCMjAyNHjoSbmxusrKzQqFEjnDhxQrG/MK8rEZEmSBaAbdu2xfXr19GpUyckJycjKSkJnTp1wrVr19C2bVvtRmkghBBIy8zW+o8QQiP5RERE4ODBg7h79y7wqlDy9fVFzZo13/jYGzdu4JdffsHmzZuVxrlFRUWhU6dOOH/+PPr164etW7di1KhRGDduHC5cuIDBgwejb9++ShOW8nuclMWLF6Nhw4Y4ffo02rVrh4iICPTu3Ru9evXCqVOn4Ofnh969eyt+ZzExMejWrRs++OADnD9/HtOnT8fUqVOVuvkiIyNx+/Zt/P3339i0aROWL1+uVEQJIdCuXTskJCRgx44diImJQc2aNdGiRYsiX1Hnk08+wcKFC3Hy5EmYmZkp5dqzZ0+UKVMGJ06cQExMDCZNmgRzc3M0aNAAS5YsgYODA+Lj4xEfH4/x48cDrwqdTz/9FGfPnsWvv/6K2NhYREZGqjzvhAkTsGDBApw4cQJubm7o2LHjW7W0fvHFF6hfv75imEF8fDzKli2Lfv36qYwtjo6ORuPGjVGxYsV8z7VixQoMGzYMgwYNwvnz57Ft2zb4+fkpHVPQ30VBv5e+ffvi8OHD2LBhA86dO4f3338frVu3xvXr1xWPT0tLw9y5c/Hdd9/h4sWL+baCTpw4EZs3b8aaNWsUf1/h4eEqr3tBrysRkSYU2L9UpkwZzJ49u8ATDB06FDNnzoSLi4u6YzM4L7NyUGXan1p/3kszw2Fj8fZdib6+vvkWkW5ubmjTpg1Wr16NadOmITo6utAfXJmZmfjxxx/h6uqqtL1Hjx5K5+jRowciIyMxdOhQAMDYsWNx9OhRLFiwAGFhYZKPk9K2bVsMHjwYADBt2jSsWLECtWvXxvvvvw8A+Oijj1C/fn08fPgQHh4eWLRoEVq0aIGpU6cCr66RfenSJXz++eeIjIzEtWvXsHPnThw9ehR169YFAHz//fcIDAxUPOe+fftw/vx5JCYmwtLSEgCwYMEC/Prrr9i0aRMGDRpUqN8ZAMyePRtNmzYFAEyaNAnt2rVDeno6rKyscPfuXUyYMAGVK1cGAFSqVEnxOEdHR8hkMnh4eCidL+/vrEKFCli6dCnq1KmD58+fw87OTrEvKioKLVu2BACsWbMGZcqUwdatW9GtW7dCx54bh4WFhWKYQa6+ffti2rRpOH78OOrUqYOsrCz89NNP+PzzzyXPNWvWLIwbNw6jRo1SbKtdu7bSMa//XcTGxipuW1hY5Pt7uXnzJtavX4979+7By8sLeDXubteuXVi1ahXmzJkDvOpSXr58OapVq5ZvfC9evMCKFSuwevVqtGnTBgDw7bffYs+ePfj+++8xYcIExbEFva6RkZH5FuVERMUh2QJYWD/99FORFpAmw9KvXz+sXr0at27dwpEjR9CzZ89CPa5cuXIqxR9eTTLK6/Lly2jYsKHStoYNG+Ly5csFPk5KSEiI4ra7uzsAIDg4WGVbbgue1PNfv34dOTk5uHz5MszMzJSev3LlykpdojExMXj+/DlKly4NOzs7xU9sbGyB3Ztvij/3El+5sY4dOxYDBgzAO++8g3nz5hXq3KdPn0anTp1Qrlw52NvbKyYb5Lbq5spdCgqvrhkeEBCg8hoUh6enJ9q1a6cYUrB9+3akp6crCvPXJSYm4sGDB2jRokWB5y3s30Vep06dghAC/v7+Sq/XgQMHlH6nFhYWSq/H627evImsrCylvx9zc3PUqVNH5XdX0OtKRKQJxR5hrqnuRUNkbW6KSzPDdfK8mpLbota/f3906NCh0Je9s7W1LfT21wfs544XLMz5Xmdubq5y3vy2yeVyyefK+zefe7ugSQVyuRyenp5KY8hyFXXCQEGxTp8+HT169MAff/yBnTt3IioqChs2bMC7776b77levHiBVq1aoVWrVvjpp5/g6uqKu3fvIjw8HJmZmW+MRd0TKQYMGICIiAgsXrwYq1atQvfu3SWvpJN7qbM3KezfRV5yuRympqaIiYlRmWiRt1XU2tq6wN+B1N9Gfn9TBb2uRESaUPKnGBoQmUxWrK7YksjU1BQRERGYP38+du7cqfbzBwYG4tChQ+jdu7di27///qvUxapJVapUwaFDh5S2/fvvv/D394epqSkCAwORnZ2NkydPok6dOgCAq1evKq0zV7NmTSQkJMDMzAy+vr4ajdff3x/+/v4YM2YMPvzwQ6xatQrvvvsuLCwsVK7ec+XKFTx+/Bjz5s1TzOx/fQJMrqNHjypmdz958gTXrl1TdDUXVX6x4NWXCVtbW6xYsQI7d+7EP//8I3kOe3t7+Pr6Yu/evUpDAdQRS40aNZCTk4PExEQ0btz4rc/t5+cHCwsLHDp0CD169ABedRufPHlSZ8vgEBHlMqxqhHTi008/xYQJEwrd+lcUEyZMQLdu3RSTJn7//Xds2bIFf/31l9qfKz/jxo1D7dq18emnn6J79+44cuQIli1bhuXLlwMAAgIC0Lp1awwcOBDffPMNzMzMMHr0aKUWqnfeeQf169dH586d8dlnnyEgIAAPHjzAjh070Llz57fqpnzdy5cvMWHCBLz33nsoX7487t27hxMnTqBr167Aq3Gcz58/x969e1GtWjXY2NigbNmysLCwwJdffokhQ4bgwoUL+PTTT/M9/8yZM1G6dGm4u7vjk08+gYuLy1sveOzr64tjx47h9u3bsLOzQ6lSpWBiYgJTU1NERkZi8uTJ8PPzU+p2zs/06dMxZMgQxVjUZ8+e4fDhwxgxYkSRYnn99+Lv74+ePXuid+/eWLhwIWrUqIHHjx/j77//RnBwcKEnwdna2uL//u//MGHCBJQqVQply5bF/PnzkZaWhv79+xc6xpKmXKmit6oSUclT7DGARBYWFnBxcdHItUc7d+6ML774Ap9//jmqVq2Kr7/+GqtWrVJaGFeTatasiV9++QUbNmxAUFAQpk2bhpkzZyoNyl+1ahV8fHzQtGlTdOnSRbHcSy6ZTIYdO3agSZMm6NevH/z9/fHBBx/g9u3bijGHxWVqaoqkpCT07t0b/v7+6NatG9q0aYMZM2YAr2a8DhkyBN27d4erqyvmz58PV1dXrF69Ghs3bkSVKlUwb948laVYcs2bNw+jRo1CrVq1EB8fj23btsHCwuKtYh0/fjxMTU1RpUoVRbdzrv79+yMzM7NQE3r69OmDJUuWYPny5ahatSrat2+vNEu3MPL7veDVa9q7d2+MGzcOAQEB6NixI44dO1bkNVDnzZuHrl27IiIiAjVr1sSNGzfw559/wtnZuUjnKQl2jmqMdQPqomzp/LvliUi/yEQxB/HZ29vj7NmzqFChgvqi0hOpqalwdHRESkoKHBwclPalp6cjNjYW5cuXh5WVlc5iLKpmzZqhevXqvMQU6cThw4fRrFkz3Lt3T23FsaHy9fXF6NGjJbuT9fU9iEgbCvr8NhbFbgHs1auX0f7yDNXy5cthZ2entNgvkSZlZGTgxo0bmDp1Krp168birwBz5syBnZ2dykxtIqKikGwB9PX1Rb9+/RAZGalyaS/6jyG2AN6/fx8vX74EAMUYMSJNW716Nfr374/q1atj27Zt8Pb21nVIJVZycrJiIWlXV1fFdY1fp6/vQUTawBbAAloAx40bh99++w0VKlRAy5YtsWHDBmRkZGg3OtI6b29v+Pn5KWYwEmlDZGQkcnJyEBMTw+LvDUqVKqX4PypV/BERvYlkAThixAjExMQgJiYGVapUwciRIxXX2jx16pR2oyQiIiIitXnjGMBq1arhiy++wP379xEVFYXvvvsOtWvXRrVq1RAdHc2FoN+Avx8i0gW+9xBRQd64DmBWVha2bt2KVatWYc+ePahXrx769++PBw8e4JNPPsFff/2FdevWaSdaPZK7sn9aWlqhr1pARKQuaWlpwGtXGSEiyiVZAJ46dQqrVq3C+vXrFVd7WLx4sdLq/61atUKTJk20FateMTU1hZOTk+J6njY2NhpZJ4+IKC8hBNLS0pCYmAgnJyeVy9kREaGgArB27dpo2bIlVqxYgc6dO+f7LbJKlSr44IMPNB2j3vLw8AB4UXci0gEnJyfFexAR0eskl4G5c+cOypUrp/2I9Ehhp5Hn5OQgKytLq7ERkfEyNzdnyx9RAbgMTAEtgCz+1MfU1JRvxkRERFRiKBWAzs7OhR6nlrsQKRERERHpF6UCkNd/JSIiIjJ8SgVgnz59dBcJEREREWlFgesA5uTkYOvWrbh8+TJkMhkCAwPRqVMnmJm9cflAo5A7fyY1NVXXoRAREVEh5X5uG/OC6ZKV3IULF9CpUyckJCQgICAAAHDt2jW4urpi27ZtCA4O1macJdKzZ88AAD4+ProOhYiIiIro2bNnRntNbcllYOrVqwc3NzesWbMGzs7OAIAnT54gMjISiYmJOHLkiLZjLXHkcjkePHgAe3t7tS/ynJqaCh8fH8TFxRnkFHXmp/8MPUfmp/8MPUfm9/aEEHj27Bm8vLxgYvLGq+IaJMkWwLNnz+LkyZOK4g+vZgnPnj0btWvX1lZ8JZqJiQnKlCmj0edwcHAwyP/YuZif/jP0HJmf/jP0HJnf2zHWlr9ckmVvQEAAHj58qLI9MTERfn5+mo6LiIiIiDREsgCcM2cORo4ciU2bNuHevXu4d+8eNm3ahNGjR+Ozzz5Damqq4oeIiIiI9IdkF3D79u0BAN26dVOMb8sdLtihQwfFfZlMhpycHO1Ea0QsLS0RFRUFS0tLXYeiEcxP/xl6jsxP/xl6jsyPikNyEsiBAwcKfZKmTZuqMyYiIiIi0iDJApCIiIiIDJNSF/C5c+cQFBQEExMTnDt3rsAHhoSEaDo2IiIiItIApRZAExMTJCQkwM3NDSYmJpDJZPmuks1xf0RERET6S6kFMDY2Fq6urorbRERERGR4lArAcuXKKW7fuXMHDRo0ULnub3Z2Nv7991+lY6nkyp2pbcgMPUfmR6Rb/BvVf3wNVUmuAxgWFobk5GSV7SkpKQgLC9N0XFQMebvnDfUPPjs7W3FbJpNBLpfrNB51u379OpYuXQoY6GuY9/UyxPzyY2h/ozdu3MDWrVt1HYbGPH/+HImJicjKyjLov1FDHs5lDJ+FxSG5DqBUtZyUlARbW1tNx0Vv6cqVK5g/fz7S0tJgY2ODKVOmoEyZMrCwsNB1aGpz9epVfPbZZ0hKSoKlpSV++eUXmJiYGMw3vHPnzqFNmzZo27Ytrl27Bn9/f8CAvsFeuXIFS5cuRXp6OszMzPDJJ5/Ay8sL5ubmug5NbW7duoX169fj8ePH8PT0xMSJEw3qb/TMmTOoXbs2vvrqK12HohEXL17EsGHDkJycjIyMDHz33Xdo3Lixwbx+AHD37l04OTnBwcEBOTk5MDU11XVIamUMn4XFpdIC2KVLF3Tp0gUymQyRkZGK+126dEGnTp0QHh6OBg0a6CZaKtCVK1dQp04dvHz5EnZ2djh79iwaNWqEb7/9Fo8fP9Z1eGpx/vx5NGzYEDk5OahYsSJOnTqF999/H3j1DU/fVzW6f/8+OnbsiA8++ADffvutoviDgXyDvXTpEho0aIAXL17AyckJZ86cQWhoKL799lskJSXpOjy1OH/+POrXr49Tp07hxIkTWL16NXr16gUYyGuY+74yatQoDBo0SGW/vv8fvHLlCpo2bYrg4GDMnj0bFSpUwNChQwEDeY/Bqy/RFStWRI0aNZCcnAxTU1ODagk0hs9CdVBpAcy9OLIQAvb29rC2tlbss7CwQL169TBw4EDtRklvJJfLsXDhQrRp0wbr169XbB86dCi++OILvHz5EgMGDICTk5NO4yyOmzdvomvXrhg0aBDmzJkDIQQ8PT1x584dxTH6/gF75coVBAQEYOHChcjOzsaYMWNw//59PH78GCNGjEDTpk3h5uam6zDfSnp6OiZOnIgePXpg2bJlwKv3mYCAAMybNw8vX77EwIED9fqi9nFxcejWrRv69u2LefPmIS0tDVu3bsWiRYtw/fp1VKpUSdchFsuVK1fQpEkT9O/fHwsWLEBOTg5++eUX3Lt3DwAwYsQIWFlZ6TrMt5aZmYmoqCh07twZX375JQDAy8sLn332GR4/fgxra2uYm5vrdSvSw4cPMXToUDRt2hRpaWlo1qwZ9u/fj1KlShlES6AxfBaqi0oBuGrVKgCAr68vxo8fz+5ePWFiYoJnz57BxsYGePVGZmFhgeXLl2PMmDH48ssv4e/vj44dO+ptN8bOnTsRFBSEjz76CHhV7N29exd79+5Fo0aNYGpqilmzZul1V01cXBwePnyItLQ0tGvXDiYmJmjcuDFOnz6NqVOn4sMPP8TYsWNhb2+v61CLTC6X4/Hjx4rWsOfPn8POzg61atXC1atXsXjxYtSpU0dvXz8hBP7880+UKVMGY8aMgRACNjY2aNasGUaOHIm7d+/qfQG4atUqPHv2DC1atMDTp0/x3nvvIS0tDcnJyXjx4gVWrlyJzZs3o3r16nr5GpqYmODRo0cIDg5WbNu4cSP27NmDRo0aQS6Xo0ePHhg+fDhcXFx0GuvbOnfuHEqVKoWhQ4fCwsICEydOVCoCs7OzVSZ/6hNj+CxUG0EGY+DAgSI4OFhxPz09XXH7/fffFwEBASIrK0tH0anHiRMnFLcXLVokZDKZmDNnjvj5559Fu3btRJkyZcSjR490GmNxbN++XQQGBoq9e/eKzp07i8TERMW+qKgoUa5cOXH9+nWdxvi2cnJyRI0aNUSvXr0U2xISEkTFihXFqVOnRMuWLUXDhg11GmNxHTx4UKxYsUJxPysrS2RlZYlKlSqJ7du3qxyfk5Oj5QiLr2/fvqJSpUoiMDBQtGvXTty6dUs8ffpUPH36VISFhYng4GCRnZ2t6zDfWvv27YW/v79Yvny5GDVqlLCyshKrVq0SFy9eFDNmzBB+fn5iz549ug6zWPbt26e4ffDgQdGgQQMRFBSkeO/Mff3kcrnOYiwOY/gsVAelArB69eqiRo0ahfqhkic2NlZ4e3uLvn37KralpaUJIYS4deuWcHV1FXv37tVhhOrz6NEjMXr0aKV80tPThZ2dnYiOjtZpbMUVFBQknJ2dRXBwsHj8+LHSvjJlyoiFCxfqLLbi+uGHH0T58uVFw4YNxSeffCIcHR1F//79hRBC7NixQ1StWlU8fPhQ12G+tYyMDMXtvB+eNWvWFBs3blTc//HHH7UeW3Hl/cDs27evqF27trh69arSMcePHxf29vZi//79OoiweHKL8RcvXoj27duLfv36iYCAALFgwQKl4ypVqiRGjhypoyjVTy6Xi0OHDimKwNz3nDlz5ojDhw/rOry3YkyfhcWh1M7buXNn3TVFUpHExsbiwIEDePToEcLCwlClShX4+voiKioK8+fPx/Dhw7Fs2TLFGM7s7Gw4ODjoVddhfjnmNuu7uLhg9uzZivt4NfPS398fgYGBOoy68F7Pr1KlSnB0dMRXX32F4cOH4+HDh7hz5w5Kly4NAHj58iUqVqwIX19fXYdeKHnza9asGWrWrInu3bvD1dUVS5YswY0bNzBt2jSMHTsWAPD48WNkZ2fr1fiqFy9ewNbWVtGVlDf23AkDMpkML1++VGyPiorCp59+inr16sHPz09HkRdO3vzMzMwUY8Sio6Oxd+9exXqwuXmmpqbC09MTZcqU0XXohZI3PxMTE8jlctjY2OD3338HALRs2RKVK1cGAGRkZAAAKlSoUOJft7zi4+Nx7tw5PH36FA0bNkSZMmUUr2Pu61a/fn3Mnz8fEydORIsWLRAaGoro6GhcunRJ1+G/kTF8FmqMritQKrpz584JFxcXERoaKgIDA4W5ubkYNmyYOHPmjMjOzhaLFy8W5cqVE+3btxfXr18Xly5dEtOmTRPlypUT9+/f13X4hZJfjsOHDxdHjhxRHPN6N9OUKVNErVq1REJCgg4iLhqp/HK7uLdu3SrKlSsnAgMDxZYtW8Q///wjpk2bJry8vMStW7d0Hf4bSf2Nnjx5UnFM7jfyXMOHDxcdOnRQ2V5SXbhwQbi6uoq1a9dKHpOdnS3S09NFhQoVxO7du8XChQuFjY2NiImJ0Wqsb0Mqv4K6dydOnCiaNm0qkpOTtRBh8byeX26Lbd6W2xYtWoh27doJ8Wq4wqxZs4SHh4feDMM4d+6cqFy5sggJCRFlypQRLi4u4sqVK0Lkk69cLhf79+8XNjY2olSpUuLUqVM6jb0wjOGzUJNYAOqZlJQU0ahRIzF+/Hjx8uVLIV51J9WqVUt07NhRnDhxQsjlcrFjxw4RFBQkSpcuLfz8/ET58uX14kNHFJBjaGio6Ny5szhw4IDS8WfOnBETJ04Ujo6O4syZMzqKuvAKyq9jx46KIvfy5csiLCxM+Pr6igoVKoiQkBC9eFMuKL9OnTopjT8SQoh///1XjB8/Xtjb24uzZ8/qKOqiuXv3rggJCRFeXl7C0tKywCJQCCEaNWokAgIChLW1tdI41pKqqPmdOnVKTJgwQTg6OurFa/im/HK7g/fs2SPKlCkjHBwcRM2aNUX58uX14v+gEEJcv35deHl5icmTJ4sHDx6Iq1eviq5du4p27dop/l++7v/+7/+ElZWVuHDhgtbjLSpj+CzUNMkCUCaTCRMTE8kf0o3U1FRRuXJlsXr1aqXtO3bsEA0bNhTdu3dXtBDJ5XKxb98+ERMTIx48eKCjiIvuTTl+8MEHirFHN2/eFB999JGoVq2aXnzwiEK+hteuXVNsv3Llirh165bKeMCS6k35ffjhh0pjx1auXCmaNm2qF8W7eDUWbunSpaJLly7iyJEj4uOPPxYmJiaSRVJ2draoWrWqMDMzE+fOndN6vEVV1PyuX78uRo4cKQICAvTiNSxMfrmtYtnZ2eLu3bti4cKFYtOmTeL27ds6jLzw0tPTxdChQ0VERIRSi+2KFStE1apV823F/ffff0Xt2rX1pjgyhs9CTZOc6/36JX6ysrJw+vRprFmzBjNmzNBG7zS9RgiBFy9ewMTEBKmpqcCrcSmWlpZo06YNsrKyMGzYMGzbtg2jRo2CTCZDs2bNdB12kRQ2x127dsHf3x+enp4YNGgQxowZA3d3d12H/0aFzW/nzp2KJUMCAgJ0HHXhFfX1A4DBgweje/fuerMul5mZGerVqwdPT0/Uq1dPseRJREQEAKBHjx6KY+VyOUxNTTFjxgxUr14dFStW1GHkhVOU/ACgfPnyGDZsGCZNmgRPT08dRV14hckvd/ymqakpfHx8FONU9YWFhQV8fX0VOeRq3rw55s6di8TERLi6uiot91K/fn3s2rULpUqV0lHUhWcMn4VaUdSKce3ataJjx46aKUepUCZMmCAcHBwU30bzzjyMiooSZcuWFc+fP9dhhMVXmByfPXumwwiLx9Bfw8Lmp8/LheT17NkzMXnyZKWWpPT0dLFjxw69XpYoV0H5GUKLilR+O3fuFHFxcboO7608ffpUZdvly5dFmTJlRGJioqKVMzY2VgfRqYehv49qWpFXe6xbty6vBKIjuTO2Pv74Y8TExKBZs2Y4cuQIPDw8IJfLYWJiAn9/f7i4uOjtau5FyVEfFys19NewqPnpY4555eZrZ2eHyZMnQyaTISIiAnK5HCdOnMC6detw/vx5XYf51t6U3/r163Hu3Dldh/nWDDm/vFf1kslkkMvlSE9PhxACVlZWkMlkmDBhAtavX48rV67Azs5O1yEXmqG/j2pNUarFtLQ0MWrUKOHv76+5kpQK5eTJk6JRo0bC29tbnDhxQtEaNmrUKNGgQQO9bh3LZeg5Mj/9zi8/uS1JMplMODo6iuPHj+s6JLVifvrt6tWrwsvLSzx79kx88sknws7OThw9elTXYRWLMb7PqItkAejk5CScnZ0VP05OTsLU1FTY29uL3377TbtRUr4uXrwounbtKqysrES1atVEkyZNhKOjozh9+rSuQ1MbQ8+R+RmWrKwsERkZKZycnMSlS5d0HY7aMT/9dvPmTVGlShXRu3dvYWFhobQskz6QujKJsb3PqItMCCHyaxlcs2aN0n0TExO4urqibt26cHZ21lYDpVEr7IW5N27ciDt37sDU1BQdOnTQq0VKDT1H5vcffc0PRchRCIENGzZg2LBh2LNnD2rVqqWV+IqL+f1HX/NDIXMUQiAmJgZ16tSBnZ0d/vnnH1SvXl1rMb6NtLQ0ZGVlQS6XK+qO3C7e/Ojz+4xO6LoCJWVXrlwR48aNU1qGQIq+XqfR0HNkfv+jj/mJIuaY1/nz58WdO3c0HF3xMb/86Ut+4i1zTE9PF8OGDRPnz5/XQoTFc+HCBdGhQwdRuXJl0bRpUzFv3jzJY/X1fUbXChxFn56ejnPnziExMRFyuVxpX8eOHTVdmxqdmzdvIiwsDElJSUhOTsb3338PU1NTyW93MplM6X7uwNiSzNBzZH7K9C0/vEWOeQUFBWktzrfF/KTpQ34oRo6WlpZYtGhRib/c4qVLl9CkSRNERESgQ4cOuHTpEjZt2oTQ0FC0aNFC5Xh9fJ8pCSS7gHft2oWIiAgkJSWpPkgmQ05OjjbiMxqpqakYPHgwsrKyEBoail9++QVBQUFYs2aN4vet77OZDD1H5qff+cEIcmR++p0fjCDH5ORkvP/++wgJCcHixYsBAImJiWjZsiW6deuGTz75RNchGoz8O9IBDB8+HN26dUN8fDzkcrnSD4s/9bO3t0eFChXw3nvvYdSoURg8eDAuXLiAPn36KBbz1Pffu6HnyPz0Oz8YQY7MT7/zgxHk+PDhQ7i4uKBt27bAqzF/bm5uaNWqFeLi4oBXYx5zSbRhUWFI9Q3b29uLGzduaLU/2ljljt3IO47h2bNnYsWKFaJGjRoiIiJCsS89PV2kpaXpLNa3Zeg5Mj/9zk8YQY7MT7/zE0aS48OHD8X27dsV93PzGTNmjOjWrZvSNioeyRbA9957D/v379duNWpkMjMzgVczrJFnHENWVhbs7OzQu3dvDBw4UPHtLj09HSNHjkS/fv1UxmSWVIaeI/PT7/xgBDkyP/3OD0aQY25+ua197dq1A1619OXmamZmhuzsbOBV/lFRUZgyZYoOo9Z/kpNAli1bhvfffx8HDx5EcHAwzM3NlfaPHDlSG/EZrEuXLiEqKgrPnz+HTCbDlClTUK1aNdja2sLc3BxyuRw2Njbo06cPZDIZvv/+e1SqVAlJSUnYt2+f5DT4ksTQc2R++p0fjCBH5qff+cEIcsybHwBMnToVISEhsLOzg6mpqWJCh7Ozs2Lyyscff4zFixfj0KFDOo5ev0n+Zaxbtw5//vknNm/ejC+//BKLFy9W/CxZskS7URqYa9euoX79+nBwcEBISAjMzMwQHh6OBQsW4Pbt28Crb3q5/7G7d+8OU1NTpKWl4fjx46hbt66uU3gjQ8+R+el3fjCCHJmffucHI8jx9fzMzc0RHh6OhQsXKvLLbQF8+fIlLCwsMGfOHCxatAiHDx/Wq7UaSySpvmF3d3cxe/ZskZOTo91OaSMwbtw40b59e6Vt8+bNE76+vmLSpEni/v37iu2ZmZni448/FlZWVuLs2bM6iPbtGHqOzE+/8xNGkCPz0+/8hBHkWJT8Pv74YyGTyYSDg4PeXcGkpJLsAs7MzET37t1LfPOxPnrx4oXi95qVlQVzc3N89NFHMDc3x/z581GhQgUMHDgQcrkcZmZmyMnJwdGjRxESEqLr0AvN0HNkfvqdH4wgR+an3/nBCHIsbH4AULlyZfj7+2PTpk16s15jiSdVGY4ePVrMnj1bu+WokZg+fbpwd3cXL168EEIIkZGRodg3YcIEUbp0aZGYmKjDCIvP0HNkfvqdnzCCHJmffucnjCDHouR37949ER8fr7NYDZFkAThixAjh6OgomjRpIoYPHy7GjBmj9ENFlzt1PTk5WVSvXl20atVK8QefO13/8ePHwtvbW2zZskWnsb4tQ8+R+el3fsIIcmR++p2fMIIci5Lf5s2bdRqrIZPs3z1//jxq1KgBExMTXLhwAadPn1b8nDlzRrvNlAYidzCrk5MTJk+ejMTERHz44YfIycmBtbU1ACA7Oxt2dnawtbXVcbRvx9BzZH76nR+MIEfmp9/5wQhyLEp+dnZ2Oo7WcEmOAdy3b592IzESueMc3n33XeTk5GD27NkIDQ3FihUrIITAzp078ezZMwQGBuo61Ldm6DkyP/3OD0aQI/PT7/xgBDkaen56QddNkMYkdxX369evi5kzZwohhDhy5Iho3bq1KFWqlPDz8xOVK1cWMTExOo707Rl6jsxPv/MTRpAj89Pv/IQR5Gjo+ekLmeCF9LQidzHLO3fuoGHDhmjUqBE2bNig2H/69Gk4ODjAwcEBrq6uOo31bRl6jsxPv/ODEeTI/PQ7PxhBjoaen17RdQVqaC5fviyWLFkiXr58qbLvyZMnokaNGmLQoEGKQbD6uM6ioefI/PQ7P2EEOTI//c5PGEGOhp6fIWABqEbXr18XpUqVEjKZTHzyySciKytLaf+TJ0/E5s2b9fpC1oaeI/PT7/yEEeTI/PQ7P2EEORp6foaCBaCaPHv2TAwYMEB8+OGHYsWKFcLMzExMnDhR8YdvCH/ohp4j89Pv/IQR5Mj89Ds/YQQ5Gnp+hkRyFjAVTXp6OgICAuDr64v33nsPpUuXRo8ePSCTyTBr1iyYman+qnPHQugLQ8+R+el3fjCCHJmffucHI8jR0PMzKLquQA3J66uUb9iwQeXbT3Z2trh9+7aOIiw+Q8+R+el3fsIIcmR++p2fMIIcDT0/Q8ECsBiysrJEZmam0ja5XK7UxJ37h//RRx+J58+fi5EjR4rIyEjFpW9KOkPPkfnpd37CCHJkfvqdnzCCHA09P0PFAvAtXbx4UXTr1k00atRIREZGinXr1in2vT7gdcOGDcLa2loEBgYKExMTcerUKR1EXHSGniPz+x99zE8YQY7M73/0MT9hBDkaen6GjOsAvoVr166hTp066NChAypVqoS9e/fi2bNnqFatGlatWgUAyMnJgampqeIxLVq0wJkzZ7B//34EBwfrMPrCMfQcmZ9+5wcjyJH56Xd+MIIcDT0/Q/f/sdtd4iEg7dEAAAAASUVORK5CYII=", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "fe6fe4b0633c42c0882da0803614f477" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 13 + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": "", + "id": "7558642fc08c06cc" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/speasy/products/variable.py b/speasy/products/variable.py index ab797dc5..65ac1f9b 100644 --- a/speasy/products/variable.py +++ b/speasy/products/variable.py @@ -217,6 +217,11 @@ def __mul__(self, other): def __rmul__(self, other): return self.__mul__(other) + def __pow__(self, power, modulo=None): + res = self.copy() + np.power(self.__values_container.values, power, out=res.__values_container.values) + return res + def __add__(self, other): if type(other) is SpeasyVariable: if self.__values_container.shape != other.__values_container.shape: @@ -263,9 +268,24 @@ def __truediv__(self, other): res = self.copy() np.divide(self.__values_container.values, float(other), out=res.__values_container.values) return res + if type(other) is SpeasyVariable: + return np.divide(self, other) raise TypeError( f"Can't divide SpeasyVariable by {type(other)}") + def __rtruediv__(self, other): + if type(other) in (int, float): + res = self.copy() + np.divide(float(other), self.__values_container.values, out=res.__values_container.values) + return res + if type(other) is SpeasyVariable: + res = self.copy() + np.divide(other.__values_container.values, self.__values_container.values, + out=res.__values_container.values) + return np.divide(other, self) + raise TypeError( + f"Can't divide {type(other)} by SpeasyVariable") + def __array_function__(self, func, types, args, kwargs): if func.__name__ in SpeasyVariable.__LIKE_NP_FUNCTIONS__: return SpeasyVariable.__dict__[func.__name__](self) diff --git a/speasy/signal/filtering/__init__.py b/speasy/signal/filtering/__init__.py index 683df044..d6032114 100644 --- a/speasy/signal/filtering/__init__.py +++ b/speasy/signal/filtering/__init__.py @@ -1,28 +1,65 @@ from scipy import signal -from typing import Callable +from typing import Callable, Union, Collection from speasy.products import SpeasyVariable import numpy as np -def apply_sos_filter(sos: np.ndarray, filter_function: Callable, var: SpeasyVariable) -> SpeasyVariable: - res = np.empty_like(var) +def _apply_filter(filter_function: Callable, sos: np.ndarray, var: SpeasyVariable) -> SpeasyVariable: + res = SpeasyVariable.reserve_like(var) res.values[:] = filter_function(sos, var.values, axis=0) return res -def sosfiltfilt(sos: np.ndarray, var: SpeasyVariable) -> SpeasyVariable: - """Apply an IIR filter to the data using :func:`scipy.signal.sosfiltfilt`. +def apply_sos_filter(sos: np.ndarray, filter_function: Callable, + var: Union[SpeasyVariable, Collection[SpeasyVariable]]) -> Union[ + SpeasyVariable, Collection[SpeasyVariable]]: + """Apply an IIR filter to the variable(s) using the given filter function. This function just applies the filter to the + values of the variable without any resampling, it assumes that the variable has a regular time axis. + + Parameters + ---------- + sos: np.ndarray + Second-order sections representation of the filter, as returned by :func:`scipy.signal.iirfilter` with `output='sos'` for example. + filter_function: Callable + The filter function to use (e.g. :func:`scipy.signal.sosfiltfilt`) + var: SpeasyVariable or Collection[SpeasyVariable] + The variable(s) to filter + + Returns + ------- + SpeasyVariable or Collection[SpeasyVariable] + The filtered variable(s) + + Notes + ----- + It only supports 1D variables. + """ + + if isinstance(var, SpeasyVariable): + return _apply_filter(filter_function, sos, var) + else: + return [_apply_filter(filter_function, sos, v) for v in var] + + +def sosfiltfilt(sos: np.ndarray, var: Union[SpeasyVariable, Collection[SpeasyVariable]]) -> Union[ + SpeasyVariable, Collection[SpeasyVariable]]: + """Apply an IIR filter to the variable(s) using :func:`scipy.signal.sosfiltfilt`. This function just applies the filter to + the values of the variable without any resampling, it assumes that the variable has a regular time axis. Parameters ---------- sos: np.ndarray Second-order sections representation of the filter - var: SpeasyVariable - The variable to filter + var: SpeasyVariable or Collection[SpeasyVariable] + The variable(s) to filter Returns ------- - SpeasyVariable - The filtered variable + SpeasyVariable or Collection[SpeasyVariable] + The filtered variable(s) + + Notes + ----- + It only supports 1D variables. """ return apply_sos_filter(sos, signal.sosfiltfilt, var) diff --git a/speasy/signal/resampling/__init__.py b/speasy/signal/resampling/__init__.py index f761af57..a99fc9af 100644 --- a/speasy/signal/resampling/__init__.py +++ b/speasy/signal/resampling/__init__.py @@ -60,13 +60,13 @@ def _interpolate(ref_time: np.ndarray, var: SpeasyVariable, interpolate_callback def resample(var: Union[SpeasyVariable, Collection[SpeasyVariable]], new_dt: Union[float, np.timedelta64], interpolate_callback: Optional[Callable] = None, *args, **kwargs) -> Union[SpeasyVariable, Collection[SpeasyVariable]]: - """Resample a variable to a new time step. The time vector will be generated from the start and stop times of the + """Resample a variable(s) to a new time step. The time vector will be generated from the start and stop times of the input variable. Uses :func:`numpy.interp` to do the resampling by default. Parameters ---------- var: SpeasyVariable or Collection[SpeasyVariable] - The variable or a collection of variables to resample + The variable(s) to resample new_dt: float or np.timedelta64 The new time step in seconds or as a numpy timedelta64 interpolate_callback: Callable or None @@ -75,7 +75,11 @@ def resample(var: Union[SpeasyVariable, Collection[SpeasyVariable]], new_dt: Uni Returns ------- SpeasyVariable or Collection[SpeasyVariable] - The resampled variable or a collection of resampled variables + The resampled variable(s) with all metadata preserved except for the new time axis + + Notes + ----- + It only supports 1D variables. """ if type(var) in (list, tuple): return [resample(v, new_dt, interpolate_callback, *args, **kwargs) for v in var] @@ -87,21 +91,26 @@ def resample(var: Union[SpeasyVariable, Collection[SpeasyVariable]], new_dt: Uni def interpolate(ref: Union[np.ndarray, SpeasyVariable], var: Union[SpeasyVariable, Collection[SpeasyVariable]], interpolate_callback: Optional[Callable] = None, *args, **kwargs) -> Union[SpeasyVariable, Collection[SpeasyVariable]]: - """Interpolate a variable to a new time vector. The time vector will be taken from the reference variable. Uses :func:`numpy.interp` to do the resampling by default. + """Interpolate a variable(s) to a new time vector. The time vector will be taken from the reference variable. + Uses :func:`numpy.interp` to do the resampling by default. Parameters ---------- ref: np.ndarray or SpeasyVariable The reference time vector var: SpeasyVariable or Collection[SpeasyVariable] - The variable or a collection of variables to interpolate + The variable(s) to interpolate interpolate_callback: Callable or None The interpolation function to use, defaults to :func:`numpy.interp` (Optional) Returns ------- SpeasyVariable or Collection[SpeasyVariable] - The interpolated variable or a collection of interpolated variables + The interpolated variable(s) with all metadata preserved except for the new time axis + + Notes + ----- + It only supports 1D variables. """ if isinstance(ref, SpeasyVariable): ref_time = ref.time