{ "cells": [ { "cell_type": "markdown", "id": "d68537ba", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "(bart_categorical)=\n", "# Categorical regression\n", "\n", ":::{post} May, 2024\n", ":tags: BART, regression\n", ":category: beginner, reference\n", ":author: Pablo Garay, Osvaldo Martin\n", ":::" ] }, { "cell_type": "markdown", "id": "0cf4f392-fdc7-4175-9e72-c8a334abea84", "metadata": {}, "source": [ "In this example, we will model outcomes with more than two categories. \n", ":::{include} ../extra_installs.md\n", ":::" ] }, { "cell_type": "code", "execution_count": 1, "id": "7c087cca", "metadata": {}, "outputs": [], "source": [ "import os\n", "import warnings\n", "\n", "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pymc as pm\n", "import pymc_bart as pmb\n", "import seaborn as sns\n", "\n", "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" ] }, { "cell_type": "code", "execution_count": 2, "id": "25cf7b45", "metadata": {}, "outputs": [], "source": [ "# set formats\n", "RANDOM_SEED = 8457\n", "az.style.use(\"arviz-darkgrid\")" ] }, { "cell_type": "markdown", "id": "e73740d8-8e70-48b4-b6f9-eb0c1f7ce72f", "metadata": {}, "source": [ "## Hawks dataset \n", "\n", "Here we will use a dataset that contains information about 3 species of hawks (*CH*=Cooper's, *RT*=Red-tailed, *SS*=Sharp-Shinned). This dataset has information for 908 individuals in total, each one containing 16 variables, in addition to the species. To simplify the example, we will use the following 5 covariables: \n", "- *Wing*: Length (in mm) of primary wing feather from tip to wrist it attaches to. \n", "- *Weight*: Body weight (in gr). \n", "- *Culmen*: Length (in mm) of the upper bill from the tip to where it bumps into the fleshy part of the bird. \n", "- *Hallux*: Length (in mm) of the killing talon. \n", "- *Tail*: Measurement (in mm) related to the length of the tail. \n", "\n", "Also we are going to eliminate the NaNs in the dataset. With these we will predict the \"Species\" of hawks, in other words, these are our dependent variables, the classes we want to predict. " ] }, { "cell_type": "code", "execution_count": 3, "id": "71f3a9bc-979f-44fc-8227-133349e4dfb1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | Wing | \n", "Weight | \n", "Culmen | \n", "Hallux | \n", "Tail | \n", "Species | \n", "
---|---|---|---|---|---|---|
0 | \n", "385.0 | \n", "920.0 | \n", "25.7 | \n", "30.1 | \n", "219 | \n", "RT | \n", "
2 | \n", "381.0 | \n", "990.0 | \n", "26.7 | \n", "31.3 | \n", "235 | \n", "RT | \n", "
3 | \n", "265.0 | \n", "470.0 | \n", "18.7 | \n", "23.5 | \n", "220 | \n", "CH | \n", "
4 | \n", "205.0 | \n", "170.0 | \n", "12.5 | \n", "14.3 | \n", "157 | \n", "SS | \n", "
5 | \n", "412.0 | \n", "1090.0 | \n", "28.5 | \n", "32.2 | \n", "230 | \n", "RT | \n", "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 155 seconds.\n", "Sampling: [y]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a9c75e927fb440f6ae09ae520f115714", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with model_hawks:\n", " idata = pm.sample(chains=4, compute_convergence_checks=False, random_seed=123)\n", " pm.sample_posterior_predictive(idata, extend_inferencedata=True)" ] }, { "cell_type": "markdown", "id": "fb2e357e-502e-4ac5-9d53-928437bd2a4e", "metadata": {}, "source": [ "## Results \n", "\n", "### Variable Importance \n", "\n", "It may be that some of the input variables are not informative for classifying by species, so in the interest of parsimony and in reducing the computational cost of model estimation, it is useful to quantify the importance of each variable in the dataset. PyMC-BART provides the function {func}`~pymc_bart.plot_variable_importance()`, which generates a plot that shows on his x-axis the number of covariables and on the y-axis the R$^2$ (the square of the Pearson correlation coefficient) between the predictions made for the full model (all variables included) and the restricted models, those with only a subset of the variables. The error bars represent the 94 % HDI from the posterior predictive distribution. " ] }, { "cell_type": "code", "execution_count": 10, "id": "a9d1c616-8c1f-4907-ad5a-adffb290c0c2", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAE3CAYAAACq3N6VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9JUlEQVR4nO3de1yUZcL/8e/AAIKDiorWGolFS5ZhHvDIRpo+5anWSvNnmlt2UNO0Hmtr17ae1q3snLrt5lp2EDfNbD1vbplbuAHmsVozMzmWiIAoicDM3L8/xhlBZgBhkBv8vF8vXsJ9X3PNNXpxO9+5DrfFMAxDAAAAAGAyAY3dAAAAAADwhrACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJSsjd2AmqxevVrbt2/X119/re+++07l5eV65plndPPNN59VPU6nU8uWLdPy5cuVkZGhsLAw9e3bVw8++KCio6MbpvEAAAAA6sz0YeXVV19VTk6OIiIi1KFDB+Xk5NSpnieeeEIrVqxQTEyMJkyYoPz8fG3YsEFbt27Ve++9p5iYGD+3HAAAAEB9mH4a2Ny5c7V582alpKRo3LhxdaojJSVFK1asUO/evfXhhx/qkUce0bx587Ro0SIVFxfrySef9G+jAQAAANSb6cPKgAED1KlTp3rV8f7770uSZs2apeDgYM/x/v37KyEhQdu2bdPBgwfr9RwAAAAA/Mv0YcUfUlNTFRYWpp49e1Y5l5CQIEnatm3buW4WAAAAgGo0+7By4sQJ5eXl6aKLLlJgYGCV8+7F9enp6ee2YQAAAACq1ezDyvHjxyVJNpvN63n38eLi4mrrcTqd/m0YAAAAgGqZfjcwsygqKmrsJpheRESECgsLG7sZaMLoQ/AH+hHqiz4Ef6Af1SwiIqLGMs1+ZCU8PFyS75ET93FfIy8AAAAAGkezDythYWGKjIxUdna2HA5HlfPutSrcGBIAAAAwl2YfViSpT58+OnHihHbs2FHlXHJysiQpPj7+XDcLAAAAQDWaVVgpKCjQgQMHVFBQUOn42LFjJUmvvPKKysrKPMe/+OILJScnKz4+Xl26dDmnbQUAAABQPdMvsH///fe1fft2SdJ3333nOZaWliZJGjJkiIYMGSJJSkpK0sKFCzV9+nTNmDHDU0e/fv00ZswYvf/++xo9erQSExOVn5+vDRs2yGazcQd7AAAA+EVJiaGhwwxJ+frXRotCQy2N3aQmzfRhZfv27frwww8rHduxY4dnSlenTp08YaU6Tz31lGJjY7V8+XK9++67CgsL06BBg/Tggw8yqgIAAHiTCZiQxTAMo7Eb0RSw9VzN2KIP9UUfgj/Qj1BXp8OKCCuoM/pR7bF1MQAAAIAmi7ACAGgWSkoMJVzr1JXd81VSwqQBAGgOCCsAAAAATImwAqDR8Yk4AADwhrACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibACAAAAwJQIKwAAAABMibCCeispMZRwrVNXds9XSYnR2M0BAABAM0FYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApmRt7AbUxp49e7RgwQLt2rVL5eXliomJ0aRJkzRq1Kha13Hs2DEtWbJEH3/8sbKzsxUcHKyLLrpIo0eP1pgxYxQSEtKArwAAAADA2TJ9WElNTdXkyZMVFBSkESNGKDw8XJs2bdLs2bOVk5OjKVOm1FjHsWPHdPPNNysrK0u9evXSuHHjVFZWps8++0x//OMf9a9//UtLlixRQAADTQAAAIBZmDqs2O12zZkzRxaLRUlJSbriiiskSffff7/GjRunBQsW6IYbblB0dHS19SxfvlxZWVn6zW9+o8cee8xzvKysTOPHj1dKSoq2b9+u+Pj4hnw5AAAAAM6CqYcSUlJSlJmZqZEjR3qCiiTZbDZNmzZNdrtdq1atqrGerKwsSVJiYmKl48HBwRo4cKAkKT8/348tBwAAAFBfph5ZSUtLkyQlJCRUOecOGe4y1bnsssskSZ9//rkGDBjgOV5eXq7//Oc/atGihXr06OGPJgMAgCaopMRQVrbh+Tl5qyFbuBQSLAUHS0FBrj+Dg08fcx8PDLQ0YsuB5s3UYSU9PV2S1Llz5yrnWrdurYiICGVkZNRYz5gxY7R69Wq9+eab+vrrr9WtWzeVl5fr888/V1FRkV588UV17NjR380HAAAmcOKEobw86fCpr7w8Kfew65jre6m4uPJj/m+uJBneqqsiMNDwGmKCg6sGnKDgysdd5yxVyrsf47OeM45bLAQmNE+mDivFp64c4eHhXs/bbDYdOnSoxnpatGihd999V3/4wx+0Zs0az2hMQECAbr/9dvXs2bPGOlq3bs0CfB9CQgxJBZKkNm0iFBbGBRNnhz4Ef6AfnZ9+/tnQT4ccys11KjfXqUOnvnJznTp0yKncw04dP1670NGypUU//+wq2z3OKrvdUGmpVFpmqKxUKis3VFpqqKxMsttPP87hkEpKXF91U7v2VSc42FBIiEXBQVJwiMUVakIsCgm2KCRECg6yKDhECjnjnCvsuMqEBFtOPzb4VHn344MtrvqD3eXc51zHrFYCkxvXIv8ydVjxl4KCAk2bNk0FBQVatGiRevbsqdLSUm3evFnPPvustmzZog8++ECtW7f2WUdRUdE5bHHTUlJy+iJ79GihSkv5pcTZoQ/BH+hHzc/PPxvKPawKIyBVR0h+/rl2ddlaSpGRrq+OHaTISIs6REodOriOdYiUAgKkocNc5V+Y51Bo6Jl9yHLqS7LbDZWXS2VlUnm5VFrm+r7SV7mXY54vw/txX4/xctyokHHcdbrUP/ycLYuldiNDni+vxy3ey/oof2b9ZpmOx7Wo9iIiImosY+qwYrPZJEnHjx/3er64uNjnqEtFzz77rHbu3KnVq1fr8ssvl+QarRk7dqwcDoeefPJJvf3223rggQd81lFeXu7znMVikdVqbdSykhQUFFSnsna7XYbh+8JWU9nyckOBge5jwX6rtyKr1er5xKahyjocDjmdTr+UDQwM9IzEmaGs0+mUw+HwWTYgIECBgYGNVra83FCAxSKn4SprGIbsFT+2rKbemspW/D1qqLKu18A1ojHLVvx7CAhwqLzckNXq/Q0C14iqZc/1771hSD+fkI4cCVB+foAniOQediovz+IKJ0ekEydq9ybPZjMU2d4dRgx1iLSoY8cAdewgtW9vqG2EQ2Fh3h9b8ffoxAmnAgNdvxvl5ZYqfahi2cBAyTDssvp4J9XQv/eua5R04oTda6hxBSmLHI5Az7mSEofKyqVybwGoXLLbAzw/l540XOfd9Z0RlMrLLSqr0EzDkGsUqtRn02uhfiErMNDpfRpekDwjQMHBUpDVUFCwUSkEVSwf2iJAwafKBwY4XGXPDEunfg4NC/SMUAUGOmUYzkrviyr2I64RVcvWhqnDintL4oyMDHXr1q3SuaKiIhUWFtZqYfy///1vtWnTxhNUKurXr58k6Ztvvqm2jr/85S/VtvPGG2/0/Py3v/3N55ucTp066ZZbbvH8vGTJEp08edJr2Q4dOmjcuHGen5cuXeozuLVt21YTJkzw/Lx8+XIVFBR4LRseHq4777zT8/PKlSt1+PBhr2VbtGihe++91/Pz6tWrlZOTU6XcLTdJdrtV0lTPsQ0bNnjWHXlTMRxu2rRJ33//vc+yU6dO9bxx+fTTT7V3716fZe+++26Fnfpf6fPPP9dXX33ls+xvfvMbtWrVSpL0xRdfaMeOHT7L3n777WrXrp0kadu2bdVu7nDbbbd51kHt2rVLW7du9Vn25ptv1kUXXSRJ+vrrr/Xvf//bZ9lRo0apS5cukqR9+/bp448/9ll22LBhns0lDhw4oI0bN/osO2TIEM+OexkZGVq7dq3PsomJierevbsk6ccff6x2R76BAweqV69ekqS8vDwtX77cZ9mul8frm719JblGQ5OSknyW7dmzp2fjjePHj+utt97yWfaqq67SoEGDJEklJSVavHix7zZ07aqhQ4dKcr2Jre73PiYmRsOHD/f8zDXCxdc1QnL9Rz1t2jTPz/6+Rrj/S+vV41O99da3PstyjXBpqGvENdck6pJL45SXJ/13b74+/niPSkpa6kSJTSUlNs+fdrv7w62Kb1CrTrcODjqp0NBiXXihVb/8ZRt1iLQoNLRYe3b/S6GhxQoN/VlB1spv7mMvO32NOHas9teIkydP6pabXNcIbw8x0zXCYnGtc1m+vPbXiCVL3vV6jQgKlDpGVr5GLF26tMZrhNPpGl16771VOpRbKKcjUA5n4Kk/rXI4AhUYGKohQ4Z7AlFKyg7l5x+X0xkoh8NVxul0Pc4wgnTppV09gejHH/P088/lrnLueh3WCuWDVfE9scNhqeV0vNMjZL65+2VtlgC4y1oUEOAKOFarQ9Gd9+qtt77wlOIa4VLxfURtmDqsxMfH6/XXX1dycrJGjBhR6Zz7H61Pnz411lNWVub5Cg4OrnTO/Yt45nEAAFCZYRgqLg7Q0aJ2ruBxwlYlhPxjXasK05HaSRrks76wMLsuvNCqDpFSuK1EP/64W6GhxQoLLa4SRPr06eP5gDE/v0y5P2U38KtFTQICTq11CSlXaIsTXsu0aNFCv0o4HQxOFKfX8IHGlZ6f16xJrfEDDfd0vPUbPtaB7zPlcAa6AtCpgOMOQ0OHjpDDaVVZmbR7917l5OR5ApPzVDmH0yqnI1Cdo38pp9M1InX48FEdO3by1PnTQcwTzJxWGcbp1+d0up5PdunIkV+c5d8ovLEY1Y1/NzK73a4bbrhBubm5WrFihbp27SrJNf1r3LhxOnjwoNatW+dJiAUFBSosLFRERITatm3rqWfy5MlKTk7W1KlTNWvWLM/xsrIyTZ06VcnJyXr88ccrfaJwJl+fKkpM8SgpMTRqtOvYP9cFe+b4Mg2s+Q/f+qtsSYmhkTe5poH9a6NFLVqIaWB+Lis1/2lgJ09KQ4cZCghwaN0/DC/rDU6X5RpR9Rphtzt07Lg8U7DyDltcf+ZJeUcsrq88yceH+FW0buWaltW+vaEOka5pWu0j5fo+UmrfTmrZsm7Xk4aaKnrihFPDRrl+N9Z+aKnSh7hGeC/bVK4RDfE+IjAwUA6HReXlUkmJUydPOnXsuKH7prkes+4fp/tRU79GNMR7gya/ZsVqtWru3Lm6++67NX78eI0cOVI2m02bNm1Sdna2Zs2a5QkqkpSUlKSFCxdq+vTpmjFjhuf47NmztXPnTv3lL3/R1q1bPQvsk5OTlZWVpSuvvFJjxoypti0VfyFq0tTKVryQ1aWs3W7I4aj6i1zfes912cDAwFrPoWxqZQMCAmq9m11jlLXbDTkr/Gfgmt5Quz5shrKSOX6XzXqNOHdlXX3I6QxUUJBFQUE1r3c4X64RhmGoqEiVFqcfzjOUl+fU4cPuxeqBtV5v0Kb16cXqHTpIHSItFRauS5HtpRYtzm5R8dlcTxry997hcJWtqQ9xjTit6VwjGqZsQIBrzUtYWKCkQLUtMeR0uq5HvvqR2a4RNWmo9wa1YeqwIrnWlCxbtkzz58/Xxo0bVV5erpiYGM2cObPS/M7qdO3aVatWrdLrr7+ulJQUJSUlKTAwUBdffLFmzJihyZMnKyQkpIFfCQAA/mcYho66g4gneBgVvnf9WVZWu/ratNEZu2S5ds1y75gVGelarAwA54Lpw4okxcXFVbso1m3GjBmVRlQqio6O1jPPPOPvpgEA0GCcTteIyOFqgkhenirtylSdiIhTW/dGukdGLJ7RkA6RUvv2BBEA5tIkwgoAAM2N02no6NHTQSTv1NSsiqMheXmuLWhro21ExWlZp+4j0uH0aEhke9d9LACgKSGsAADgZ06nocLCyjcvPHy4chA5cqR2QcRiqTmItG9HEAHQPBFWAAA4Cw6HocKjFUdDzggih6Uj+VI1G0V5WCxS27bVrxFp31612iwAAJojwgoAAKc4HK4RkdwzgkjFXbSOHJGq2ZXTw2KR2rWreY3ImXdJBwCcRlgBAJwXHA5D+fmVd8hybd/r+jn3sJR/RHL4vu2AR0CA1K5tzVOzCCIAUD+EFQBAk2MYhsrKXDcoPFkqnSyRjhadvlfPylWGjh41Ki1ez8+vfRBp3676INKuLUEEAM4FwgoAwO/cYaK0VCo56QoTJ0td4aKkpHLIKDl5qlyJ4eO49/LV3HhZr//N+/HAAKld+wo3L/SyRqQtQQQATIOwAgDnoZrChPfjxlmVry5M+JPVKrVoIYWEuEZPJOnaa6QLLpA6dDgdRDp2cN1nJDCQIAIATQVhBQBMqDmGCfdXaIvKP1c8HtJCCm1hqbF8xWPuUZCSEkNDh7mmgv3+MYtCQwklANDUEVYAoA4Mw1B5+akwcFIqPfXnyQpfZ/588qTh47j3eppbmAAA4GwRVgA0Sw0dJkpP1m6xtj8QJgAA5yvCCgBT+fprQ07DV6AwqgkaVUMFYQIAgKaNsALgnDEMQwUFUnqG+8tQerrre7cHH5Ykw0cNdUOYAACgaSKsAPA7w3Dd38IdRNLTDU9AOX68+sd26iSFhRImAAAAYQVAPTgchg4dkg5mSBlnhJKSEu+PCQiQfnGhFB0tdb5Yio626MILDE2f6Tr/1mJ2cQIAAC6EFQA1stsN5eS4QsjBdNf0rYwMKSNTKivz/pjAQCnqIim6syuYRHe2qHNn6eIoKSSkchhxBRv/Tv0CAABNH2EFgEdpqaGsbNf0rYxM41QwkbKzJbvd+2OCg6SLLz4dSNzh5KJOTL8CAAD1Q1gBzkMlJYYyMquuJ/nxR9/39ghtcWrqVudToSRair5YuvBC7ggOAAAaBmEFaMaOHz8VStJP7byV4fr+UK7vx9hsUpdo9/Qty6l1JVKHSCkggFACAADOHcIK0AwUHjVOLXB3hZKD6a71JEeO+H5MRIQrkHTuLHVxj5R0ltq2lSwWQgkAAGh8hBWgiTAMQ/n58qwjyagwUnK0yPfjItvLE0Sio11rSjpfLLVpQyABAADmRlgBTMbpNJSbW/XGiRkZUvHPvh934YVSl1MjJRVDic1GKAEAAE0TYQVoJA6HoR9/qnDjRHcoyZROnvT+mMAA6RedTm8H3KXCdsDcmwQAADQ3hBWggZWXu7YDzsiovPtWVpZUVu79MVarK4BUvHFil87SRRdJwcGEEgAAcH4grAB+UlpqKDPTdTf39PTTu3BlZ0sOH9sBh4ScCiMVpm5Fd5Z+8QvuUQIAAEBYAc7SiROnF7anZxqeaVw//SQZPm7CHhZW+U7u7u8v6Mh2wAAAAL4QVgAfjh07fbPEijdOPHzY92NatXLdo8R948Qu0a6Rk8hItgMGAAA4W4QVnNcMw1BhoSuEHEyvvB1wQaHvx7VrW3k74M4Xu0JKmzaEEgAAAH8hrOC8YBiGDufJc+PEgxV23jp2zPfjOnaUZx1J9KkbJ3buLLUKJ5AAAAA0NMIKmhWn09BPh05vB5yRYehghiuknDjh/TEWi/SLC73fODEsjFACAADQWAgraJLsdkM5OVJ6pjuYuEZKMrOk0lLvjwkMlKIucgWSijdOvDhKCgkhlAAAAJgNYQWmVlbmukeJJ5CcmsaVlS3Z7d4fExwkRZ3aDriLe5Sks3RRJykoiFACAADQVBBWYAolJa77kmRUuJP7wQzpxx8lp497lIS2cI+QVNgOuLN04YVSYCChBAAAoKkjrOCcKi6uvB2w+8aJPx3y/RhbS6lLl9N3cnffo6RDJPcoAQAAaM4IK2gQR49WCCUZp2+ceOSI78e0aXP6xoldOls8oybt2rIdMAAAwPmIsIJ6czpP37b9fx8xlJll6OhR3+Uj28uzBXDFGye2aUMgAQAAwGmEFdTbd/tPf7/nq9PfX3hB1RsnRneWbDZCCQAAAGrWJMLKnj17tGDBAu3atUvl5eWKiYnRpEmTNGrUqLOqp7i4WG+++aY2bdqkrKwsBQUFKSoqStddd52mT5/eQK1v/mIuPf39bx+WfnmZRRdHSaGhhBIAAADUnenDSmpqqiZPnqygoCCNGDFC4eHh2rRpk2bPnq2cnBxNmTKlVvX8+OOPmjRpkrKysjRgwAAlJiaqrKxMmZmZ+uijjwgr9WC1WiS5poINGWwhpAAAAMAvTB1W7Ha75syZI4vFoqSkJF1xxRWSpPvvv1/jxo3TggULdMMNNyg6OrraehwOhx544AEdPnxYb731lvr161fleQAAAACYS0BjN6A6KSkpyszM1MiRIz1BRZJsNpumTZsmu92uVatW1VjPRx99pK+++kp33XVXlaAiSVarqTMbAAAAcF4y9bv0tLQ0SVJCQkKVcwMHDqxUpjobNmyQJN1www366aeftGXLFh0/flxRUVG65ppr1LJlSz+2GgAAAIA/mDqspKenS5I6d+5c5Vzr1q0VERGhjIyMGuv5+uuvJUnbt2/XM888o7KyMs+5tm3b6pVXXlHfvn3902gAAAAAfmHqsFJcXCxJCg8P93reZrPp0KFqbn1+Sn5+viRp7ty5uuuuuzRhwgQFBwdr/fr1mjdvnu6//35t2LBBHTp08FlH69atFRBg6llzjSYkxJBUIElq0yZCYWEssMfZoQ/BH+hHqC/6EPyBfuRfpg4r/mIYrp2qrr32Ws2ePdtzfOLEicrNzdXf/vY3rVy5UtOmTfNZR1FRUYO3s6kqKTl9U8ijRwtVWsovJc4OfQj+QD9CfdGH4A/0o9qLiIiosYyphwpsNpsk6fjx417PFxcX+xx18VbP4MGDq5wbNGiQpNNTxQAAAACYg6nDintLYm/rUoqKilRYWOh1PcuZunTpIklq1apVlXPuY6WlpfVoKQAAAAB/M3VYiY+PlyQlJydXObd161ZJUp8+fWqsx71d8ffff1/lnPtYp06d6txOAAAAAP5n6rDSv39/RUVFad26ddq7d6/neHFxsV577TVZrVaNHj3ac7ygoEAHDhxQQUFBpXpuvvlmBQcHa+nSpcrNza1Uz+uvvy5JGjZsWAO/GgAAAABnw9QL7K1Wq+bOnau7775b48eP18iRI2Wz2bRp0yZlZ2dr1qxZnilekpSUlKSFCxdq+vTpmjFjhud4VFSUHnnkEc2dO1c33nijhg4dquDgYG3ZskU5OTm67bbb1L9//8Z4iQAAAAB8OOuwEhsbW+lni8Wili1b6tJLL9Xw4cN1++23KygoyG8N7Nevn5YtW6b58+dr48aNKi8vV0xMjGbOnKkbb7yx1vVMnDhRnTp10htvvKH169fL4XAoJiZGU6ZM0dixY/3WXgAAAAD+UeeRFff0K4fDoZycHO3cuVO7d+/Wli1btHjxYlmtp6teu3atlixZoqKiItntdnXs2FH333+/EhMTa/VccXFxWrx4cY3lZsyYUWlE5UyDBw/2uiMYAAAAAPOpc1h59tlnK/28e/duTZw4UV988YXWr1+vm266yXMuMjJSL7zwgi655BI5nU4tWLBA999/v9auXVtpGhcAAAAAuPltgX337t09oy1n7t7Vr18/XXLJJa4nDAhQYmKiysvLlZ6e7q+nBwAAANDM+HU3sMsuu0ySquzGVVF5ebmee+45RUdHa8CAAf58egAAAADNiF93A/v5558lSW3btvV63m6366GHHtIPP/ygpKQkhYSE+PPpAQAAADQjfg0rn3/+uSTpV7/6VZVzJSUlmjVrlr799lu98847uvTSS/351AAAAACamXqHFafTqezsbL3xxhvatm2bBg8erOHDh1cqU1BQoKlTp6qsrEwrVqxQx44d6/u0AAAAAJq5OoeVM++3Ikm33nqr/vjHPyogoPJSmJdeekm7du1Sx44dNWHCBM/xiRMn6o477qhrEwAAAAA0Y/W+z0ppaan27t2rgwcPauXKlbr66qs1ZsyYSmXnzp2ruXPn1q+lAAAAAM4rfrvPyt/+9je98MILmjt3rgYMGKBOnTrVu3EAAAAAzl9+27r4nnvuUUJCgk6ePKmFCxf6q1oAAAAA5ym/3mdl9uzZslgsWrNmjXJycvxZNQAAAIDzjF/DSteuXXXdddfJbrdr8eLF/qwaAAAAwHnGr2FFkmbMmCGLxaIPPvhAeXl5/q4eAAAAwHnC72Hl8ssv19ChQ1VaWqolS5b4u3oAAAAA54mz3g1s3759NZZZsGBBnRoDAAAAAG5+H1kBAAAAAH8grAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwpSYRVvbs2aN77rlH8fHxuvrqq3Xrrbdq7dq1da6vvLxcN910k2JjY3XDDTf4saUAAAAA/MXa2A2oSWpqqiZPnqygoCCNGDFC4eHh2rRpk2bPnq2cnBxNmTLlrOt87bXXlJmZ2QCtBQAAAOAvph5ZsdvtmjNnjiwWi5KSkjR37lz99re/1erVq3XZZZdpwYIFSk9PP6s6v/nmGy1atEgPPfRQwzQaAAA0SaGhFiVvCdA3u9spNNTS2M0BIJOHlZSUFGVmZmrkyJG64oorPMdtNpumTZsmu92uVatW1bq+srIyPfroo+revbsmTJjQEE0GAAAA4CemngaWlpYmSUpISKhybuDAgZXK1MbChQuVkZGh1atXy2LhExMAAADAzEw9suKe4tW5c+cq51q3bq2IiAhlZGTUqq49e/Zo8eLFmjFjhrp06eLPZgIAAABoAKYeWSkuLpYkhYeHez1vs9l06NChGuspKyvTY489pq5du+quu+6qU1tat26tgABTZ7tGExJiSCqQJLVpE6GwMEatcHboQ/AH+hH8KSIiorGbgCaKa5F/mTqs+Msrr7yijIwMffDBBwoMDKxTHUVFRX5uVfNRUmJ4vj96tFClpfxS4uzQh+AP9CP4S0REhAoLCxu7GWiiuBbVXm0+FDD1UIHNZpMkHT9+3Ov54uJin6Mubt98843eeustTZkyRbGxsX5vIwAAAICGYeqwEh0dLUle16UUFRWpsLDQ63qWivbt2yeHw6EFCxYoNja20pckHTx4ULGxserdu7ff2w8AAACg7kw9DSw+Pl6vv/66kpOTNWLEiErntm7dKknq06dPtXVER0fr1ltv9Xpu5cqVCg8P1/XXX6/Q0FD/NBoAAACAX5g6rPTv319RUVFat26d7rjjDnXt2lWSa/rXa6+9JqvVqtGjR3vKFxQUqLCwUBEREWrbtq0kqWfPnurZs6fX+leuXKn27dvrT3/6U8O/GAAAAABnxdTTwKxWq+bOnSvDMDR+/Hg9/vjjmjdvnm666Sbt379f06dPr7QNcVJSkoYPH66kpKRGbDUAAAAAfzD1yIok9evXT8uWLdP8+fO1ceNGlZeXKyYmRjNnztSNN97Y2M0DAAAA0EBMH1YkKS4uTosXL66x3IwZMzRjxoxa17tv3776NAsAAABAAzL1NDAAAAAA5y/CCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMKUmcZ8VAAAAoCkIDbUoeYtFERERKiwsbOzmNHmMrAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFPippAAgGaBG7EBQPNDWAHQ6HiTCQAAvGEaGAAAAABTIqwAAAAAMCXCCgAAAABTYs0K6o31BgAAAGgIjKwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCVrYzegNvbs2aMFCxZo165dKi8vV0xMjCZNmqRRo0bV6vFffvmlPv74Y6WlpSknJ0cnTpxQp06ddN111+m+++5Tq1atGvgVAAAAADhbpg8rqampmjx5soKCgjRixAiFh4dr06ZNmj17tnJycjRlypQa65g5c6YKCwvVq1cv3XTTTbJYLEpLS9PixYu1adMmvffee2rXrt05eDUAAAAAastiGIbR2I3wxW63a9iwYTp06JCWL1+uK664QpJUXFyscePG6eDBg1q/fr2io6OrrWfRokX69a9/rQ4dOniOGYah//u//9Pf//53jR8/Xk888US1dRQWFtb79TR3ERER/D2hXuhD8Af6EeqLPgR/oB/VLCIiosYypl6zkpKSoszMTI0cOdITVCTJZrNp2rRpstvtWrVqVY313HvvvZWCiiRZLBZNmzZNkrRt2zb/NhwAAABAvZk6rKSlpUmSEhISqpwbOHBgpTJ1YbW6ZsEFBgbWuQ4AAAAADcPUYSU9PV2S1Llz5yrnWrdurYiICGVkZNS5/g8++EDS6eADAAAAwDxMvcC+uLhYkhQeHu71vM1m06FDh+pU9969e/XnP/9Z7dq10913311j+datWysgwNTZzhRqM/cQqA59CP5AP0J90YfgD/Sj+jN1WGkoWVlZuu++++RwOPTSSy+pbdu2NT6mqKjoHLSsaWMhGeqLPgR/oB+hvuhD8Af6Uc1qE+ZMHVZsNpsk6fjx417PFxcX+xx18SUnJ0eTJk1SQUGBFixYoH79+tW7nQAAAAD8z9TzmtxbEntbl1JUVKTCwkKv61l8yc7O1sSJE3X48GG98sorGjRokL+aCgAAAMDPTB1W4uPjJUnJyclVzm3dulWS1KdPn1rVlZ2drTvuuEOHDx/Wyy+/rCFDhvivoQAAAAD8ztRhpX///oqKitK6deu0d+9ez/Hi4mK99tprslqtGj16tOd4QUGBDhw4oIKCgkr1uINKbm6uXnrpJQ0dOvScvQYAAAAAdWPqNStWq1Vz587V3XffrfHjx2vkyJGy2WzatGmTsrOzNWvWLHXp0sVTPikpSQsXLtT06dM1Y8YMz/E77rhDOTk5uvrqq7Vv3z7t27evynNVLA8AAACg8Zk6rEhSv379tGzZMs2fP18bN25UeXm5YmJiNHPmTN144421qiMnJ0eStGvXLu3atctrGcIKAAAAYC4WwzCMxm5EU8DWczVjiz7UF30I/kA/Qn3Rh+AP9KOa1WbrYlOvWQEAAABw/iKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAU2oSYWXPnj265557FB8fr6uvvlq33nqr1q5de1Z1OJ1OLV26VKNGjVJcXJz69eunmTNnKj09vWEaDQAAAKBerI3dgJqkpqZq8uTJCgoK0ogRIxQeHq5NmzZp9uzZysnJ0ZQpU2pVzxNPPKEVK1YoJiZGEyZMUH5+vjZs2KCtW7fqvffeU0xMTAO/EgAAAABnw2IYhtHYjfDFbrdr2LBhOnTokJYvX64rrrhCklRcXKxx48bp4MGDWr9+vaKjo6utJyUlRZMmTVLv3r21ZMkSBQcHS5K++OIL3Xnnnerdu7eWLl1abR2FhYV+eU3NWUREBH9PqBf6EPyBfoT6og/BH+hHNYuIiKixjKmngaWkpCgzM1MjR470BBVJstlsmjZtmux2u1atWlVjPe+//74kadasWZ6gIkn9+/dXQkKCtm3bpoMHD/r/BQAAAACoM1OHlbS0NElSQkJClXMDBw6sVKY6qampCgsLU8+ePaucc9e9bdu2+jQVAAAAgJ+ZOqy4F7937ty5yrnWrVsrIiJCGRkZ1dZx4sQJ5eXl6aKLLlJgYGCV8+4pZCy0BwAAAMzF1Avsi4uLJUnh4eFez9tsNh06dKjaOo4fP+4p66uOis/lS23m1IG/J9QffQj+QD9CfdGH4A/0o/oz9cgKAAAAgPOXqcOKe9TDPTpypuLiYp+jLm7u875GTtzHfY28AAAAAGgcpg4r7vUk3talFBUVqbCw0Ot6lorCwsIUGRmp7OxsORyOKufda1Vq2v4YAAAAwLll6rASHx8vSUpOTq5ybuvWrZKkPn361FhPnz59dOLECe3YsaPKOXfd7ucCAAAAYA6mDiv9+/dXVFSU1q1bp71793qOFxcX67XXXpPVatXo0aM9xwsKCnTgwAEVFBRUqmfs2LGSpFdeeUVlZWWe41988YWSk5MVHx+vLl26NPCrAQAAAHA2TB1WrFar5s6dK8MwNH78eD3++OOaN2+ebrrpJu3fv1/Tp0+vFDKSkpI0fPhwJSUlVaqnX79+GjNmjL788kuNHj1azz33nH7729/q3nvvlc1m05NPPnmOX1nzMXHiRMXGxlY6lpqaqtjYWC1YsKCRWoWmLjY2VhMnTmzsZqAZom/BX7z9/1cX9ElUx9d7Kn/1v6bA1GFFcgWNZcuWqVevXtq4caOWLVumNm3a6Pnnn9fUqVNrXc9TTz2lOXPmyGKx6N1339WWLVs0aNAgvf/++4qJiWnAV9A4srOzFRsbq8mTJ/sss2vXLsXGxurRRx89hy1Dc5Cbm6sXX3xRo0ePVu/evdWtWzclJCTo3nvv1apVqyqNYAJ1VVZWpquvvlq9e/f2uuZw3bp1io2NVdeuXXX06NEq57/88kvFxsbqnnvuOQetxbl2vvaP8+lNalPi/nep7VdqampjN7nJMPV9Vtzi4uK0ePHiGsvNmDFDM2bM8HouICBAEydO5NMLoJ7WrVun3//+9zp58qSuvPJK3XjjjQoPD1deXp5SUlL02GOPafXq1Xr77bcbu6lo4oKDg9WzZ09t3bpV//3vf3XVVVdVOp+WliaLxSKn06kvv/xSQ4YMqXTe/Wagb9++nmMbNmxQaGhowzceDa4h+sfZmDdvnkpKSurWeDQ7o0ePrrKO+sMPP1ROTo7uuOMOtWrVqtK5Tp061areuLg4bdiw4by+X0uTCCsA6m/ixInKycnR5s2b61zHZ599pocfflitWrXSa6+9poEDB1Y6bxiGPv74Y73//vv1bS6aifr2u759+2rr1q1KTU2t8mY0NTVVAwYM0K5du5SSklLlzWhaWpqnDrdLL720Tu1AwzBb/zgbv/jFL+r0OJhXffrjzTffXOVYWlqacnJyNGnSJF100UV1alNoaOh5f90y/TQwnFtff/21nnrqKY0cOVK9evVSXFycRo0apUWLFqm8vLxedVc3L3fw4MEaPHiw5+cffvhBPXr00LXXXquioqJKZQ8cOKDu3btr8ODBPu/BA/9zOBx66qmn5HQ69corr1QJKpJksVg0dOhQLVy40HNswYIFPoe8V61apdjYWK1atarG53/00UcVGxurrKwsvfHGG7r++usVFxen4cOHa/369ZKk8vJyvfrqqxo8eLCuuuoqjRo1Sp9//rnX+oqLizV//nyNGDFCcXFx6t27tyZPnqwvv/yySln38L7dbtef//xnDR48WN26ddP1119fZY0c/Mv9RtL9xtLt8OHDSk9PV//+/dWzZ88q58vKyrRr1y6Fh4friiuu8Bz3dh1y963s7GwlJSVp2LBhuuqqqzRo0CAtXLhQTqezSrtKSkr03HPPKTExUVdddZVGjhypFStWsGbvHPN3/8jPz9fTTz+toUOHqlu3burbt69mzJih7777rspz+5qOVde+UVBQoMcee0z9+/dXXFycxo4dW+W6GRsb63ktFacUMZ276SgrK9O7776ryZMnKzExUd26dVP//v01ffp0/fe//61SnmsKIys4w4oVK/Tpp58qPj5e11xzjU6ePKm0tDS9+OKL+uqrr87ZL8sll1yi3/3ud5ozZ44ef/xxzZ8/X5Lrl/yhhx5SeXm5nn/++RpvCgr/SU1NVVZWlnr06KH+/ftXWzY4OLjB2vHMM89oz549GjRokAICArRhwwb97//+r1q1aqWkpCTt379fiYmJKi0t1bp16zR16lRt3LhRUVFRnjqOHj2qCRMmaP/+/erdu7cSEhJ0/PhxffLJJ5o0aZJeffXVKp/CStJDDz2kPXv26JprrlFAQIA2btyop556SkFBQZ5dB+Ff3bp1U8uWLfXll1/K4XAoMDBQ0ukpPH369JFhGHrppZdUWFjomSqxe/dunTx5Uv379/c8pibPPfec0tLSNGjQIA0cOFCffPKJFixYoPLycj344IOecg6HQ/fdd59SU1N1+eWXa+TIkSoqKtKzzz5bq+304T/+7B+ZmZmaOHGicnNzNXDgQA0ZMkT5+fnatGmTkpOT9dZbb6l79+7VtqeufePYsWP6f//v/6lly5YaNWqU8vPztXHjRk2ePFmrVq3SL3/5S0nS9OnTPVOLpk+f7nl8165d6/6XiHOqqKhITz/9tHr37q3ExES1atVKWVlZ2rx5sz777DMtXbpUcXFxjd1MUyGsNHOZmZk+A8ahQ4eqHLvvvvv0xBNPVPrP3TAM/f73v9cHH3yg7du3q1evXg3W3orGjBmjzz//XB999JFWrFihsWPH6oUXXtC3336r6dOnn7N2wMV9n6J+/fo1ajsOHDigNWvWqG3btpJcQ+9jx47VQw89pMsuu0xr165VWFiYJCkhIUEPPvig3n77bc2ZM8dTxx//+Eft379fTz/9tG655RbP8SNHjujWW2/V448/rl/96lcKCQmp9NyHDh3SunXrZLPZJEl33HGHRo0apTfffJOw0kCsVqt69eqlzz77TN98843nP/G0tDSFhYXpyiuvlNPplGEY2rZtm/7nf/5H0uk3q2fTX7/55hutWbNGHTp0kCRNmzZN119/vd59913df//9nhC+atUqpaamavDgwfrzn/+sgADXJIU777yz0nb6aHj+7B+PPPKIjhw5ojfeeKPSyPHUqVN1yy23aM6cOVq7dm217alr3/j22289u566H9OvXz/NmTNHS5cu1VNPPSXJtTbXPbXI1xpdmFvr1q21ZcsWdezYsdLx/fv3a+zYsXr55Ze1ZMmSRmqdOTENrJnLzMzUwoULvX6tXLmySvlOnTpV+RTSYrHo9ttvl+S6N825NHfuXF144YV6+umn9c477+idd95Rjx49NG3atHPaDkh5eXmSpAsuuKBR2zFlyhRPUJGk7t27KyoqSseOHdODDz7oCSqSdP311ysoKEj79u3zHCsoKNDGjRvVv3//SkFFktq3b6/JkyeroKBA//nPf6o890MPPeQJKpJrBLBnz546ePCgiouL/fkyUYF7qk9KSornWFpamnr27Cmr1apu3bopNDS00pSZuqxHmDZtmieoSFLbtm113XXX6eeff9bBgwc9x9esWSNJmjlzpueNpeRaD/PrX//67F4c6s0f/eO///2vdu7cqV//+tdVprh26dJFY8eO1Xfffed1OlhFde0bYWFhmj17dqXHjB49WlarVV9//XW1z4mmJTg4uEpQkaTLLrtMffv21bZt2+o97b65YWSlmUtISNAbb7zh9dyuXbt02223VTpWVlampKQkrV+/Xj/88INOnDghwzA85w8fPtyg7T1Tq1at9MILL+iOO+7Qn/70J4WHh+uFF16o9bSO81V121p6O/fJJ5/UefHfueZtukNkZKSysrKqnAsMDFTbtm2Vm5vrOfbVV1/J4XCotLTU66hjenq6JNe6qUGDBlU6d+WVV1Yp7/5P5/jx45WCzPmoofqd+w1lamqq7r33XuXm5io9Pd0TNoOCgtSjRw/Pm9GysjLt3r1bbdq00eWXX17r9tf07+u2b98+hYWFea27Z8+eWr58ea2f83xi5v6xa9cuSa7RVW/XhR9++MHzp3tKljd17RudO3dWy5YtKx2zWq1q166djh07Vt3LRx015v+Te/fu1eLFi7V9+3YdOXKkSjgpLCys9MHJ+Y6wgkoeeOABffrpp4qOjtbw4cPVrl07Wa1WHTt2TO+8806j3D+jW7duuuCCC5STk6PExMQm86a6MVWcy+z24Ycf6tixY5o0aVKVc2duqehNZGSkJFV6498YvAUCq9Va7Tm73e752b1hw44dOzxT27zxtiWptzVS7uf2dp+H801D9DvJFSJatWqlHTt2yG63e950xsfHe8rEx8dr/vz5Kigo0Pfff6+TJ0/qmmuukcViqXX7q+tbFf99i4uLfY4wtmvXrtbPd74xc/9wXxe2bNmiLVu2+HyumrYqrmvf8LX+0mq1et3gAfXXUP2xJjt27PDUP3DgQEVHRyssLEwWi0Uff/yxvv32W+5VdgbCCjz27NmjTz/9VAkJCVq0aFGl0Ytdu3bpnXfeqVf9Foul0pvGio4fP+7zYv3ss88qJydHbdq00fr16zV69GglJCTUqy3Nnbe5zO5pD3Wd59yzZ09JrqmAM2fOrPXj3G8GvL2Zb4zd3NxvSO+66y799re/PefP35w1RL+TXPfJ6t27tzZv3qyvvvrKsx6hW7dunjLuhdSpqak6cOCApLpvSVsTm82mwsJCr+fy8/Mb5DmbAzP3D/d14fHHH9eECRPq3Bb6RtPRUP2xJn/9619VVlbmueF5Re4RPlTGmhV4ZGVlSZKuvfbaKtOsvG3nerZat27t9VP57Oxsn8Pcn3zyif7+97+rb9++WrlypWw2mx599FEVFBTUuz04O3379lVUVJR27txZaW64NxU/FWrdurUk7yMye/fu9W8ja+Gqq66SxWLRzp07z/lzo+4qblGblpamHj16KCgoyHM+Li5OLVq0UFpaWr1v9leT2NhYnThxQt9++22Vc9WN1qHh1Ld/uHf5qu914Vz0Dfe6FkZzm6bMzEy1adOmSlApKSnxunUxCCuowH2Dq+3bt1c6vn//fi1atKje9Xfr1k05OTmVFjmWlZXp2Wef9Vr+8OHD+v3vf682bdro+eefV1RUlJ588knl5eXpd7/7Xb3bg7MTGBioP/zhDwoICNCsWbN8brawefNmPfDAA56f3Z9u/uMf/6g0nWHnzp017qzTECIjIzVs2DDt3LlTixcvrrQmy2337t3cmdpk3Ls2rVu3ThkZGVW2gQ0ODlb37t2VnJys3bt3q127drrssssapC033nijJOnVV1+t1KcPHDigf/zjHw3ynKhefftHXFycunfvrvXr12vDhg1V6nc6nVXu1eLNuegb7g+AvO3oCfPr1KmTioqKtH//fs8xh8OhefPm8UGsD0wDg0dcXJzi4uK0ceNG5eXlqXv37vrpp5+0efNmJSYm6qOPPqpX/ZMmTVJycrLuu+8+jRgxQqGhodq6datatWrlWQ/hZhiGHn30URUWFmrBggWeRa4jR47UZ599ptWrV2vp0qX1Gq7H2bvmmmv03HPPac6cOfrNb36jbt26qUePHmrZsqWOHDmitLQ0ZWZmasCAAZ7HXH311erRo4dSUlJ02223qXfv3vrxxx+1efNmDRo0SP/617/O+et44okndPDgQT3//PNavXq1evToIZvNpkOHDumbb75Renq6kpOTFRoaes7bBu9iY2PVpk0bz25MFdcjuPXp08ezOLriTWb97eabb9bq1au1efNm3XzzzUpISFBRUZHWr1+vAQMG6NNPPz2rtTKoP3/0jxdffFGTJk3ybHd+5ZVXKiQkRD/++KN27dqlgoICffXVV9W241z0jX79+umjjz7SzJkzlZiYqJCQEP3yl7/UtddeW696cW5MmDBBycnJGj9+vIYNG6bg4GClpaUpNzdXffr0qVUoPt8wsgKPwMBAvf7667rllluUmZmppUuX6vvvv9cjjzyihx9+uN71X3PNNXr55ZcVFRWl1atX65///KcGDhyoN998s8pNBN98801t3bpVY8aM8eyL7/aHP/xBUVFReu6552rcRhL+N2rUKG3atEn33HOPnE6nPvzwQy1evFhbtmxR586d9fTTT1caibNYLHrttdd00003KTMzU8uWLdOhQ4f0l7/8pUHfUFanTZs2eu+99/Twww8rKChIa9euVVJSknbv3q2YmBjNmzfPc/M4mIPFYvF8Wh4aGur1pmkV36A21BQwyXWtXLRoke666y7l5+fr7bff1o4dO/Too4967qVxvu8Md675o39ERUXpww8/1NSpU3XixAl98MEHeu+99/Ttt9+qd+/eeumll2psx7noG2PHjtXdd9+t/Px8/fWvf9WLL76of/7zn/WqE+fOoEGDNH/+fEVFRWnNmjVat26dLrnkEq1cuVKdOnVq7OaZksXwNgcCAACctZdffll//etftWjRIiUmJjZ2c2Ai9A2gbhhZAQDgLHm759T333+vd999V61ataqyZgLnD/oG4F+sWQEA4Cw9+eSTysnJUVxcnFq1aqWsrCxt3rxZdrtdf/rTn1jvdB6jbwD+xTQwAADO0po1a/Tee+/pwIEDKi4uVlhYmK666irdeeed+tWvftXYzUMjom8A/kVYAQAAAGBKrFkBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACm9P8Bxyg34nKUVG0AAAAASUVORK5CYII=", "text/plain": [ "
<xarray.DataArray 'y' ()> Size: 8B\n", "array(96.34186308)
<xarray.DataArray 'y' ()> Size: 8B\n", "array(96.34186308)
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 423 seconds.\n", "Sampling: [y]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "fd3bc21fab6642fa80eb59bfd081901f", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n" ], "text/plain": [ "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with pm.Model(coords=coords) as model_t:\n", " μ_t = pmb.BART(\"μ\", x_0, y_0, m=50, separate_trees=True, dims=[\"species\", \"n_obs\"])\n", " θ_t = pm.Deterministic(\"θ\", pm.math.softmax(μ_t, axis=0))\n", " y_t = pm.Categorical(\"y\", p=θ_t.T, observed=y_0)\n", " idata_t = pm.sample(chains=4, compute_convergence_checks=False, random_seed=123)\n", " pm.sample_posterior_predictive(idata_t, extend_inferencedata=True)" ] }, { "cell_type": "markdown", "id": "60dc23a9-2351-4502-824a-944e0f454c4c", "metadata": {}, "source": [ "Now we are going to reproduce the same analyses as before. " ] }, { "cell_type": "code", "execution_count": 16, "id": "a05a3d39-307a-4c08-93ec-3a0503ea6c25", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAE3CAYAAACq3N6VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA87UlEQVR4nO3deXxU1cH/8e9MJoGELASKtMUU0GAEQpAlQCAVQXkUWSwoyIMgVVwAQdAftdpi9bHUirtAsSLuhAoqlr2iItVQkyAQWQQEBBJiWSQxJBJIZub+/hhmyDCThcwkuYTP+/XiRebekzNn4HC535zlWgzDMAQAAAAAJmOt7wYAAAAAgD+EFQAAAACmRFgBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACmRFgBAAAAYEqEFQAAAACmZKvvBlRl2bJl2rRpk7Zv365vv/1WZWVl+utf/6rhw4efVz1Op1OLFi3S4sWLdfDgQUVERKhnz5564IEH1KZNm9ppPAAAAIAaM31Yeemll5SXl6fY2FhdcsklysvLq1E9jz32mJYsWaL4+HiNGTNGx48f1+rVq7Vhwwa9++67io+PD3LLAQAAAATC9NPAZs6cqXXr1ikjI0OjRo2qUR0ZGRlasmSJunfvrg8//FAPPfSQZs2apfnz56u4uFiPP/54cBsNAAAAIGCmDyu9e/dWq1atAqrjvffekyRNmzZNYWFhnuMpKSlKTU3Vxo0btX///oDeAwAAAEBwmT6sBENmZqYiIiLUtWtXn3OpqamSpI0bN9Z1swAAAABUosGHlZMnT+rYsWO69NJLFRIS4nPevbj+wIEDddswAAAAAJVq8GGlqKhIkhQZGen3vPt4cXFxpfU4nc7gNgwAAABApUy/G5hZFBYW1ncTTC82NlYFBQX13QxcwOhDCAb6EQJFH0Iw0I+qFhsbW2WZBj+yEhUVJanikRP38YpGXgAAAADUjwYfViIiItSiRQsdOnRIDofD57x7rQoPhgQAAADMpcGHFUnq0aOHTp48qc2bN/ucS09PlyQlJyfXdbMAAAAAVKJBhZX8/Hzt27dP+fn5XsdHjhwpSXrxxRdVWlrqOf7ll18qPT1dycnJatu2bZ22FQAAAA1PSYmh1Guc6tj5uEpKjPpuzgXP9Avs33vvPW3atEmS9O2333qOZWVlSZKuu+46XXfddZKktLQ0zZ07V5MnT9aUKVM8dfTq1UsjRozQe++9p2HDhqlv3746fvy4Vq9ercjISJ5gDwAAAJiQ6cPKpk2b9OGHH3od27x5s2dKV6tWrTxhpTJPPPGEEhIStHjxYr3zzjuKiIhQv3799MADDzCqAgAAVFJiaMBAQ9JxfbzGovBwS303CbjoWQzDYHyqGth6rmps0YdA0YcQDPQj1NTZsCLCCmqMflR9bF0MALhoME8cABoewgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAAAAADAlwgoAAAAAUyKsAKh37OIEAAD8IawAAAAAMCXCCgAAAABTstV3Ay4UZWVlFZ6zWCyy2Wz1WlaSQkNDa1TWbrfLMCqeelPdsqWlpV6vg1WvJNlsNlksllot63A45HQ6g1I2JCREVqvVNGWdTqccDkeFZa1Wq0JCQuqtbFmZIavFIqfhKmsYhux2e7Xqraps+X9HtVXW9Rm4RtRn2fJ/DlarQ2Vlhmw2/0+N5hrhW9bs14hzy9bmNSIkxPVvo6zM4tOHuEb4L3uhXCPO59+9q5x06pRDdrshu11yOCS7w3XcYXd9LYXI4bDIbpdKS50qK3PqZIlktRqSDK9+xDXCt2x1EFaq6eWXX67wXJs2bTR06FDP61dffbXCC1irVq108803e16/8cYbOnXqlN+yl1xyiUaNGuV5vXDhQhUVFfkt26xZM40ZM8bzevHixcrPz/dbNioqSnfccYfn9fvvv6+jR4/6Ldu4cWPdc889ntfLli1TXl6e37KhoaGaOHGi5/Xq1at14MABv2Ul6f777/d8vXbtWu3du7fCshMnTvRclD777DPt3LmzwrJ33XWXIiIiJElffPGFtm3bVmHZ3/72t4qOjpYkffnll9q8eXOFZW+77TY1b95ckrRx40ZlZWVVWPbWW29Vy5YtJUnZ2dnasGFDhWWHDx+uSy+9VJK0fft2/fvf/66w7JAhQ9S2bVtJ0u7du/XJJ59UWHbgwIFq166dJGnfvn1as2ZNhWWvu+46dejQQZJ08OBBrVixosKyffv2VefOnSVJ33//vZYuXVph2T59+qhbt26SpGPHjmnx4sUVlm1/ZbJ27OwpScrPz1daWlqFZbt27arU1FRJUlFRkd58880Ky3bq1En9+vWTJJWUlGjBggUVt6F9ew0YMECS6z+zyv7dx8fH68Ybb/S85hrhUtk1wmazadKkSZ7Xwb5GuP9L69blM7355q4Ky3KNcLnQrhE9evRQr169JNXeNeLUqVO6+SbXNcLft3CNcCl/jTAMQ0uWfKAjR36Q02mV07DKcFrldFplGFaFhobr5ltu9dzgf/Lp5zpy5PiZ8yFyOi0yDFd5iyVU/foNcIUCu7Ql+xv98MOPZ8pa5HSGeMo6DauSOl0l+5nQcPDgf3Wi8Cc5z5w33G0587pFi1/I4bTIYZcKT5ToVEnZ2baW+93ptEpyyul0B9XqTEJyBx+LpLM34U1jjunNN8/2aa4RLuWvEdVBWAEAAKgDhiHfG2SnVSdONFZenuG5Sc8vaK6yMsNzo+25QXdadbKkhT5dZ3hu0r/dm6DS044zN9oh5b7Hosb7o/RDvlN2h2tU4JsdKTp9puy5N+khIY2UsdHpacPRYze42nBOW52GVTJClLbY6WqDXZJGVPq5/7my/CjGNZWW/U9m+bKJlZbd4/Xzi1aVlv3hePlXTSotWxGLnLJYnbJanbJaXF9HNmms0FCrbDaptPSkTpeedJ2zONX6V9/W6H3gzWJUNg4Gj4p+qigxfOvWtGlT/fTTT0GvV2IaWKBlzTBto7KyJSWGBt/kmgb28RqLGjcW08CCXFaq/2tEbZa12Ww6dUoaMNCQ1erQyn8aCg9nGphZrxGGIbmaZpVklcPhmg5aWuqQ0+m6sS7/y+nUmZtk65kbb1dZTxln+bIWz0/rHQ7JYTdUZnf61Omu1+GwyDAscjik06elfy53tXfAde4+IjnsFs8Nv9Np8dykl5UZnqlBjjPH3OVc04QsZ6cP2Y1yP61v+Gw2QyEhki1EstksstmkEJsUYjU8X7vOSSEhrl82mxQaavF8HWI1FGIzZDtzvvz3uMpaPXVbLA6verzqthlqFBYiW6hFthBJcigkxPDU41W/TWrUKEShNlc73GVDQiSrn0EWf//uS0oMDRnmuq6s+NDiuRZxjfAtGxsbW2E5N0ZWqqn8f5oNrWz5m51AyoaFhXmFlWDVW1dlQ0JCqj2H8kIra7VaPRccM5a12w05y90wWiyWavdhM5SVzPFv2ezXiNov6+pDTmeIQkMtCg2t+sawLq8RhmF4bsRPn3bdcLteW+VwWMvdDPveVDscxpnfrWd+VVLW6VvWu5xx9r0c5W/iQ7zK2svV6TrmW9bhL1g4/JT1U67cn0z5P7kq/mRrWvZ89hNy1bv24+qUPZ/w4b+s1XrmZtlW7kb7nJvt8jf3XmX83Gjb/JQ7+70Wr+89GybOfn/5Y+eWC6nge73aeea9LZa63sPpfG5pa6es+9+93W7I4XD1v4quRWa4NzDDfUR1EFYAABcku93QiRNS4QmpsFA6euzszen7Sw1ZreVvtg2vm2/vG+tzb/Yr+gl8NctWUg5VC7GeveG1hpz9uvwvzw3+Ocd8yvo7VlG9VslpSO+eWWJwxzgpPNzidZN+7k/tz71JP/dYRYGj/DGr9eIZbQFqgrACAKh3DoehoiLpxAnpx0JX+CgsdH9tlPv67PHi4orre+VVyfsn6+Z2vjfVFd3AV6esvzptIa6b5orqtFX1PpWECv/vVXFZ9zSZ+lBSYujdxa5+M3qUpcKphADqDmEFABBUhmGouFh+A4Y7eJx77kSRax3D+bJYpKgoKSbG9fs337iOX9tfCgutPAC4p8Wc90/la3CzXtnNvtVavzfoAGBmhBUAQIUMw1BJydlg4R0+zgaP8udOFLqmPdVEZBNX8IiJkZrGnP06Jsbi9dr9dVSUK3BIrp+KDxjoSjwP/46figNAQ0BYAYCLyKlThp/RDteIR/nj7nMnTkhVbBxWofBw+Q0YMTEWnzDSNEaKjlaFD3IEAFycCCsAcIE6fdqoYI2Hd/goX+b06Zq9V1iY1LSpv/Bh8Qkj7uARFkbwAAAEhrACACZgt/tfy3HuAvMTJ84eLymp2XuFhvqfatXUz6hHdLTr68aNCR4AgLpHWAGAIHM4DJ0o8jfa4VrnccLP8eKfqq7XnxBr+XUd54YPi9c59/nwcBZ0AwAuDIQVAKiE0+na2crvVKsT/heYFwWws1VMdEXhw+KZXlX+eGQkwQMA0HARVgBcNAzD0E8/+Y5quB4qWPECc2cNd7Zyb6lb4ToPP8HDvbMVAAAgrAC4QLm31PUd7fD/EEH365o+RbxJBVvqNo2xeNZ1lD8fFcXOVgAABIqwAsBU9u41dOq0/2d5nBs+Smu6pW7jiqda+X3GR7QUGkrwAACgrhFWANSL06cN7dwlbd0mZX99doHHxCmSVP0FH2Ghri11Y84JF+cuMC8fPho1IngAAHAhIKwAqBMFPxratk3aus3Qtu3S7m8lu923XLNmUmxT3+BR0RPMGzdmgTkAAA0VYQVA0BmGodxcaev2s+EkN9e3XPPmUlIn6coE6eVXXMcWp1kUHk74AABcmMLDLUpfb1FsbKwKCgrquzkXPMIKgICVlhr6do9rSte2M+Hkx0Lfcm3buMJJp04WJXWSfvFz16hISYmhl1+pwV6/AACgQSOsADhvJ4oMbS83arJzp+9i97BQqX37M+Ek0aLEjlJ0NCMmAACg+ggrACplGIa+/6+0zR1Otkn7D/iWaxojdUo8O2pyRTspLIxwAgAAao6wAsCL3W5o717v9SbHj/uWi4tzhZOkThYlJbpes9AdAAAEE2EFASspMTRgoCHpuD5ew+LoC81PPxna8Y20bbuhrdukb76RSk55l7HZpIQrzk7p6pQoxcby9wwAAGoXYQW4yBw96gol7nCy7zvJ6fQuExl5dtSkU6LU/kqeTQKg4WMXJ8B8CCtAA+ZwGNp/QF7h5MgR33K/+IWU5F5vkii1aSNZrYQTAABQvwgrQANy6pShb3aeDSfbd0g//eRdxmqV2sVLnTqdXW/ys58RTAAAgPlcEGFl69atmjNnjrKzs1VWVqb4+HiNGzdOQ4YMqXYdJ06c0BtvvKFPPvlEhw4dUlhYmC699FINGzZMI0aMUKNGjWrxEwC1Iz/f8OzStXW79O23ksPhXSY8XErs6FprktRJ6tBeioggnAAAAPMzfVjJzMzU+PHjFRoaqkGDBikqKkpr167V9OnTlZeXpwkTJlRZx4kTJzR8+HDl5uaqW7duGjVqlEpLS/X555/rz3/+sz7++GO98cYbslqtdfCJgJoxDEMHc6Rt287u0nUoz7dci595P3jxsraSzUY4QcPHegMAaHhMHVbsdrtmzJghi8WitLQ0dejQQZJ03333adSoUZozZ45uuOEGtWnTptJ6Fi9erNzcXP32t7/VI4884jleWlqq0aNHKyMjQ5s2bVJycnJtfhzgvJSWGtq1++zzTbZvlwpPeJexWFxhpFMnKenMyEnLlmwhDAAAGgZTh5WMjAzl5ORo+PDhnqAiSZGRkZo0aZIeeOABLV26VA8++GCl9eTm5kqS+vbt63U8LCxMffr00bZt23Tc34MkgDpUWGho2w5p2zbXQvjdu32fCt+okWtnLvfISWIHKSqKYAIAABomU4eVrKwsSVJqaqrPuT59+niVqUy7du0kSV988YV69+7tOV5WVqb//Oc/aty4sbp06RKMJgPVYhiG8vLKPRV+u3TgoG+52FjvLYSvaCeFhhJOAADAxcHUYeXAgQOSpNatW/uci4mJUWxsrA4e9HOHd44RI0Zo2bJlev3117V9+3YlJiaqrKxMX3zxhQoLC/Xcc8+pZcuWwW4+4GG3G/p2T7lwsk3K9zOlvvWvyoWTTtKlrZjSBQAALl6mDivFxcWSpKioKL/nIyMjdfjw4Srrady4sd555x396U9/0vLlyz2jMVarVbfddpu6du1aZR0xMTEswK9Ao0aGpHxJUtOmsew0JamoyKnsrXZt2WLX5i1l2r7d7vNU+NBQqWNHm7pcZVPXLqHq0tmm2NiLs4/RhxBssbGx9d0EXODoQwgG+lHgTB1WgiU/P1+TJk1Sfn6+5s+fr65du+r06dNat26dnnrqKa1fv14ffPCBYmJiKqyjsLCwDlt8YSkpMTxf//hjgU6fvvhuNA8f8X4q/HffSYbhXSYqyjVq4t5C+MoEqVEjp6TSM7+ki3UDI/oQgondwBAo+hCCgX5UteqEOVOHlcjISElSUVGR3/PFxcUVjrqU99RTT2nLli1atmyZrrzySkmu0ZqRI0fK4XDo8ccf11tvvaX7778/eI1Hg+VwGPruO9eDF7dud03pOnrMt9wvf3lmIfyZcNL6VzwVHgAA4HyYOqy4tyQ+ePCgEhMTvc4VFhaqoKCgWgvj//3vf6tp06aeoFJer169JEk7duwIvMFokEpKfJ8Kf/Kkd5kQq9Su3dlw0qmT9LPmBBMAAIBAmDqsJCcn65VXXlF6eroGDRrkdW7Dhg2SpB49elRZT2lpqedXWFiY17n8fNc8+XOP4+L1w3HD68GLe/ZIDqd3mYgI11Ph3bt0dWjveiAdAAAAgsfUYSUlJUVxcXFauXKlbr/9drVv316Sa/rXvHnzZLPZNGzYME/5/Px8FRQUKDY2Vs2aNfMc79q1q9LT0zVv3jxNmzbNc7y0tFTz5s2TJPXs2bNuPhRMxek0dOCg9xbC33/vW+6SS87u0pWUKF12mRQSQjgBAACoTaYOKzabTTNnztRdd92l0aNHa/DgwYqMjNTatWt16NAhTZs2TW3btvWUT0tL09y5czV58mRNmTLFc3z69OnasmWLXn75ZW3YsMGzwD49PV25ubnq2LGjRowYUR8fEXXs9GnXU+G3bnM9fHHbDuncJVEWi3T5ZWeeCn9m5OTnLQkmAAAAdc3UYUVyrSlZtGiRZs+erTVr1qisrEzx8fGaOnWqhg4dWq062rdvr6VLl+qVV15RRkaG0tLSFBISol/96leaMmWKxo8fr0aNGtXyJ0F9+PFH12iJe5eu3d9KZec8Fb5xY9c0LvfISccOUmQk4QQAAKC+WQzj3A1W4Q9bz1WspMTQgIGubvTxGku9rd0wDEOH8s6OmmzdJuXk+pZrFntmIXwn1y5d7eIlm41wUp/M0ofQMLBdKAJFH0Iw0I+qdsFvXQxUpqzM0O5vvdeb/Pijb7k2bc6MmpzZQviXv+Sp8AAAABcCwgouGEVFrm2D3c82+WanVFrqXSY01PWwxaQz600SO0oxMQQTAACACxFhBaZkGIYOH5bXU+H3H/B9KnxMtGshfKdE10L4hCukRo0IJwAAAA0BYQWmYLcb2veddzj54Qffcpe28n4q/K9+xZQuAACAhoqwgnpx8qShHd+cXWuy4xuppMS7TEiIdMUVZ6Z0nRk5adaMYNIQhYdblL7ewmJEAADghbCCOnHsmKGt28/s0rVd2rtXcp7zVPgmTXyfCt+4MeEEAADgYkVYQdA5nYb2H5C2bTs7cvLfw77lft7S+8GLbdvwVHgAAACcRVhBUP3hUUM7d0nFxd7HrVYp/vKzD17slChdcgnBBAAAABUjrCBg+/ef3aJr41eu38MbSx06nN1CuEN7qUkTwgkAAACqj7CCgF1yydmvJ94rdetqUfzlPBUeAAAAgSGsIGCuERPX6Mrw31gUHk5IAQAAQOCs9d0AAAAAAPCHsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlAgrAAAAAEyJsAIAAADAlGz13YDq2Lp1q+bMmaPs7GyVlZUpPj5e48aN05AhQ86rnuLiYr3++utau3atcnNzFRoaqri4OF177bWaPHlyLbUeAAAAQE2YPqxkZmZq/PjxCg0N1aBBgxQVFaW1a9dq+vTpysvL04QJE6pVz/fff69x48YpNzdXvXv3Vt++fVVaWqqcnBx99NFHhJUAhIdblL7eotjYWBUUFNR3cwAAANBAmDqs2O12zZgxQxaLRWlpaerQoYMk6b777tOoUaM0Z84c3XDDDWrTpk2l9TgcDt1///06evSo3nzzTfXq1cvnfQAAAACYi6nXrGRkZCgnJ0eDBw/2BBVJioyM1KRJk2S327V06dIq6/noo4+0bds23XnnnT5BRZJsNlNnNgAAAOCiZOq79KysLElSamqqz7k+ffp4lanM6tWrJUk33HCD/vvf/2r9+vUqKipSXFycrr76ajVp0iSIrQYAAAAQDKYOKwcOHJAktW7d2udcTEyMYmNjdfDgwSrr2b59uyRp06ZN+utf/6rS0lLPuWbNmunFF19Uz549g9NoAAAAAEFh6rBSXFwsSYqKivJ7PjIyUocPH66ynuPHj0uSZs6cqTvvvFNjxoxRWFiYVq1apVmzZum+++7T6tWrdckll1RYR0xMjKxWU8+aM4XY2Nj6bgIucPQhBAP9CIGiDyEY6EeBM3VYCRbDMCRJ11xzjaZPn+45PnbsWB05ckSvvvqq3n//fU2aNKnCOgoLC2u9nRc6dgNDoOhDCAb6EQJFH0Iw0I+qVp0wZ+qhgsjISElSUVGR3/PFxcUVjrr4q6d///4+5/r16yfp7FQxAAAAAOZg6rDi3pLY37qUwsJCFRQU+F3Pcq62bdtKkqKjo33OuY+dPn06gJYCAAAACDZTh5Xk5GRJUnp6us+5DRs2SJJ69OhRZT3u7Yr37t3rc859rFWrVjVuJwAAAIDgM3VYSUlJUVxcnFauXKmdO3d6jhcXF2vevHmy2WwaNmyY53h+fr727dun/Px8r3qGDx+usLAwLVy4UEeOHPGq55VXXpEkDRw4sJY/DQAAAIDzYeoF9jabTTNnztRdd92l0aNHa/DgwYqMjNTatWt16NAhTZs2zTPFS5LS0tI0d+5cTZ48WVOmTPEcj4uL00MPPaSZM2dq6NChGjBggMLCwrR+/Xrl5eXp1ltvVUpKSn18RAAAAAAVOO+wkpCQ4PXaYrGoSZMmuvzyy3XjjTfqtttuU2hoaNAa2KtXLy1atEizZ8/WmjVrVFZWpvj4eE2dOlVDhw6tdj1jx45Vq1at9Nprr2nVqlVyOByKj4/XhAkTNHLkyKC1FwAAAEBwWAz3vr7V5A4r7ulXDodDeXl52rJli5xOp1JSUrRgwQLZbGdz0IoVK/TGG2+osLBQdrtdLVu21H333ae+ffsG8aPULraeqxpb9CFQ9CEEA/0IgaIPIRjoR1WrztbFNZ4G9tRTT3m9/vrrrzV27Fh9+eWXWrVqlW666SbPuRYtWujZZ5/VZZddJqfTqTlz5ui+++7TihUrvKZxAQAAAIBb0BbYd+7c2TPacu7uXb169dJll13mekOrVX379lVZWZkOHDgQrLcHAAAA0MAEdTewdu3aSZLPblzllZWV6emnn1abNm3Uu3fvYL49AAAAgAYkqLuB/fTTT5KkZs2a+T1vt9v14IMP6rvvvlNaWpoaNWoUzLcHAAAA0IAENax88cUXkqRf//rXPudKSko0bdo07dq1S2+//bYuv/zyYL41AAAAgAYm4LDidDp16NAhvfbaa9q4caP69++vG2+80atMfn6+Jk6cqNLSUi1ZskQtW7YM9G0BAAAANHA1DivnPm9Fkm655Rb9+c9/ltXqvRTm+eefV3Z2tlq2bKkxY8Z4jo8dO1a33357TZsAAAAAoAEL+Dkrp0+f1s6dO7V//35J0syZMzVixIggN7P+sU921dhPHIGiDyEY6EcIFH0IwUA/qlqdPmfl1Vdf1bPPPquZM2eqd+/eatWqVU2rBgAAAIDgbV189913KzU1VadOndLcuXODVS0AAACAi1RQn7Myffp0WSwWLV++XHl5ecGsGgAAAMBFJqhhpX379rr22mtlt9u1YMGCYFYNAAAA4CIT1LAiSVOmTJHFYtEHH3ygY8eOBbt6AAAAABeJoIeVK6+8UgMGDNDp06f1xhtvBLt6AAAAABeJ894NbPfu3VWWmTNnTo0aAwAAAABuQR9ZAQAAAIBgIKwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMCXCCgAAAABTIqwAAAAAMKULIqxs3bpVd999t5KTk3XVVVfplltu0YoVK2pcX1lZmW666SYlJCTohhtuCGJLAQAAAASLrb4bUJXMzEyNHz9eoaGhGjRokKKiorR27VpNnz5deXl5mjBhwnnXOW/ePOXk5NRCawEAAAAEi6lHVux2u2bMmCGLxaK0tDTNnDlTv//977Vs2TK1a9dOc+bM0YEDB86rzh07dmj+/Pl68MEHa6fRAAAAAILC1GElIyNDOTk5Gjx4sDp06OA5HhkZqUmTJslut2vp0qXVrq+0tFQPP/ywOnfurDFjxtRGkwEAAAAEiamngWVlZUmSUlNTfc716dPHq0x1zJ07VwcPHtSyZctksViC00gAAAAAtcLUIyvuKV6tW7f2ORcTE6PY2FgdPHiwWnVt3bpVCxYs0JQpU9S2bdtgNhMAAABALTD1yEpxcbEkKSoqyu/5yMhIHT58uMp6SktL9cgjj6h9+/a68847a9SWmJgYWa2mznamEBsbW99NwAWOPoRgoB8hUPQhBAP9KHCmDivB8uKLL+rgwYP64IMPFBISUqM6CgsLg9yqhic2NlYFBQX13QxcwOhDCAb6EQJFH0Iw0I+qVp0wZ+qhgsjISElSUVGR3/PFxcUVjrq47dixQ2+++aYmTJighISEoLcRAAAAQO0wdVhp06aNJPldl1JYWKiCggK/61nK2717txwOh+bMmaOEhASvX5K0f/9+JSQkqHv37kFvPwAAAICaM/U0sOTkZL3yyitKT0/XoEGDvM5t2LBBktSjR49K62jTpo1uueUWv+fef/99RUVF6frrr1d4eHhwGg0AAAAgKCyGYRj13YiK2O123XDDDTpy5IiWLFmi9u3bS3JN/xo1apT279+vlStXenb3ys/PV0FBgWJjY9WsWbMq609ISFDbtm31r3/9q8qyzDmsGnMzESj6EIKBfoRA0YcQDPSjql3wa1ZsNptmzpwpwzA0evRoPfroo5o1a5Zuuukm7dmzR5MnT/bahjgtLU033nij0tLS6rHVAAAAAILB1NPAJKlXr15atGiRZs+erTVr1qisrEzx8fGaOnWqhg4dWt/NAwAAAFBLTD0NzEwYxqsaw50IFH0IwUA/QqDoQwgG+lHVLvhpYAAAAAAuXoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSoQVAAAAAKZEWAEAAABgSrb6bkB1bN26VXPmzFF2drbKysoUHx+vcePGaciQIdX6/q+++kqffPKJsrKylJeXp5MnT6pVq1a69tprde+99yo6OrqWPwEAAACA82X6sJKZmanx48crNDRUgwYNUlRUlNauXavp06crLy9PEyZMqLKOqVOnqqCgQN26ddNNN90ki8WirKwsLViwQGvXrtW7776r5s2b18GnAQAAAFBdFsMwjPpuREXsdrsGDhyow4cPa/HixerQoYMkqbi4WKNGjdL+/fu1atUqtWnTptJ65s+fr9/85je65JJLPMcMw9D//d//6R//+IdGjx6txx57rNI6CgoKAv48DV1sbCx/TggIfQjBQD9CoOhDCAb6UdViY2OrLGPqNSsZGRnKycnR4MGDPUFFkiIjIzVp0iTZ7XYtXbq0ynruuecer6AiSRaLRZMmTZIkbdy4MbgNBwAAABAwU4eVrKwsSVJqaqrPuT59+niVqQmbzTULLiQkpMZ1AAAAAKgdpg4rBw4ckCS1bt3a51xMTIxiY2N18ODBGtf/wQcfSDobfAAAAACYh6kX2BcXF0uSoqKi/J6PjIzU4cOHa1T3zp079be//U3NmzfXXXfdVWX5mJgYWa2mznamUJ25h0Bl6EMIBvoRAkUfQjDQjwJn6rBSW3Jzc3XvvffK4XDo+eefV7Nmzar8nsLCwjpo2YWNhWQIFH0IwUA/QqDoQwgG+lHVqhPmTB1WIiMjJUlFRUV+zxcXF1c46lKRvLw8jRs3Tvn5+ZozZ4569eoVcDsBAAAABJ+p5zW5tyT2ty6lsLBQBQUFftezVOTQoUMaO3asjh49qhdffFH9+vULVlMBAAAABJmpw0pycrIkKT093efchg0bJEk9evSoVl2HDh3S7bffrqNHj+qFF17QddddF7yGAgAAAAg6U4eVlJQUxcXFaeXKldq5c6fneHFxsebNmyebzaZhw4Z5jufn52vfvn3Kz8/3qscdVI4cOaLnn39eAwYMqLPPAAAAAKBmTL1mxWazaebMmbrrrrs0evRoDR48WJGRkVq7dq0OHTqkadOmqW3btp7yaWlpmjt3riZPnqwpU6Z4jt9+++3Ky8vTVVddpd27d2v37t0+71W+PAAAAID6Z+qwIkm9evXSokWLNHv2bK1Zs0ZlZWWKj4/X1KlTNXTo0GrVkZeXJ0nKzs5Wdna23zKEFQAAAMBcLIZhGPXdiAsBW89VjS36ECj6EIKBfoRA0YcQDPSjqlVn62JTr1kBAAAAcPEirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFMirAAAAAAwJcIKAAAAAFO6IMLK1q1bdffddys5OVlXXXWVbrnlFq1YseK86nA6nVq4cKGGDBmipKQk9erVS1OnTtWBAwdqp9EAAAAAAmKr7wZUJTMzU+PHj1doaKgGDRqkqKgorV27VtOnT1deXp4mTJhQrXoee+wxLVmyRPHx8RozZoyOHz+u1atXa8OGDXr33XcVHx9fy58EAAAAwPmwGIZh1HcjKmK32zVw4EAdPnxYixcvVocOHSRJxcXFGjVqlPbv369Vq1apTZs2ldaTkZGhcePGqXv37nrjjTcUFhYmSfryyy91xx13qHv37lq4cGGldRQUFATlMzVksbGx/DkhIPQhBAP9CIGiDyEY6EdVi42NrbKMqaeBZWRkKCcnR4MHD/YEFUmKjIzUpEmTZLfbtXTp0irree+99yRJ06ZN8wQVSUpJSVFqaqo2btyo/fv3B/8DAAAAAKgxU4eVrKwsSVJqaqrPuT59+niVqUxmZqYiIiLUtWtXn3Puujdu3BhIUwEAAAAEmanDinvxe+vWrX3OxcTEKDY2VgcPHqy0jpMnT+rYsWO69NJLFRIS4nPePYWMhfYAAACAuZh6gX1xcbEkKSoqyu/5yMhIHT58uNI6ioqKPGUrqqP8e1WkOnPqwJ8TAkcfQjDQjxAo+hCCgX4UOFOPrAAAAAC4eJk6rLhHPdyjI+cqLi6ucNTFzX2+opET9/GKRl4AAAAA1A9ThxX3ehJ/61IKCwtVUFDgdz1LeREREWrRooUOHTokh8Phc969VqWq7Y8BAAAA1C1Th5Xk5GRJUnp6us+5DRs2SJJ69OhRZT09evTQyZMntXnzZp9z7rrd7wUAAADAHEwdVlJSUhQXF6eVK1dq586dnuPFxcWaN2+ebDabhg0b5jmen5+vffv2KT8/36uekSNHSpJefPFFlZaWeo5/+eWXSk9PV3Jystq2bVvLnwYAAADA+TB1WLHZbJo5c6YMw9Do0aP16KOPatasWbrpppu0Z88eTZ482StkpKWl6cYbb1RaWppXPb169dKIESP01VdfadiwYXr66af1+9//Xvfcc48iIyP1+OOP1/EnazjGjh2rhIQEr2OZmZlKSEjQnDlz6qlVuNAlJCRo7Nix9d0MNED0LQSLv///aoI+icpUdE8VrP53ITB1WJFcQWPRokXq1q2b1qxZo0WLFqlp06Z65plnNHHixGrX88QTT2jGjBmyWCx65513tH79evXr10/vvfee4uPja/ET1I9Dhw4pISFB48ePr7BMdna2EhIS9PDDD9dhy9AQHDlyRM8995yGDRum7t27KzExUampqbrnnnu0dOlSrxFMoKZKS0t11VVXqXv37n7XHK5cuVIJCQlq3769fvzxR5/zX331lRISEnT33XfXQWtR1y7W/nEx3aReSNx/L9X9lZmZWd9NvmCY+jkrbklJSVqwYEGV5aZMmaIpU6b4PWe1WjV27Fh+egEEaOXKlfrjH/+oU6dOqWPHjho6dKiioqJ07NgxZWRk6JFHHtGyZcv01ltv1XdTcYELCwtT165dtWHDBn3zzTfq1KmT1/msrCxZLBY5nU599dVXuu6667zOu28Gevbs6Tm2evVqhYeH137jUetqo3+cj1mzZqmkpKRmjUeDM2zYMJ911B9++KHy8vJ0++23Kzo62utcq1atqlVvUlKSVq9efVE/r+WCCCsAAjd27Fjl5eVp3bp1Na7j888/1+9+9ztFR0dr3rx56tOnj9d5wzD0ySef6L333gu0uWggAu13PXv21IYNG5SZmelzM5qZmanevXsrOztbGRkZPjejWVlZnjrcLr/88hq1A7XDbP3jfPzyl7+s0ffBvALpj8OHD/c5lpWVpby8PI0bN06XXnppjdoUHh5+0V+3TD8NDHVr+/bteuKJJzR48GB169ZNSUlJGjJkiObPn6+ysrKA6q5sXm7//v3Vv39/z+vvvvtOXbp00TXXXKPCwkKvsvv27VPnzp3Vv3//Cp/Bg+BzOBx64okn5HQ69eKLL/oEFUmyWCwaMGCA5s6d6zk2Z86cCoe8ly5dqoSEBC1durTK93/44YeVkJCg3Nxcvfbaa7r++uuVlJSkG2+8UatWrZIklZWV6aWXXlL//v3VqVMnDRkyRF988YXf+oqLizV79mwNGjRISUlJ6t69u8aPH6+vvvrKp6x7eN9ut+tvf/ub+vfvr8TERF1//fU+a+QQXO4bSfeNpdvRo0d14MABpaSkqGvXrj7nS0tLlZ2draioKHXo0MFz3N91yN23Dh06pLS0NA0cOFCdOnVSv379NHfuXDmdTp92lZSU6Omnn1bfvn3VqVMnDR48WEuWLGHNXh0Ldv84fvy4nnzySQ0YMECJiYnq2bOnpkyZom+//dbnvSuajlXTvpGfn69HHnlEKSkpSkpK0siRI32umwkJCZ7PUn5KEdO5LxylpaV65513NH78ePXt21eJiYlKSUnR5MmT9c033/iU55rCyArOsWTJEn322WdKTk7W1VdfrVOnTikrK0vPPfectm3bVmf/WC677DL94Q9/0IwZM/Too49q9uzZklz/yB988EGVlZXpmWeeqfKhoAiezMxM5ebmqkuXLkpJSam0bFhYWK21469//au2bt2qfv36yWq1avXq1fp//+//KTo6WmlpadqzZ4/69u2r06dPa+XKlZo4caLWrFmjuLg4Tx0//vijxowZoz179qh79+5KTU1VUVGRPv30U40bN04vvfSSz09hJenBBx/U1q1bdfXVV8tqtWrNmjV64oknFBoa6tl1EMGVmJioJk2a6KuvvpLD4VBISIiks1N4evToIcMw9Pzzz6ugoMAzVeLrr7/WqVOnlJKS4vmeqjz99NPKyspSv3791KdPH3366aeaM2eOysrK9MADD3jKORwO3XvvvcrMzNSVV16pwYMHq7CwUE899VS1ttNH8ASzf+Tk5Gjs2LE6cuSI+vTpo+uuu07Hjx/X2rVrlZ6erjfffFOdO3eutD017RsnTpzQ//7v/6pJkyYaMmSIjh8/rjVr1mj8+PFaunSprrjiCknS5MmTPVOLJk+e7Pn+9u3b1/wPEXWqsLBQTz75pLp3766+ffsqOjpaubm5WrdunT7//HMtXLhQSUlJ9d1MUyGsNHA5OTkVBozDhw/7HLv33nv12GOPef3nbhiG/vjHP+qDDz7Qpk2b1K1bt1prb3kjRozQF198oY8++khLlizRyJEj9eyzz2rXrl2aPHlynbUDLu7nFPXq1ate27Fv3z4tX75czZo1k+Qaeh85cqQefPBBtWvXTitWrFBERIQkKTU1VQ888IDeeustzZgxw1PHn//8Z+3Zs0dPPvmkbr75Zs/xH374QbfccoseffRR/frXv1ajRo283vvw4cNauXKlIiMjJUm33367hgwZotdff52wUktsNpu6deumzz//XDt27PD8J56VlaWIiAh17NhRTqdThmFo48aN+p//+R9JZ29Wz6e/7tixQ8uXL9cll1wiSZo0aZKuv/56vfPOO7rvvvs8IXzp0qXKzMxU//799be//U1Wq2uSwh133OG1nT5qXzD7x0MPPaQffvhBr732mtfI8cSJE3XzzTdrxowZWrFiRaXtqWnf2LVrl2fXU/f39OrVSzNmzNDChQv1xBNPSHKtzXVPLapojS7MLSYmRuvXr1fLli29ju/Zs0cjR47UCy+8oDfeeKOeWmdOTANr4HJycjR37ly/v95//32f8q1atfL5KaTFYtFtt90myfVsmro0c+ZM/eIXv9CTTz6pt99+W2+//ba6dOmiSZMm1Wk7IB07dkyS9POf/7xe2zFhwgRPUJGkzp07Ky4uTidOnNADDzzgCSqSdP311ys0NFS7d+/2HMvPz9eaNWuUkpLiFVQk6Wc/+5nGjx+v/Px8/ec///F57wcffNATVCTXCGDXrl21f/9+FRcXB/Njohz3VJ+MjAzPsaysLHXt2lU2m02JiYkKDw/3mjJTk/UIkyZN8gQVSWrWrJmuvfZa/fTTT9q/f7/n+PLlyyVJU6dO9dxYSq71ML/5zW/O78MhYMHoH9988422bNmi3/zmNz5TXNu2bauRI0fq22+/9TsdrLya9o2IiAhNnz7d63uGDRsmm82m7du3V/qeuLCEhYX5BBVJateunXr27KmNGzcGPO2+oWFkpYFLTU3Va6+95vdcdna2br31Vq9jpaWlSktL06pVq/Tdd9/p5MmTMgzDc/7o0aO12t5zRUdH69lnn9Xtt9+uv/zlL4qKitKzzz5b7WkdF6vKtrX0d+7TTz+t8eK/uuZvukOLFi2Um5vrcy4kJETNmjXTkSNHPMe2bdsmh8Oh06dP+x11PHDggCTXuql+/fp5nevYsaNPefd/OkVFRV5B5mJUW/3OfUOZmZmpe+65R0eOHNGBAwc8YTM0NFRdunTx3IyWlpbq66+/VtOmTXXllVdWu/1V/f267d69WxEREX7r7tq1qxYvXlzt97yYmLl/ZGdnS3KNrvq7Lnz33Xee391Tsvypad9o3bq1mjRp4nXMZrOpefPmOnHiRGUfHzVUn/9P7ty5UwsWLNCmTZv0ww8/+ISTgoICrx+cXOwIK/By//3367PPPlObNm104403qnnz5rLZbDpx4oTefvvtenl+RmJion7+858rLy9Pffv2vWBuqutT+bnMbh9++KFOnDihcePG+Zw7d0tFf1q0aCFJXjf+9cFfILDZbJWes9vtntfuDRs2b97smdrmj78tSf2tkXK/t7/nPFxsaqPfSa4QER0drc2bN8tut3tuOpOTkz1lkpOTNXv2bOXn52vv3r06deqUrr76alkslmq3v7K+Vf7vt7i4uMIRxubNm1f7/S42Zu4f7uvC+vXrtX79+grfq6qtimvaNypaf2mz2fxu8IDA1VZ/rMrmzZs99ffp00dt2rRRRESELBaLPvnkE+3atYtnlZ2DsAKPrVu36rPPPlNqaqrmz5/vNXqRnZ2tt99+O6D6LRaL101jeUVFRRVerJ966inl5eWpadOmWrVqlYYNG6bU1NSA2tLQ+ZvL7J72UNN5zl27dpXkmgo4derUan+f+2bA3818fezm5r4hvfPOO/X73/++zt+/IauNfie5npPVvXt3rVu3Ttu2bfOsR0hMTPSUcS+kzszM1L59+yTVfEvaqkRGRqqgoMDvuePHj9fKezYEZu4f7uvCo48+qjFjxtS4LfSNC0dt9ceq/P3vf1dpaanngefluUf44I01K/DIzc2VJF1zzTU+06z8bed6vmJiYvz+VP7QoUMVDnN/+umn+sc//qGePXvq/fffV2RkpB5++GHl5+cH3B6cn549eyouLk5btmzxmhvuT/mfCsXExEjyPyKzc+fO4DayGjp16iSLxaItW7bU+Xuj5spvUZuVlaUuXbooNDTUcz4pKUmNGzdWVlZWwA/7q0pCQoJOnjypXbt2+ZyrbLQOtSfQ/uHe5SvQ60Jd9A33uhZGcy9MOTk5atq0qU9QKSkp8bt1MQgrKMf9gKtNmzZ5Hd+zZ4/mz58fcP2JiYnKy8vzWuRYWlqqp556ym/5o0eP6o9//KOaNm2qZ555RnFxcXr88cd17Ngx/eEPfwi4PTg/ISEh+tOf/iSr1app06ZVuNnCunXrdP/993teu3+6+c9//tNrOsOWLVuq3FmnNrRo0UIDBw7Uli1btGDBAq81WW5ff/01T6Y2GfeuTStXrtTBgwd9toENCwtT586dlZ6erq+//lrNmzdXu3btaqUtQ4cOlSS99NJLXn163759+uc//1kr74nKBdo/kpKS1LlzZ61atUqrV6/2qd/pdPo8q8Wfuugb7h8A+dvRE+bXqlUrFRYWas+ePZ5jDodDs2bN4gexFWAaGDySkpKUlJSkNWvW6NixY+rcubP++9//at26derbt68++uijgOofN26c0tPTde+992rQoEEKDw/Xhg0bFB0d7VkP4WYYhh5++GEVFBRozpw5nkWugwcP1ueff65ly5Zp4cKFAQ3X4/xdffXVevrppzVjxgz99re/VWJiorp06aImTZrohx9+UFZWlnJyctS7d2/P91x11VXq0qWLMjIydOutt6p79+76/vvvtW7dOvXr108ff/xxnX+Oxx57TPv379czzzyjZcuWqUuXLoqMjNThw4e1Y8cOHThwQOnp6QoPD6/ztsG/hIQENW3a1LMbU/n1CG49evTwLI4u/5DZYBs+fLiWLVumdevWafjw4UpNTVVhYaFWrVql3r1767PPPjuvtTIIXDD6x3PPPadx48Z5tjvv2LGjGjVqpO+//17Z2dnKz8/Xtm3bKm1HXfSNXr166aOPPtLUqVPVt29fNWrUSFdccYWuueaagOpF3RgzZozS09M1evRoDRw4UGFhYcrKytKRI0fUo0ePaoXiiw0jK/AICQnRK6+8optvvlk5OTlauHCh9u7dq4ceeki/+93vAq7/6quv1gsvvKC4uDgtW7ZM//rXv9SnTx+9/vrrPg8RfP3117VhwwaNGDHCsy++25/+9CfFxcXp6aefrnIbSQTfkCFDtHbtWt19991yOp368MMPtWDBAq1fv16tW7fWk08+6TUSZ7FYNG/ePN10003KycnRokWLdPjwYb388su1ekNZmaZNm+rdd9/V7373O4WGhmrFihVKS0vT119/rfj4eM2aNcvz8DiYg8Vi8fy0PDw83O9D08rfoNbWFDDJda2cP3++7rzzTh0/flxvvfWWNm/erIcfftjzLI2LfWe4uhaM/hEXF6cPP/xQEydO1MmTJ/XBBx/o3Xff1a5du9S9e3c9//zzVbajLvrGyJEjddddd+n48eP6+9//rueee07/+te/AqoTdadfv36aPXu24uLitHz5cq1cuVKXXXaZ3n//fbVq1aq+m2dKFsPfHAgAAHDeXnjhBf3973/X/Pnz1bdv3/puDkyEvgHUDCMrAACcJ3/PnNq7d6/eeecdRUdH+6yZwMWDvgEEF2tWAAA4T48//rjy8vKUlJSk6Oho5ebmat26dbLb7frLX/7CeqeLGH0DCC6mgQEAcJ6WL1+ud999V/v27VNxcbEiIiLUqVMn3XHHHfr1r39d381DPaJvAMFFWAEAAABgSqxZAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApkRYAQAAAGBKhBUAAAAApvT/AR4UCydVTohCAAAAAElFTkSuQmCC", "text/plain": [ "
<xarray.DataArray 'y' ()> Size: 8B\n", "array(97.26565657)
<xarray.DataArray 'y' ()> Size: 8B\n", "array(97.26565657)