diff --git a/docs/analysis.rst b/docs/analysis.rst index 4dcbd2ec..7f784499 100644 --- a/docs/analysis.rst +++ b/docs/analysis.rst @@ -11,7 +11,7 @@ The analysis module provides tools to characterize the type of holes. The MNAR case is the trickiest, the user must first consider whether their missing data mechanism is MNAR. In the meantime, we make assume that the missing-data mechanism is ignorable (ie., it is not MNAR). If an MNAR mechanism is suspected, please see this article :ref:`An approach to test for MNAR [1]` for relevant actions. -Then Qolmat proposes a test to determine whether the missing data mechanism is MCAR or MAR. +Then Qolmat proposes two tests to determine whether the missing data mechanism is MCAR or MAR. 2. How to use the results ------------------------- @@ -50,7 +50,11 @@ The best-known MCAR test is the :ref:`Little [2]` test, and it h b. PKLM Test ^^^^^^^^^^^^ -The :ref:`PKLM [2]` (Projected Kullback-Leibler MCAR) test compares the distributions of different missing patterns on random projections in the variable space of the data. This recent test applies to mixed-type data. It is not implemented yet in Qolmat. +The :ref:`PKLM [2]` (Projected Kullback-Leibler MCAR) test compares the distributions of different missing patterns on random projections in the variable space of the data. This recent test applies to mixed-type data. The :class:`PKLMTest` is now implemented in Qolmat. +To carry out this test, we perform random projections in the variable space of the data. These random projections allow us to construct a fully observed sub-matrix and an associated number of missing patterns. +The idea is then to compare the distributions of the missing patterns through the Kullback-Leibler distance. +To do this, the distributions for each pattern are estimated using Random Forests. + References ---------- diff --git a/examples/pklm/p_value_validity/p_value_validity.ipynb b/examples/pklm/p_value_validity/p_value_validity.ipynb new file mode 100644 index 00000000..6dc082cf --- /dev/null +++ b/examples/pklm/p_value_validity/p_value_validity.ipynb @@ -0,0 +1,154 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "5cffe020", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from qolmat.analysis.holes_characterization import PKLMTest" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d347e25d", + "metadata": {}, + "outputs": [], + "source": [ + "def generate_randunif(N=500):\n", + " Z = np.random.uniform(0,1, N)\n", + "\n", + " X = np.sort(Z)\n", + " F = np.array(range(N))/float(N)\n", + " \n", + " return X, F" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "16b5f189", + "metadata": {}, + "outputs": [], + "source": [ + "def create_df_with_nan(n_rows: int, n_cols: int, nan_ratio: float) -> pd.DataFrame:\n", + " data = {f\"Colonne_{i}\": np.random.normal(size=n_rows).astype(float) for i in range(n_cols)}\n", + " df = pd.DataFrame(data)\n", + " nb_valeurs_manquantes = int(nan_ratio * df.size)\n", + " indices_valeurs_manquantes = np.random.choice(df.size, nb_valeurs_manquantes, replace=False)\n", + " df.values.flat[indices_valeurs_manquantes] = np.nan\n", + " return df" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "797a1534", + "metadata": {}, + "outputs": [], + "source": [ + "N_sim = 500\n", + "list_res_exact = []\n", + "list_res_approximation = []\n", + "\n", + "for _ in range(N_sim):\n", + " df = create_df_with_nan(500, 10, 0.35).to_numpy()\n", + " pklm_test_exact = PKLMTest(\n", + " nb_projections=100,\n", + " nb_permutation=30,\n", + " nb_trees_per_proj=200,\n", + " exact_p_value=True\n", + " )\n", + " pklm_test_approximation = PKLMTest(\n", + " nb_projections=100,\n", + " nb_permutation=30,\n", + " nb_trees_per_proj=200,\n", + " exact_p_value=False\n", + " )\n", + " list_res_exact.append(pklm_test_exact.test(df))\n", + " list_res_approximation.append(pklm_test_approximation.test(df))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "5a3a0cd3", + "metadata": {}, + "outputs": [], + "source": [ + "Z_pv_exact = np.array(list_res_exact, dtype=np.float32)\n", + "X_pv_exact = np.sort(Z_pv_exact)\n", + "F_pv = np.array(range(N_sim))/float(N_sim)\n", + "\n", + "Z_pv_approximation = np.array(list_res_approximation, dtype=np.float32)\n", + "X_pv_approximation = np.sort(Z_pv_approximation)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "d94f9f27", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAHHCAYAAAC1G/yyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACoVElEQVR4nOzdd3hT1RvA8W+S7payKXtvQQpl741MEaFlyJIhIgoyBYEyBEQFGT8Q2YiMliEoIBsEZI8iey9B2VAo3bm/P45ddNCM7vfzPHmanNx77sntTfLmTJ2maRpCCCGEEAIAfUoXQAghhBAiNZHgSAghhBAiGgmOhBBCCCGikeBICCGEECIaCY6EEEIIIaKR4EgIIYQQIhoJjoQQQgghopHgSAghhBAiGgmOhBBCCCGikeAoBfXo0YPChQtbNc+lS5ei0+m4efOmVfM117hx49DpdDHSChcuTI8ePZL82Ddv3kSn07F06dLItB49euDi4pLkx46g0+kYN25csh0vuq1bt+Lu7o6DgwM6nY5nz56lSDneJCXPkTXFda2nRt9++y1FixbFYDDg7u5u8v579+5Fp9Oxdu1a6xcuFYt43Xv37k3poqS49PKeTUiaD46uXbvGRx99RNGiRXFwcMDV1ZVatWoxc+ZMAgMDU7p4SWby5Mls2LAhpYuRbLZs2ZJq34ypsWyPHz/G09MTR0dH5syZw/Lly3F2dk6x8qTGc5QRbd++neHDh1OrVi2WLFnC5MmT49125cqVzJgxI/kKJ0QcdDodAwYMiPO5iMqA48ePx0h/9uwZffv2JWfOnDg7O9OgQQNOnjxp0nFtzC5xKrB582Y6dOiAvb093bp1o1y5coSEhHDgwAGGDRvGuXPnmD9/fkoXM0lMnjyZ9u3b07Zt2xjpXbt2pWPHjtjb26dMwRLh0qVL6PWmxeVbtmxhzpw5Jn3BFipUiMDAQGxtbU0soWkSKltgYCA2Nsn/Njt27BgvXrxg4sSJNG7cONmP/7rUeI4yot27d6PX61m0aBF2dnYJbrty5UrOnj3LoEGDkqdwQliB0WikZcuWnD59mmHDhpEjRw7mzp1L/fr1OXHiBCVKlEhUPmn2E+nGjRt07NiRQoUKsXv3bvLkyRP53CeffMLVq1fZvHlzCpYwZRgMBgwGQ0oXI0FJHbiFhYVhNBqxs7PDwcEhSY/1Jil1/AcPHgCQJUuWFDm+KVL6f5SRPHjwAEdHxzcGRkIkl6CgIOzs7Ez+wRyftWvXcvDgQdasWUP79u0B8PT0pGTJknh7e7Ny5cpE5ZNmm9W++eYbXr58yaJFi2IERhGKFy/OwIEDgbj7nkR4ve00ot/A5cuX+eCDD8icOTM5c+ZkzJgxaJrGnTt3ePfdd3F1dSV37txMmzYtRn7x9flJbHv1d999R82aNcmePTuOjo54eHjEatvX6XQEBASwbNkydDodOp0usg/P68dv1aoVRYsWjfNYNWrUoHLlyjHSfv75Zzw8PHB0dCRbtmx07NiRO3fuJFjmCAcOHKBKlSo4ODhQrFgxfvzxxzi3e73PUWhoKOPHj6dEiRI4ODiQPXt2ateuzY4dOwDVT2jOnDmRrz3iBlH/2++++44ZM2ZQrFgx7O3tOX/+fIL/9+vXr9OsWTOcnZ3JmzcvEyZMQNO0yOfj+3+9nmdCZYtIe7225NSpUzRv3hxXV1dcXFxo1KgRhw8fjrFNxP/xzz//ZPDgwZHVw++99x4PHz6M+x/wn/r169O9e3cAqlSpEuP6iK+/V/369alfv36s1+/r68ukSZPInz8/Dg4ONGrUiKtXr8ba/8iRI7Ro0YKsWbPi7OzM22+/zcyZM1PtOfruu+/Q6XTcunUr1nMjR47Ezs6Op0+fArB//346dOhAwYIFsbe3p0CBAnz++edvbLY35XMH4O7du3z44Ye4ublhb2/PW2+9xeLFixM8RoSwsDAmTpwYef0XLlyYUaNGERwcHOOYS5YsISAgIPJ/EFfZQF0Pmzdv5tatW5Hbvt4/0mg0JvraeOedd8icOTNOTk7Uq1ePP//8842vKeIa9PHxYdSoUeTOnRtnZ2fatGnzxs+k48ePo9PpWLZsWazntm3bhk6nY9OmTQDcunWL/v37U6pUKRwdHcmePTsdOnRIVL/NxL6fAIKDg/H29qZ48eKR19Hw4cNj/I8AduzYQe3atcmSJQsuLi6UKlWKUaNGJVgOc77jrl69So8ePciSJQuZM2emZ8+evHr1KlaZP//8c3LmzEmmTJlo06YNf//9d5xlSMz1G/E/Xb16NaNHjyZfvnw4OTnh7++f4Oszxdq1a3Fzc6Ndu3aRaTlz5sTT05ONGzfGOt/xSbM1R7/99htFixalZs2aSZK/l5cXZcqU4euvv2bz5s189dVXZMuWjR9//JGGDRsydepUVqxYwdChQ6lSpQp169a1ynFnzpxJmzZt6NKlCyEhIaxevZoOHTqwadMmWrZsCcDy5cvp3bs3VatWpW/fvgAUK1Ys3tfRrVs3jh07RpUqVSLTb926xeHDh/n2228j0yZNmsSYMWPw9PSkd+/ePHz4kNmzZ1O3bl1OnTqVYC3EmTNnaNq0KTlz5mTcuHGEhYXh7e2Nm5vbG1/zuHHjmDJlSuRr8vf35/jx45w8eZImTZrw0Ucfce/ePXbs2MHy5cvjzGPJkiUEBQXRt29f7O3tyZYtG0ajMc5tw8PDeeedd6hevTrffPMNW7duxdvbm7CwMCZMmPDG8kaXmLJFd+7cOerUqYOrqyvDhw/H1taWH3/8kfr16/PHH39QrVq1GNt/+umnZM2aFW9vb27evMmMGTMYMGAAPj4+8R7jyy+/pFSpUsyfP58JEyZQpEiReK+PN/n666/R6/UMHTqU58+f880339ClSxeOHDkSuc2OHTto1aoVefLkYeDAgeTOnZsLFy6wadMmBg4cmCrPkaenJ8OHD8fX15dhw4bFeM7X15emTZuSNWtWANasWcOrV6/4+OOPyZ49O0ePHmX27Nn8/fffrFmzxpTTGa/79+9TvXr1yP4VOXPm5Pfff6dXr174+/u/sWmrd+/eLFu2jPbt2zNkyBCOHDnClClTuHDhAr/88gugPjfmz5/P0aNHWbhwIUC8n59ffvklz58/5++//+b7778HiDWQITHXxu7du2nevDkeHh54e3uj1+tZsmQJDRs2ZP/+/VStWvWN52bSpEnodDpGjBjBgwcPmDFjBo0bN8bPzw9HR8c496lcuTJFixbF19c38odCBB8fH7JmzUqzZs0A1QR98OBBOnbsSP78+bl58yY//PAD9evX5/z58zg5Ob2xjG9iNBpp06YNBw4coG/fvpQpU4YzZ87w/fffc/ny5cj+o+fOnaNVq1a8/fbbTJgwAXt7e65evZqoYNJUnp6eFClShClTpnDy5EkWLlxIrly5mDp1auQ2vXv35ueff6Zz587UrFmT3bt3R34PRWfq9Ttx4kTs7OwYOnQowcHBb6zJDAoK4tGjR7HSX758GSvt1KlTVKpUKVZNVNWqVZk/fz6XL1+mfPnyCR4PAC0Nev78uQZo7777bqK2v3HjhgZoS5YsifUcoHl7e0c+9vb21gCtb9++kWlhYWFa/vz5NZ1Op3399deR6U+fPtUcHR217t27R6YtWbJEA7QbN27EOM6ePXs0QNuzZ09kWvfu3bVChQrF2O7Vq1cxHoeEhGjlypXTGjZsGCPd2dk5xnHjO/7z5881e3t7bciQITG2++abbzSdTqfdunVL0zRNu3nzpmYwGLRJkybF2O7MmTOajY1NrPTXtW3bVnNwcIjMT9M07fz585rBYNBev8wKFSoUo+wVKlTQWrZsmWD+n3zySax8NC3qf+vq6qo9ePAgzuei/9+7d++uAdqnn34amWY0GrWWLVtqdnZ22sOHDzVNi/v/FV+e8ZVN02JfX23bttXs7Oy0a9euRabdu3dPy5Qpk1a3bt3ItIj/Y+PGjTWj0RiZ/vnnn2sGg0F79uxZnMd7ff9jx47FSH/93EeoV6+eVq9evcjHEa+/TJkyWnBwcGT6zJkzNUA7c+aMpmnqvVGkSBGtUKFC2tOnT2PkGb3cqfEc1ahRQ/Pw8IiRdvToUQ3Qfvrpp8i019+TmqZpU6ZMifH+0bSoz44Ipnzu9OrVS8uTJ4/26NGjGNt17NhRy5w5c5xliODn56cBWu/evWOkDx06VAO03bt3R6Z1795dc3Z2jjev6Fq2bBnr80nTEn9tGI1GrUSJElqzZs1i/H9evXqlFSlSRGvSpEmCx484Tr58+TR/f//IdF9fXw3QZs6cmeD+I0eO1GxtbbUnT55EpgUHB2tZsmTRPvzwwxjled2hQ4diXQdxfSYk9v20fPlyTa/Xa/v374+x3bx58zRA+/PPPzVN07Tvv/9eAyI/hxLLnO+46OdA0zTtvffe07Jnzx75OOK66t+/f4ztOnfubPb1G3EOixYtmuA1/Xr533SL/jnn7Owc67VpmqZt3rxZA7StW7cm6rhpslktogouU6ZMSXaM3r17R943GAxUrlwZTdPo1atXZHqWLFkoVaoU169ft9pxo/8Sevr0Kc+fP6dOnTom97SP4OrqSvPmzfH19Y3RbOTj40P16tUpWLAgAOvXr8doNOLp6cmjR48ib7lz56ZEiRLs2bMn3mOEh4ezbds22rZtG5kfQJkyZSJ/nSUkS5YsnDt3jitXrpj1GgHef/99cubMmejto49+iPi1ExISws6dO80uw5uEh4ezfft22rZtG6OpM0+ePHTu3JkDBw7Eql7u27dvjCaoOnXqEB4eHmdzUFLo2bNnjF91derUAYi85k+dOsWNGzcYNGhQrJpFc4a1J+c58vLy4sSJE1y7di0yzcfHB3t7e959993ItOjvyYCAAB49ekTNmjXRNI1Tp06Z/Bpfp2ka69ato3Xr1miaFuP916xZM54/f57g+3/Lli0ADB48OEb6kCFDAJKs7+Wbrg0/Pz+uXLlC586defz4ceRrCggIoFGjRuzbty/e2t3ounXrFuOzvn379uTJkyfydcfHy8uL0NBQ1q9fH5m2fft2nj17hpeXV2Ra9P9vaGgojx8/pnjx4mTJksXsz93XrVmzhjJlylC6dOkY/9+GDRsCRH6+RryHNm7cmKhzY4l+/frFeFynTh0eP34c+f6KOL+fffZZjO1erwUy5/rt3r17vLV+cXn33XfZsWNHrNvrtb6gBnjE1a81om9jYkexp8ngyNXVFYAXL14k2TGif8kDZM6cGQcHB3LkyBErPaJvgjVs2rSJ6tWr4+DgQLZs2ciZMyc//PADz58/NztPLy8v7ty5w6FDhwA1/cGJEydifEBcuXIFTdMoUaIEOXPmjHG7cOFCZAffuDx8+JDAwMA4RwGUKlXqjeWbMGECz549o2TJkpQvX55hw4bx119/mfQaixQpkuht9Xp9rH5YJUuWBEjS+aEePnzIq1ev4jwnZcqUwWg0xupL8fp1GNHUY81rLiFvOn5EYFGuXDmrHC85z1GHDh3Q6/WRzW+aprFmzZrIvk4Rbt++TY8ePciWLRsuLi7kzJmTevXqAVj0vozw8OFDnj17xvz582O993r27AmQ4Pvv1q1b6PV6ihcvHiM9d+7cZMmSJckC6Ted94gfO927d4/1uhYuXEhwcHCizt/rnys6nY7ixYtHvldfvnzJv//+G3mL6G9WoUIFSpcuHaN51cfHhxw5ckQGJaC+LMeOHUuBAgWwt7cnR44c5MyZk2fPnlnl/wvqXJw7dy7WeYj43In4/3p5eVGrVi169+6Nm5sbHTt2xNfXN0kCpTf9/yKuq9eb5F9/b5pz/ZryeQ2QP39+GjduHOtWtmzZWNs6OjrG2a8oKCgo8vnESJN9jlxdXcmbNy9nz55N1Pbx/YINDw+Pd5+4RnzFNwoseo2MOceKsH//ftq0aUPdunWZO3cuefLkwdbWliVLliS6h31cWrdujZOTE76+vtSsWRNfX1/0ej0dOnSI3MZoNKLT6fj999/jfJ1JOXFi3bp1uXbtGhs3bmT79u0sXLiQ77//nnnz5sWowUuIKb9CEsOS/6M1JeaaM0VCr8vcaz6lmVvGvHnzUqdOHXx9fRk1ahSHDx/m9u3bMfpchIeH06RJE548ecKIESMoXbo0zs7O3L17lx49eiT4pZXYaygijw8++CBW/5gIb7/9doKvJaHjJZU3nfeI1/Xtt9/GO9mkNT5XvvvuO8aPHx/5uFChQpGBk5eXF5MmTeLRo0dkypSJX3/9lU6dOsWYOuLTTz9lyZIlDBo0iBo1apA5c2Z0Oh0dO3Z8Y1CS2PeT0WikfPnyTJ8+Pc7tCxQoAKjPsX379rFnzx42b97M1q1b8fHxoWHDhmzfvj3ec26t7zgw/b1tzvVr7c/r6PLkycM///wTKz0iLW/evInKJ00GR6BGYc2fP59Dhw5Ro0aNBLeNiIhfnyE4KX5RWXKsdevW4eDgwLZt22JUCy5ZsiTWtqZ8EDo7O9OqVSvWrFnD9OnT8fHxoU6dOjEukmLFiqFpGkWKFIn8NZNYOXPmxNHRMc5msUuXLiUqj2zZstGzZ0969uzJy5cvqVu3LuPGjYsMjqz5wW80Grl+/XqM13n58mWAyBE5pvwfE1u2nDlz4uTkFOc5uXjxInq9PvJDMqlkzZo1zpmyb926Fe+oxoRE/Ko8e/ZsgvMppdZz5OXlRf/+/bl06RI+Pj44OTnRunXryOfPnDnD5cuXWbZsGd26dYtMjxhJmZDEXkMRI4HCw8PNmpOqUKFCGI1Grly5QpkyZSLT79+/z7NnzyhUqJDJeYLl77mIa8PV1dWiubZe/1zRNI2rV69GfuF269aN2rVrRz4f/YvXy8uL8ePHs27dOtzc3PD396djx44x8lu7di3du3ePMfI4KCgoUTPKJ/b9VKxYMU6fPk2jRo3eeF71ej2NGjWiUaNGTJ8+ncmTJ/Pll1+yZ8+eeM9jUnzHRVxX165di1Fb9Pp709Lr19rc3d3Zv38/RqMxRqfsI0eO4OTklOjvtzTZrAYwfPhwnJ2d6d27N/fv34/1/LVr1yKHEru6upIjRw727dsXY5u5c+davVwRHwjRjxUeHp6oySgNBgM6nS5GtH/z5s04Z8J2dnY2aTkILy8v7t27x8KFCzl9+nSMJjWAdu3aYTAYGD9+fKxfDpqm8fjx4wTL3axZMzZs2MDt27cj0y9cuMC2bdveWLbX83ZxcaF48eIxqkYjZne21hIY//vf/yLva5rG//73P2xtbWnUqBGgPhgMBkOirpnEls1gMNC0aVM2btwYo/nu/v37rFy5ktq1a8dozkkKxYoV4/Dhw4SEhESmbdq0KdHTNbyuUqVKFClShBkzZsR6/dGvo9R6jt5//30MBgOrVq1izZo1tGrVKsZM4hG/rqO/Fk3TIj9bEpLYzx2DwcD777/PunXr4qwNf9O0BC1atACINZt1RC1FXKOLEsPZ2dmiZiUPDw+KFSvGd999F+eooje9rgg//fRTjC4Ua9eu5Z9//qF58+YAFC1aNEZTS61atSK3LVOmDOXLl8fHxwcfHx/y5MkTa2SxwWCI9Zk3e/bsRNUSJ/b95Onpyd27d1mwYEGsPAIDAwkICADgyZMnsZ6PqHVLaAh6UnzHRZzfWbNmxUh//Tqz9Pq1tvbt23P//v0Yfc0ePXrEmjVraN26daLn2UuzNUfFihVj5cqVkUPuo8+QHTEBVPT5J3r37s3XX39N7969qVy5Mvv27YusLbCmt956i+rVqzNy5EiePHlCtmzZWL16NWFhYW/ct2XLlkyfPp133nmHzp078+DBA+bMmUPx4sVj9cHx8PBg586dTJ8+nbx581KkSJFYQ5yja9GiBZkyZWLo0KGRF3N0xYoV46uvvmLkyJHcvHmTtm3bkilTJm7cuMEvv/xC3759GTp0aLz5jx8/nq1bt1KnTh369+9PWFgYs2fP5q233npj/6GyZctSv359PDw8yJYtG8ePH2ft2rUxOk17eHgAqnNgs2bNMBgMsX4BJpaDgwNbt26le/fuVKtWjd9//53NmzczatSoyE7dmTNnpkOHDsyePRudTkexYsXYtGlTnH0/TCnbV199FTmPSf/+/bGxseHHH38kODiYb775xqzXY4revXuzdu1a3nnnHTw9Pbl27Ro///yz2UP99Xo9P/zwA61bt8bd3Z2ePXuSJ08eLl68yLlz5yKD49R6jnLlykWDBg2YPn06L168iPWjoXTp0hQrVoyhQ4dy9+5dXF1dWbduXaL7fCX2c+frr79mz549VKtWjT59+lC2bFmePHnCyZMn2blzZ5xfmhEqVKhA9+7dmT9/Ps+ePaNevXocPXqUZcuW0bZtWxo0aGDaSfmPh4cHPj4+DB48mCpVquDi4hKjVu1N9Ho9CxcupHnz5rz11lv07NmTfPnycffuXfbs2YOrqyu//fbbG/PJli0btWvXpmfPnty/f58ZM2ZQvHhx+vTpk6hyeHl5MXbsWBwcHOjVq1esId6tWrVi+fLlZM6cmbJly3Lo0CF27txJ9uzZ35h3Yt9PXbt2xdfXl379+rFnzx5q1apFeHg4Fy9exNfXl23btlG5cmUmTJjAvn37aNmyJYUKFeLBgwfMnTuX/Pnzx6gdi68s1vyOc3d3p1OnTsydO5fnz59Ts2ZNdu3aFedcVpZcv9bWvn17qlevTs+ePTl//nzkDNnh4eExml/fKFFj2lKxy5cva3369NEKFy6s2dnZaZkyZdJq1aqlzZ49WwsKCorc7tWrV1qvXr20zJkza5kyZdI8PT21Bw8exDvM8fWhlPENga1Xr5721ltvxUi7du2a1rhxY83e3l5zc3PTRo0ape3YsSNRQ/kXLVqklShRQrO3t9dKly6tLVmyJNYQYU3TtIsXL2p169bVHB0dNSByOGl8UwlomqZ16dIlcuhzfNatW6fVrl1bc3Z21pydnbXSpUtrn3zyiXbp0qV494nwxx9/aB4eHpqdnZ1WtGhRbd68eXGW/fXhr1999ZVWtWpVLUuWLJqjo6NWunRpbdKkSVpISEjkNmFhYdqnn36q5cyZU9PpdJF5Rgxh/fbbb2OVJ76h/M7Oztq1a9e0pk2bak5OTpqbm5vm7e2thYeHx9j/4cOH2vvvv685OTlpWbNm1T766CPt7NmzsfKMr2yaFnsYraZp2smTJ7VmzZppLi4umpOTk9agQQPt4MGDMbaJbyh+fFMMvC6+/TVN06ZNm6bly5dPs7e312rVqqUdP3483qH8a9asibFvfEOGDxw4oDVp0kTLlCmT5uzsrL399tva7NmzU/U5irBgwQIN0DJlyqQFBgbGev78+fNa48aNNRcXFy1Hjhxanz59tNOnT8c6D3Fd64n93NE0Tbt//772ySefaAUKFNBsbW213Llza40aNdLmz5//xtcQGhqqjR8/XitSpIhma2urFShQQBs5cmSMz0BNM20o/8uXL7XOnTtrWbJk0YDIzypTr41Tp05p7dq107Jnz67Z29trhQoV0jw9PbVdu3YlePyI46xatUobOXKklitXLs3R0VFr2bJljCkU3uTKlSuRQ74PHDgQ6/mnT59qPXv21HLkyKG5uLhozZo10y5evBjrcyq+6yox7ydNU9OyTJ06VXvrrbc0e3t7LWvWrJqHh4c2fvx47fnz55qmadquXbu0d999V8ubN69mZ2en5c2bV+vUqZN2+fLlN75OS7/j4vruCAwM1D777DMte/bsmrOzs9a6dWvtzp07Zl+/8V07CQG0Tz75JM7n4vsMePLkidarVy8te/bsmpOTk1avXr04PwsTovvv4EIIIUSqsXfvXho0aBBjGQghkkua7XMkhBBCCJEUJDgSQgghhIhGgiMhhBBCiGikz5EQQgghRDRScySEEEIIEY0ER0IIIYQQ0aTZSSDNZTQauXfvHpkyZUr2tYiEEEIIYR5N03jx4gV58+aNNZmntWW44OjevXtJvn6VEEIIIZLGnTt3yJ8/f5IeI8MFR5kyZQLUybX2OlahoaFs376dpk2bYmtra9W8RRQ5z8lDznPykPOcfORcJw+rnufr16FnTzh3Dv/Royng7R35PZ6UMlxwFNGU5urqmiTBkZOTE66urvLGS0JynpOHnOfkIec5+ci5Th5WO89r1kDv3pAjBxw6BCVKgLd3snSJkQ7ZQgghhEg9goKgf3/w9IR33oGTJ+G/xauTS4arORJCCCFEKnXligqKLlyAH36Ajz6CFBg8JcHRa8LDwwkNDTVr39DQUGxsbAgKCiI8PNzKJRMR5DzHZjAYsLGxkRGYQoi0a9Uq6NsX8uSBw4fB3T3FiiLBUTQvX77k77//xtxJwzVNI3fu3Ny5c0e+pJKQnOe4OTk5kSdPHuzs7FK6KEIIkXiBgfDZZ7BwIXTuDPPmQTJ0uk6IBEf/CQ8P5++//8bJyYmcOXOa9aVrNBp5+fIlLi4uST4HQ0Ym5zkmTdMICQnh4cOH3LhxgxIlSsh5EUKkDRcuqGa0a9dUcPThhynSjPY6CY7+ExoaiqZp5MyZE0dHR7PyMBqNhISE4ODgIF9OSUjOc2yOjo7Y2tpy69atyHMjhBCp2k8/wccfQ6FCcPQolCuX0iWKJN8sr5FmGpFWSaAohEgTAgLU3EXdu6tao2PHUlVgBFJzJIQQQojkcvasCohu3YJly6Bbt5QuUZzkp2YqVrhwYUqVKoW7u3vk7cyZM1Y/jp+fH6tXr7Z6vqZaunQpbdu2TeliCCGEsDZNg0WLoEoVMBjg+PFUGxiB1BzFSdM0XoW+Mnk/o9FIQGgAhhBDopo4nGyd3tiM5+Pjg3sSD2f08/Njw4YNdOzYMUmPI4QQIgN68UL1LVqxAvr0gZkzwcy+vcklRYOjffv28e2333LixAn++ecffvnllzfWHOzdu5fBgwdz7tw5ChQowOjRo+nRo4dVy/Uq9BUuU1ysmmdcXo58ibOds8n7Xbp0iUaNGrFv3z6KFi3Kd999x86dO9myZQvnzp3j448/5tWrVwQFBdG5c2dGjx4NQEhICF9++SW///47BoOBPHny8NNPPzF27FieP3+Ou7s71atXZ968eTGON27cOM6cOcPTp0+5d+8eJUqUYOnSpWTPnj3O8q1YsYJVq1axadMmQAWbxYoV45dffsHNzY1OnTrh7+9PUFAQDRo0YNasWbGCyb179zJo0CD8/PwAOHv2LK1ateLmzZsA7Nq1i++//57AwEAMBgNTp06lQYMGJp9LIYQQSej0adWMdu+eCo46d07pEiVKijarBQQEUKFCBebMmZOo7W/cuEHLli1p0KABfn5+DBo0iN69e7Nt27YkLmnK8fLyitGsFhgYSKlSpfj222/x9PRk7969zJkzh+XLl6PX6ylcuDC7du3i5MmTnDhxgnXr1nH48GEApkyZwuXLlzlx4gSnT59m+fLl5MqViwkTJkSe09cDowj79+9n5cqVXLx4kQIFCjBy5Mh4y9yuXTsOHz7Mv//+C6hAJ2vWrFSoUIEsWbLw22+/ceLECf766y9u3ryJr6+vSefk+vXrTJ06lU2bNnHixAlWrlxJ586dCQ4ONikfIYQQSUTT0M+fD9WqqVqiEyfSTGAEKVxz1Lx5c5o3b57o7efNm0eRIkWYNm0aAGXKlOHAgQN8//33NGvWzGrlcrJ14uXIlybvZzQa8X/hj2sm10Q3q71JfM1qnTp1Ys+ePTRr1oxdu3aRM2dOAAIDA+nfvz9+fn7o9Xru3LmDn58f1atXZ9OmTUydOhV7e3uAyH0So2XLluTOnRuAvn370q5du3i3dXR05P3332f58uUMGzaMpUuX0rNnT0CdoxEjRnDgwAE0TePBgweUK1fOpCa9bdu2cf36derXrx+ZptfruX37NiVKlEh0PkIIIawjYu7kV68g5OFzKn/3HYY//1RrpE2bBmlsepE01efo0KFDNG7cOEZas2bNGDRoULz7BAcHx6hR8Pf3B9S8RtGXCYmY58hoNKJpGo42preHappGuG14ovoSRWz/ptm4jUYjRqMxVnpYWBhnz54lW7Zs3LlzJ3KbkSNHkj17dk6cOIGNjQ3vv/8+gYGBkc/HlV/Ea47rONHLGT0PnU4X7/YAPXr0oFevXnz00Uds2rSJadOmYTQamTZtGvfv3+fQoUM4ODgwZMiQyPJFL4deryc8PDzyGK9evYpR/vr16+Pj4xPrPCdUpvQu4vyFhoZiMBgszi/i/WHucjoiceQ8Jx8519alaXD8uI6ZM3X88ouO0FAdjap8TbvO3jjUMtKl33L0Xl5qYyuc8+T8v6Wp4Ojff//Fzc0tRpqbmxv+/v4EBgbGOXnjlClTGD9+fKz07du34+QUVXNjY2ND7ty5efnyJSEhIRaV88WLFxbtH8FoNBIQEBAZ0EU3ZswYihQpwv/+9z9at25NqVKlKFq0KA8fPqRIkSK8evWKK1eusHPnTqpWrYq/vz9NmzZl+vTplC9fHnt7ex49ekSOHDmwtbXlyZMncR4HVIC5ZcsWrl69Sq5cufjhhx+oU6dOvNuDqtUzGo0MHDiQevXqYWNjg7+/P/fv3ydbtmyEhIRw584dfH19adOmTWQfpLCwMPz9/cmZMye3bt3i+vXr5MiRg8WLF6uaOX9/atWqxYQJEzh06BDl/psb48SJE3gk86rNqU1ISAiBgYHs27ePsLAwq+W7Y8cOq+Ul4ifnOfnIubaM0QgbNhRn/frivHxp/1+qxgD+R5Mio3F1N3I+nzOLH2Uj75YtVjtuxI/k5JCmgiNzjBw5ksGDB0c+9vf3p0CBAjRt2hRXV9fI9KCgIO7cuYOLi4vZswtrmsaLFy/IlCmTVSaT1Ov19O7dO0bQN23aNAICAtizZw+HDx/GycmJ6dOn07t3bw4cOIC3tzfdu3fH19eXokWL0qBBAxwcHHB1dWXMmDGMHj2ahg0bYmtrS548edi8eTOtWrXihx9+oG7dutSoUYMffvghRjns7e2pU6cOH3/8MXfv3qV48eIsWbIkxvmLy4cffsiIESPYvHlz5LZDhw7F09OTWrVqkSdPHpo0aYKdnR2urq44ODhgY2ODq6srrq6uDBs2jCZNmuDm5sY777yDXq/H1dWVChUqsGDBAoYOHUpgYCAhISG4u7tn+A7ZQUFBODo6UrduXavMkB0aGsqOHTto0qQJtra2ViihiIuc5+Qj59pymga1axs4diyq60gWnrKIXrTjF7ypTAOOY0RHj671sbe33nlO6Ae51WmpBKD98ssvCW5Tp04dbeDAgTHSFi9erLm6uib6OM+fP9cA7fnz5zHSAwMDtfPnz2uBgYGJzut14eHh2tOnT7Xw8HCz80iNvL29Y533lJRez7OlrHENRxcSEqJt2LBBCwkJsUp+Im5ynpOPnGvLnD+vacWKaZoKkdStCke06xTWnpBFe5dftFqeI7U9e9Dm/eJk9fMc3/d3UkhTk0DWqFGDXbt2xUjbsWMHNWrUSKESCSGEEOmX0aiWQMuaFcqWVevDKhqfM50/qcV93KjIKTbSNnI/gyHh/rSpXYo2q718+ZKrV69GPr5x4wZ+fn5ky5aNggULMnLkSO7evctPP/0EQL9+/fjf//7H8OHD+fDDD9m9eze+vr5s3rw5pV5ChjBu3Lg40ydMmMD69etjpa9bt45ixYolcamEEEIkpbNnoVkzNUVRdNl4zFJ60JpNfMcQRjGZUOxibKPTSXBktuPHj8foJxLRN6h79+4sXbqUf/75h9u3b0c+X6RIETZv3sznn3/OzJkzyZ8/PwsXLrTqMH6ReGPHjmXs2LEpXQwhhBBWNn48xPW7uAYHWU1HnAmgFb+xmVbJXrbkkKLBUf369RMcyr506dI49zl16lQSlkoIIYTImMLDVWA0cWLMdB1GhvEtk/iSw1SnE6v4mwKx9m/UKJkKmsTS/Wg1IYQQQrzZ3Lnw2WcqQIouBw/5iW40ZyuTGYk34wkj9ii0YsWgQgUgMHnKm5QkOBJCCCEysDt31MoeBw7Efq4O+1hFJ+wIoRlb2U783ViWLYNn9qSL4ChNjVYTQgghhGWMRpg0CUqUUKt6FCwYOzDSE86XfMUeGnCFElTgdJyBUfbs0L49zJwJNWsm0wtIBlJzJIQQQmQQf/0FI0dCQhNX5+I+P/MBjdjFRMYwkTGExxEu7NsHdeokYWFTkNQcCQoXLoyfn1+yHvPevXvUifauGjduHEFBQZGPe/TowYwZM6x6zA0bNnD48GGr5BUSEkKrVq0oX748n3zyCfPmzePbb78FwM/Pj9WrV8fY/vXXN3bsWFasWGGVsgghRGLMn6/6BCUUGDVgN6epQHnO0IQdjGN8nIHRypXpNzACqTmKk6ZpGI2mr+FiNBoJDw8gPNxAYubX1OsTt0BtepQ3b172798f+Xj8+PEMGjTIKstexGfDhg24u7tTvXr1N24bFhaGjU38b49Tp05x5coVLl26FOs5Pz8/NmzYQMeOHSPTXn99EyZMMOMVCCGE6fz9VdNXQkvK6QlnLBMYw0R205AP+Jn75I5z2wsXoHTpJCpsKiHBURyMxlfs3++S5MepU+clBoNzvM/rdDq+/PJLNm/eTEBAAN7e3nTp0iXe7UuWLMnKlSupXLkyoKZC2LhxI7/88gvTp09n1apVhIaGYmtry6xZs+KcWbx+/foMGjSItm3bAtC+fXtatWpFjx49ePHiBYMHD+b06dMEBQVRvXp1/ve//2FnF3Pyr86dO9OqVSs6d+7M3LlzGTRoEE+fPsXZ2ZmGDRsybtw4ChYsiLu7O8+ePaNfv37/nY86GAwGtm/fDsCFCxdo1KgRd+7coVy5cqxevRo7OztevnzJgAEDImu7OnTogLe3d4Llz5UrF7/++is7duxg6dKlDBgwgN69e8d67W+//TbHjh3D0dGR3bt389133+Hr60tYWBi5cuXixx9/JCAggC5dunD37l3c3d0ZPHgw169f59mzZ4waNYqxY8fy/PnzWIFY9Nc3fPhw3N3dGTRoEOPGjePChQu8evWKa9eukTt3btauXUu2bNkIDQ1l4MCB7Ny5k2zZslGrVi1OnDjB3r17470OhBAiwrZt4OUFz5/Hv00e7rGCLtRlH96MZzKjMGKIc9tff03/gRFIs1qqp9PpOHXqFFu3buXTTz/l5s2b8W7bo0ePGHNDLVmyhA8//BCArl27cuzYMfz8/Jg9ezY9e/Y0uSxDhgyhTp06HD16lNOnT2M0Gpk5c2as7Ro3bszOnTsBtbxL5cqV+eOPP3j16hWnT5+OFZTNmzcPgP379+Pn50euXLkAVQPz22+/ceHCBe7fv8+6desA+OqrrwgODsbPz48jR46wYcMGfHx8Eix7ixYtaNOmDcOGDcPPzy9WYBTh8uXL7Nu3j927d7Ny5UouXbrEoUOHOHnyJF26dKF///6ULVuWhQsXUqpUKfz8/OjWrVvk/rly5WLChAk0aNAAPz8/5s2bF+/ri+7IkSMsXbqU8+fPRwZhAPPnz+fKlSucO3eO/fv389dffyX4OoUQAtQItA8/hHfeSTgwaso2/HCnJJdpyG6+Yky8gVH16tC6dRIVOJWRmqM46PVO1Knz0uT9jEYj/v7+uLq6otcnrlntTSK+xIsWLUrdunXZt28fhQsXjnPbbt26UbFiRaZNm8bdu3e5fPkyzZs3B1Qz0KRJk3j8+DE2NjZcunSJwMBAHB0dE/36NmzYwKFDh5g+fToAgYGBGAyx30SNGzdm/PjxhIeHc/78eSZNmsTOnTsxGAxUrVo10athv/feezg5qXNUtWpVrv23qM+uXbsYP348er0eZ2dnunXrxo4dO/Dy8kr0a4nPBx98EFm+DRs2cOzYMTw8PAAIf33yDyt65513yJ49O6DWEDxz5gygXmv0MnXv3p2FCxcmWTmEEGlbSIjqcD1zZuz5iqIzEMYExjKKKWylGV1ZziNyxrt9mzYwa1YSFDiVkuAoDjqdLsHmrvj3M2IwhGMwOCcqODJHQn2U8ufPT+XKldm4cSPnzp3jgw8+wMbGhpCQENq1a8eePXuoUqUK/v7+ZM6cmeDg4FjBkY2NTYwgIHonYk3TWLduHSVLlkywjAULFsTe3p4VK1bg4eFBo0aNmDRpEgaDgUYmTJ8avf+RwWAgLCwszu2in5OEyp8YLi5RzamapjFy5Ej69u1rUh7mMOe1CiHE68aPh/9+v8YrH3+zik7U4BBfMIVvGI6WQEPSmDGQ0bpJSrNaKrdkyRIAbt68yf79+2OM8IpLz549Wbx4MT/99FNkk1pQUBAhISEULFgQgNmzZ8e7f/HixTly5AigFgI+EG3yi7Zt2zJ16tTIL+6nT5/GWDg4usaNGzN27FgaN25M1qxZsbW1Zc2aNTRu3DjO7TNlysTzhOp+o2nUqBHLly9H0zQCAgJYvnw5TZs2fWP5XV1dE32MiNc7b948njx5AkBoaGiilq6J6zimvL7oGjZsyMqVKwkNDSU0NDRyEWYhhIjw77/g7Q21asHkyQlv24LN+OFOYW5Sjz+YyhcJBkYATZpYsbBphARHqVx4eDgVK1akadOmzJo1K94mtQjvvvsux44dw83NjTJlygDqy/qrr76iatWqeHh4xOpAHd3w4cPZs2cP5cuXZ+TIkVSrVi3yue+//x5HR0fc3d15++23adSoUbx9oBo3bsytW7cig6HGjRsTEBBAhQoV4tx+yJAhNGnSBHd3dx48eJDgaxw9ejS2trZUqFCBatWq0aZNGzw9Pd9Y/q5du+Lr60vFihUT1TTVpUsXevToQYMGDahQoQLu7u7s3r37jfs1atSI4OBg3n777cjO5qa8vug++ugjChcuTNmyZalVqxbFihUjS5Ysid5fCJH+tWqlanYOHox/GxtC+YZhbKYVh6iBO34cpNYb87azg3LlrFjYNEKnJbTyazoU0aT0/PlzXF1dI9ODgoK4ceMGRYoUMXs4ual9jt5Ep9Px9OlT+TJ8jbXPc2r34sULMmXKRGhoKF26dMHDw4MRI0bE2s4a13B0oaGhbNmyhRYtWiS6n5gwnZzn5JNezvWNG7B+PezfD7//rvoZJaQgt1hNRypznC/4mukMBt7cRF+tGgwYAB98kPiybfYbhfOzKVx76UC3Zv5WPc/xfX8nBelzJEQq17hxY4KDgwkKCqJ27dp89tlnKV0kIUQyu31bDcvfvFkNp09stUYbNrKUHjwnM3XYzxHePM9b7dqwbh3EMbA2w5DgKBWLq1LvwYMHkf1romvSpEnkDM0ifYnoQyWEyJgmTVKdok1p57ElhG8YziBm8gtt+ZDFPCPrG/ebMgVGjICMPvZDgqM0JleuXMm+1IcQQoiU8d13MHq0uq/TJS5AKsJ1fPCiAqf5jJnM5lMS04zWrRt88YVl5U0v0n+HDRNlsC5YIh0xGo0pXQQhhJWEhalanGHD1GMHh8QFRu+zllNUJBtPqMlBZvMZiQmMOneGBAYyZzhSc/QfW1tbdDodDx8+JGfOnGbNJ2M0GgkJCSEoKChDdBROKXKeY9I0jZCQEB4+fIher09wNKIQIvXTNNUJOvrE/2+ass2eIKYxhE+Yiy8d6MMC/MmcqOMdPAhxrCaVoUlw9B+DwUD+/Pn5+++/E1yiIyGapkXOOi2T9SUdOc9xc3JyomDBghIwCpFGPXsGCxeqWxxrWserOFfwxZMyXOBj5jKPfsRXW/R605yNDbi7W1Lq9EmCo2hcXFwoUaIEoaGhZu0fGhrKvn37qFu3bpoeJprayXmOzWAwYGNjI8GiEGnY559DtOUxE6Ujq5hPX/4hD9U5zGncE9w+emDk5ASbNoEJq0hlGBIcvcZgMMS5Xlhi9w0LC8PBwUG+tJOQnGchRHpz7x78/HPit3cgkJkMpC8LWEFn+jGPl2RK9P56PezZA1WrmlHYDECCIyGEECIFXbgAlSurTtiJUYqL+OJJCa7QmwUsoheJ6XQd3YgREhglRDonCCGEEMnsxQtYtAhKl4ayZeHVq8Tt15WfOIEHtoRSlaMsojemBEZFi8KCBWruJBE/qTkSQgghklhAgOpkfeGCuv30E9y5k/j9nQjgfwygJ0tZSnc+YQ6vcE70/rlzw9q1ULOmTPCYGBIcCSGEEEnk5k01X9HChWDuVGRlOYcvnhTmJt1Zyk90N2l/Bwe1HEjNmuYdPyOS4EgIIYSwkhcvYO9e2L5drYV25UrUc87OapHYxA+I1ujJEv7HAK5RjCoc4wJlTSpPpUqqGa1SJZN2y/AkOBJCCCHMdP067N4Nx46p25kzMTtWGwxQoYJK++uvxOfrwgt+4GM+YAUL6M1AZhKIU6L3d3WF+fOhQwc1Mk2YRoIjIYQQwgSaBj/+CHPmwNmzsZ8vWhSaNYOmTcHWFt59F8LDE5//25zGF0/yco/OrGAVnU0u4//+B15eJu8m/iPBkRBCCJFIYWFqvbMZM9RjgwFq1VLLb1SpoobkFyqkAqhWrWDLFlNy1+jLfGYykIuUxoMTXKGkSeUrVUp1vC5XzqTdxGskOBJCCCHe4MYNWL8efvgBrl1TaV99BR9/DNmyqceapkaivfsu/P67KX2LIBP+LKAPXvgyl48ZzHSCcUjUvra26ljt2qkarRw5THxxIhYJjoQQQog43LsHS5aokV6nTkWl58gBs2ZBp05RaefOwXvvxeyAnVgVOYkvnuTiAZ74sAbPRO1na6sCs/v3VYdrX19VkyUsJ8GREEII8Rp/fzWD9N276rFeD/XqqQ7O3bqpkWcAx4/Dl1+q0Wmm0/iEOUxjCGcoTzO2cZ1iidozWzZ48kQFRgaDmi5AAiPrkeBICCGEiEbToH17FRjZ26vOze++Czlzxt52wAA4csT0Y2TmGYvoxfusZxafMoxvCcE+UfvmzavK9eSJmrto8WLV10hYjwRHQgghBHD7Nixfrpb1uHEDbGzU4qw1asTeduNGGDIkqv+RKapwFB+8yMIz3mM9G3gv0fs6OMCjR2q+pDx5YPNmyJLF9DKIhMnsB0IIITK8YcOgcGEYPVoFRs7O8P33sQOjM2egTh1o29acwEhjEN9zgNo8IBcVOWVSYGRvD0FBKjCqWxcOHZLAKKlIzZEQQogM69o16NEDDhxQj+vWhZ49VbOai0vUdsHB0LUrrFlj3nGy8oSl9KANv/EdQxjFZEKxMymP4GBVm/XTT9Cxo6yRlpQkOBJCCJFhffyxCozs7NTQ/GHD4t6udm3V+docNTjIajriTACt+ZVNtDa7vH/9BWXKmL27SCRpVhNCCJHh/PwzeHjAjh3q8YkTcQdGT59C+fLmBUY6jAzjG/ZRlzsUwB0/swOjDh1UGSQwSh5ScySEECLDuHvXmeHD9TFmuB4yJO4ZpbdsURMrBgebfpwcPGQZ3WnB70zhC8YygTBszSrzL7+oPk4i+UhwJIQQIkM4dw4+/7w+ISFqQqB8+dTkjnEN0Q8NVcP3oy8im1h12McqOmFHCO/wO9t4x+wyL1kigVFKkGY1IYQQ6VZoqJqgcdw46NrVhpAQG95+W+Pnn+HSpbgDI1CzX5saGOkwMopJ7KEBVymOO34mB0bRO1nPnas6i4vkJzVHQggh0q0PP1T9ixQVecycGU79+nF//T18qJrSIkavJVYu7rOcrjRmJ5P4kvF4E27iV+w770CuXGo0WocOqrO4SBkSHAkhhEiXVq6MCoxatoTmzcPR6fZSq1bdyG3Cw9V8QcePw8mTaqh+UJBpx2nAblbQBR0aTdnOLhqbXNZs2eC336KWJYm+bptIfhIcCSGESHfu3YMuXdT9IUPgu+8gNNTIli0vAdVktm8fjBwJR4+adww94YxhImOZwB4a0IUV3Ce3yfm4u0ODBmrx2JAQNQt2ixbmlUlYh/Q5EkIIka4MGaI6W0PUoqyg1kw7cyY7n36qJ18+aNRIBUbmTKaYm3/YSWPGMBFvxtOU7SYHRlWrwqxZcOeOmo37zBlV3q++UrNhi5QjNUdCCCHSjf37Yfp0dd/DA4YPB1tbFRgNG6Zn1qzakdsaDKpZTdNMO0YTtvMzHxCGDY3YxR/UN7mcJUvC33/DZ5+pxy4u4O2tZuF2czM5O2FlEhwJIYRIF8LCooKNPn1g/nx1/+lT6NsX1q5VQ/gLF9a4eVNHeLhp+RsIYzzejGQKO2hCV5bzkFxmlfXy5aj7JUrA779DsWJmZSWSgARHQggh0oWZM8HPT3VuHj8+Kr1dO9i7F/R6jTx5XnDzZiaT887H36yiEzU4xCgm8w3D0SzomTJlCtSvr2a8zpzZ7GxEEpHgSAghRJp29y4MGAAbNqjHfftCnjzqvq+vCoxANZ/dvetqcv7N2cJPdCMQR+qzlz+p/ead4mBrq+ZdeustGDpULSIrUifpkC2EECLNunQJKlZUgZFOp2qJRo5Uz61cCV5eUdtqmmk9r20IZSrD2UJLDlOdipwyOzACFRiVLw+7dklglNrJv0cIIUSaNXy4mrixXDmYOBHs7GDxYti2DbZuNT/fgtxiNR2pzHGG8i3TGWxWM5qbm6rVKlZMNffVq6eG6ovUTYIjIYQQaYbRqDovb9kCBw/C6dMq3c4O3nvPOsdow0aW0BN/XKnDfo5Q3ax8mjRRS5eItEeCIyGEEGnClStqMdgLF2KmZ8miZre2lC0hTGUEnzODDbxLT5bwjKwm52Njo+Yq+uQTy8skUoYER0IIIdKE1atVYOTqCh98ANeuqeazZ88sz7swN/DBC3f8GMgMZvEZEWuxmaJCBbUum4uL5WUSKUc6ZAshhEj1goOj1kkbPBj++UcFRtbQjnWcoiLZeUxNDjKLgZgTGNWpA3v2SGCUHkhwJIQQItXbskVNnGhvDwsWwC+/WJ6nPUHMZgDraM8OmlCJk5ygssn5GAywdCn88QdkNb0VTqRCKR4czZkzh8KFC+Pg4EC1atU4+oYVAGfMmEGpUqVwdHSkQIECfP755wSZuoSyEEKINGPt2qhV6oOD1bxGlirGVQ5Sk94s5GPm4okv/pg+G2OuXKoWq3t389ZoE6lTigZHPj4+DB48GG9vb06ePEmFChVo1qwZDx48iHP7lStX8sUXX+Dt7c2FCxdYtGgRPj4+jBo1KplLLoQQIilpGvj4QKtW0KGDCoqsxYvVnKQSLrykOoeZx8eY04zWoIGqwcqZ03plE6lDigZH06dPp0+fPvTs2ZOyZcsyb948nJycWLx4cZzbHzx4kFq1atG5c2cKFy5M06ZN6dSp0xtrm4QQQqQte/dCx46webP18nQgkHl8xGo6sYlWeHCC07iblVfVqrB7N9Ssab3yidQjxUarhYSEcOLECUZGTGUK6PV6GjduzKFDh+Lcp2bNmvz8888cPXqUqlWrcv36dbZs2ULXrl3jPU5wcDDB0X5y+Pv7AxAaGkpoaKiVXg2ReUb/K5KGnOfkIec5ech5ju3+fejWzQZzanPiU4qL+OJJCa7Qh/kspHci89dibGdvrxEcrKN+/XBCQ41WK196EW6MWs03qb5jk0OKBUePHj0iPDwcNze3GOlubm5cvHgxzn06d+7Mo0ePqF27NpqmERYWRr9+/RJsVpsyZQrjo69A+J/t27fj5ORk2YuIx44dO5IkXxGTnOfkIec5ech5hvBwHatWlWTt2lJYMzD6gOX8wMfcoQBVOcpZypuwd0Q5VJAUHKyjYEF/3N33s2VLmNXKmF5cD7lGxf+aGa19Tb969cqq+SUkTc1ztHfvXiZPnszcuXOpVq0aV69eZeDAgUycOJExY8bEuc/IkSMZPHhw5GN/f38KFChA06ZNcXU1fQHChISGhrJjxw6aNGmCra2tVfMWUeQ8Jw85z8lDzrOiadCqlYEdO6zX28OJAGbzKR+yhGV04xPmEIDp4+yzZdN48kQFSVOmhPPpp47Y2TW1WjnTk9/PHATVQGP1azqi5Sc5pFhwlCNHDgwGA/fv34+Rfv/+fXLnzh3nPmPGjKFr16707t0bgPLlyxMQEEDfvn358ssv0etjv6ns7e2xt7ePlW5ra5tkH0RJmbeIIuc5ech5Th4Z/TyfOgXWrGgoyzl88aQwN+nBEpbRw+Q8MmUCvZ7IwGjQIPjiCwNgsF5B0xmDPurcWPuaTs73R4p1yLazs8PDw4Ndu3ZFphmNRnbt2kWNGjXi3OfVq1exAiCDQf0jNE1LusIKIYRIEkYjrFmjRqVFsGwSRY0eLOEYVdDQUYVjZgVGOXLAixfw/LmakXvAAPjmG0vKJdKSFG1WGzx4MN27d6dy5cpUrVqVGTNmEBAQQM+ePQHo1q0b+fLlY8qUKQC0bt2a6dOnU7FixchmtTFjxtC6devIIEkIIUTqFxQEx46pNcheX5z15Uvz8nTmJT/wMV35mYX04jNmEYjpfUtLloyacHL2bOjWTd0XGUeKBkdeXl48fPiQsWPH8u+//+Lu7s7WrVsjO2nfvn07Rk3R6NGj0el0jB49mrt375IzZ05at27NpEmTUuolCCGEMMHjx2puoGHDrLMmWoTy/IUvnuTnb7rwMyvpYtL+Op3q95Q9uwqMDAYZqp+RpXiH7AEDBjBgwIA4n9u7d2+MxzY2Nnh7e+Pt7Z0MJRNCCGENU6bApk1w6ZIKjqxLoy/zmclALlKaSpzkCiVNysHBAT77LGazWcOGEhhlZCm+fIgQQoj06+JFGDUKDh6MCowKFFA1NJbKhD+r6MSP9GMJPanBIZMDo5Yt1Zpo69erxxFl/Pxzy8sn0i4JjoQQQiSJLVugbl1138ZGjUi7exfy5bO8BqkiJzmBBy3Ygic+9OcHgnBM9P6OjtCvHxQqBPXqwdWrUc/VqQPNm1tWPpG2pXizmhBCiPQnOFh1ZH78GAoWVDUzxYqpZTfimec3kTQ+YQ7TGMJZytGc37lGcZNyMBg0qlbVMX++Gi0XoVAhtYDs8OGWlE+kBxIcCSGEsKpff1V9eB4/VrVEESO/atWyLDDKzDMW0pv2rGMWnzKMbwnB1GFkGuHhOv74IyrF0RGWL4f33ze/bCJ9keBICCGEVbx4AT/+qEaiAWTJAnPngq0tNG6s+h2ZqzLH8MGLbDyhHev4hXZm5qQmdNTroV07+OADqF8fMmc2v2wi/ZHgSAghhMXOnIE2beDmTfX4rbfg6FEVGFWsCH/9ZW7OGgOZyTcMxw93GrGLmxQxu5w6nYam6ejXD+bMMTsbkc5Jh2whhBAWCQ2FTz9VgZGTEyxYoOYyWrIEihQxPzDKyhM20JYZfM5sPqU2BywKjAA0TUfjxkZmzLAoG5HOSc2REEIIs5w6BaNHw759UbNat2kDPj7w0UcxOzubqjqHWE1HMvGC1vzKJlpbXF6DQaN+/dv8/HNebG2lbkDET4IjIYQQJjtzRvUjevJEPc6eXd1Wr7YsXx1GhjCNyYziKFXpxCruUNCyPHUqiOvdO4xTp/zIli2vZYUU6Z6EzkIIIRLt0CFo0ADeflsFRpkyQfXqqubo8mXL8s7OI36jNd8ynGkMoT57zQqMXlufnIkTYcIEyJPHsvKJjENqjoQQQsTrxQvYsUM1oZ08qSZ2BFUb07KlGq5/6JDlx6nNflbRCXuCac4WtmL+LIzRm/OqVYOPP7a8fCJjkeBICCFEnIKC1OzRp07FTK9fXy0J0q0b/PuvZcfQYWQkU5jAWP6kFp1YxT3yWZYpan6llSujZugWwhQSHAkhhIjl+XPVXBYxaWPnzlCjhppF+uBBaN1azYJtiVzcZzldacxOJvEl4/Em3ApfS+XLq1ouG/mGM1tQWBBH7x7lRfALAkIDeBnyMsYtIOS/tNCY6fltbvF50ZQuveXk0hFCCBHDy5fw3nsqMNLr1SKsL1/CrFlw5Yp1jlGfPaykM3qMNGMbO2liUX42NhAWBrlywe7dEhhZqu3qtmy7ts3k/Rz/W1DYXmfqzOWpi1w+QgghANVXZ9cuGDsWDh+OSps2zXrH0BPOaL5iLBP4g3p0YQX/YllP6datYft2FRx9+SXkyGGlwmZg155eA6B0jtK4ObvhYucS4+Zs6xw7zc4Z51A/tAdfkts+dwq/AstIcCSEEIIFC+Cbb2KuTh9BpwNNs/wYufmHFXShHn8wHm8m8SVGDGbnp9erfkWLFqkmvoYNoU8fy8spoixsvZBaBWslevtHj0I5+yAJC5RMJDgSQogM7vJl6Ns3dnq5cnD2rHUCo8bs4Gc+IBwDjdjFH9Q3O68iRaB3b/D0VOu3de+u0sePV4vICmEpmedICCEysNBQGDQoZlr16irt7NmY6QUKmJ6/gTAmMpptNMMPd9zxsygwKllSjZ4bNUoFQqNGqVqjwoVVuYWwBqk5EkKIDCosDJo3V/2MIkyeDM+eqSa2CHq9qj26c8e0/PPxNyvpTE0O8iWTmMoINDN/k2fODCNHqjmLrl5V0wicOxf1fKdO0glbWI9cSkIIkQGdO6dmjY4eGDVsCMePw/r1UWk5csCjR6bn35wt/EQ3gnCgPnv5k9om56HXqw7hBQqoUXK2tvDjj/DFF+Dvr7bx8IAWLVQNkhDWIsGREEJkIPfvq5mtT5yImW5np2a6DgyMSsufH/7+27T8bQhlEl8ynG/ZTAu6s4zHmDd8zGhUAVLfvupvkyZqmD6oJrT162VJkNTCaAwlPPwloaFPUrooViHBkRBCZCCrVqnAyNYWihaFGzcgJETdojMnMCrAbVbTkSocYyjfMp3BZjej5cunRp716qXKsmWLCoycnODrr6F/fzCYP9BNvEE++2DqFIbwR7O4FLiM8PCXcdwCIu9rWsgb80xLJDgSQogM4sULmDFD3Xd0hEuX4t/W1MCoNb+ylB68IBN12cdhaphVRr0eNm5UfaEigp8lS6BfP3W/Qwf49FOzshYm6JH/AYWdwOjvyz/+id9Pp7PFYHAlMLBe0hUuGUhwJIQQGUBwMMycCbduqVojfxO+8BJiSwhf8wWD+Z6NtKEnS3hKNrPzO31aTSEQ4dkzGDZM1WxVqACTJlleZvFmjga1eq/OpRmFctTGYHB57eYcZ5peb0doaChbIlYoTqMkOBJCiHTsyRMYOBB8faOazkJD4962SBHVzJZYhbmBD16448dAZjCLzwCdWeX08ICFC1Vg9O+/sG+fWsNt40Z4/Fj1LTp8GBwczMpemEnv+i6FC3+c0sVIdhIcCSFEOvXwoVooNnon67hErEtmSmD0HutZzIc8IRu1+JPjVDG7nN9/Dz17qhFpY8eqaQSiL2pboAD88osERiL5yCSQQgiRDl2+rEZ3RQRGjo4wZEjsTswGgwqMEsuOYGbxKet5n500phInLQqM3NxUB+ssWaBKFZg4UQVG2bPDgAGwYoWajNLDw+xDCGEyqTkSQoh0JCwMpk6FMWOilv2oWlV1xo5rAdnw8MTnXYyr+OBFOc7Snzn8wMeY04zWpQs8fw6bNqmpBSK4uakZsKtXVzVILi4mZy2EVUhwJIQQaVxgIPz5J+zdq+b+uXAh6rn69eGvv1TfI0t44sMC+nAfN6pzGD8qmpyHXg9Ll6qJHP/8U9VaDRyoZrcuWRJcXS0roxDWIsGREEKkYQ8eQI0acP167OeyZVO1L5YERg4E8j2f048fWUVHPuJHXmBeFLNpE2zbpgIjR0dYs0ZNSClEaiPBkRBCpGF9+6rAKEcOsLeHu3dVeq9eUKYMDB1qft4luYQvnpTkMn2Yz0J6Y2ozmk6nmvdq1lQzXs+cqdJnz5bASKRe0iFbCCHSqKVL1VB3UJMj3r2rRp59+incu2dZYNSFnzmBB/YEU40jLKQP5vQv0jRVs1W2rBqRBmp26169zC+bEElNao6EECINOncOPvxQ3ff0hD171H17e1UrYy5HXjGbT+nFYn6iK/2ZSwCm94zW68HLS91fs0at2wZqagFZJFakdlJzJIQQacy1a2o5DU2DevVUrczFi+q5gADz8y3DeY5RhU6soieL6c4yswIjgA0b1Ii0VavUCLqWLVWfo0uX1LppQqRmUnMkhBBpyJEj0KJFVCfrP/5QN8to9GApc/iEGxShMse5QFmzcrK1hUWLoGJFtVisTgc+PtC+vbovRFogwZEQQqQRd+9CgwZvnvHaFM68ZC796cZyFvEhnzKbQJzMyqtgQbh6NSpAAqhWTfWHEiItkWY1IYRI5e7dU81oRYvGDowaNVL9e8xRnr84RhXasZ4PWE5vFpkdGNWvD35+qi/UwoUwcqRKr1vXvLIJkZKk5kgIIVKhY8dg3Tr4/Xc1iePrMmdWi7Hu2mVO7hp9WMBMBnKZknhwgsuUsqi8w4apWqIrV2Kmt21rUbZCpAgJjoQQIhW5c0eNPjt8OO7n3dzU8hobN6oOz6bKhD8/8hGdWM08PuJzvicIR7PKWq+eWg+tcWMVpEUERg0bQuXKKr1GDbOyFiJFSXAkhBCpwI0bsGCBWqE+ofXO7t+PmtvIVO6cwhdP3LiPF6vxxcu8jFCdrD091Ui0lSth3jyVvmwZdOtmdrZCpAoSHAkhRAratk3H4MH1uH7dNkZ64cJQqZLq4By9Wc3R0ZwO2Rr9mct0BnOWcjTnd65R3Kzy6vVqbTRPTzVMf+RIuHVLPVetGnTsaFa2QqQq0iFbCCFSyN9/g6engevXswBaZHqHDqqv0fbtsfsbmRoYZeYZvngyhwHMpy81OWh2YAQwerSaV6lzZ3W7dUstXTJlippSwM7O7KyFSDWk5kgIIVLIlCkQGKhDp9PQNB05c8KsWWrBWHd3CA62LP/KHMMHL7LxhPdZy3retyi/7Nnh669hwoSotOHDYdw4VaMlRHohNUdCCJHMbt+Gjz6CuXPVY03TkT+/xqFDcPo0tGplaWCkMZAZ/EktHpGDipyyODACePwYQkJUTZGXF5w8CVOnSmAk0h+pORJCiGS0cqWas+jFi5jpbdoYqV7dwKNHluWflScs5kPaspHpfM4XfE0o1mnrql9f9TcqUUJmuxbpmwRHQgiRhB4/Vn2HDhxQt7jmLAKYO9dg8bGqcRgfvMjEC9qwkd9oY1F+ej0Yjer+smXwwQfmTzgpRFoiwZEQQiSBY8fgiy9UJ+WEhuar/kYAuhjBiCl0GBnCNCYzKnLh2NsUMrfokSLKMmaMDM8XGYv8BhBCiCQwZAjs3q0Co7JlY47icnOD9etV52tN0wGqU7Y5gVF2HvErbfiW4UxnMPX4wyqBUYQqVWJ2wBYiI5DgSAghrOzVKzh4UN2fMkV1wA4JUY+LFFEzST95ogIoRfsvSDJNLQ7ghzvVOUwLNvMFUwnD9s07vkG2bPDOO6rG6JdfLM5OiDRHmtWEEMKKHj+GoUNVjVGmTFELsAIUL66W2YjolB3FtMBIh5ERTGUiYzhITTqzkrvkt7jstWqpofq1akmHa5GxSc2REEJYgaap2qKqVWHpUpUWFhb1fIsWqsbop59eD4xMk5MH/E5zJvElX/MFDdltlcBo9mzVYbx2bQmMhJCaIyGEsFBwMLz3nprVGiBLFrUobPTZrIcOVYHRmDHmH6cee1lJZwyE04xt7KSJReUGFQg1aQL9+1uclRDphtQcCSGEBc6ehdatVWCk16vO18+e8d8INKheHSZOhBEjoHt3846hJ5yxjGcXjbhIadzxszgwMhigd2949Ai2bZMh+kJEJzVHQghhhuBgtcjqhg1RaUYjnD+v7tvZwaefwvTpcPiw+cdx419W0IUG7GE83nzFaIxYNidSo0aq3C4uFmUjRLolwZEQQiRSeDhs3AgrVsCWLRAUFPP5iHmKnJzUaK9p0yw7XiN2soIuGNHTiF3spYFlGaKaz2bNUjVHQoi4SXAkhBAJeP5cBUJ79qiZrm/divl87txqKZCAABUYubiojtjr15t/TANhjGMco5jMThrTleU8wM2yFwJMnqwmppQO10IkLMVbmefMmUPhwoVxcHCgWrVqHD16NMHtnz17xieffEKePHmwt7enZMmSbNmyJZlKK4TIaBo1gs6dYcECFRhlyQKF/ptjsXBh9TggQD12cYGXL2PXKJkiL3fZRSNGMoXRfMU7bLVKYLR1q5pWQAIjId4sRWuOfHx8GDx4MPPmzaNatWrMmDGDZs2acenSJXLlyhVr+5CQEJo0aUKuXLlYu3Yt+fLl49atW2TJkiX5Cy+ESPcCAuDECXX/449VJ+utW+HmTZWWPXvU86ACI0s0YyvL6Uow9tRnLweoY1mG/9HroXFjq2QlRIaQosHR9OnT6dOnDz179gRg3rx5bN68mcWLF/PFF1/E2n7x4sU8efKEgwcPYmurZoEtXLhwchZZCJFBhIXBoEHqft68cPcu/Pqrepw5s2puix4YWcKGUCYyhi+Yyhaa042feEwO62QOTJokfYyEMEWKBUchISGcOHGCkdGmj9Xr9TRu3JhDhw7Fuc+vv/5KjRo1+OSTT9i4cSM5c+akc+fOjBgxAkM87/zg4GCCg4MjH/v7+wMQGhpKaGioFV8RkflZO18Rk5zn5JGRz/PDh9CkiQ3nz6s2qNGjw+jfX31c6vUaz59br22qALdZRSeqcYRhfMM0hqBZqcdDkSJGfvrJSLVqGhnw3xhLRr6mzRUeHm7y+Uqq85yc/7cUC44ePXpEeHg4bm4x29Ld3Ny4ePFinPtcv36d3bt306VLF7Zs2cLVq1fp378/oaGheHt7x7nPlClTGD9+fKz07du34+TkZPkLicOOHTuSJF8Rk5zn5JERz/Pq1aU4f740mTKF0K7dZebOdQNyAmA0xh8Y2diEERaW+I/VVvzGMrrzgkzUYT+HqWFyWQ2GcMLD9by+BEnx4k+YOPEQjx+HId0yY8qI17SpNDu1CvLZs2d5ftm8C8ja5/nVq1dWzS8hOk2LmKosed27d498+fJx8OBBatSI+kAYPnw4f/zxB0eOHIm1T8mSJQkKCuLGjRuRNUXTp0/n22+/5Z9//onzOHHVHBUoUIBHjx7h6upq1dcUGhrKjh07aNKkSWSzn7A+Oc/JIyOe55AQ+O03HZ06qQDn7beNnDunIzxcBR52dtp/C8hGD0S0/x5rJHaNNFtCmMJIhjCdjbShJ0t4SrZE7Pn6MaI/Vvf1eg2jUcfRo6G4uyeqOBlGRrymzbV+tzNu9qHgNptaxT8yad+kOs/+/v7kyJGD58+fW/37+3UpVnOUI0cODAYD9+/fj5F+//59cufOHec+efLkwdbWNkYTWpkyZfj3338JCQnBzs4u1j729vbY29vHSre1tU2yN0dS5i2iyHlOHhnpPHftCj4+UY//+iuqectggJCQuIIf3Wt/E1aYG6ymIxU5xSC+ZyYDE71v7O10se4bjToaNoQqVTLG/8wcGematpTBYDD7XFn7PCfn/yzFhvLb2dnh4eHBrl27ItOMRiO7du2KUZMUXa1atbh69SpGozEy7fLly+TJkyfOwEgIIRIrPBx2746dHvHbKjzc8mO8x3pOUZFcPKAWfzKTQSQ+MHqzwoVh5syYs3YLIUyXovMcDR48mAULFrBs2TIuXLjAxx9/TEBAQOTotW7dusXosP3xxx/z5MkTBg4cyOXLl9m8eTOTJ0/mk08+SamXIIRIB+7cgfffVx2xXxetVd5sdgQzi09Zz/vsohEVOcVxqliecTT16sHly/DZZ5Apk1WzFiLDSdGh/F5eXjx8+JCxY8fy77//4u7uztatWyM7ad++fRt9tNUQCxQowLZt2/j88895++23yZcvHwMHDmTEiBEp9RKEEGncihXQq5d1gqC4FOMqPnhRjrN8wv+YS3+sWVvk5KSG6vftC9JSJIR1pPjyIQMGDGDAgAFxPrd3795YaTVq1OCwJas4CiFENNOnJ11g1AFfFtKbB+SiBoc4RSWr5Js3L9y7B0WKqOAunp4IQggzpfjyIUIIkRK2bVPBxcmTsZ+ztAbGgUDm8jG+eLGFFlTipFUCo4hy3bsHbm5w5IgERkIkBQmOhBAZjqbBgAFRy4C8zpK55kpyicNUpydL6MuPdGIVL7Bs2HFE74LQULU2WrNmsHMn5MxpUbZCiHikeLOaEEIkJU2DK1fU7epV9XfHDnXf2jqzgh/5iL/JT1WOcoa3rZKv0aj6FvXuDUOHQoECVslWCBEPCY6EEOlav34wf37cz+n1KniydCpcR14xi8/ozSKW8wEf8wMBuFiW6X88PGDRInjrLbCRT2whkoW81YQQ6dalS7BkibpfqBC4uMC5c+qxwWCduYvKcB5fPCnKdXqymKX0wFqj0Wxs4OjRqGY1IUTykLecECLd0TQ103W9eqqfTrNm0LZtVGAE1gmMurOUY1RBh0YVjrGUnlgrMCpfHubNk8BIiJQgNUdCiHRF08DLC9asUY9dXCBfPjVztLU485I5fEJ3fmIxPfmU2bzC2eJ8s2SBRo3g88+hVi3LyymEMI8ER0KIdCMgAAYNigqMRoxQHZm9va13jHKcwRdPCnCHrvzEz3S1OE97+zAOHtSoVElmcRQiNZDgSAiRLuzbB+3bRy0B4uAAU6da8wgavVnILD7jCiWozHEuUdoqOY8YcYzy5StbJS8hhOWkNVsIkaZpGmzcCE2bxlwbLSjozftmzpy4Y7jwghV0YQF9+YluVOOI1QKjSpWMVKz4wCp5CSGsQ4IjIUSaFBgICxZA2bKqs3XEEiDFi7+5E3OBAmpx1ufP33wcd05xkkq05jc6sop+/EgQjhaXP0KPHho66y21JoSwAgmOhBBpTmAglCqlFlu9eDEq3c5OTe5oNMa/b9eucOcOvHjxpqNofMxcDlGDF2SiEifxoaM1ih+paFH44IMECiuESBESHAkh0pSQEGjeXAU4QIxal5AQyJo1/n179gQ/vzcfw5Xn+OLJXD5hIb2pyUGuUsKicr9Or4dp09RoOiFE6iIdsoUQacqvv8Iff0Q9jpjdunBhFTA9fRr3fsOHw6ZNcP58wvl7cBwfvMjOY95nLet53yrlBlWzFRKimv527lQTU1qyjpsQImlIzZEQIk0wGmH0aOjRIyrNYFB/dToVGMU3saOHB6xf/6bASOMzZnKQmjwmO5U4adXAKE+eqJqtXbtUYCSESJ2k5kgIkSZMngyTJsVMi2hS07SEZ7w+cSLhvLPwlMV8yHts4HsGMYKphGJnWYFf888/6u/o0VCwoFWzFkJYmUXBUWhoKP/++y+vXr0iZ86cZMuWzVrlEkIIAMLC4JNP4l48NizM8vyrcZjVdCQzz3mXDfzKu5Zn+ppKldRUA+XKQUfr9ukWQiQBk4OjFy9e8PPPP7N69WqOHj1KSEgImqah0+nInz8/TZs2pW/fvlSpUiUpyiuEyGAOH44KjPT6hEeimUKHkcFMZwojOU5l6vEHt7F+W1fTprBqFchvRyHSDpP6HE2fPp3ChQuzZMkSGjduzIYNG/Dz8+Py5cscOnQIb29vwsLCaNq0Ke+88w5XrlxJqnILITKIbdvU35w5rRcYZeMxv9KG7xjG93xOXfZZJTCKGCkX0dw3e7YqvwRGQqQtJtUcHTt2jH379vHWW2/F+XzVqlX58MMP+eGHH1i6dCn79++nRAnrDn8VQmQcly9HrZMWffZrML8WqRYHWEUnHAmkJZvYQkvLCwqULKnKC6oPVK5c0oQmRFplUnC0atWqRG3n4OBAv379zCqQEEKA6sBcsyY8fhz386YGRjqMjGAqExnDIWrQiVXcJb/lBQU6dYK1a6MdS6eG6ufIYZXshRDJzOSh/A8evHkNoP3795tVGCGEePUKxoyBMmXiD4xMlZMHbKEFk/iSqYygAXusFhhVr65Gw0XMV5Qrl6rtKl/eKtkLIVKAyR2yy5Urx9y5c2nfvn2s5wIDAxkxYgTz5s0jJCTEKgUUQqRvoaGqb86uXXDpEly/HjWxozXU5Q9W0QkbwniHreygqdXytrFRHcYjdO4My5apdCFE2mVyzdGIESPo1q0bnTp14mm0qWj3799P+fLl2bp1K3v27LFqIYUQ6c+pU/Dll/DWWzBkCGzZAteuWS8w0hPOGCawm4ZcohQVOG1xYJQ1K9jaRj2OPpWAmxvMmyeBkRDpgcnB0ZAhQzh+/DhXr17lrbfeYu3atQwcOJCGDRvSokULTp8+Ta1atZKirEKINC4sTPXNqV1bzf0zeTJcuYLVV6V341+205RxjGMiY2jMTv4lj8X5Pn0a93Ifnp4q2MuUyeJDCCFSAbN+45QtW5bDhw/TpUsXvLy8cHJyYufOndSrV8/a5RNCpAMvX6paldmz4fbt2M9bsxmtETv5mQ/Q0NGYneyhofUyf02WLPDVV2qSSiFE+mFWcBQaGoq3tzfr16/Hy8uLrVu3MnnyZIoVK0b+/Nbp5CiESD969IB162Kn63RxB0bxpSfEQBjejOdLJrGLRnzAzzzAzazyvkmTJvDhh9C2LTg4JMkhhLAqo2YkICSAlyEvY9wCQmOnvQx5SVljAuvxZAAmB0d+fn507dqVgIAAtm3bRoMGDbh79y59+vShXLlyTJs2jV69eiVFWYUQadTr3RCzZ1cj0TQNChSA+/fVoqwRTA2M8nKXlXSmNgcYw0SmMBLNiutqRw/Wxo0Db2+rZS1Ekrny+AptVrfh9vPbvAp9ZdK+y/9b5MLFziUJSpb6mRwcVatWje7duzN9+nRcXNRJy5cvH1u2bGHhwoUMHjyYdevWsWXLFqsXVgiRthw/Du+9B0+eqMft2kGDBvDpp+pxs2YqcLJkcGsztrKcroRgRwP2sJ+6lhf8NRGBUZ06qhO5EGnBnpt7uPjoYow0vU6Ps60zLnYusW7Odv+l27qQ1XEF8JQiWYqkTOFTmMnB0YYNG2jevHmcz/Xu3ZsmTZrQu3dviwsmhEi7QkPBx0c1p4X/VztfqhQMHgwRFctOTlFLg5jDhlAmMJaRfM3vvEM3fuIROc3Pz0aVNa5aK3t7WLQIOnSQ0Wgi7WlarCk/v/czLnYuONg4oEvECIgjR7YSGPg0UdumRya/zeMLjCIUKlSIHTt2mF0gIUTatX49LFgA+/dDQEBU+vLlULgw1K0bNbP1K9Nq+WPIzx1W05FqHGE4U/mOoRY1o9nZxV97Vbo0LF4MNWqYnb0QKcrJ1omczub/cMiITAqObt++TcGCBRO9/d27d8mXL5/JhRJCpC1GI8yZA599Fvs5Nzfw8oL8+a2zcGxLNrGM7gTgTF32cYiaFuf5emBUuzZUrgwtW0KjRtafakAIkbqZ9FOrSpUqfPTRRxw7dizebZ4/f86CBQsoV64c6+IaniKESHe+/joqMKpSJSq9SRNYsUI1nyVi5aEE2RLCdwxhE635k1q442eVwCiCTqdqiS5fVjVf338PjRtLYCRERmRSzdH58+eZNGkSTZo0wcHBAQ8PD/LmzYuDgwNPnz7l/PnznDt3jkqVKvHNN9/QokWLpCq3ECIV+Ocf8PWN6qTs6goRv53eegu2b4e9e+Hddy07TiFuspqOVOIknzOdGQwCrBe1VKkCf/wBjo5Wy1IIkYaZFBxlz56d6dOnM2nSJDZv3syBAwe4desWgYGB5MiRgy5dutCsWTPKlSuXVOUVQqQSly+rRVejrSKEv78KMN55B4YNgx071FxAljSnteUXFvMhz8hCbQ5wjKoWl/11o0ZJYCSEiGJyh+zr169TpEgR2rdvH+fis0KI9O/331U/ohcvYqY3agSrV6tamK++UuulmcuOYL5lGJ8xm3W0oxeLeE4Wi8oNKggKDIx6/MUXltdsCSHSF5OHd5QoUYKHDx9GPvby8uL+/ftWLZQQIvXauRPatIkZGNWqpYa6r14NrVtD+/aWBUZFucaf1OIjfmQAs2nPWqsERjY2UYGRgwNMm6bWd5N+RUKI6EyuOdJemwRky5YtTJkyxWoFEkKkLq9ewYkT6nbsGGzYELUafbFiqumsyH/zxL3zDhw+bNnx2rOGhfTmITmpyUFO4mFZhtFElLtSJVW75ZIxJ/8VGUh2O3CzecrTp3sID38Z7Rbw2uOXGI1RaUFBcSyCmIHIdGZCiDjduAFLlrxF9+42PH8e8zlbWzXR4wcfRAVGc+ZYNqmjPUFMZzD9+QEfPOnDAl7gan6G8ejUCVautHq2QqQ6jmFX8KkOBt0fnD5tzgLMeuztC1i9XGmBycGRTqeLNWNmRp1BU4j06rPPYM4cG4zG4gDkywdvvw2nT8O9eyowypUL3n9fbb92LQwYYP7xSnAZXzwpzUU+Yh7z6Ys1R6NFyJcPxo+3erZCpEp2xn8w6CBUM5DZuRQGg8t/N+do913iTXd0LImDQ+LnNkxPzGpW69GjB/b29gAEBQXRr18/nJ2dY2y3fv1665RQCJHkNA1WrYJNm2D3brUQLOhwd3/A+PHZqFjRhtKlo2a1HjsWMmVSK9OfPq2CJXN1YiU/8hH3yEs1jvAXFSx6LdEXiY1uwwbVH0pvvfVohUgTbgXnoneDcyldjDTF5OCoe/fuMR5/8MEHViuMECJlTJkSc0FVvR769AmnefNDNG/eghkzogKjbdvUpI979lh2TEdeMZOB9GEhy/mAj/mBACzvBBRXYDR8uIxIE0IknsnB0ZIlS5KiHEKIFOLjExUYffghdO4MVauCg4ORLVtg7lw9X3yhnp82DX75xfLAqDQX8MWTYlzjQxaxhJ4kRTNavXqqlquhOd0thBAZlnTIFiIDu3ULunRR92vVgh9/jFp1PjQUtmwpwvz5BgCyZoXRo2POEWSObixjLv25RSGqcIzzvGVZhvFo0kRNJ2Ajn3JCCBNJ67sQGdQ//0CrVhAerhaF3bMnZiDx22865s9/O/Lx06eWBUZOBLCU7iyjBz54JVlgZGOj1kbbtk0CIyGEeeSjQ4gMqk8fOHtW3e/XTw3Pj27lyqjfTgaDur2+en1ileMMvnhSkNt05Sd+pquZpY5iMKgyBwXFTFu6FGrXtjh7IUQGJjVHQmRAp07B5s3q/oYNqs+RpsGVK7Bsmeqrs25dVB+g8HBzAyONXizkKFUJxRYPTlgcGNnZwZAhaobr6IFRpUpq5FxEM6EQQphLao6EyED+/VcFQ199pR7Xrg0tW0KPHmoY/+PH0be2rIO0Cy+YRz+6sJL59GEgMwnC8tVdf/wRFiyAgICotC5dYP58cHKyOHshhJDgSIiM4sgRqFs3qgbI1VUt/1GokJrYEdQQfqPR8mNVwA9fPMnDP3RiJavpZHmmwK+/qjIePKgeOziopsFixaySvRBCANKsJkSGEB4OfftGBUZffw01aqgmtIjAyNHRGoGRRj9+4DDVCcCZSpy0SmBkMEDTpmqR2FatotKXLpXASAhhfRIcCZHO/f03DBwIf/2lHh85Ap9/rhaSBZg4ERo1snyIvivP8cGLH+jPInpRg0NcpYRlmf4nPBy2b4+5qG316uDpaZXshRAiBmlWEyId0jT45hvVN+fataj0IUOgeHE12eOjR2Bvr4Ijc0ehRfDgOD54kYNHtGcN62hvUX5vat5zcVF9jGRZRyFEUpDgSIh06Kuv1MzQoAINDw/o3VutLVa+fFRTWnCwpUfS+JTZfMdQTlOBJuzgBkUtzTTewChbNpg9W72OTJksPowQQsRJgiMh0hFNg8GDYcYM9XjwYBUkZc4Me/dC2bLw7Jl1jpWFpyyiF+34he8ZxBd8TQj21sn8P0WLQpkyULIklCgB778PuXJZ9RBCCBGLBEdCpBOhoTB0KMyapR536KAWlP3tN9WR+eRJ6x2rKkfwwYvMPOddNvAr1l/VtWdPWLzY6tkKIcQbSYdsIdKJWbOiAqMBA8DXV/XLad/emoGRxmCmcYDa/EMe3PFLksCoZEnVX0oIIVJCqgiO5syZQ+HChXFwcKBatWocPXo0UfutXr0anU5H27Ztk7aAQqRiYWEqkJg4UT3u2FE1q/XpA599Zr3jZOMxv9KGaQxlBoOoyz5uU8h6B0BNSHnwIFy8qIbvCyFESkjxZjUfHx8GDx7MvHnzqFatGjNmzKBZs2ZcunSJXAl0Lrh58yZDhw6lTp06yVhaIVKX8+fVcPZz59TjSpXUDNKtWsHWrdY7Tk3+ZDUdcSSQlmxiCy2tlzlqjqUFC2TpDyFE6pDiNUfTp0+nT58+9OzZk7JlyzJv3jycnJxYnEBng/DwcLp06cL48eMpWtTykTFCpEU3b8J776nAKGtWmD4dxoyBihWtFxjpMDKCr/mDetyiEO74WT0wevttNXpOAiMhRGqRojVHISEhnDhxgpEjR0am6fV6GjduzKFDh+Ldb8KECeTKlYtevXqxf//+BI8RHBxMcLTxyv7+/gCEhoYSGhpq4SuIKSI/a+crYpLzDIcP63jvPQOPH+twc9PYuDGMSZMMDB78+u8dDXPXSMvBQ36iG83ZymRGMpYJhFv8kRGzPDY2GrNmhePsrJFR/51yPSefjHaujZr23z0tWV9zUp3n5HwNKRocPXr0iPDwcNzc3GKku7m5cfHixTj3OXDgAIsWLcLPzy9Rx5gyZQrjx4+Plb59+3ackmiVyh07diRJviKmjHqew8NhwIBGPH7sQpEizyhZ8ik1ahRG06IHQRFBiHmBUV3+YCWdsSWUZmxlO82sUfTI8uj1Rjw9L/HOO7d49iyYLVuslH0allGv55SQUc713eC75MkFwUHBbEmBN5m1z/OrV6+sml9CUrzPkSlevHhB165dWbBgATly5EjUPiNHjmTw4MGRj/39/SlQoABNmzbF1dXVquULDQ1lx44dNGnSBFtbW6vmLaJk9PP85586/vnHBldXDU3LzLZtWaI9a1lQpCecUUxmHOPYR126sIJ/yGthiWPXXm3YYOSdd4oDxS3MO+3L6Ndzcspo53rtMV8IAXsHe1o0apFsx02q8xzR8pMcUjQ4ypEjBwaDgfv378dIv3//Prlz5461/bVr17h58yatW7eOTDP+N5WujY0Nly5dothrq1Da29tjbx97YjpbW9ske3MkZd4iSkY9z76+6m9IiI6bN19/1vz1NNz4l5/5gIbsZiJjmMBYjFhjyFhUmXLmhJUroXHjNPW7LFlk1Os5JWSUc62PXF9HlyKv19rnOTlfQ4p+QtnZ2eHh4cGuXbsih+MbjUZ27drFgAEDYm1funRpzpw5EyNt9OjRvHjxgpkzZ1KgQIHkKLYQyU7T1MKxAwbAgQMqLSjIevk3ZBcrUD2iG7OTPTS0Xub/yZVLja7Lnt3qWQshhFWl+M+3wYMH0717dypXrkzVqlWZMWMGAQEB9OzZE4Bu3bqRL18+pkyZgoODA+XKlYuxf5YsWQBipQuRXjx+DM2awYkTUWkGg+p7ZCk94XgzntF8xS4a8QE/8wC3N+9ogrfegl69oEcPNapOCCFSuxQPjry8vHj48CFjx47l33//xd3dna1bt0Z20r59+zZ6fYrPOCBEipk4MWZglCULvHxpeb55uMdKOlOH/YxlAlMYaaVmNGXAAOjfX62NJoQQaUmKB0cAAwYMiLMZDWDv3r0J7rt06VLrF0iIVOLPP2H9+php1lg4tinb+JkPCMGOhuxmH/Usz/Q/VavCokUglblCiLRKqmSESIX8/KBJE6hdG+7csV6+BsKYzEi28Q7HqYw7flYLjAwG2LgRDh+WwEgIkbalipojIUSUZcvUivSR87dZSX7usIpOVOcwI/iabxmGZqXfRzqdmpW7cWOrZCeEEClKao6ESCVevoQfflAdl60dGLVgM364U4hb1OMPvmGE1QKjBg3g0CEJjIQQ6YcER0KksGPHoFMnNcS9f3/r5m1DKN8ylM204iA1ccePg9SyWv4lSsCOHVCtmtWyFEKIFCfNakKkkLAw6N1bNaO9zsZGPW+JQtxkNR3x4ASDmcb3fI4lk0S+7r33YNo01ddICCHSEwmOhEgB16+Dp2fMIfqZM8Pz5+q+pYHRu2xgCT15TmZqc4CjmFa1o9Ml3LTXoUPUTN1CCJHeSLOaEMnswgUoWzZmYFS6dFRgZAk7gvmeQWzgPfbQgIqcMjkwgrgDo8yZ1d86deKu7RJCiPRCgiMhktHly/D22xAcrB7Xrw+TJ8OlS5bnXZRr/EktPuYHPmUW77OOZ1hjSmoVKT1/Dvb2ag4jR0crZCuEEKmUNKsJkUw0DT78MKrJbMIEcHGBwYMtz7s9a1hIbx6Rg5oc5CQelmcaSUfRolCjBvTpozphCyFEeibBkRDJ5LPP1IzXAJUrw86dsG+fZXnaE8R0BtOfH/DBk77Mx5/Mlhf2P1myBHH6tIGCBdP/CuZCCBFBgiMhksG//8K8eVGPjx+3PM/iXMEXT8pwgX78wI98hDVHowFMnbqPPHkaWDVPIYRI7aTPkRBJ6O+/oV07yJPH8hFo0XVkFSephDMBVOcwP9IPawdGnToZcXMLtGqeQgiRFkhwJEQSCQ2FunXhl1+sl6cDgcynD6vozK+0wYMTnMbdegf4z6xZsGxZuNXzFUKItECa1YRIAnfvqg7XN26ox2+aNygxSnMBXzwpxjV6sZDFfIi1a4sitGqVJNkKIUSaIMGREFa2fz80awaB0VqkLA2MuvITP/AxtyhEVY5yDustex89cCtUCD7/HIoUUTVfQgiREUlwJISVDRgQMzCyhBMB/I8B9GQpS+jBAP7HK5ytk/l/IgKj6tVVYGcjnwpCiAxO+hwJYSW3b0Pt2vDXX9bJ7y3OcowqeOJLN5bxIUusHhgB5MqllgI5dEgCIyGEAKk5EsIqNA1at7ZWYKTxIYv5HwO4SnEqc5yLlLFGxpGcnSEgQAVzu3eDrUxjJIQQkSQ4EsIKPvnEOoGRCy/4gY/5gBXMpw8DmUkQpq/VodeD0RgzzcYG3npLre0WEKBqjFavlsBIiLRI0zSCw4N5GfIyzltASAA3nlzFzfqVzRmCBEdCmOHWLTh6FI4dg02bVMBhqbc5jS+e5OUenVnBKjqblY+9fdTabRH691druJUoASEhaoqBJUsgXz7Lyy2ESFrTD01nxZkVsYKfcC3h6TZa5YHaJcGglx40ppLgSAgTvHwJw4bFnO06IXZ2KhhJmMZH/MgMBnGR0nhwgiuUNLuM0QOjGjXgf/+DSpWgbVt4+FA1qfn6gpub2YcQQiSjifsm8izoWax0gw4c9JDV3oHsDo5ks3cki709me3syWRrQ2nn58AdSmYz//Mko5LgSIhE2rcP3n8fHj1Sj+NqunrdmwIjV54zn7544csc+jOEaQTjYJXybtgA776r7u/aBRs3qvsLFkhgJERa0jxXIHWyQ/Gs+bDThYMWBMZANC3il1DQf7ence6f1Une8KaS4EiIRHjxAnr2VIFR1qyq/87Dh5blWYkT+OBFTh7SAV/W0sE6hUX1LYoIjMaOhYkT1X13d+jY0WqHEUIkg84FgnGxAcLvEteUaTqdDQZDJgwGl2g3ZwwGF2xsspA//5DkLnKaJ8GREG+gafDpp3D9Ojg5QbVqsHWrRTkygP/xHUM5Q3masY3rFLNKWUuWhI8+UgHQzZswfjwsXaqe691bBUm6pJlUWwiRRCJ6DGUrMI3COevECH4MBhf0ersULV96JMGREAkIC4MRI2DZMvX41SvLAqMsPGURvWjHL8xgICOYSgj2FpfT3h6GD1fBkE4HmzerqQU0TT2eMAFGj7b4MEKIFGTv5I6ra5WULkaGIMGREHH4919VA7Nrlxr2bg1VOIoPXmThGW35hY20tUq+mTOraQQKFlRNfd26RQVwmTOrQKlWLascSgghMgQZ3ydENJoGly6pZqlff7VWYKTxOdP5k1rcx42KnLJKYGRrC336qJm5CxZU/aHatYsKjGrVUp3IJTASQgjTSM2REKhRZePGwYoVKtiwlmw8Zik9aM0mvmMIo5hMKJb1DyhaFPr1g169IFu2qPR+/eDAATV9wPr10LKlhYUXQogMSoIjIYBOnVRAAapGJiwsakFWc9XgIKvpiDMBtOI3NtPKovzKlFELw2bPHvu5zZth3Tp1/7ffoGlTiw4lhBAZmjSriQxL02DLFqhfPyowGj1azV9kSWCkw8hwprKPutymIO74WRwY2drCqVNxB0bnzsHQoep+p04SGAkhhKWk5khkKJqmalY2boTt2+Hvv1W6jY2aF+jrr1Wtkbly8JCf6EZztjKZkXgznjAsW7zM1lZNJfDjj3DvXszbP//As2dqu6xZwdvbokMJIYRAgiORwWzeHDU5IkCmTNC3rxqZVr68ZYFRHfaxik7YEUIztrKdZhaVVa+Hd96B8HCYPj3hbd9/H6ZMUWunCSGEsIwERyLDWLUKBg1S911dYc4cNWR/3z6oUCH2Yq2JpSeckUxhPN7spw6dWck/5LW4vEajavaL7rPPIG/emLd8+dTrEUIIYR0SHIl0LTwcTp5UTVKLFqm0EiVU81rXrpbnn4v7/MwHNGIXXzGaCYwl3My3VblycPasup83LxQooIboFyigbu3aqcdCCCGSlgRHIt169UrVCF29GpVWsyYcPGid/Buwm5V0BqAp29lFY7PysbNT8xXNmaMef/QRzJtnnTIKIYQwnYxWE+nW9OlRgVHTpuDsbJ3ASE844/BmJ405Sznc8TM7MALw8IgKjAYOhJkzLS+jEEII80nNkUiXLl1S64kBNGwIO3ZYPm8RQB7usYIu1GUf3oxnMqMwYjA7Pzc3OHRI3R88GKZNs7yMQgghLCPBkUh3nj5VI9BCQ8HBAXbvtk6+TdnGcroSii0N2c0+6r1xHzs7Nfv26xwdoXhxOHNGPf76a7XArRBCiJQnzWoiXTl8GIoVUyPQAIKCYm9jMLGix0AYkxjFNt7hJJVwxy9RgRHEHRh16gR+flFzLA0bBsOHm1YmIYQQSUdqjkSa9+QJrFkDv/yims+MxoS3Dw9PfN75+JtVdKIGh/iCKXzDcDQLflM4OangqH59VcNVpQpMnAg6ndlZCiGEsDIJjkSa5usLH3ygmtDeRKczrd9RCzazjO4E4kg9/uAgpi9vbzBEBWP29vD772oSymfPoFQpWLtWpQshhEg9pFlNpFlHjujo1UsFRpkzx79dgwaq71FiAyMbQvmGYWymFYeogTt+ZgVGhQrFrKXasgVOnFCBUaFCav4lmbdICCFSH6k5EmnWmDF6Xr6E0qXhypW4t8mVC/buTXxgVJBbrKYjlTnOEL5jOoMB89q8bt2Kul+xolreY+dO9fiTT1QTmxBCiNRHgiORJl24kI2DB1XQcvFi3NsYDPDgQeLzbMNGltKD52SmDvs5QnWzylawoFqWJHpn7FOnou737Quff25W1kIIIZKBNKuJNOfCBRg7tiYhIQnX6CS247UtIXzPIDbSlr3UpyKnzAqM8uZVgZq/f+xRanXrwqxZcO6cWsrERn6WCCFEqiUf0SLNWbhQT2io+RMvRleE6/jgRQVO8xkzmc2nmNOMZmcH778PZcvGHC3n4aGWBunbV0akCSFEWiHBkUj1nj+HXbvUHEbbt8Pp09YJjN5nLYvoxSNyUJODnKCyWfnodKpT+OzZMdOrVIGjR61QUCGEEMlKgiORqm3YAL17w+PH1svTniCmMYRPmIsvHejDAvxJYLjbG7ze2dvJSdUiDRtmYUGFEEKkCAmORKoVHg7du6s+PNHnC1I0zGn+Ks4VfPGkDBf4mLnMo59Z+bxOr1dl7dULqlWTPkVCCJGWyUe4SLW2bVOBkU4XMzBycNAICjI9oOnIKubTl3/IQ3UOcxp3i8oXMalklSrw559ga2tRdkIIIVIJGa0mUq0ff1R/X2+2MjUwciCQH+nLKjrzK23w4ITZgZGzM7RoAd7eqlx2dmrZEgmMhBAi/ZDgSKQ627ZBpUrw66+W51WKixyhGl1ZTm8W8AE/85JMZuVlMMDx46pc//yj0mrVgnz5LC+nEEKI1EOa1USqcPEi7N+vFpDdscM6eXblJ37gY25TkKoc5Szlzcqndm1o106tiVaoEDRpAnv2qOeGDrVOWYUQQqQeEhyJFKVpalmNL7+0Xp5OBPA/BtCTpSylO58wh1c4m5SHoyMEBqoRZ998E5V+5IgKjBwc1KSOLVpYr9xCCCFSBwmORIq5fBlatoSrV62XZ1nO4YsnhblJd5byE91N2t/JCcqVi5qfyMsr6rngYJg/X92vXFlN7iiEECL9kT5HIlm9fAnjx6sApHRpawZGGj1ZzDGqYERPFY6ZHBgBvHoVFRj166dmuI7Qpw8sXhz1nBBCiPQpVQRHc+bMoXDhwjg4OFCtWjWOJjCt8IIFC6hTpw5Zs2Yla9asNG7cOMHtRcp7+RL27lXNUG+9BePGqTXGIkahxbWshilLbbjwguV0ZTG9WEEXqnGEC5RN9P45c0LDhjB4MEyfDqtWwbFjMHeuev7uXfjiC1i+XD1euBC6dEl8+YQQQqQtKd6s5uPjw+DBg5k3bx7VqlVjxowZNGvWjEuXLpErV65Y2+/du5dOnTpRs2ZNHBwcmDp1Kk2bNuXcuXPkk2FDqYqmwcaNcc9wbTCojs5//glhYVHpjo5qv6CgxB3jbU7jiyd5uUdnVrCKziaV8YsvVJ+nuJw/rwK59euj5lnq109N9CiEECL9SvGao+nTp9OnTx969uxJ2bJlmTdvHk5OTiyOaL94zYoVK+jfvz/u7u6ULl2ahQsXYjQa2bVrVzKXXCQkNBTee0/dogdGer3qzBweDn/8ERUY6XSqBicwMLGBkUZffuQI1QjEEQ9OmBwYeXnFHxj9/TfUqaNGz4WHq/s+PjBnjkmHEEIIkQalaM1RSEgIJ06cYOTIkZFper2exo0bc+jQoUTl8erVK0JDQ8mWLVtSFVOYKCwMPvtM1RpF5+4Ofn6xgx9XVzUT9sOHics/E/4soA9e+DKXjxnMdIJxMKmMxYur5rO4aBqsWAFPnkDBgvDbb/D22yZlL4QQiWI0hhIe/jLaLeC1xy+xSfFqjIwnRYOjR48eER4ejpubW4x0Nzc3Ll68mKg8RowYQd68eWncuHGczwcHBxMcHBz52N/fH4DQ0FBCQ0PNLHncIvKzdr5pyatX0LKlgT//VO9mg0Gjc2eNzJk1/vc/AwBt2oTz6696QIeTk8aLF5DY9c0qchJfPMnFAzzxYQ2eJpexcmUjBw+Gx2jOAzhzBtau1bNunZ7Ll1V5WrQIp0wZIxnxXyrXc/KQ85x8UvJcBwZe5tq1zwgJ+QejMSoA0rSQN+5r919wFGY0pInrJKnOc3K+9hTvc2SJr7/+mtWrV7N3714cHOKuOZgyZQrjx4+Plb59+3acnJySpFw7rDWLYRrz8KEjM2ZU4ty5HEQsDNuly3maN79Jr15NAQNNm97g118LExEMhYaGoWm2vHkhWY1PmMM0hnCG8jRjG9cpZnIZHR1DGTHid7ZsibkmybZthfjhB/fIxzY24VSq9IAKFc6xZUuAycdJTzLq9Zzc5Dwnn5Q41/b2a3Bw2B3v85pmg6Y5AvZomsN/9x3QNAf+fHqSs/7hNH5wgav2z5KryBaz9nl+9eqVVfNLSIoGRzly5MBgMHD//v0Y6ffv3yd37twJ7vvdd9/x9ddfs3PnTt5OoM1j5MiRDB48OPKxv78/BQoUoGnTpri6ulr2Al4TGhrKjh07aNKkCbYZbLEtoxFq1jRw7pwevV7DaNSh12vkzFmadevK8OqVnmzZNP78szDRg6DQ0IjzFH9glJlnLKIX77OeWXzKML4lBHuzytm7t553320eI+3pU5gwQdVq1aljpFcvI61aabi65gDqmXWc9CAjX8/JSc5z8knJc3379knu3IFs2VqRP/8oDAZnDAYXDAYX9Hpn9Hq7ePd9/7vsvAh5wVdt6lMsq+k/CpNbUp3niJaf5JCiwZGdnR0eHh7s2rWLtm3bAkR2rh4wYEC8+33zzTdMmjSJbdu2Ubly5QSPYW9vj7197C9SW1vbJHtzJGXeqdWqVXDyJNjYQFiYCnSMRh3Tphkit/H318VqynqTKhzFBy+y8Iz3WM8G3jO5bJkyweefQ7duUKyYAVBl0jSYNw9Gj1b9i2xtYfFiPcWLSwN/dBnxek4Jcp6TT0qca4NBfe7Y2+cjW7YaZuVha5O2rhFrn+fkfO0p3qw2ePBgunfvTuXKlalatSozZswgICCAnj17AtCtWzfy5cvHlP+GFU2dOpWxY8eycuVKChcuzL///guAi4sLLi4uKfY6MrJ79+Cjj9T96MFP/vxQowaUKgWLFkUt1po4GoOYwVRGcIqKNGAPtyhsctnatVMjzOKqiDx8GPr3V/ffekvNw1S8uMmHEEIIkc6keHDk5eXFw4cPGTt2LP/++y/u7u5s3bo1spP27du30eujfsn/8MMPhISE0L59+xj5eHt7M27cuOQsukCNPOvenf86VauZr8+eVbU1vXvD9evw44+JH4kGkJUnLKUHbfiN7xjCKCYTSvxVzvHp0gV+/jl2+o0bsGCBKheogMjPT9V6CSGEEKni62DAgAHxNqPt3bs3xuObN28mfYFEoty6BVWqRAU++fJB4cIqOHrxQk2gaKoaHGQ1HXEmgNb8yiZam1W2qVNh+PCYaf7+4O2tapIiBj24uqoASgIjIYQQEeQrQZglPBzat48KjPLkgVy5YNMm8/LTYWQo3zGZURyhGh1Zzd8UMDmfnDnVkPzXZofgt9/U7Nb37qnHjRqptdJatgRpjRVCCBGd9DwVZunSBY4fV/fbtIGBA+HUKfPyysFDNtGKbxjBtwyjPntNDoz0etV/6PbtmIHRq1cwc6Yq4717qglt61bYuVPNkC2BkRBCiNdJzZEwyfnz8O23aikNUDNHr18PHTqYl18d9rGKTtgRwjv8zjbeMTkPe3vVxBc9KDp1CmbPVst/vHyp0qpVgz171PptQgiRWmiaRnB4MC9DXkbeAkICIu+HGlP/xI/pjQRHItHOnIGaNaOCDYBt22DHDvjlF9Py0mFkJFOYwFgOUJvOrOQepi4crFGzpo7ly2MGRqGhUL++6mMEUKSIGsY/aJAERkKIlPEs6BkAv17+ldV/HIsVCIVr4W/Mw9aQdobxp3USHIk3evpUjUj77beY6cOHQ3AwtGhhWn65uM9yutKYnUziS8bjTbgZl+KXXxr56itDjDR/f9XB2t9fjZjbsgVq1VIL2wohREo5c/8vCgD/vPiHkwnMa+Jo44iLnQvOds642LlE3qrlq0bBzAWTr8AZnARHIkGaBj16xAyM3n9f9TGqWROyZlXbJFYDdrOCLujQaMp2dhH3mnhv8t57V/D2LgwY8POD6dPh0CG4ejVqm9atoXZts7IXQgirCteMAOR3zceWzgtiBD4RgZCzrTMGveENOYnkIMGRSND8+fDrr1GPGzaEsWPVrNJbtqhRa4mhJ5wxTGQsE9hDA7qwgvskvERMXFxcVLNemTJP0LTCjBwJX38dc5sCBVSz2owZJmcvhBBJytXeleYlmr95Q5GiJDgS8QoNVSO9ItSrp0aDVahgWj65+YeVdKYu+xjHOCbxJUZM/3WUIwc8ehRxP5AlS3SRgVHHjtCzJ1SqpLYTQgghzCXBkYjl6FE1yuv33+HCBZWm08Eff5ieVxO28zMfEIYNjdjFH9Q3aX8HBzULN0RN1FimjMayZWU5fVoljB0L48ebXjYhhBAiLhIciRi2bYPmzWP3IzKlXxGAgTDG481IprCDJnRlOQ/JZVIeNjZRgRHAf8voceGCDv7Lq3dv1cQnhBBCWItMAiki/f23mjVa01S/HXPl42/20IARTGUUk2nO7yYHRhBzEVuA/xa1xt5eo3PnC5w5E8qCBZCGFqkWQgiRBkjNkQBUjZGnpxoC7+ysmrPM0Zwt/EQ3AnGkPnv5E8uGi+l0UbVW4eFqhutFi8J5+vQypUoVtyhvIYQQIi5ScyTQNBgwQAVGxYtD1apw5YppedgQylSGs4WWHKY6FTllcWAUUTZbWyhbFkaOhNOnoUYNE9v4hBBCCBNIzZFg1iw1P5BOBx9/DEOGmLZ/QW6xmo5U5jhD+ZbpDEazQtw9dCj07atmuLaJdqWGykz6QgghkpAERxnc5ctqWQ1QtTSmBkZt2MgSeuKPK3XYzxGqW1Se0qVhwgQ1eWOePBZlJYQQQphFmtUyKE1TC8iWKmXe/raEMJ3P2Uhb9lGXipwyOzAqVEjNwv3jj/DXX2oRWwmMhBBCpBSpOcqA/P3h889h8WLz9i/MDXzwwh0/BjKDWXwGmL54mV6vFrMtW9a8cgghRHILCbnPixfHCQ9/+dotII60qFt+bqd00YUJJDjKQMLC1FIgQ4fCjRvm5dGOdSyiF4/JTk0OcoLKZpenXz8JjIQQaYemhXPsWAVCQ++bvG/El+3TcJnCPy2Q4CiD+OcfqFgR7pv+ngbAniC+YygDmMMa2tObhfiT2ay8bG1h1Ch1E0KItMJoDIkMjDJnro2NTRYMBhf0emcMBpcEb4v8VvDN4Xm8V05+EaYFEhylc3fvwqefwi+/mJ9HMa7iiydlOc/HzGUe/TCnGQ2gSxe11EexYuaX5//t3Xl4U1X6wPFvkrbpQkthgLZAHdldQBCQgsgPxWIdEEFH2oGCgLLIMiB1YdXWYZXFYRFFcYDKIAUXGJV90A677IqyFspOK4tAF9qkzfn9EQkN3ZLaJE14P8+TR+7Nufe8eQPmfc659x4hhHC1Jk3W4uVVyeb2WZrv+DXXgQGJciXFkQdLSzOvTp+SUvZzxJDExwwkjVBas5MfaVbmc50/DzVrlj0WIYQQwhnkbjUPlZoKTzxhXRg9/7ztx/tyk/kMIokefMsztGCv3YVRQID5v1qt+VonKYyEEEK4Axk58kCrVpmXArn1sEQvL/Mt+199ZdvxjTjCCqJpwHEG8DGf0J+yTKNlZYG/P0ycCF262H24EEII4RJSHHmgqVOtnyJduzb88ottx/ZiCR8ymLOE04pd/EyTMsWg0cCyZfDss+DnV6ZTCCFEhaaU4mbeTTINmVavLENWoX3Jp5JdHa6wgxRHHubgQfjhh9vbej2cOlX6cf5kMZe/8xKLSORFhjKPLGy/2PBO48dDTEyZDxdCiAotbGYoV3KyUdi31mOgT6CDIhLlSYojD6IUjBp1e9vLC3JtuDviAX5hBdHcyyn6sohE+v6hOFasMD/lWgghPEluXo7lz5mGLKuyKMA7gEo+lajkU4kAn9t/tuzzDqCafzUGtxzs/MCF3aQ48iCJibB27e3tvLzSjlD0ZTHzGMpJ6vIIuzlM2Z/B4eMD48ZJYSSE8ExK3S6HdvbfQWhgHQJ8AvD39kerkfubPIkURx7i2jVzYWKrADL5kMH05t98wssMZw438S9T3506mZ+63aYN+PqW6RRCCOFW7g2+l8p+Ia4OQziIFEceol8/uHDBtrZN+IkVRFObc8Tybz4jtkx9ajTw5pvmC8CFEEIITyHjgG7u/Hl44QXz7fulUwzkI3bRihx8ac6+MhVGfn7Qp495tEoKIyGEEJ5GiiM3pRQkJMA998CXX5bePpAbLKMHH/EKi+hHG3ZwnIZ29ztnjvn5RYsXQ1CQ3YcLIYTbUcpEfn4WRuMlV4cinESm1dzQypUwciScPm1b+4fZx3JiCCGdaJbzOdF29afRmIuxmBjzOm1CCOHOjMYrXLy4BKPxMvn5mXe8sgrtM5myXB2ycDIpjtzMO++YR4xsoxjKPGbyGj/TmL+wlhPUt7mvBx6AEyfMjwOIiIBFi8oSsRBCVCwXLvyTc+emlenY7DzYdRVa6sr+HDhR8Ulx5Cby8yE+HiZNsq19Za7xCf15gS+Zw995g+kY0Nt0rL+/eV20Q4fM2x06mJ92LU+6FkJ4gry8GwAEBT1KlSqR6HSV0OkCyEeP0aTFoHTkmLTczIfsfEWWMZ+svHyu5mTx93UjAHi9mws/gHA4KY7cxCuvwCef2Na2JbtZTgxVucrzfMlK7FhxFsjONr80Gujb19yvVq5OE0J4iF+zfgVgacoRks6dsizxYVImm4730nrhpZWfT08m364bmDTJ1sJIMYLZTONNDtCMJ9nEKerY1dfTT5vXQ3v4YWjSxDyCJIQQnuTY1WPU94Yr2Ve5kFH4fT8vv2KfdF3JpxKRdSLx9ZKHunkyKY4quJQU83RaaapwlUX0oytfM5M4xjAFIz429xMSAh98AM89Zx4xEkIIT3XrSdcPhzZjeNS/Ci3zodPqXByhcDUpjiqon36CF1+EH38svW1rdpDE3wgkgy58zbd0sauvdu3gm2+gcuUyBiuEEG6osm9lmoc1d3UYogKSK0kqoLw881IcpRVGGky8znQ283+cozbNOGB3YfTEE7BpkxRGQgghxC1SHFVArVubL4guyZ+4zDd0YTpvMpPXeJxkznKPTeevWhWSkszLjXz3HXh7l0PQQgghhIeQabUKZuFC2Lu35DaPsYVl9EBPLn9hDev4i83nr10b9uwxX2MkhBBCiMJk5KgCycyEgQOLf1+DibFMIpnHOUldmnHArsJo6FA4c0YKIyGEEKIkMnJUQUybBqNGFf9+DdJZQm8i+S+TGMc7xJNvx9fXrBnMni13ogkhPIfJlEde3pUSl/4ouN9ovIGf3zHCvM64OnRRwUlx5GK//QbDh8O//118m8f5ns/oiRYTUaznv3S0q4969czPSdLJ3alCCA9hMhnYvbsxN28et+s4Hx8sDzkxIEuAiKJJceQio0eb1yr79dfi22jJZzwTeZt/8D/aE8tS0giz6fyBgRAbCx07mh/q6CXftBDCg+TmnrcURjkmLbn5GnJM5rXPsvJMZOcrbuZDTj7cLOJ11QhdHopw8acQFZX8ZLrAlCnw7rsltwnlIkuJpT3/4x3imcQ4TJQ+9KPRwGOPyXOLhBCeLS3zImAudDptLX7ZDx+dj9UDHvOy86hdvTZ1qocR27SPs8IVbkaKIyfbvh3Gji25TSQb+Te9MKHlSTbxPx636dz33w8HD8r0mRDC8xnz8wDQAN/2+LbI5T4CvAPw1t1+VonRaGTNmjV06tQJb3mGiSiBFEdOdOWKeZqrODrySCCBsUxmIx3pzRIuUaPU81avDu+/Dy+8IAvECiHuPp0bdnZ1CMLDyE+pkyQlwZ//XPzDHWtxju/owGimMo5J/IW1pRZGderArFlw+jRER0thJIQQQpQHGTlygqQk6NGj+PefZi1L6E0OvjxOMtt4rNRzdu5svq5Ibs0XQrizLEMWa46v4bec38gyZJFpyLz9MmaSZcgk13gdY94N8vIyyM/PxGTKprIuk7gGro5eeCopjhxs8+biCyMvjExkPKOYxmo60YdErlCt1HM2bw7fflvOgQohhAtM2fQyuozlVPICPx3U1IGvDvy8wc/XvK8kGo38jInyJ3+rHOytt4reH84Zkvgbj7Cb15nOe8ShbJjljIgwT6UJIYQnCNdsp1Hpl1ai0IDGD43WD63WH52uEl5elagf+qLjgxR3HSmOHOjYMfPI0Z268DWL6UsGgfwfm9lJG5vO9/zz8OWX5RykEEK4kA7zbfjXvFrTttGb6HSV7ngFoNNVQqv1QyPXEQgnkeLIQfbsKTyd5o2BqYwmjn/yH56lH4v4jao2nW/OHBgyxAGBCiFEOcjMPEhW1i9Wy3aYTEUt52G9r44+DYCb2nupXv05F38KIcykOHKAr77S8Le/We+7l1SWE0MzDjCCWcxhOOYndJTMyws+/bTkC7qFEMKVcnJOs2fPQ2U6VqcBgwlytbXLOSohyk6KIwdISLC+gvA5vmIhL3GVqrRlG3t4xKbzhIbC99/Dffc5IkohhCgfubnmp1Vrtb5UqdKx0JTYrZfS+GIw6TAoLTn5GnJMGt7dPoeNp3YzPaquiz+FELdJcVTOfvyxCkeOmEeEfMhlBq/zd97nC/5Kfz7hOsGlnqNOHUhMhHbtHBysEEIUQan836e/ilvl3nqK7ObNFAB+M+qY/rMi03DR6pb8W7foG03GYvvUaeXR/qLiqBDF0bx585g+fTppaWk0bdqUuXPn0qpVq2Lbf/7557z11lucOnWKBg0a8O6779KpUycnRly8+HjzM4rqkcJyYmjMzwxhHh8yGFum0UaPNq+9JoQQznbu3BxSU8eRn59ZpuPTsrP49ljpzxnR6/RWy3zUDqpN5wbylGtRcbi8OFq+fDlxcXHMnz+fiIgIZs2aRVRUFEePHqVGjcL3d27fvp0ePXowZcoUnnnmGT777DO6devGvn37aNy4sQs+wW39+2sADdEsZwEDSCeE1uzkAA+XeFxwMHToAM88A337OiNSIYQo7NdfV9xRGGnRaP1B6w8aX0z4kIcPecobo9JhUDpyTBou3cxkb9ohfs6ozCddZloVPneueXbnemdCVEQuL47ee+89BgwYQL9+/QCYP38+q1evZuHChYwePbpQ+9mzZ/P000/zxhtvADBhwgQ2btzI+++/z/z5850ae0EZGbDiUwMfEscrfMQy/sYgPiKDoCLb+/lBs2bmp2ffc49zYxVCiKJczLxIADD1CHx/CQwmE5D5+6t0zcPq8XLzlx0ZohBO4dLiyGAwsHfvXsaMGWPZp9VqiYyMZMeOHUUes2PHDuLi4qz2RUVFsWrVqiLb5+bmkpuba9m+ceMGYF6d2Wgsfv7bXi2DjrOTHjTkGAP4mE/oT1HTaNWrK0aNMjFokAm9nt9jKbcw7gq3vrfy/P5EYZJn56hIeb6UdYkAP8jON99BBhDgHWAZ8bGMAHnfHg2y7PeuRLdG3SrE5yhORcq1J3NUnp35vbm0OLp8+TL5+fmEhIRY7Q8JCeHIkSNFHpOWllZk+7S0tCLbT5kyhXfeeafQ/g0bNuDv71/GyK3VTk5mLws4R20i2MlBbt3SqgDz+mdNm6bTq9dh6tW7gUYDmzaVS9d3tY0bN7o6hLuC5Nk5KkKer5v8OW3K4umq3RgS1h29Vo9WU8qT+01Ajvl1Zs8ZznDGGaH+IRUh13eD8s5zdnErtzuAy6fVHG3MmDFWI003btwgPDycp556iqCgoqe8bJadje7VV9EuXowpNpbEKu9x5CMN7068SYsW3kREwOXLEBQElSr9CY2m9AVlRemMRiMbN26kY8eOeHvLtQuOInl2joqV57Mu7t+xKlauPZej8nxr5scZXFocVatWDZ1OR3p6utX+9PR0QkNDizwmNDTUrvZ6vR79rfmrAry9vf/Yl3boEHTvDqmpsHAh2r59+UdeHq0j19CpUyfLuStVKnsXomR/+DsUNpE8O4fk2Xkk185R3nl25ndW+kqnDuTj40OLFi3YVGCOyWQysWnTJtq0KXq9sTZt2li1B/PQXXHty51SsGgRtGxp3t6zB/r1M8+dCSGEEMLtubQ4AoiLi2PBggUkJiZy+PBhBg8eTFZWluXutRdffNHqgu0RI0awbt06Zs6cyZEjR0hISGDPnj0MGzbM8cFmZkKfPvDSS+b1PHbvhgcecHy/QgghhHAal19zFBMTw6VLl3j77bdJS0ujWbNmrFu3znLR9ZkzZ9Bqb9dwjz76KJ999hnjx49n7NixNGjQgFWrVjn+GUc//QQxMXD2LCxZAr16ObY/IYQQQriEy4sjgGHDhhU78pOcnFxoX/fu3enevbuDo/qdUrBgAYwYAQ0bmqfRZLEzIYQQwmO5fFqtQrtxA3r2hEGDzNNpO3dKYSSEEEJ4uAoxclQh7d8P0dGQnm5+jHVMjKsjEkIIIYQTyMjRnZSCDz6A1q3NDyjat08KIyGEEOIuIsVRQdeumUeLhg6FgQNh+3aoX9/VUQkhhBDCiWRa7Zbdu80jRFevwhdfwF//6uqIhBBCCOECMnKkFMyaBW3bQrVq5muNpDASQggh7lp3d3F09So89xyMHAnDhsHWrVCnjqujEkIIIYQL3b3Tart3m590nZEB//kPPPusqyMSQgghRAVw944cRUVBrVrmaTQpjIQQQgjxu7tu5EgpBcCNAQNg4kTw9jY/7LEcGI1GsrOzuXHjhqz47ECSZ+eQPDuH5Nl5JNfO4ag83/j9t/rW77gjaZQzeqlAzp07R3h4uKvDEEIIIUQZnD17ltq1azu0j7uuODKZTFy4cIHAwEA0Gk25nvvGjRuEh4dz9uxZgoKCyvXc4jbJs3NInp1D8uw8kmvncFSelVJkZGRQs2ZNqwXpHeGum1bTarUOrziDgoLkH54TSJ6dQ/LsHJJn55FcO4cj8ly5cuVyPV9x7t4LsoUQQgghiiDFkRBCCCFEAVIclSO9Xk98fDx6vd7VoXg0ybNzSJ6dQ/LsPJJr5/CEPN91F2QLIYQQQpRERo6EEEIIIQqQ4kgIIYQQogApjoQQQgghCpDiSAghhBCiACmO7DRv3jzuvfdefH19iYiIYNeuXSW2//zzz7nvvvvw9fWlSZMmrFmzxkmRujd78rxgwQLatWtHlSpVqFKlCpGRkaV+L8LM3r/PtyQlJaHRaOjWrZtjA/QQ9ub52rVrDB06lLCwMPR6PQ0bNpT/d9jA3jzPmjWLRo0a4efnR3h4OCNHjiQnJ8dJ0bqnzZs306VLF2rWrIlGo2HVqlWlHpOcnEzz5s3R6/XUr1+fxYsXOzzOP0wJmyUlJSkfHx+1cOFC9csvv6gBAwao4OBglZ6eXmT7bdu2KZ1Op6ZNm6YOHTqkxo8fr7y9vdXBgwedHLl7sTfPPXv2VPPmzVP79+9Xhw8fVn379lWVK1dW586dc3Lk7sXePN+SmpqqatWqpdq1a6e6du3qnGDdmL15zs3NVS1btlSdOnVSW7duVampqSo5OVkdOHDAyZG7F3vzvHTpUqXX69XSpUtVamqqWr9+vQoLC1MjR450cuTuZc2aNWrcuHHqq6++UoBauXJlie1Pnjyp/P39VVxcnDp06JCaO3eu0ul0at26dc4JuIykOLJDq1at1NChQy3b+fn5qmbNmmrKlClFto+OjladO3e22hcREaEGDRrk0Djdnb15vlNeXp4KDAxUiYmJjgrRI5Qlz3l5eerRRx9Vn3zyierTp48URzawN88ffvihqlu3rjIYDM4K0SPYm+ehQ4eqDh06WO2Li4tTbdu2dWicnsSW4ujNN99UDz74oNW+mJgYFRUV5cDI/jiZVrORwWBg7969REZGWvZptVoiIyPZsWNHkcfs2LHDqj1AVFRUse1F2fJ8p+zsbIxGI1WrVnVUmG6vrHn+xz/+QY0aNXj55ZedEabbK0uev/76a9q0acPQoUMJCQmhcePGTJ48mfz8fGeF7XbKkudHH32UvXv3WqbeTp48yZo1a+jUqZNTYr5buOvv4F238GxZXb58mfz8fEJCQqz2h4SEcOTIkSKPSUtLK7J9Wlqaw+J0d2XJ851GjRpFzZo1C/2DFLeVJc9bt27lX//6FwcOHHBChJ6hLHk+efIk3333HbGxsaxZs4aUlBSGDBmC0WgkPj7eGWG7nbLkuWfPnly+fJnHHnsMpRR5eXm88sorjB071hkh3zWK+x28ceMGN2/exM/Pz0WRlUxGjoRHmTp1KklJSaxcuRJfX19Xh+MxMjIy6N27NwsWLKBatWquDsejmUwmatSowccff0yLFi2IiYlh3LhxzJ8/39WheZTk5GQmT57MBx98wL59+/jqq69YvXo1EyZMcHVoogKQkSMbVatWDZ1OR3p6utX+9PR0QkNDizwmNDTUrvaibHm+ZcaMGUydOpX//ve/PPTQQ44M0+3Zm+cTJ05w6tQpunTpYtlnMpkA8PLy4ujRo9SrV8+xQbuhsvx9DgsLw9vbG51OZ9l3//33k5aWhsFgwMfHx6Exu6Oy5Pmtt96id+/e9O/fH4AmTZqQlZXFwIEDGTduHFqtjB2Uh+J+B4OCgirsqBHIyJHNfHx8aNGiBZs2bbLsM5lMbNq0iTZt2hR5TJs2bazaA2zcuLHY9qJseQaYNm0aEyZMYN26dbRs2dIZobo1e/N83333cfDgQQ4cOGB5PfvsszzxxBMcOHCA8PBwZ4bvNsry97lt27akpKRYik+AY8eOERYWJoVRMcqS5+zs7EIF0K2CVMmSo+XGbX8HXX1FuDtJSkpSer1eLV68WB06dEgNHDhQBQcHq7S0NKWUUr1791ajR4+2tN+2bZvy8vJSM2bMUIcPH1bx8fFyK78N7M3z1KlTlY+Pj/riiy/UxYsXLa+MjAxXfQS3YG+e7yR3q9nG3jyfOXNGBQYGqmHDhqmjR4+qb7/9VtWoUUNNnDjRVR/BLdib5/j4eBUYGKiWLVumTp48qTZs2KDq1aunoqOjXfUR3EJGRobav3+/2r9/vwLUe++9p/bv369Onz6tlFJq9OjRqnfv3pb2t27lf+ONN9Thw4fVvHnz5FZ+TzR37lx1zz33KB8fH9WqVSu1c+dOy3vt27dXffr0sWq/YsUK1bBhQ+Xj46MefPBBtXr1aidH7J7syfOf//xnBRR6xcfHOz9wN2Pv3+eCpDiynb153r59u4qIiFB6vV7VrVtXTZo0SeXl5Tk5avdjT56NRqNKSEhQ9erVU76+vio8PFwNGTJE/fbbb84P3I18//33Rf7/9lZu+/Tpo9q3b1/omGbNmikfHx9Vt25dtWjRIqfHbS+NUjJ+KIQQQghxi1xzJIQQQghRgBRHQgghhBAFSHEkhBBCCFGAFEdCCCGEEAVIcSSEEEIIUYAUR0IIIYQQBUhxJIQQQghRgBRHQgiXePzxx3n11VddHYZD9O3bl27durk6DCFEGUlxJIQQbkSj0bBq1apC+4sqyObNm8e9996Lr68vERER7Nq1yzlBCuHmpDgSQogKRilFXl7eHzrH8uXLiYuLIz4+nn379tG0aVOioqL49ddfyylKITyXFEdCuJFLly4RGhrK5MmTLfu2b9+Oj49PoZWvS5KQkECzZs346KOPCA8Px9/fn+joaK5fv17qsRs2bMDX15dr165Z7R8xYgQdOnQA4MqVK/To0YNatWrh7+9PkyZNWLZsWYnnLWpEJDg4mMWLF1u2z549S3R0NMHBwVStWpWuXbty6tQpy/vJycm0atWKgIAAgoODadu2LadPny6yv+TkZDQajdXnOHDgABqNxnLOxYsXExwczPr167n//vupVKkSTz/9NBcvXrQck5+fT1xcHMHBwfzpT3/izTffLLSqu8lkYsqUKdSpUwc/Pz+aNm3KF198USiWtWvX0qJFC/R6PVu3bi0xX6V57733GDBgAP369eOBBx5g/vz5+Pv7s3Dhwj90XiHuBlIcCeFGqlevzsKFC0lISGDPnj1kZGTQu3dvhg0bxpNPPgnAqVOn0Gg0JCcnl3iulJQUVqxYwTfffMO6devYv38/Q4YMKTWGJ598kuDgYL788kvLvvz8fJYvX05sbCwAOTk5tGjRgtWrV/Pzzz8zcOBAevfu/YemdYxGI1FRUQQGBrJlyxa2bdtmKVYMBgN5eXl069aN9u3b89NPP7Fjxw4GDhyIRqMpc58A2dnZzJgxgyVLlrB582bOnDnD66+/bnl/5syZLF68mIULF7J161auXr3KypUrrc4xZcoUPv30U+bPn88vv/zCyJEj6dWrF//73/+s2o0ePZqpU6dy+PBhHnrooTLHbDAY2Lt3L5GRkZZ9Wq2WyMhIduzYUebzCnHXcO26t0KIshgyZIhq2LCh6tmzp2rSpInKycmxvHfu3DnVqFEj9cMPPxR7fHx8vNLpdOrcuXOWfWvXrlVarVZdvHix1P5HjBihOnToYNlev3690uv1Ja5o3rlzZ/Xaa69Zttu3b69GjBhh2QbUypUrrY6pXLmyZQXvJUuWqEaNGimTyWR5Pzc3V/n5+an169erK1euKEAlJyeXGr9St1cXLxjz/v37FaBSU1OVUkotWrRIASolJcXSZt68eSokJMSyHRYWpqZNm2bZNhqNqnbt2qpr165KKaVycnKUv7+/2r59u1X/L7/8surRo4dVLKtWrSo1bkD5+vqqgIAAq5eXl5elz/PnzyugUJ9vvPGGatWqVal9CHG383JdWSaEKKsZM2bQuHFjPv/8c/bu3Yter7e8V6tWLY4cOVLqOe655x5q1apl2W7Tpg0mk4mjR48SGhpa4rGxsbG0bt2aCxcuULNmTZYuXUrnzp0JDg4GzCNJkydPZsWKFZw/fx6DwUBubi7+/v5l+8DAjz/+SEpKCoGBgVb7c3JyOHHiBE899RR9+/YlKiqKjh07EhkZSXR0NGFhYWXuE8Df35969epZtsPCwizX7Vy/fp2LFy8SERFhed/Ly4uWLVtaptZSUlLIzs6mY8eOVuc1GAw8/PDDVvtatmxpU0z//Oc/rUaFAEaNGkV+fr7tH0wIUSwpjoRwQydOnODChQuYTCZOnTpFkyZNnNr/I488Qr169UhKSmLw4MGsXLnS6tqg6dOnM3v2bGbNmkWTJk0ICAjg1VdfxWAwFHtOjUZT6Fodo9Fo+XNmZiYtWrRg6dKlhY6tXr06AIsWLWL48OGsW7eO5cuXM378eDZu3Ejr1q0LHaPVmq8qKNhnwf5u8fb2LjXOkmRmZgKwevVqq2IUsCpqAQICAmw6Z2hoKPXr17faFxgYaLl+qlq1auh0OtLT063apKenl1r4CiHkmiMh3I7BYKBXr17ExMQwYcIE+vfvX6Y7kM6cOcOFCxcs2zt37kSr1dKoUSObjo+NjWXp0qV88803aLVaOnfubHlv27ZtdO3alV69etG0aVPq1q3LsWPHSjxf9erVrS50Pn78ONnZ2Zbt5s2bc/z4cWrUqEH9+vWtXpUrV7a0e/jhhxkzZgzbt2+ncePGfPbZZ8X2B1j1eeDAAZs++y2VK1cmLCyMH374wbIvLy+PvXv3WrYfeOAB9Ho9Z86cKRR3eHi4Xf3ZysfHhxYtWlhdpG8ymdi0aRNt2rRxSJ9CeBIpjoRwM+PGjeP69evMmTOHUaNG0bBhQ1566SXL++fPn+e+++4r9eJnX19f+vTpw48//siWLVsYPnw40dHRNo8sxMbGsm/fPiZNmsQLL7xgNQrSoEEDNm7cyPbt2zl8+DCDBg0qNIpxpw4dOvD++++zf/9+9uzZwyuvvGI1ahMbG0u1atXo2rUrW7ZsITU1leTkZIYPH865c+dITU1lzJgx7Nixg9OnT7NhwwaOHz/O/fffX2R/t4qThIQEjh8/zurVq5k5c6ZNn72gESNGMHXqVFatWsWRI0cYMmSI1R1wgYGBvP7664wcOZLExEROnDjBvn37mDt3LomJiXb3Z6u4uDgWLFhAYmIihw8fZvDgwWRlZdGvXz+H9SmEp5BpNSHcSHJyMrNmzeL7778nKCgIgCVLltC0aVM+/PBDBg8ejNFo5OjRo1ajLkWpX78+zz//PJ06deLq1as888wzfPDBBzbHUr9+fVq1asWuXbuYNWuW1Xvjx4/n5MmTREVF4e/vz8CBA+nWrVuJjwqYOXMm/fr1o127dtSsWZPZs2dbjcD4+/uzefNmRo0axfPPP09GRga1atXiySefJCgoiJs3b3LkyBESExO5cuUKYWFhDB06lEGDBhXZn7e3N8uWLWPw4ME89NBDPPLII0ycOJHu3bvbnAOA1157jYsXL9KnTx+0Wi0vvfQSzz33nNVnnTBhAtWrV2fKlCmcPHmS4OBgmjdvztixY+3qyx4xMTFcunSJt99+m7S0NJo1a8a6desICQlxWJ9CeAqNsmfyXAjhERISEli1apXd00hCCHE3kGk1IYQQQogCpDgSQhRSqVKlYl9btmxxdXhCCOFQMq0mhCgkJSWl2Pdq1aqFn5+fE6MRQgjnkuJICCGEEKIAmVYTQgghhChAiiMhhBBCiAKkOBJCCCGEKECKIyGEEEKIAqQ4EkIIIYQoQIojIYQQQogCpDgSQgghhChAiiMhhBBCiAL+H+dS3/6a6Q2IAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "for _ in range(100):\n", + " X, F = generate_randunif()\n", + " plt.plot(X, F, 'b')\n", + " \n", + "ax.axline((0, 0), (1, 1), linewidth=1, color='r')\n", + "plt.plot(X_pv_exact, F_pv, 'g', label=\"Exact p_value\")\n", + "plt.plot(X_pv_approximation, F_pv, 'y', label=\"p_value without refitting\")\n", + "\n", + "plt.legend(\n", + " loc='upper left',\n", + " fontsize=8,\n", + ")\n", + "plt.xlabel(\"x: p_values under H0\")\n", + "plt.ylabel(\"F(X)\")\n", + "plt.title(\"Cumulative distribution function value of the p-values under H0\")\n", + "plt.grid()\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.8.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/tutorials/plot_tuto_mcar.py b/examples/tutorials/plot_tuto_mcar.py index c43d1217..cbd2d2fd 100644 --- a/examples/tutorials/plot_tuto_mcar.py +++ b/examples/tutorials/plot_tuto_mcar.py @@ -3,7 +3,7 @@ Tutorial for Testing the MCAR Case ============================================ -In this tutorial, we show how to test the MCAR case using the Little's test. +In this tutorial, we show how to test the MCAR case using the Little and the PKLM tests. """ # %% @@ -14,7 +14,7 @@ import pandas as pd from scipy.stats import norm -from qolmat.analysis.holes_characterization import LittleTest +from qolmat.analysis.holes_characterization import LittleTest, PKLMTest from qolmat.benchmark.missing_patterns import UniformHoleGenerator plt.rcParams.update({"font.size": 12}) @@ -31,22 +31,32 @@ q975 = norm.ppf(0.975) # %% +# 1. Testing the MCAR case with the Little's test and the PKLM test. +# ------------------------------------------------------------------ +# # The Little's test -# --------------------------------------------------------------- +# ================= +# # First, we need to introduce the concept of a missing pattern. A missing pattern, also called a # pattern, is the structure of observed and missing values in a dataset. For example, in a # dataset with two columns, the possible patterns are: (0, 0), (1, 0), (0, 1), (1, 1). The value 1 # (0) indicates that the column value is missing (observed). # # The null hypothesis, H0, is: "The means of observations within each pattern are similar.". + +# %% +# The PKLM test +# ============= +# The test compares distributions of different missing patterns. # +# The null hypothesis, H0, is: "Distributions within each pattern are similar.". # We choose to use the classic threshold of 5%. If the test p-value is below this threshold, # we reject the null hypothesis. -# -# This notebook shows how the Little's test performs on a simplistic case and its limitations. We -# instanciate a test object with a random state for reproducibility. +# This notebook shows how the Little and PKLM tests perform on a simplistic case and their +# limitations. We instantiate a test object with a random state for reproducibility. -test_mcar = LittleTest(random_state=rng) +little_test_mcar = LittleTest(random_state=rng) +pklm_test_mcar = PKLMTest(random_state=rng) # %% # Case 1: MCAR holes (True negative) @@ -77,11 +87,13 @@ plt.show() # %% -result = test_mcar.test(df_nan) -print(f"Test p-value: {result:.2%}") +little_result = little_test_mcar.test(df_nan) +pklm_result = pklm_test_mcar.test(df_nan) +print(f"The p-value of the Little's test is: {little_result:.2%}") +print(f"The p-value of the PKLM test is: {pklm_result:.2%}") # %% -# The p-value is larger than 0.05, therefore we don't reject the HO MCAR assumption. In this case -# this is a true negative. +# The two p-values are larger than 0.05, therefore we don't reject the H0 MCAR assumption. +# In this case this is a true negative. # %% # Case 2: MAR holes with mean bias (True positive) @@ -110,11 +122,13 @@ # %% -result = test_mcar.test(df_nan) -print(f"Test p-value: {result:.2%}") +little_result = little_test_mcar.test(df_nan) +pklm_result = pklm_test_mcar.test(df_nan) +print(f"The p-value of the Little's test is: {little_result:.2%}") +print(f"The p-value of the PKLM test is: {pklm_result:.2%}") # %% -# The p-value is smaller than 0.05, therefore we reject the HO MCAR assumption. In this case -# this is a true positive. +# The two p-values are smaller than 0.05, therefore we reject the H0 MCAR assumption. +# In this case this is a true positive. # %% # Case 3: MAR holes with any mean bias (False negative) @@ -149,17 +163,221 @@ # %% -result = test_mcar.test(df_nan) -print(f"Test p-value: {result:.2%}") +little_result = little_test_mcar.test(df_nan) +pklm_result = pklm_test_mcar.test(df_nan) +print(f"The p-value of the Little's test is: {little_result:.2%}") +print(f"The p-value of the PKLM test is: {pklm_result:.2%}") # %% -# The p-value is larger than 0.05, therefore we don't reject the HO MCAR assumption. In this case -# this is a false negative since the missingness mechanism is MAR. +# The Little's p-value is larger than 0.05, therefore, using this test we don't reject the H0 MCAR +# assumption. In this case this is a false negative since the missingness mechanism is MAR. +# +# However the PKLM test p-value is smaller than 0.05 therefore we don't reject the H0 MCAR +# assumption. In this case this is a true negative. # %% -# Limitations -# ----------- +# Limitations and conclusion +# ========================== # In this tutoriel, we can see that Little's test fails to detect covariance heterogeneity between # patterns. # # We also note that the Little's test does not handle categorical data or temporally # correlated data. +# +# This is why we have implemented the PKLM test, which makes up for the shortcomings of the Little +# test. We present this test in more detail in the next section. + +# %% +# 2. The PKLM test +# ------------------------------------------------------------------ +# +# The PKLM test is very powerful for several reasons. Firstly, it covers the concerns that Little's +# test may have (covariance heterogeneity). Secondly, it is currently the only MCAR test applicable +# to mixed data. Finally, it proposes a concept of partial p-value which enables us to carry out a +# variable-by-variable diagnosis to identify the potential causes of a MAR mechanism. +# +# There is a parameter in the paper called size.res.set. The authors of the paper recommend setting +# this parameter to 2. We have chosen to follow this advice and not leave the possibility of +# increasing this parameter. The results are satisfactory and the code is simpler. +# +# It does have one disadvantage, however: its calculation time. +# + +# %% + +""" +Calculation time +================ + ++------------+------------+----------------------+ +| **n_rows** | **n_cols** | **Calculation_time** | ++============+============+======================+ +| 200 | 2 | 2"12 | ++------------+------------+----------------------+ +| 500 | 2 | 2"24 | ++------------+------------+----------------------+ +| 500 | 4 | 2"18 | ++------------+------------+----------------------+ +| 1000 | 4 | 2"48 | ++------------+------------+----------------------+ +| 1000 | 6 | 2"42 | ++------------+------------+----------------------+ +| 10000 | 6 | 20"54 | ++------------+------------+----------------------+ +| 10000 | 10 | 14"48 | ++------------+------------+----------------------+ +| 100000 | 10 | 4'51" | ++------------+------------+----------------------+ +| 100000 | 15 | 3'06" | ++------------+------------+----------------------+ +""" + +# %% +# 2.1 Parameters and Hyperparmaters +# ================================================ +# +# To use the PKLM test properly, it may be necessary to understand the use of hyper-parameters. +# +# * ``nb_projections``: Number of projections on which the test statistic is calculated. This +# parameter has the greatest influence on test calculation time. Its defaut value +# ``nb_projections=100``. +# Est-ce qu'on donne des ordres de grandeurs utiles ? J'avais un peu fait ce travail. +# +# * ``nb_permutation`` : Number of permutations of the projected targets. The higher is better. This +# parameter has little impact on calculation time. +# Its default value ``nb_permutation=30``. +# +# * ``nb_trees_per_proj`` : The number of subtrees in each random forest fitted. In order to +# estimate the Kullback-Leibler divergence, we need to obtain probabilities of belonging to +# certain missing patterns. Random Forests are used to estimate these probabilities. This +# hyperparameter has a significant impact on test calculation time. Its default +# value is ``nb_trees_per_proj=200`` +# +# * ``compute_partial_p_values``: Boolean that indicates if you want to compute the partial +# p-values. Those partial p-values could help the user to identify the variables responsible for +# the MAR missing-data mechanism. Please see the section 2.3 for examples. Its default value is +# ``compute_partial_p_values=False``. +# +# * ``encoder``: Scikit-Learn encoder to encode non-numerical values. +# Its default value ``encoder=sklearn.preprocessing.OneHotEncoder()`` +# +# * ``random_state``: Controls the randomness. Pass an int for reproducible output across +# multiple function calls. Its default value ``random_state=None`` + +# %% +# 2.2 Application on mixed data types +# ================================================ +# +# As we have seen, Little's test only applies to quantitative data. In real life, however, it is +# common to have to deal with mixed data. Here's an example of how to use the PKLM test on a dataset +# with mixed data types. + +# %% +n_rows = 100 + +col1 = rng.rand(n_rows) * 100 +col2 = rng.randint(1, 100, n_rows) +col3 = rng.choice([True, False], n_rows) +modalities = ['A', 'B', 'C', 'D'] +col4 = rng.choice(modalities, n_rows) + +df = pd.DataFrame({ + 'Numeric1': col1, + 'Numeric2': col2, + 'Boolean': col3, + 'Object': col4 +}) + +hole_gen = UniformHoleGenerator( + n_splits=1, + ratio_masked=0.2, + subset=['Numeric1', 'Numeric2', 'Boolean', 'Object'], + random_state=rng +) +df_mask = hole_gen.generate_mask(df) +df_nan = df.where(~df_mask, np.nan) +df_nan.dtypes + +# %% +pklm_result = pklm_test_mcar.test(df_nan) +print(f"The p-value of the PKLM test is: {pklm_result:.2%}") + +# %% +# To perform the PKLM test over mixed data types, non numerical features need to be encoded. The +# default encoder in the :class:`~qolmat.analysis.holes_characterization.PKLMTest` class is the +# default OneHotEncoder from scikit-learn. If you wish to use an encoder adapted to your data, you +# can perform this encoding step beforehand, and then use the PKLM test. +# Currently, we do not support the following types : +# +# - datetimes +# +# - timedeltas +# +# - Pandas datetimetz + +# %% +# 2.3 Partial p-values +# ================================================ +# +# In addition, the PKLM test can be used to calculate partial p-values. We denote as many partial +# p-values as there are columns in the input dataframe. This “partial” p-value corresponds to the +# effect of removing the patterns induced by variable k. +# +# Let's take a look at an example of how to use this feature + +# %% +data = rng.multivariate_normal( + mean=[0, 0, 0, 0], + cov=[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], + size=400 +) +df = pd.DataFrame(data=data, columns=["Column 1", "Column 2", "Column 3", "Column 4"]) + +df_mask = pd.DataFrame( + { + "Column 1": False, + "Column 2": df["Column 1"] > q975, + "Column 3": False, + "Column 4": False, + }, + index=df.index +) +df_nan = df.where(~df_mask, np.nan) + +# %% +# The missing-data mechanism is clearly MAR. Intuitively, if we remove the second column from the +# matrix, the missing-data mechanism is MCAR. Let's see how the PKLM test can help us identify the +# variable responsible for the MAR mechanism. + +# %% +pklm_test = PKLMTest(random_state=rng, compute_partial_p_values=True) +p_value, partial_p_values = pklm_test.test(df_nan) +print(f"The p-value of the PKLM test is: {p_value:.2%}") + +# %% +# The test result confirms that we can reject the null hypothesis and therefore assume that the +# missing-data mechanism is MAR. +# Let's now take a look at what partial p-values can tell us. + +# %% +for col_index, partial_p_v in enumerate(partial_p_values): + print(f"The partial p-value for the column index {col_index + 1} is: {partial_p_v:.2%}") + +# %% +# As a result, by removing the missing patterns induced by variable 2, the p-value rises +# above the significance threshold set beforehand. Thus in this sense, the test detects that the +# main culprit of the MAR mechanism lies in the second variable. + + +# %% +# Calculation time -> TO BE DELETED +# | **n_rows** | **n_cols** | **Calculation_time** | +# |------------|------------|----------------------| +# | 200 | 2 | 2"12 | +# | 500 | 2 | 2"24 | +# | 500 | 4 | 2"18 | +# | 1000 | 4 | 2"48 | +# | 1000 | 6 | 2"42 | +# | 10000 | 6 | 20"54 | +# | 10000 | 10 | 14"48 | +# | 100000 | 10 | 4'51" | +# | 100000 | 15 | 3'06" | \ No newline at end of file diff --git a/qolmat/analysis/holes_characterization.py b/qolmat/analysis/holes_characterization.py index 5669ac7b..7e585ee4 100644 --- a/qolmat/analysis/holes_characterization.py +++ b/qolmat/analysis/holes_characterization.py @@ -1,21 +1,61 @@ from abc import ABC, abstractmethod -from typing import Optional, Union +from itertools import combinations +from typing import List, Optional, Tuple, Union import numpy as np import pandas as pd +from category_encoders.one_hot import OneHotEncoder +from joblib import Parallel, delayed from scipy.stats import chi2 +from sklearn import utils as sku +from sklearn.ensemble import RandomForestClassifier from qolmat.imputations.imputers import ImputerEM +from qolmat.utils.input_check import check_pd_df_dtypes class McarTest(ABC): """ - Astract class for MCAR tests. + Abstract base class for performing MCAR (Missing Completely At Random) tests. + + Parameters + ---------- + random_state : int or np.random.RandomState, optional + Seed or random state for reproducibility. + + Methods + ------- + test(df) + Abstract method to perform the MCAR test on the given DataFrame or NumPy array. """ + def __init__(self, random_state: Union[None, int, np.random.RandomState] = None): + """ + Initializes the McarTest class with a random state. + + Parameters + ---------- + random_state : int or np.random.RandomState, optional + Seed or random state for reproducibility. + """ + self.rng = sku.check_random_state(random_state) + @abstractmethod - def test(self, df: pd.DataFrame) -> float: - pass + def test(self, df: Union[pd.DataFrame, np.ndarray]) -> Union[float, Tuple[float, List[float]]]: + """ + Perform the MCAR test on the input data. + + Parameters + ---------- + df : pd.DataFrame or np.ndarray + Data to be tested for MCAR. Can be provided as a pandas DataFrame or a NumPy array. + + Returns + ------- + float or tuple of float and list of float + Test statistic, or a tuple with the test statistic and additional details if applicable. + """ + raise NotImplemented class LittleTest(McarTest): @@ -27,7 +67,7 @@ class LittleTest(McarTest): References ---------- Little. "A Test of Missing Completely at Random for Multivariate Data with Missing Values." - Journal of the American Statistical Association, Volume 83, 1988 - Issue 404 + Journal of the American Statistical Association, Volume 83, 1988 - no 404, p. 1198-1202 Parameters ---------- @@ -44,17 +84,16 @@ def __init__( imputer: Optional[ImputerEM] = None, random_state: Union[None, int, np.random.RandomState] = None, ): - super().__init__() + super().__init__(random_state=random_state) if imputer and imputer.model != "multinormal": raise AttributeError( "The ImputerEM model must be 'multinormal' to use the Little's test" ) self.imputer = imputer - self.random_state = random_state def test(self, df: pd.DataFrame) -> float: """ - Apply the Little's test over a real dataframe. + Apply the Little's test to a real dataframe. Parameters @@ -67,7 +106,7 @@ def test(self, df: pd.DataFrame) -> float: float The p-value of the test. """ - imputer = self.imputer or ImputerEM(random_state=self.random_state) + imputer = self.imputer or ImputerEM(random_state=self.rng) imputer = imputer._fit_element(df) d0 = 0 @@ -86,9 +125,571 @@ def test(self, df: pd.DataFrame) -> float: obs_mean = df_pattern.mean().to_numpy() diff_means = obs_mean - ml_means[list(tup_pattern)] - inv_sigma_pattern = np.linalg.inv(ml_cov[:, tup_pattern][tup_pattern, :]) + inv_sigma_pattern = np.linalg.solve( + ml_cov[:, tup_pattern][tup_pattern, :], + np.eye(len(tup_pattern)) + ) - d0 += n_rows_pattern * np.dot(np.dot(diff_means, inv_sigma_pattern), diff_means.T) + d0 += n_rows_pattern * np.dot( + np.dot(diff_means, inv_sigma_pattern), diff_means.T + ) degree_f += tup_pattern.count(True) return 1 - float(chi2.cdf(d0, degree_f)) + + +class PKLMTest(McarTest): + """ + This class implements the PKLM test, a fully non-parametric, easy-to-use, and powerful test + for the missing completely at random (MCAR) assumption on the missingness mechanism of a + dataset. The null hypothesis is "The missing data mechanism is MCAR". + + This test is applicable to mixed data (quantitative and categoricals) types. + + + If you're familiar with the paper, this implementation of the PKLM test was made for the + parameter size.resp.set=2 only. + + References + ---------- + Spohn, M. L., Näf, J., Michel, L., & Meinshausen, N. (2021). PKLM: A flexible MCAR test using + Classification. arXiv preprint arXiv:2109.10150. + + Parameters + ----------- + nb_projections : int + Number of projections. + nb_projections_threshold : int + If the maximum number of possible permutations is less than this threshold, then all + projections are used. Otherwise, nb_projections random projections are drawn. + nb_permutation : int + Number of permutations. + nb_trees_per_proj : int + Number of trees per projection. + compute_partial_p_values : bool + If true, compute the partial p-values. + exact_p_value : bool + If True, compute exact p-value. + encoder : OneHotEncoder or None, default=None + Encoder to convert non numeric pandas dataframe values to numeric values. + random_state : int, RandomState instance or None, default=None + Controls the randomness. + Pass an int for reproducible output across multiple function calls. + """ + + def __init__( + self, + nb_projections: int = 100, + nb_projections_threshold: int = 200, + nb_permutation: int = 30, + nb_trees_per_proj: int = 200, + compute_partial_p_values: bool = False, + exact_p_value: bool = False, + encoder: Union[None, OneHotEncoder] = None, + random_state: Union[None, int, np.random.RandomState] = None, + ): + super().__init__(random_state=random_state) + self.nb_projections = nb_projections + self.nb_projections_threshold = nb_projections_threshold + self.nb_permutation = nb_permutation + self.nb_trees_per_proj = nb_trees_per_proj + self.compute_partial_p_values = compute_partial_p_values + self.exact_p_value = exact_p_value + self.encoder = encoder + + if self.exact_p_value: + self.process_permutation = self._parallel_process_permutation_exact + else: + self.process_permutation = self._parallel_process_permutation + + def _encode_dataframe(self, df: pd.DataFrame) -> np.ndarray: + """ + Encodes the DataFrame by converting numeric columns to a numpy array + and applying one-hot encoding to objects and boolean columns. + + Parameters: + ----------- + df : pd.DataFrame + The DataFrame to be encoded. + + Returns: + ------- + np.ndarray + The encoded DataFrame as a numpy ndarray, with numeric data concatenated + with one-hot encoded categorical and boolean data. + """ + if not df.select_dtypes(include=["object", "bool"]).columns.to_list(): + return df.to_numpy() + + if not self.encoder: + self.encoder = OneHotEncoder( + cols=df.select_dtypes(include=["object", "bool"]).columns, + return_df=False, + handle_missing='return_nan' + ) + + return self.encoder.fit_transform(df) + + def _pklm_preprocessing(self, X: Union[pd.DataFrame, np.ndarray]) -> np.ndarray: + """ + Preprocesses the input DataFrame or ndarray for further processing. + + Parameters: + ----------- + X : Union[pd.DataFrame, np.ndarray] + The input data to be preprocessed. Can be a pandas DataFrame or a numpy ndarray. + + Returns: + ------- + np.ndarray + The preprocessed data as a numpy ndarray. + + Raises: + ------- + TypeNotHandled + If the DataFrame contains columns with data types that are not numeric, string, or + boolean. + """ + if isinstance(X, np.ndarray): + return X + + check_pd_df_dtypes( + X, + [ + pd.api.types.is_numeric_dtype, + pd.api.types.is_string_dtype, + pd.api.types.is_bool_dtype + ] + ) + return self._encode_dataframe(X) + + @staticmethod + def _get_max_draw(p: int) -> int: + """ + Calculates the number of possible projections. + + Parameters: + ----------- + p : int + The number of columns of the input matrix. + + Returns: + -------- + int + The number of possible projections. + """ + return p*(2**(p-1) - 1) + + def _draw_features_and_target_indexes(self, X: np.ndarray) -> Tuple[List[int], int]: + """ + Randomly selects features and a target from the dataframe. + This corresponds to the Ai and Bi projections of the paper. + + Parameters: + ----------- + X : np.ndarray + The input dataframe. + + Returns: + -------- + Tuple[np.ndarray, int] + Indices of selected features and the target. + """ + _, p = X.shape + nb_features = self.rng.randint(1, p) + features_idx = self.rng.choice(p, size=nb_features, replace=False) + target_idx = self.rng.choice(np.setdiff1d(np.arange(p), features_idx)) + return features_idx.tolist(), target_idx + + @staticmethod + def _check_draw(X: np.ndarray, features_idx: List[int], target_idx: int) -> np.bool_: + """ + Checks if the drawn features and target are valid. Here we check that the number of induced + classes is equal to 2. Using the notation from the paper, we want |G(Ai, Bi)| = 2. + + Parameters: + ----------- + X : np.ndarray + The input dataframe. + features_idx : np.ndarray + Indices of the selected features. + target_idx : int + Index of the target. + + Returns: + -------- + bool + True if the draw is valid, False otherwise. + """ + target_values = X[~np.isnan(X[:, features_idx]).any(axis=1)][:, target_idx] + is_nan = np.isnan(target_values).any() + is_distinct_values = (~np.isnan(target_values)).any() + return is_nan and is_distinct_values + + def _generate_label_feature_combinations(self, X: np.ndarray) -> List[Tuple[int, List[int]]]: + """ + Generates all valid combinations of features and labels for projection if + nb_projections_threshold > _get_max_draw(X.shape[1]). + + Parameters: + ----------- + X : np.ndarray + The input data array. + + Returns: + -------- + List[Tuple[int, List[int]]] + A list of tuples where each tuple contains a label and a list of selected features that + can be used for projection. + """ + _, p = X.shape + indices = list(range(p)) + result = [] + + for label in indices: + feature_candidates = [i for i in indices if i != label] + + for r in range(1, len(feature_candidates) + 1): + for feature_set in combinations(feature_candidates, r): + if self._check_draw(X, list(feature_set), label): + result.append((list(feature_set), label)) + + return result + + def _draw_projection(self, X: np.ndarray) -> Tuple[List[int], int]: + """ + Draws a valid projection of features and a target. + If nb_projections_threshold < _get_max_draw(X.shape[1]). + + Parameters: + ----------- + X : np.ndarray + The input dataframe. + + Returns: + -------- + Tuple[np.ndarray, int] + Indices of selected features and the target. + """ + is_checked = False + while not is_checked: + features_idx, target_idx = self._draw_features_and_target_indexes(X) + is_checked = self._check_draw(X, features_idx, target_idx) + return features_idx, target_idx + + @staticmethod + def _build_dataset( + X: np.ndarray, + features_idx: np.ndarray, + target_idx: int + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Builds a dataset by selecting specified features and target from a NumPy array, excluding + rows with NaN values in the feature columns. + For the label, we create a binary classification problem where yi =1 if target_idx_i is + missing. + + Parameters: + ----------- + X: np.ndarray + Input data array. + features_idx: np.ndarray + Indices of the feature columns. + target_idx: int + Index of the target column. + + Returns: + -------- + Tuple[np.ndarray, np.ndarray]: A tuple containing: + - X (np.ndarray): Full observed array of selected features. + - y (np.ndarray): Binary array indicating presence of NaN (1) in the target column. + """ + X_features = X[~np.isnan(X[:, features_idx]).any(axis=1)][:, features_idx] + y = np.where( + np.isnan(X[~np.isnan(X[:, features_idx]).any(axis=1)][:, target_idx]), + 1, + 0, + ) + return X_features, y + + @staticmethod + def _build_label( + X: np.ndarray, + perm: np.ndarray, + features_idx: np.ndarray, + target_idx: int + ) -> np.ndarray: + """ + Builds a label array by selecting target values from a permutation array, + excluding rows with NaN values in the specified feature columns. + For the label, we create a binary classification problem where yi =1 if target_idx_i is + missing. + + Parameters: + ----------- + X: np.ndarray + Input data array. + perm: np.ndarray + Permutation array from which labels are selected. + features_idx: np.ndarray + Indices of the feature columns. + target_idx: int + Index of the target column in the permutation array. + + Returns: + -------- + np.ndarray: Binary array indicating presence of NaN (1) in the target column. + """ + return perm[~np.isnan(X[:, features_idx]).any(axis=1), target_idx] + + def _get_oob_probabilities(self, X: np.ndarray, y: np.ndarray) -> np.ndarray: + """ + Trains a RandomForestClassifier and retrieves out-of-bag (OOB) probabilities. + + Parameters: + ----------- + X: np.ndarray + Feature array for training. + y: np.ndarray + Target array for training. + + Returns: + -------- + np.ndarray: Out-of-bag probabilities for each class. + """ + clf = RandomForestClassifier( + n_estimators=self.nb_trees_per_proj, + min_samples_split=10, + bootstrap=True, + oob_score=True, + random_state=self.rng, + max_features=1., + ) + clf.fit(X, y) + return clf.oob_decision_function_ + + @staticmethod + def _U_hat(oob_probabilities: np.ndarray, labels: np.ndarray) -> float: + """ + Computes the U_hat statistic, a measure of classifier performance, using + out-of-bag probabilities and true labels. + + Parameters: + ----------- + oob_probabilities: np.ndarray + Out-of-bag probabilities for each class. + labels: np.ndarray + True labels for the data. + + Returns: + -------- + float: The computed U_hat statistic. + """ + if oob_probabilities.shape[1] == 1: + return 0. + + oob_probabilities = np.clip(oob_probabilities, 1e-9, 1 - 1e-9) + + unique_labels = np.unique(labels) + label_matrix = (labels[:, None] == unique_labels).astype(int) + p_true = oob_probabilities * label_matrix + p_false = oob_probabilities * (1 - label_matrix) + + p0_0 = p_true[:, 0][np.where(p_true[:, 0] != 0.0)] + p0_1 = p_false[:, 0][np.where(p_false[:, 0] != 0.0)] + p1_1 = p_true[:, 1][np.where(p_true[:, 1] != 0.0)] + p1_0 = p_false[:, 1][np.where(p_false[:, 1] != 0.0)] + + if unique_labels.shape[0] == 1: + if unique_labels[0] == 0: + n0 = labels.shape[0] + return ( + np.log(p0_0 / (1 - p0_0)).sum() / n0 + - np.log(p1_0 / (1 - p1_0)).sum() / n0 + ) + else: + n1 = labels.shape[0] + return ( + np.log(p1_1 / (1 - p1_1)).sum() / n1 + - np.log(p0_1 / (1 - p0_1)).sum() / n1 + ) + + n0, n1 = label_matrix.sum(axis=0) + u_0 = ( + np.log(p0_0 / (1 - p0_0)).sum() / n0 - np.log(p0_1 / (1 - p0_1)).sum() / n1 + ) + u_1 = ( + np.log(p1_1 / (1 - p1_1)).sum() / n1 - np.log(p1_0 / (1 - p1_0)).sum() / n0 + ) + + return u_0 + u_1 + + def _parallel_process_permutation( + self, + X: np.ndarray, + M_perm: np.ndarray, + features_idx: np.ndarray, + target_idx: int, + oob_probabilities: np.ndarray, + ) -> float: + y = self._build_label(X, M_perm, features_idx, target_idx) + return self._U_hat(oob_probabilities, y) + + def _parallel_process_permutation_exact( + self, + X: np.ndarray, + M_perm: np.ndarray, + features_idx: np.ndarray, + target_idx: int, + oob_probabilites_unused: np.ndarray, + ) -> float: + X_features, _ = self._build_dataset(X, features_idx, target_idx) + y = self._build_label(X, M_perm, features_idx, target_idx) + # In this case, we fit the classifier in each permutation. It takes much more longer. + oob_probabilities = self._get_oob_probabilities(X_features, y) + return self._U_hat(oob_probabilities, y) + + def _parallel_process_projection( + self, + X: np.ndarray, + list_permutations: List[np.ndarray], + features_idx: np.ndarray, + target_idx: int, + ) -> Tuple[float, List[float]]: + X_features, y = self._build_dataset(X, features_idx, target_idx) + oob_probabilities = self._get_oob_probabilities(X_features, y) + u_hat = self._U_hat(oob_probabilities, y) + # We iterate over the permutation because for a given projection, we fit only one classifier + # to get oob probabilities and compute u_hat nb_permutations times. + result_u_permutations = Parallel(n_jobs=-1)( + delayed(self.process_permutation)( + X, M_perm, features_idx, target_idx, oob_probabilities + ) + for M_perm in list_permutations + ) + return u_hat, result_u_permutations + + + @staticmethod + def _build_B(list_proj: List, n_cols: int) -> np.ndarray: + """ + Constructs a binary matrix B based on the given projections. + + Parameters: + ----------- + list_proj : List + A list of tuples where each tuple represents a projection, and the + second element of each tuple is an index used to build the target. + n_cols : int + The number of columns in the resulting matrix B. + + Returns: + -------- + np.ndarray + A binary matrix of shape (n_cols, len(list_proj)) where each column corresponds to a + projection, and the entries are 0 or 1 based on the projections. + """ + list_bi = [projection[1] for projection in list_proj] + B = np.ones((len(list_proj), n_cols), dtype=int) + + for j in range(len(list_bi)): + B[j, list_bi[j]] = 0 + + return B.transpose() + + def _compute_partial_p_value( + self, + B: np.ndarray, + U: np.ndarray, + U_sigma: np.ndarray, + k: int + ) -> float: + """ + Computes the partial p-value for a statistical test based on a given permutation. + + Parameters: + ----------- + B : np.ndarray + Pass matrix indicating the column used to create the target in each projection. + U : np.ndarray + A vector of shape (nb_permutations,) representing the observed test statistics. + U_sigma : np.ndarray + A matrix of shape (nb_permutations, nb_observations) where each row represents the test + statistics for a given projection and all the permutations. + k : int + The index of the column on which to compute the partial p_value. + + Returns: + -------- + float + The partial p-value. + """ + U_k = B[k, :]@U + p_v_k = 1 + + for u_sigma_k in (B[k, :]@U_sigma).tolist(): + if u_sigma_k >= U_k: + p_v_k += 1 + + return p_v_k / (self.nb_permutation + 1) + + def test(self, X: Union[pd.DataFrame, np.ndarray]) -> Union[float, Tuple[float, List[float]]]: + """ + Apply the PKLM test over a real dataset. + + Parameters + ---------- + X : np.ndarray + The input dataset with missing values. + + Returns + ------- + float + If compute_partial_p_values=False. Returns the p-value of the test. + Tuple[float, List[float]] + If compute_partial_p_values=True. Returns the p-value of the test and the list of all + the partial p-values. + """ + X = self._pklm_preprocessing(X) + _, n_cols = X.shape + + if self._get_max_draw(n_cols) <= self.nb_projections_threshold: + list_proj = self._generate_label_feature_combinations(X) + else: + list_proj = [self._draw_projection(X) for _ in range(self.nb_projections)] + + M = np.isnan(X).astype(int) + list_perm = [self.rng.permutation(M) for _ in range(self.nb_permutation)] + U = 0.0 + list_U_sigma = [0.0 for _ in range(self.nb_permutation)] + + parallel_results = Parallel(n_jobs=-1)( + delayed(self._parallel_process_projection)( + X, list_perm, features_idx, target_idx + ) + for features_idx, target_idx in list_proj + ) + + for U_projection, results in parallel_results: + U += U_projection + list_U_sigma = [x + y for x, y in zip(list_U_sigma, results)] + + # Je suggère d'alléger le code de cette manipulation même si théoriquement ça a de la valeur + U = U / self.nb_projections + list_U_sigma = [x / self.nb_permutation for x in list_U_sigma] + + p_value = 1 + for u_sigma in list_U_sigma: + if u_sigma >= U: + p_value += 1 + + p_value = p_value / (self.nb_permutation + 1) + + if not self.compute_partial_p_values: + return p_value + else: + B = self._build_B(list_proj, n_cols) + U = np.array([item[0] for item in parallel_results]) + U_sigma = np.array([item[1] for item in parallel_results]) + p_values = [self._compute_partial_p_value(B, U, U_sigma, k) for k in range(n_cols)] + return p_value, p_values diff --git a/qolmat/utils/exceptions.py b/qolmat/utils/exceptions.py index 513e843b..b03375db 100644 --- a/qolmat/utils/exceptions.py +++ b/qolmat/utils/exceptions.py @@ -37,7 +37,7 @@ def __init__(self, shape: Tuple[int, ...]): class NotDataFrame(Exception): def __init__(self, X_type: Type[Any]): - super().__init__(f"Input musr be a dataframe, not a {X_type}") + super().__init__(f"Input must be a dataframe, not a {X_type}") class NotEnoughSamples(Exception): diff --git a/qolmat/utils/input_check.py b/qolmat/utils/input_check.py new file mode 100644 index 00000000..1aeec922 --- /dev/null +++ b/qolmat/utils/input_check.py @@ -0,0 +1,29 @@ +import pandas as pd + +from qolmat.utils.exceptions import TypeNotHandled + +def check_pd_df_dtypes(df: pd.DataFrame, allowed_types: list): + """ + Validates that the columns of the DataFrame have allowed data types. + + Parameters: + ----------- + df : pd.DataFrame + DataFrame whose columns' data types are to be checked. + + Raises: + ------- + TypeNotHandled + If any column has a data type that is not numeric, string, or boolean. + """ + def is_allowed_type(dtype): + return any(check(dtype) for check in allowed_types) + + invalid_columns = [ + (col, dtype) + for col, dtype in df.dtypes.items() + if not is_allowed_type(dtype) + ] + if invalid_columns: + for column_name, dtype in invalid_columns: + raise TypeNotHandled(col=str(column_name), type_col=dtype) diff --git a/tests/analysis/test_holes_characterization.py b/tests/analysis/test_holes_characterization.py index c794b94e..6c9b63f6 100644 --- a/tests/analysis/test_holes_characterization.py +++ b/tests/analysis/test_holes_characterization.py @@ -3,9 +3,12 @@ import pytest from scipy.stats import norm -from qolmat.analysis.holes_characterization import LittleTest +from qolmat.analysis.holes_characterization import LittleTest, PKLMTest from qolmat.benchmark.missing_patterns import UniformHoleGenerator from qolmat.imputations.imputers import ImputerEM +from qolmat.utils.exceptions import TypeNotHandled + +### Tests for the LittleTest class @pytest.fixture @@ -60,3 +63,180 @@ def test_little_mcar_test(df_input: pd.DataFrame, expected: bool, request): def test_attribute_error(): with pytest.raises(AttributeError): LittleTest(random_state=42, imputer=ImputerEM(model="VAR")) + + +### Tests for the PKLMTest class + + +@pytest.fixture +def supported_multitypes_dataframe() -> pd.DataFrame: + return pd.DataFrame({ + 'int_col': [1, 2, 3], + 'float_col': [1.1, 2.2, 3.3], + 'str_col': ['a', 'b', 'c'], + 'bool_col': [True, False, True] + }) + + +@pytest.fixture +def np_matrix_with_nan_mcar() -> np.ndarray: + rng = np.random.default_rng(42) + n_rows, n_cols = 10, 4 + matrix = rng.normal(size=(n_rows, n_cols)) + num_nan = int(n_rows * n_cols * 0.40) + nan_indices = rng.choice(n_rows * n_cols, num_nan, replace=False) + matrix.flat[nan_indices] = np.nan + return matrix + + +@pytest.fixture +def missingness_matrix_mcar(np_matrix_with_nan_mcar): + return np.isnan(np_matrix_with_nan_mcar).astype(int) + + +@pytest.fixture +def missingness_matrix_mcar_perm(missingness_matrix_mcar): + rng = np.random.default_rng(42) + return rng.permutation(missingness_matrix_mcar) + + +@pytest.fixture +def oob_probabilities() -> np.ndarray: + return np.matrix([[0.5, 0.5], [0, 1], [1, 0], [1, 0]]).A + +def test__encode_dataframe(supported_multitypes_dataframe): + mcar_test_pklm = PKLMTest(random_state=42) + np_dataframe = mcar_test_pklm._encode_dataframe(supported_multitypes_dataframe) + n_rows, n_cols = np_dataframe.shape + assert n_rows == 3 + assert n_cols == 7 + + +def test__draw_features_and_target_indexes(np_matrix_with_nan_mcar): + mcar_test_pklm = PKLMTest(random_state=42) + _, p = np_matrix_with_nan_mcar.shape + features_idx, target_idx = mcar_test_pklm._draw_features_and_target_indexes(np_matrix_with_nan_mcar) + assert isinstance(target_idx, np.integer) + assert isinstance(features_idx, list) + assert target_idx not in features_idx + assert 0 <= target_idx <= (p-1) + for feature_index in features_idx: + assert 0 <= feature_index <= (p-1) + + +@pytest.mark.parametrize("dataframe_fixture, features_idx, target_idx, expected", + [ + ("np_matrix_with_nan_mcar", np.array([1, 0]), 2, True), + ("np_matrix_with_nan_mcar", np.array([1, 0, 2]), 3, False) + ] +) +def test__check_draw(request, dataframe_fixture, features_idx, target_idx, expected): + dataframe = request.getfixturevalue(dataframe_fixture) + mcar_test_pklm = PKLMTest() + result = mcar_test_pklm._check_draw(dataframe, features_idx, target_idx) + assert result == expected + + +@pytest.mark.parametrize("matrix_fixture", [("np_matrix_with_nan_mcar")]) +def test__generate_label_feature_combinations(request, matrix_fixture): + X = request.getfixturevalue(matrix_fixture) + _, n_cols = X.shape + mcar_test_pklm = PKLMTest() + result = mcar_test_pklm._generate_label_feature_combinations(X) + # Check that number of projections is smaller than possible + assert len(result) <= mcar_test_pklm._get_max_draw(n_cols) + # Check there are no duplicates + assert len(set([x for x in result if result.count(x) > 1])) == 0 + for features, label in result: + assert isinstance(label, int) + assert isinstance(features, list) + assert label not in features + for feature_index in features: + assert 0 <= feature_index <= (n_cols-1) + + + +@pytest.mark.parametrize("dataframe_fixture, features_idx, target_idx", + [ + ("np_matrix_with_nan_mcar", np.array([1, 0]), 2), + ] +) +def test__build_dataset(request, dataframe_fixture, features_idx, target_idx): + dataframe = request.getfixturevalue(dataframe_fixture) + mcar_test_pklm = PKLMTest() + X, y = mcar_test_pklm._build_dataset(dataframe, features_idx, target_idx) + assert X.shape[0] == len(y) + assert not np.any(np.isnan(X)) + assert not np.any(np.isnan(y)) + assert np.all(np.unique(y) == [0, 1]) + assert X.shape[1] == len(features_idx) + assert len(y.shape) == 1 + + +@pytest.mark.parametrize("dataframe_fixture, permutation_fixture, features_idx, target_idx", + [ + ("np_matrix_with_nan_mcar", "missingness_matrix_mcar_perm", np.array([1, 0]), 2), + ] +) +def test__build_label( + request, + dataframe_fixture, + permutation_fixture, + features_idx, + target_idx +): + dataframe = request.getfixturevalue(dataframe_fixture) + m_perm = request.getfixturevalue(permutation_fixture) + mcar_test_pklm = PKLMTest() + label = mcar_test_pklm._build_label(dataframe, m_perm, features_idx, target_idx) + assert not np.any(np.isnan(label)) + assert len(label.shape) == 1 + assert np.isin(label, [0, 1]).all() + + +@pytest.mark.parametrize( + "oob_fixture, label", + [ + ("oob_probabilities", np.array([1, 1, 1, 1])), + ("oob_probabilities", np.array([0, 0, 0, 0])), + ] +) +def test__U_hat_unique_label(request, oob_fixture, label): + oob_prob = request.getfixturevalue(oob_fixture) + mcar_test_pklm = PKLMTest() + mcar_test_pklm._U_hat(oob_prob, label) + + +@pytest.mark.parametrize( + "oob_fixture, label, expected", + [ + ("oob_probabilities", np.array([1, 0, 0, 0]), 2/3*(np.log(1 - 1e-9) - np.log(1e-9))), + ] +) +def test__U_hat_computation(request, oob_fixture, label, expected): + oob_prob = request.getfixturevalue(oob_fixture) + mcar_test_pklm = PKLMTest() + u_hat = mcar_test_pklm._U_hat(oob_prob, label) + assert round(u_hat, 2) == round(expected, 2) + +@pytest.mark.parametrize( + "list_proj, n_cols", + [ + ( + [ + (np.array([3, 1]), 0), + (np.array([0]), 1), + (np.array([3]), 0), + (np.array([1, 2]), 3), + (np.array([3, 0]), 2), + (np.array([0, 1]), 2) + ], + 4 + ) + ] +) +def test__build_B(list_proj, n_cols): + mcar_test_pklm = PKLMTest() + B = mcar_test_pklm._build_B(list_proj, n_cols) + column_sums = np.sum(B, axis=0) + assert np.all(column_sums == 3) diff --git a/tests/utils/test_input_check.py b/tests/utils/test_input_check.py new file mode 100644 index 00000000..25c903c7 --- /dev/null +++ b/tests/utils/test_input_check.py @@ -0,0 +1,48 @@ +import pandas as pd +import pytest + +from qolmat.utils.exceptions import TypeNotHandled +from qolmat.utils.input_check import check_pd_df_dtypes + + +@pytest.fixture +def multitypes_dataframe() -> pd.DataFrame: + return pd.DataFrame({ + 'int_col': [1, 2, 3], + 'float_col': [1.1, 2.2, 3.3], + 'str_col': ['a', 'b', 'c'], + 'bool_col': [True, False, True], + 'datetime_col': pd.to_datetime(['2021-01-01', '2021-01-02', '2021-01-03']) + }) + + +@pytest.fixture +def supported_multitypes_dataframe() -> pd.DataFrame: + return pd.DataFrame({ + 'int_col': [1, 2, 3], + 'float_col': [1.1, 2.2, 3.3], + 'str_col': ['a', 'b', 'c'], + 'bool_col': [True, False, True] + }) + + +def test__check_pd_df_dtypes_raise_error(multitypes_dataframe): + with pytest.raises(TypeNotHandled): + check_pd_df_dtypes( + multitypes_dataframe, [ + pd.api.types.is_numeric_dtype, + pd.api.types.is_string_dtype, + pd.api.types.is_bool_dtype + ] + ) + + +def test__check_pd_df_dtypes(supported_multitypes_dataframe): + check_pd_df_dtypes( + supported_multitypes_dataframe, + [ + pd.api.types.is_numeric_dtype, + pd.api.types.is_string_dtype, + pd.api.types.is_bool_dtype + ] + )