{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ARD (Automatic Relevance Determination) Tutorial\n",
    "\n",
    "PHYSBO uses a Gaussian process (GP) as the surrogate model with a Gaussian kernel:\n",
    "\n",
    "$$k(\\vec{x}, \\vec{x}') = \\lambda^2 \\exp\\left[ -\\frac{1}{2} \\frac{\\|\\vec{x} - \\vec{x}'\\|^2}{\\eta^2} \\right] $$\n",
    "\n",
    "Here $\\lambda$ and $\\eta$ are hyperparameters representing the kernel scale and the input length scale, respectively.\n",
    "\n",
    "With the **ARD** kernel, a separate length scale $\\eta_i$ is assigned to each input dimension $i$. Learning these independently allows the model to automatically identify which dimensions matter for the objective and to optimize more efficiently.\n",
    "\n",
    "$$k(\\vec{x}, \\vec{x}') = \\lambda^2 \\exp\\left[ -\\frac{1}{2} \\sum_i \\frac{(x_i - x_i')^2}{\\eta_i^2} \\right] $$\n",
    "\n",
    "In PHYSBO, you can enable the ARD kernel by passing the keyword argument `ard=True` to `policy.bayes_search` when running Bayesian optimization.\n",
    "\n",
    "In this tutorial, we compare runs with and without ARD on a problem where only some dimensions affect the objective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "\n",
    "Import the packages used in this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import physbo\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem setup for this tutorial\n",
    "\n",
    "We take as the objective the weighted sum of squared inputs:\n",
    "\n",
    "$$ f(\\vec{x}; \\vec{w}) = - \\sum_i w_i x_i^2 $$\n",
    "\n",
    "With weight vector $\\vec{w} = [5, 1, 0, 0]$, the search space has $D=4$ dimensions: $x_0$ contributes most to the objective, then $x_1$, while $x_2$ and $x_3$ do not contribute at all.\n",
    "\n",
    "We prepare the candidate set as the point that attains the maximum, $\\vec{x} = \\vec{0}$, plus $N-1 = 999$ points drawn from a $D$-dimensional normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Search space shape: (1000, 4)\n"
     ]
    }
   ],
   "source": [
    "N = 1000\n",
    "\n",
    "weights = np.array([5.0, 1.0, 0.0, 0.0])\n",
    "D = len(weights)\n",
    "\n",
    "np.random.seed(137)\n",
    "test_X = np.random.randn(N, D)\n",
    "test_X[0, :] = 0.0\n",
    "\n",
    "\n",
    "def simulator(actions: np.ndarray) -> np.ndarray:\n",
    "    X2 = test_X[actions, :] ** 2\n",
    "    return - np.einsum(\"ai,i -> a\", X2, weights)\n",
    "\n",
    "\n",
    "print(\"Search space shape:\", test_X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common parameters\n",
    "\n",
    "We set the same seed and the same number of steps for both runs so that the effect of ARD can be compared fairly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_initial = 20\n",
    "n_bayes = 30\n",
    "n_perm = 20\n",
    "score = \"EI\"\n",
    "seed = 31415"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case: ard=True\n",
    "\n",
    "First, we run Bayesian optimization with `ard=True`.\n",
    "Because a length scale is learned per dimension, relevant and irrelevant dimensions can be distinguished more easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================================================\n",
      "Running with ard=True\n",
      "============================================================\n",
      "\n",
      "Best value (ard=True): -0.0\n",
      "Best point: [0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "print(\"=\" * 60)\n",
    "print(\"Running with ard=True\")\n",
    "print(\"=\" * 60)\n",
    "policy_ard = physbo.search.discrete.Policy(test_X)\n",
    "policy_ard.set_seed(seed)\n",
    "policy_ard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False)\n",
    "policy_ard.bayes_search(\n",
    "    max_num_probes=n_bayes,\n",
    "    simulator=simulator,\n",
    "    score=score,\n",
    "    ard=True,\n",
    "    is_disp=False,\n",
    ")\n",
    "\n",
    "best_fx_ard, best_actions_ard = policy_ard.history.export_sequence_best_fx()\n",
    "best_x_ard = test_X[best_actions_ard[-1], :]\n",
    "print(\"\\nBest value (ard=True):\", best_fx_ard[-1])\n",
    "print(\"Best point:\", best_x_ard)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case: ard=False\n",
    "\n",
    "Next, we run Bayesian optimization with the same settings but `ard=False`.\n",
    "The isotropic kernel (a single length scale for all dimensions) gives equal weight to irrelevant dimensions as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "============================================================\n",
      "Running with ard=False\n",
      "============================================================\n",
      "\n",
      "Best value (ard=False): -0.011203757603284015\n",
      "Best point: [-0.04345543 -0.04197482 -0.72416357 -0.13521863]\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n\" + \"=\" * 60)\n",
    "print(\"Running with ard=False\")\n",
    "print(\"=\" * 60)\n",
    "policy_noard = physbo.search.discrete.Policy(test_X)\n",
    "policy_noard.set_seed(seed)\n",
    "policy_noard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False)\n",
    "policy_noard.bayes_search(\n",
    "    max_num_probes=n_bayes,\n",
    "    simulator=simulator,\n",
    "    score=score,\n",
    "    ard=False,\n",
    "    is_disp=False,\n",
    ")\n",
    "\n",
    "best_fx_noard, best_actions_noard = policy_noard.history.export_sequence_best_fx()\n",
    "best_x_noard = test_X[best_actions_noard[-1], :]\n",
    "print(\"\\nBest value (ard=False):\", best_fx_noard[-1])\n",
    "print(\"Best point:\", best_x_noard)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing best values\n",
    "\n",
    "PHYSBO seeks the point that **maximizes** the objective.\n",
    "For this problem the maximum is 0, so the closer the best value found is to 0, the better the search has performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1119c16d0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGdCAYAAAD+JxxnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQG5JREFUeJzt3Ql4VOX1+PGTbQIhBAw7sgiFEkTcQBBEQIOC8PvXBS0oPhWloFZQllbBDZda/Ckqi7XUB4HyLwhSaVW0CEKBigERRQGBv7QoCITQYgIhkHX+z3mTGTNkksy9ySS5c7+f57nOdmfJzSVzfN9zzhvl9Xq9AgAA4ELRtf0BAAAAaguBEAAAcC0CIQAA4FoEQgAAwLUIhAAAgGsRCAEAANciEAIAAK5FIAQAAFwrtrY/QF1VVFQkR44ckYYNG0pUVFRtfxwAABAC7RN96tQpad26tURHVz7eQyBUDg2C2rZtG8oxBwAAdcyhQ4ekTZs2le5HIFQOHQnyHcikpKTq/e0AAICwOHnypBnI8H2PV4ZAqBy+6TANggiEAABwllDTWkiWBgAArkUgBAAAXItACAAAuBaBEAAAcC0CIQAA4FoEQgAAwLUIhAAAgGsRCAEAANciEAIAAK5VI4HQ73//e7ngggukXr160rt3b/n0008r3H/FihWSkpJi9u/evbt88MEHZRZUe/LJJ6VVq1ZSv359GTRokHzzzTcB+5w4cUJGjRplukI3btxYxowZI9nZ2WH5+QAAgDOFPRBavny5TJ48WaZPny6ff/65XHLJJTJ48GDJyMgIuv8nn3wit99+uwlcvvjiC7npppvMtmvXLv8+L7zwgsyZM0fmzZsnW7dulQYNGpjXPHv2rH8fDYJ2794ta9eulVWrVsmmTZtk3Lhx4f5xAQCAg0R5dXgljHQE6IorrpBXX33V3C4qKjKLoU2YMEGmTp1aZv8RI0bI6dOnTfDic+WVV8qll15qAh/9uK1bt5YpU6bIr3/9a/N4VlaWtGjRQhYtWiQjR46UPXv2yIUXXijbtm2Tnj17mn1Wr14tQ4cOle+//948P5RF2xo1amRem7XGAABwBqvf32FddDUvL0+2b98u06ZN898XHR1tprLS0tKCPkfv1xGk0nS0529/+5u5fuDAAUlPTzev4aM/sAZc+lwNhPRSp8N8QZDS/fW9dQTp5ptvLvO+ubm5Zit9IAG1/bsT8v5X6eKVsP4/AwC4zvDL28hF5zeq1c8Q1kDoP//5jxQWFprRmtL09t69e4M+R4OcYPvr/b7HffdVtE/z5s0DHo+NjZXk5GT/PueaMWOGPP3005Z/RkS246dy5Z5Fn0nWmfza/igAEHEua3deZAdCTqKjVqVHonRESKfw4G7PrvraBEGdmifK4G6BwTcAoGo6N0+U2hbWQKhp06YSExMjx44dC7hfb7ds2TLoc/T+ivb3Xep9WjVWeh/NI/Ltc24ydkFBgakkK+994+PjzQb4bNiXIe9+eUSio0Re/vklcnGbxhwcAIgwYa0a83g80qNHD1m3bp3/Pk2W1tt9+vQJ+hy9v/T+Siu/fPt36NDBBDOl99HRG8398e2jl5mZmSY/yWf9+vXmvTWXCKhMTl6BPP634krF0X07EAQBQIQK+9SYTjfdddddJnG5V69eMmvWLFMVdvfdd5vHf/GLX8j5559vcnTUQw89JAMGDJCXXnpJhg0bJsuWLZPPPvtMXn/9dfN4VFSUTJw4UX77299K586dTWD0xBNPmEowLbNXXbt2lSFDhsjYsWNNpVl+fr6MHz/eJFKHUjEGzP7oG/n+hzPSulE9mXL9TzkgABChwh4IaTn88ePHTQNETVTW6SstZfclOx88eNBUc/n07dtXli5dKo8//rg8+uijJtjRirGLLrrIv8/DDz9sgintC6QjP/369TOvqQ0YfZYsWWKCn9TUVPP6w4cPN72HgMrsPpIl8z8+YK4/c+NF0iCeVDoAiFRh7yPkVPQRcqfCIq/c/Npm+er7LBnavaW8NqpHbX8kAEAYv79ZawwoZXHatyYIahgfK9P/TzeODQBEOAIhoMSRzDMy88N95vojN6RIi6Qfp1oBAJGJQAjwLeT7zm45nVcoPdqfJ3f0asdxAQAXIBACROTD3eny0Z5jEhcTJTNu6S7R2jwIABDxCITgeifP5pvRIHVv/5/IT1s0dP0xAQC3IBCC6724ep9knMqVC5okyPhrO7n+eACAmxAIwdW2f/eD/Hnrd+b6727uLvXiYmr7IwEAahCBEFwrv7BIHl25U7ST1vDL20jfTk1r+yMBAGoYgRBc6/VN/5Z9x05JcgOPPDasa21/HABALSAQgit9+5/TMnvdN+b648O6mmAIAOA+BEJwZc+gx/62U/IKiuSqTk3k5svOr+2PBACoJQRCcJ2/fnFYNu//r8THRstzN3WXqCh6BgGAWxEIwVVOnM6TZ1d9ba4/mNpZLmjaoLY/EgCgFsXW5psDNe259/fIDzn50qVFQxnXvyO/gLrq8HaR9b8VOZtV258EQDhd85hIp1SpTQRCcI3N+/8jb3/+vehM2Izh3SUuhgHROunEAZElt4nk/Le2PwmAcDvzg9Q2AiG4wtn8QnnsrzvN9Tt7t5fL251X2x8JwZzJFFn68+IgqNWlIgOnigg5XEDEanVxbX8CAiG4w6vr98u3/82RFknx8pshXWr74yCYwnyRFaNF/vP/RJLOF7l9mUhSK44VgLBiRAghm7byK/loT4Yjj9h/s3PN5dM/6yZJ9eJq++PgXNre++8Pi/z7HyJxDQiCANQYAiGERHvuvPnpIUcfrRsuaimDu7Ws7Y+BYLb8QeSzBcXTYMPn14nhcgDuQCCEkJzOLfBfXzWhn8REOytvQz9vx6YN6BlUF+1bLfLho8XXr39WJGVobX8iAC5CIISQZJcEQtqE8KLzG3HUUD3Sd4m8PUbnxkQuv0ukz3iOLIAaRf0wQnI6rzgQSowndkY1OZUusnSESF62SIf+IsNeEtPbAABqEIEQQnI6t9BcNiAQQnXIyxF583aRk9+LNOks8vPFIjEksQOoeQRCsJQjlOCJ4YihaoqKRP52n8iRz0Xqnydyx/LiSwCoBQRCsBQIMTWGKvvHcyJfvyMSHScyYolIk59wUAHUGgIhWEqWZmoMVbLjTZF/ziy+/rM5IhdcxQEFUKsIhBCSnLziHCFGhGDbd5+IvDuh+Hq/ySKX3sHBBFDrCIRgaUSIHCHYcuLfIstGiRTli1x4o8i1T3AgAdQJBEKwlCPE1BjsLaQ6QuTMCZHWl4ncNE8kmj89AOoG/hohJCRLw/ZCqm/9InAhVU8CBxNAnUEghJBk00cIdhZS/eDXIgc2Fi+kqmXyDVnrDUDdQptghCTH31maPkK21tLKcvaCtbb85xuR7YuKF1K99Q2Rlt1r+xMBQBkEQrCYLM0pY8nRL0XeHOHus+z634p0uaG2PwUABMW3GkJCsrRNWd8XX9ZPLl5Py206DhDpcXdtfwoAKBeBECytNUYfIYtyTxVftrpE5Od/4mwDgDqGZGlY7CxNjpCtQCi+IWcaALgtEDpx4oSMGjVKkpKSpHHjxjJmzBjJzs6u8Dlnz56VBx54QJo0aSKJiYkyfPhwOXbsmP/xL7/8Um6//XZp27at1K9fX7p27SqzZ88OeI0NGzZIVFRUmS09PT1sP6t7kqUZRLQkr+R8JxACgDoprN9qGgQdPXpU1q5dK/n5+XL33XfLuHHjZOnSpeU+Z9KkSfL+++/LihUrpFGjRjJ+/Hi55ZZbZPPmzebx7du3S/PmzeXPf/6zCYY++eQT85oxMTFm39L27dtngjAffR6qNjWWQCBkDSNCAODOQGjPnj2yevVq2bZtm/Ts2dPcN3fuXBk6dKjMnDlTWrduXeY5WVlZ8sYbb5hA6dprrzX3LVy40Iz6bNmyRa688kq55557Ap7TsWNHSUtLk5UrV5YJhDTw0ZEoVE1eQZHkFRaZ64lUjVmTWzIi5EnkNAQAN02NaXCiQYgvCFKDBg2S6Oho2bp1a9Dn6GiPjhzpfj4pKSnSrl0783rl0QAqOTm5zP2XXnqptGrVSq677jr/iFJ5cnNz5eTJkwEbAivGFDlCdkeECIQAwFWBkObjnDsVFRsbawKW8nJ19H6Px1NmFKdFixblPkenxpYvX26mx3w0+Jk3b568/fbbZtMptIEDB8rnn39e7uedMWOGmYrzbfocFDtdkh8UHxstsTHk11uSR7I0ANRllr/Vpk6dGjQRufS2d+9eqQm7du2SG2+8UaZPny7XX3+9//4uXbrIvffeKz169JC+ffvKggULzOUrr7xS7mtNmzbNjCz5tkOHXNgJuJL8IBZcrcKIkIeqMQCIiByhKVOmyOjRoyvcR/N2WrZsKRkZGQH3FxQUmEoyfSwYvT8vL08yMzMDRoW0auzc53z99deSmppqRoIef/zxSj93r1695OOPPy738fj4eLOhLErnqyFHiKoxAIiMQKhZs2Zmq0yfPn1MQKN5Pzoyo9avXy9FRUXSu3fvoM/R/eLi4mTdunWmbN5X+XXw4EHzej67d+82ydR33XWXPPfccyF97h07dpgpM1ShqzSJ0lUonydHCABcVTWmlV5DhgyRsWPHmnwdTYLWqq6RI0f6K8YOHz5sRnUWL15sRmw0N0d7DU2ePNnkEmnp+4QJE0wQpBVjvukwDYIGDx5s9vPlDmn5vC9AmzVrlnTo0EG6detm+hLNnz/fBGFr1qwJ14/rikCIHkJVmRojEAIA1/URWrJkiQl+NNjRajEd5ZkzZ47/cQ2OdMQnJyfHf5/m8fj21UouDXhee+01/+N/+ctf5Pjx46aPkG4+7du3l2+//dZc1+k1ncLTQCshIUEuvvhi+eijj+Saa64J548bsU7nkSNU9amxH/tZAQDqjiiv1+ut7Q9RF2n5vI5QaeJ06aaMbvSnT76V6e/ulqHdW8pro4qnOREC/af1TLKIt0hkyj6RhsFz4wAAtff9TS00Qk+WJkfImvyc4iBIkSwNAHUSgRBCT5ZmeQ1702JR0SJxCZxpAFAHEQihUiRLV0MPoagozjQAqIMIhFApkqWr2lWaijEAqKsIhGBhaiyGo2UFpfMAUOcRCKFSJEvbRFdpAKjzCIRQKZKlbaKrNADUeQRCqFROSUNFOktblHuy+JLSeQCoswiEEPLUWAI5Qvamxlh5HgDqLAIhVIry+SomS1M1BgB1FoEQKnU6l7XGqpYj1JCzDADqKAIhVCivoEjyCouXiUhkiQ2bU2P0EQKAuopACBXKySvOD1LkCFlEsjQA1HkEQggpUdoTGy1xMZwuljA1BgB1Ht9sCCk/iNJ5G+gsDQB1HoEQQusqTem8dXSWBoA6j0AIoXWVJlHaOsrnAaDOIxBCSMnSTI1VZfX5JM4yAKijCIRQoeySHKGE+FiOlBVeL+XzAOAABEIIsat0DEfKivwzIt7iIJLO0gBQdxEIIbRkaXKE7JXOS5RIXAPOMgCoowiEEFKOUAOmxuyXzkfzzwwA6ir+QiPEdcaYGrNXMcY6YwBQlxEIIcQ+QiRL2+sqzTpjAFCXEQghxGRpAiFLGBECAEcgEEKFSJa2iZXnAcARCIRQoZw8X44QI0KWsPI8ADgCgRBCW2KDZGlrWHkeAByBQAgVIlnaJqbGAMARCIRQIZKlbSJZGgAcgUAIIfYRIkfI3oKrlM8DQF1GIIRy5RUUSV5hkbmeyBIbNjtL01ARAOoyAiFUuryGSiBZ2l6OEJ2lAaBOIxBCpYnSnthoiYvhVLGEztIA4Ah8u6HS/CC6SttAsjQAOEJYA6ETJ07IqFGjJCkpSRo3bixjxoyR7OySKYNynD17Vh544AFp0qSJJCYmyvDhw+XYsWMB+0RFRZXZli1bFrDPhg0b5PLLL5f4+Hjp1KmTLFq0KCw/YyQ77V95ngVX7ZfPkyMEAK4NhDQI2r17t6xdu1ZWrVolmzZtknHjxlX4nEmTJsl7770nK1askI0bN8qRI0fklltuKbPfwoUL5ejRo/7tpptu8j924MABGTZsmFxzzTWyY8cOmThxovzyl7+UDz/8MCw/Z8Q3UyRRugqdpakaA4C6LGw10Xv27JHVq1fLtm3bpGfPnua+uXPnytChQ2XmzJnSunXrMs/JysqSN954Q5YuXSrXXnutP+Dp2rWrbNmyRa688kr/vjrC1LJly6DvPW/ePOnQoYO89NJL5rY+/+OPP5ZXXnlFBg8eHKafOJK7SlM6b4nXS2dpAHD7iFBaWpoJVnxBkBo0aJBER0fL1q1bgz5n+/btkp+fb/bzSUlJkXbt2pnXK02nz5o2bSq9evWSBQsWiFe/fEq9d+nXUBoAnfsapeXm5srJkycDNrfLpoeQPQW5IkUlFXceRoQAoC4L2//qp6enS/PmzQPfLDZWkpOTzWPlPcfj8ZgAqrQWLVoEPOeZZ54xI0YJCQmyZs0a+dWvfmVyjx588EH/6+hzzn0NDW7OnDkj9evXL/PeM2bMkKeffrpKP3PkdpUmR8hWorQiEAKAyBoRmjp1atBk5dLb3r17JZyeeOIJueqqq+Syyy6TRx55RB5++GF58cUXq/Sa06ZNM1Nzvu3QoUPidv5kaXKE7HWV1iAomsJMAIioEaEpU6bI6NGjK9ynY8eOJn8nIyMj4P6CggJTSVZebo/en5eXJ5mZmQGjQlo1Vt5zVO/eveXZZ58101taJab7nltppre1ei3YaJDS5+mGH5EjVNWu0kyLAUDEBULNmjUzW2X69OljAhrN++nRo4e5b/369VJUVGQCl2B0v7i4OFm3bp0pm1f79u2TgwcPmtcrj1aGnXfeef5ARvf94IMPAvbRyrWKXgMVrTPG1JgldJUGAMcIW46QVmoNGTJExo4da6q4NAl6/PjxMnLkSH/F2OHDhyU1NVUWL15skp4bNWpkeg1NnjzZ5BLpCM6ECRNMAOOrGNPSeh3d0dv16tUzAc7vfvc7+fWvf+1/7/vuu09effVVM2V2zz33mADsrbfekvfffz9cP25Ed5amasxuM0VGhACgrgtrXfSSJUtM8KPBjlaL6SjPnDlz/I9rcKQjPjk5Of77tMTdt69OdWm112uvveZ/XEeMfv/735t+Q1opps0SX375ZRNw+WjpvAY9us/s2bOlTZs2Mn/+fErnbSdLUz5vb3kNmikCQF0X5S1ddw4/rTDTESpNnNaRKTf6xYJPZdP/Oy4v3XaJDO/RprY/jnN8tlBk1USRLsNEbl9a258GAFzlpMXvb0paEEKyNDlCljA1BgCOQSCEclE1ZhNTYwDgGARCKBfJ0jZRPg8AjkEghHLl5BWXz5MsbXdqjGRpAKjrCIRQ6YhQgoccIUuYGgMAxyAQQlD5hUWSV1BkrjMiZBEjQgDgGARCqDBRWtFQ0WZnaZbYAIA6j0AIFU6LeWKjJS6G08QSyucBwDH4hkNQJEpXw+rz8e5sxAkATkIghKBIlK4CpsYAwDEIhBAU64xVAcnSAOAYBEIIiq7SNhXkihTlF19n9XkAqPMIhBBUdm5xM0UqxmyOBimqxgCgziMQQlA5ecVVY4ksuGovEIprIBJNI0oAqOsIhFBJsnQsR8hWV+lEjhsAOACBEIIiWdomEqUBwFEIhBDUaX+OENM7llA6DwCOQiCEoKgasyn3ZPElK88DgCMQCCGo0yXJ0g3IEbKGlecBwFEIhBAU5fM2MTUGAI5CIIRKkqXJEbKEZGkAcBQCIQRFjpBNlM8DgKMQCKHiHKF4+gjZSpb2NOTMAgAHIBBCxeXzJEvbyxGiagwAHIFACBV2lqaPkN0cITpLA4ATEAihjPzCIskrKDLXE5kas4byeQBwFAIhlJsorcgRsojyeQBwFAIhlHE6rzg/yBMbLXExnCL2OksncWYBgAPwLYfyS+c99BCyjPJ5AHAUAiFUkChN6bztZGkPydIA4AQEQqigqzSBkCUFeSKFecXXKZ8HAEcgEEL5PYQIhOxNiylGhADAEQiEUO6IUAI5QvYSpeMSRGIYTQMAJyAQQrnLazA1ZhGl8wDgOARCKINkaZvoKg0AjhPWQOjEiRMyatQoSUpKksaNG8uYMWMkO7tUHkUQZ8+elQceeECaNGkiiYmJMnz4cDl27Jj/8UWLFklUVFTQLSMjw+yzYcOGoI+np6eH88eNGCRL20RXaQBwnLAGQhoE7d69W9auXSurVq2STZs2ybhx4yp8zqRJk+S9996TFStWyMaNG+XIkSNyyy23+B8fMWKEHD16NGAbPHiwDBgwQJo3bx7wWvv27QvY79zHUVmyNH2E7JXOs/I8ADhF2DI69+zZI6tXr5Zt27ZJz549zX1z586VoUOHysyZM6V169ZlnpOVlSVvvPGGLF26VK699lpz38KFC6Vr166yZcsWufLKK6V+/fpm8zl+/LisX7/ePO9cGvjoSBTsJkuT8GtvaoxACADE7SNCaWlpJgjxBUFq0KBBEh0dLVu3bg36nO3bt0t+fr7ZzyclJUXatWtnXi+YxYsXS0JCgtx6661lHrv00kulVatWct1118nmzZur5edyA5KlbaKrNAA4Ttj+l1/zcc6dioqNjZXk5ORyc3X0fo/HU2YUp0WLFuU+R0eC7rjjjoBRIg1+5s2bZ4Kw3NxcmT9/vgwcONAEYJdffnnQ19H9dPM5ebKkFNqFsukjZA9dpQEg8keEpk6dWm6ysm/bu3ev1AQdJdIpOE3CLq1Lly5y7733So8ePaRv376yYMECc/nKK6+U+1ozZsyQRo0a+be2bduKW/2YLE2OkCVMjQFA5I8ITZkyRUaPHl3hPh07dpSWLVv6q7h8CgoKTCWZPhaM3p+XlyeZmZkBo0JaNRbsOTrSo9NfGvBUplevXvLxxx+X+/i0adNk8uTJASNCbg2G/Iuu0lnaGqrGACDyA6FmzZqZrTJ9+vQxAY3m/fgCFU1qLioqkt69ewd9ju4XFxcn69atM2XzvsqvgwcPmtcrTcvw33rrLTOSE4odO3aYKbPyxMfHmw0/5giRLG0RI0IA4DhhyxHSSq8hQ4bI2LFjTb6OJkGPHz9eRo4c6a8YO3z4sKSmppqEZx2x0SkpnebSkRnNJdL+QxMmTDBBkFaMlbZ8+XIzwnTnnXeWee9Zs2ZJhw4dpFu3bqYvkY4caRC2Zs2acP24EVk+T2dpi+gsDQCOE9b66CVLlpjgR4MdrRbTUZ45c+b4H9fgSEd8cnJy/PdpHo9vX01e1h5Br732WtAkae0vFKw8XqfXdApPAy2tKLv44ovlo48+kmuuuSaMP20kdpYmR8gSOksDgONEeb1eb21/iLpIc4R0hEp7G+nIlFvkFxZJ58f+bq7vePI6aZzgqe2P5Bzz+omk7xS5822RTj+2gAAA1N3vb9YaQ4CckmkxRY6QRXSWBgDHIRBCgOySRGlPTLR4Yjk9bOUI0VkaAByDbzqUUzpPfpBldJYGAMchEEI5idKsM2ZJYb5Iwdni655EzioAcAgCIQTNEaJ03mZ+kGJqDAAcg0AIQUeEEjxMjdkKhGLricTEcVYBgEMQCCEAy2vYxPIaAOBIBEIIurwGU2MW0VUaAByJQAgBSJa2iXXGAMCRCIQQgGRpm/JKcoRIlAYARyEQQgCSpavaVZrSeQBwEgIhBCBZ2ia6SgOAIxEIIQDJ0jbRVRoAHIlACAFOlzRUpLO0Rbkniy/JEQIARyEQQvCpMRoq2iyfb8gZBQAOQiCEAJTPV7V8nmRpAHASAiEEzRFiaswiOksDgCMRCCFojhCdpS2iszQAOBKBEMopn2fRVXvJ0kmcUQDgIARC8CsoLJLcgiJzvYEnliNjBeXzAOBIBEIoMy2myBGyiM7SAOBIBELwyy5JlPbERIsnllPDEjpLA4Aj8W0HP/KDbCosECk4U3ydhooA4CgEQvBjnbEqrjyvWHQVAByFQAhll9cgUdretFhMvEishzMKAByEQAhBukpTOm8JXaUBwLEIhODH1JhNdJUGAMciEIJfTknVGF2lbTZTZMFVAHAcAiH4ZZfkCCWQI2QNpfMA4FgEQigzNZZIjpA1dJUGAMciEEKQZGmW17CErtIA4FgEQvAjWdompsYAwLEIhOCXk1ecI0SytN2V5xtyNgGAwxAIoczUWIKHPkKWUD4PAI5FIIQgydLkCNmaGmN5DQBwHAIh+JEsXdXO0kyNAYDThC0QOnHihIwaNUqSkpKkcePGMmbMGMnOLvk/53K8/vrrMnDgQPOcqKgoyczMtPW6X331lVx99dVSr149adu2rbzwwgvV/vNFotMlDRWpGrOI8nkAcKywBUIarOzevVvWrl0rq1atkk2bNsm4ceMqfE5OTo4MGTJEHn30Uduve/LkSbn++uulffv2sn37dnnxxRflqaeeMkEWKpZT0lCRqTGL6CwNAI4VlmSQPXv2yOrVq2Xbtm3Ss2dPc9/cuXNl6NChMnPmTGndunXQ502cONFcbtiwwfbrLlmyRPLy8mTBggXi8XikW7dusmPHDnn55ZcrDcTcjmRpmyifBwDHCsuIUFpampm28gUratCgQRIdHS1bt24N6+vqPv379zdBkM/gwYNl37598sMPP5T72rm5uWY0qfTmJgWFRZJbUGSuMyJkEVNjAOBYYQmE0tPTpXnz5gH3xcbGSnJysnksnK+rly1atAjYx3e7oveeMWOGNGrUyL9pbpGbnC6ZFlPkCFlEsjQAuCMQmjp1qklirmjbu3evONG0adMkKyvLvx06dEjcmCjtiYkWTyzFhCErKhTJzym+zurzABDZOUJTpkyR0aNHV7hPx44dpWXLlpKRkRFwf0FBgan40sfsCuV19fLYsWMB+/huV/Te8fHxZnN7D6EEFly1Nxqk4hOr95cCAKhbgVCzZs3MVpk+ffqY0net2urRo4e5b/369VJUVCS9e/e2/WFDeV3d57HHHpP8/HyJi4sz92mFWZcuXeS8886z/d6u6SHkoZmirfygGI9IrHsDaQBwqrDMgXTt2tWUwY8dO1Y+/fRT2bx5s4wfP15Gjhzprxg7fPiwpKSkmMd9NIdHK7z2799vbu/cudPc1hGfUF/3jjvuMInS2l9Iy+yXL18us2fPlsmTJ4fjR424HCESpS2iqzQAOFrYkkG0jF0DndTUVFPe3q9fv4BePjpio5Vc2jvIZ968eXLZZZeZQEdp9Zfefvfdd0N+XU10XrNmjRw4cMCMGul03pNPPknpfMhdpVlnzBISpQHA0aK8Xq+3tj9EXaTl8xpUaeK0drGOdCs//14mv/WlXN25qfzfMfanL13nX+tF/u/NIi0uErl/c21/GgBwvZMWv78pD0JAsnQDcoTsjQix4CoAOBKBEIzskhwheghZRFdpAHA0AiEEjAglkiNkM0eI0nkAcCICIZyTLE35vCV5vkCoIWcSADgQgRCMnJLO0gRCdsvnCYQAwIkIhBDQR6iBh/J5S5gaAwBHIxCCwdRYVVeeZ0QIAJyIQAjnJEuTI2QJ5fMA4GgEQjBO51E+bwudpQHA0QiEENhQkfJ5a5gaAwBHIxDCOYEQU2OWMCIEAI5GIITAZGmW2LCG1ecBwNEIhCAFhUWSW1BkjgTJ0hZRPg8AjkYgBH+itGJqzIKiIpH808XX4ytf4RgAUPcQCMGfHxQXEyWeWE4Jy4nSitXnAcCR+NYDidJVnRaLjhOJjedMAgAHIhACidJVLp1PFImK4kwCAAciEIJ/nTESpe12lWZ5DQBwKgIhyGn/yvMsuGoJPYQAwPEIhECOUHVMjQEAHIlACD8GQjRTtIYRIQBwPAIhSHZJjhA9hCyiqzQAOB6BECSnJEcokRwha3JPFl8yNQYAjkUgBH/5fAILrtrMEaKrNAA4FYEQ/DlClM9bxNQYADgegRD8fYQaeCift4RkaQBwPAIh/NhZmqkxayifBwDHIxBCqWTpWI6GnWRpOksDgGMRCMFfPk+ytM0coXiW2AAApyIQQqlkaXKE7OUI0VkaAJyKQAgssVHlHCFGhADAqQiE8GOyNEtsWEP5PAA4HoGQyxUUFkluQZG5TrK0BUVFInm+qTFGhADAqQiEXO50XnGitEogRyh0+ad/vE4gBACORSDkcr5E6biYKImPJVnacqJ0VIxIbL0w/XYAAI4NhE6cOCGjRo2SpKQkady4sYwZM0ays0uSS8vx+uuvy8CBA81zoqKiJDMzM+Dxb7/91rxOhw4dpH79+vKTn/xEpk+fLnl5eQH76HPP3bZs2RKuHzUiAiGaKVahdD4qqtp/LwCAmhG2DnoaBB09elTWrl0r+fn5cvfdd8u4ceNk6dKl5T4nJydHhgwZYrZp06aVeXzv3r1SVFQkf/zjH6VTp06ya9cuGTt2rJw+fVpmzpwZsO9HH30k3bp1899u0qRJNf+EkTU11oBEaWvIDwKAiBCWQGjPnj2yevVq2bZtm/Ts2dPcN3fuXBk6dKgJWFq3bh30eRMnTjSXGzZsCPq4L0jy6dixo+zbt0/+8Ic/lAmENPBp2bJlNf5UkT4ixLSYJawzBgARISxTY2lpaWY6zBcEqUGDBkl0dLRs3bq1Wt8rKytLkpOTy9z/s5/9TJo3by79+vWTd999t9LXyc3NlZMnTwZsbsA6YzZROg8AESEsgVB6eroJQkqLjY01AYs+Vl32799vRpruvfde/32JiYny0ksvyYoVK+T99983gdBNN91UaTA0Y8YMadSokX9r27atuKurNOuMWUJXaQBwXyA0derUoInIpTfN46kJhw8fNtNkt912m8kT8mnatKlMnjxZevfuLVdccYU8//zzcuedd8qLL75Y4etpTpKOLvm2Q4cOiaumxsgRsoau0gAQESwNA0yZMkVGjx5d4T6at6O5ORkZGQH3FxQUmEqy6sjbOXLkiFxzzTXSt29fU2lWGQ2KNGm7IvHx8WZzbbI0I0L2RoRYeR4A3BMINWvWzGyV6dOnjyl93759u/To0cPct379elPxpUFJVUeCNAjS1124cKHJO6rMjh07pFWrVlV630hFsrRNJEsDQEQIS2JI165dzbSVTlnNmzfPlM+PHz9eRo4c6a8Y04AmNTVVFi9eLL169TL3af6Qbpr7o3bu3CkNGzaUdu3amfwifY72GWrfvr2pEjt+/Lj/PX0jTX/605/E4/HIZZddZm6vXLlSFixYIPPnzw/Hj+p4JEtXdWqMlecBwMnCliG7ZMkSE/xosKOjNsOHD5c5c+b4H9fgSEvftXeQjwZNTz/9tP92//79zaWO/OiUnE5vaZCkW5s2bQLez+v1+q8/++yz8t1335kE7ZSUFFm+fLnceuut4fpRHY1k6apOjREIAYCTRXlLRxDw0/J5rR7TxGntdB2pHlj6ubz/1VF56v9cKKOv6lDbH8c5lo0S2btKZNjLIleMqe1PAwCw+f3NWmMuxxIbNlE1BgARgUDI5QiEbCJZGgAiAoGQy2XnUj5vC52lASAiEAi53I/J0qw1ZgmdpQEgIhAIuVxOnm/RVZbYsJcjFLmJ9ADgBgRCLufvI8QSG6HTQkvK5wEgIhAIuVhBYZGczS8y1xkRsiDvtEZDxdfjG4bldwMAqBkEQi7mW2dMNSBHyPq0WFS0SFz96v/FAABqDIGQi/kSpeNioiQ+lmRpWwuuRkWF6bcDAKgJBEIuRqK0TfQQAoCIQSDkYv4eQiRKW0PpPABEDAIhF/uxqzTTYpawvAYARAwCIRfzl87TQ8gaukoDQMQgEHIxX45QIoGQNbkniy8pnQcAxyMQcjFyhGxiagwAIgaBkIv5coQSyBGyhq7SABAxCIRc7McFV1lnzFaOEFNjAOB4BEIuRrJ0VafGEqvxtwEAqA0EQi6WU9JHiBEhi0iWBoCIQSDkYtklVWMJHvoI2SufZ8FVAHA6AiEX+7GhIjlCltBZGgAiBoGQi5EsbRPl8wAQMQiEXMzfR4gRIWvoLA0AEYNAyMV+7CxNjpC9qbGkMPxWAAA1iUDIxfwNFVl9PnRer0ieLxCifB4AnI5AyMV8fYQon7cgP0fEW1R83UMgBABORyDkUgWFRXI2v/gLnRwhG/lBEiXiaRCOXw0AoAYRCLlUTn5xorRqQI6QvYqxqKjq/8UAAGoUgZDL84PiYqIkPpZk6ZDRVRoAIgqBkEuRKG0TpfMAEFEIhFzeQ4hEaYvoKg0AEYVASNy+vAbTYpbQVRoAIgqBkEuxzlgVc4QonQeAiEAg5FKnS7pKN6CZor0cIbpKA0BEIBASt68zxtSYvakxmikCQCQIWyB04sQJGTVqlCQlJUnjxo1lzJgxkp3ta0YX3Ouvvy4DBw40z4mKipLMzMwy+1xwwQXmsdLb888/H7DPV199JVdffbXUq1dP2rZtKy+88EK1/3xOx9RYFZOlmRoDgIgQtkBIg6Ddu3fL2rVrZdWqVbJp0yYZN25chc/JycmRIUOGyKOPPlrhfs8884wcPXrUv02YMMH/2MmTJ+X666+X9u3by/bt2+XFF1+Up556ygRZKBsIUTVmd2qsIacTAESA2HC86J49e2T16tWybds26dmzp7lv7ty5MnToUJk5c6a0bt066PMmTpxoLjds2FDh6zds2FBatmwZ9LElS5ZIXl6eLFiwQDwej3Tr1k127NghL7/8cqWBmJuc9k+NheUUiFw0VASAiBKWEaG0tDQzHeYLgtSgQYMkOjpatm7dWuXX16mwJk2ayGWXXWZGfAoKCgLeu3///iYI8hk8eLDs27dPfvjhhyq/d8RNjXnIEbKE8nkAiChhGQ5IT0+X5s2bB75RbKwkJyebx6riwQcflMsvv9y81ieffCLTpk0z02M64uN77w4dOgQ8p0WLFv7HzjvvvKCvm5uba7bSU2yRLNtXNcaIkDV0lgYA944ITZ06tUyi8rnb3r17w/dpRWTy5Mkmofriiy+W++67T1566SUz7VY6iLFjxowZ0qhRI/+mSdaRjGTpqnaWJkcIAFw3IjRlyhQZPXp0hft07NjR5O9kZGQE3K/TV1pJVl5uj129e/c2r/3tt99Kly5dzOsfO3YsYB/f7YreW0eWNMgqPSIUycFQDkts2EP5PAC4NxBq1qyZ2SrTp08fU/quVVs9evQw961fv16KiopM4FKdNBFac498U3H63o899pjk5+dLXFycuU8r1zRIKm9aTMXHx5vNLbL9S2yQLG2vszQjQgAQCcKSLN21a1dTBj927Fj59NNPZfPmzTJ+/HgZOXKkv2Ls8OHDkpKSYh730RweDWz2799vbu/cudPc1pEkXyL0rFmz5Msvv5R///vfpkJs0qRJcuedd/qDnDvuuMMkSmvfIi3fX758ucyePTtgtAelO0uTLB0yr5fyeQCIMGHrI6RBigY6qamppmy+X79+Ab18dMRGK7m0d5DPvHnzTCWYBlBKq7/09rvvvmtu64jNsmXLZMCAAaYs/rnnnjOBUOnX1fyeNWvWyIEDB8xolE7nPfnkk5TOn4McIRsKzop4i9sO0FkaACJDlNer/5uLc2mOkAZVWVlZptN1pEl54u9yNr9I/vnwNdI2OaG2P44zZGeIzOys/2xEnjwhEs0KNQDg9O9v/pK7UGGR1wRBihwhm8trEAQBQEQgEHJxfpBKIEfIRuk8C64CQKQgEHJxflBsdJTEx3IKhIyu0gAQcfgWdHmitDbBRIjoKg0AEYdAyIWyaaZoD12lASDiEAi5UI5/RIgeQpbksbwGAEQaAiEX8nWVTvDQVdp21RgAICIQCLm4aiyR5TXs5Qix4CoARAwCIRfnCDE1ZhELrgJAxCEQcnWOEFNjthZcZUQIACIGgZCLy+eZGrNbPs/K8wAQKQiEXDw1RrK0RXSWBoCIQyDk6hEhyuctobM0AEQcAiEXyi6pGiNHyCLK5wEg4pAt60IkS1e1fD6pGn8bAJyisLBQ8vPza/tjuF5cXJzExFTfjAaBkAud9pXP01DRZmdpGioCbuL1eiU9PV0yMzNr+6OgROPGjaVly5bVsl4mgZCLO0vTR8gCr5e1xgCX8gVBzZs3l4SEBBarruWgNCcnRzIyMsztVq1aVfk1CYRciM7SNhTkihQVB5AssQG4azrMFwQ1adKktj8ORKR+/frmOGgwpL+Xqk6TkSzt4qoxkqVtJEor1hoDXMOXE6QjQag7fL+P6sjZIhBycY4QDRVt5AdpEBTNPxvAbaojFwV18/fBX3SXKSzyypl8X0NF+ghZ7ypNojQARBICIZfmBymmxux0lWZ5DQCIJARCLs0Pio2OkvhYfv0hY+V5AA6VlpZmEoqHDRsWcP+3335rpph8W3JysgwYMED++c9/Buz31FNP+feJjY2Vpk2bSv/+/WXWrFmSm5srTsc3oVt7CMXHMudtBV2lATjUG2+8IRMmTJBNmzbJkSNHyjz+0UcfydGjR83jrVu3lv/5n/+RY8eOBezTrVs3s8/BgwflH//4h9x2220yY8YM6du3r5w6VaqYxIEIhNxaMUZ+kM2pMbpKA3CO7OxsWb58udx///1mRGjRokVl9mnSpIlpTnjRRRfJo48+KidPnpStW7cG7KMjQbqPBkrdu3c3gdXGjRtl165d8r//+7/iZARCLkPpvE1MjQEo3dQvr6BWNn1vK9566y1JSUmRLl26yJ133ikLFiwo9zXOnDkjixcvNtc9Hk+lr62ve8MNN8jKlSsdfW7QUNG1XaX51VtCsjSAElp5e+GTH9bK8fj6mcGSYGF5JJ0W0wBIDRkyRLKyssxIzsCBA/379O3bV6Kjo03HZg2SevToIampqSG9vgZDa9asESdjRMhl6CptE+XzABxm37598umnn8rtt9/un94aMWKECY5KW758uXzxxRfy9ttvS6dOncz0mS5sGgoNnJzeY4lhAdcmS9NDyJLck8WXLLgKuF79uBgzMlNb7x0qDXgKCgpMXk/pwCU+Pl5effVV/31t27aVzp07m033v/nmm03uj+5XmT179kiHDh3EyRgRcm2yNDGwvRwhkqUBt9MREJ2eqo0t1NEXDWg03+ell16SHTt2+Lcvv/zSBEZvvvlm0OfdeuutZuTotddeq/Q99u7dK6tXr5bhw4eLk/Ft6DIkS9vE1BgAB1m1apX88MMPMmbMGGnUqFHAYxq46GiR5gydSwOtBx980PQOuvfee/1remlglZ6eLkVFRfLf//5XNmzYIL/97W/l0ksvld/85jfiZIwIuUx2qT5CsIBkaQAOooHOoEGDygRBvkDos88+M2Xywdx1111mMdPS02e7d++WVq1aSbt27UyitVajTZs2zTRfTEx09tJDfBu6dEQokRwhayifB+Ag7733XrmP9erVy19C7w1SSq+jQCdOnPDf1tEh3SIVI0IurRpjRMhmsrSHtcYAIJIQCLkMydJVzBFi0VUAiCgEQi5eawx2coScPRcOAKihQEjnF0eNGiVJSUnSuHFjk7mua55U5PXXXzdJWPoczVzPzMwMeFyz1EuvlFt627ZtW9DVdH3bli1bwvWjOrSzNH2EQlaQK1KUX3ydESEAiChhC4Q0CNIs87Vr15oyPl3Vdty4cRU+R9t7azmfLvoWjLYB19VvS2+//OUvTTOnnj17Bl1N17dpy3CIWatGJTIiZH1aTHkYEQKASBKW+RHtNKlNlnSUxhegzJ07V4YOHSozZ84M6HJZ2sSJE/0jP8HoInC6+q2Plve98847ZhXcc5tM+VbTRfDyeStr1bieL1E6LkEkmpE0AIgkYRkRSktLM9NhpUdptJ+BLuq2devWanufd9991zR2uvvuu8s89rOf/UyaN28u/fr1M/tVJjc31/RUKL1Fdvk8gZD10nkqxgAg0oQlENLukxqElKYtu5OTk81j1dkwavDgwdKmTRv/fdrYSVuKr1ixQt5//30TCN10002VBkMzZswwjad8m669EmkKi7xm1WRFjpCNRGmmxQDA3YHQ1KlTy01W9m269khN+P777+XDDz80SdilNW3aVCZPniy9e/eWK664Qp5//nm588475cUXX6zw9bRDZlZWln87dOiQRGoPIUXVmAWUzgNAxLI0PzJlyhQZPXp0hft07NjR5OZkZGQE3K/rlGglWXXl7SxcuNDkAekUWGU0KNKk7YroKruhrLTrZDkl+UGx0VESH0vnhJDl+UrnmRoDgEhj6duwWbNmkpKSUuGmCc19+vQxpe/bt2/3P3f9+vVmsTYNSqpKW4JrIPSLX/xC4uLiKt1fV9zVNVLczlc6n+CJCXkFY7DOGAAEo4uyxsTEmFSUc+mSHL6ZIt1H0020crz00h3qggsu8O9Xv359c/vnP/+5iRlqSliGBbp27WrK4MeOHSuffvqpbN68WcaPHy8jR470V4wdPnzYBE76uI/mD2nQsn//fnN7586d5va5B04P0IEDB0zp/Ln+9Kc/yZtvvmmm6HT73e9+JwsWLDCVZW5HorRNrDwPAGXa3Sxbtkwefvhh8x0bTLdu3Uz7moMHD5rBC60mv//++8vs98wzz5j99u3bJ4sXLzbFVlpg9dxzz0lNCFvp0JIlS0zwk5qaaqrFdLXbOXPmBJS+6w+tB9Nn3rx58vTTT/tv9+/f31zqASw9JadJ0tpTSAOpYJ599ln57rvvTIK27rN8+XK59dZbxe38y2tQMWYNXaUBOJA2KL744oulXr16Mn/+fDNjc9999wUsoHrw4EEzULBu3TrzXa2DGNrupkWLFhW+to4CXXjhhSZ3WAc4NK/23CIj/Q72pcOcf/75ctttt5nv83M1bNjQv5+ubq/f/TqL8+STT5rv7i5duogjAyGtEFu6dGm5j+vw17mr3oa6wm1Fr3vXXXeZDRV1laZ03hLK5wGUpt9d+T/+T3yN0n5mFlIbdJZEC4i0dY22ttFBhauuukquu+46k65y4403mmrrjRs3mlzeBx54QEaMGFFuP7/SAxJaiKRV1jfccIMsWrRInnjiiXL311UftMBJg7FQPPTQQ2ZQQ3sF6qhTOPGN6CI5ecXJ0vQQsls+T7I0AJ3SyBH5XfDGwGH36BERT4OQd9cRoenTp5vrnTt3lldffdWM/mggpJc7d+40qSa+0RydmtIpLW2IrJXXwXzzzTdm2aqVK1ea2xoQabD1+OOPB+Sf6mtrkFVYWChnz54197388sshD6ZoGx4NoMKN0iGXJkvDztQYgRAAZ9FAqDSdcvJVdesqEG3btg2Y0tLpLs3R0cfKozlB2sNP29UoXTVC286cm+CsU1qa56tB1SOPPGKeYyVfV2eNaqKwhxEhFyFZuqpTY6wzBqBkekpHZmrrva3sfk5ltQYWOiVmV2FhoZlu0+ImzQEqfb8GSJoX7KPTYJ06dTLXtaffsGHDTB6wTnlVRleNOH78uFlLNNwIhGrYtr/9XgqP7JDa0PpUrjwZmyspxxuK/L04kkcI0ncVX9JZGoDSUQoL01N1lVZ4Hzp0KCDR+euvvzbtb3RkKJgPPvhATp06JV988YUpi/fZtWuXWe5Kn6sjSsHo1Nm1115rKsfKW3PUZ/bs2SZ5W1eGCDcCoRoW9a91cuWpdVJr9Df+n5IN1jRkEV8AkUNL1Lt37y6jRo2SWbNmmWTpX/3qVzJgwICAtULPTZLWkZ1LLrkk4H4NnCZNmmQqxjXhOhjtMahTddrWRnOVfDSw0hEmrSbXfKU///nPpspNl77yjSiFE4FQDYtKGSppR2pvHTNPTLR0a50k9eLIE7KkURuRtlVvBgoAdYVOk73zzjsmb0dL1kuXzwdz7Ngxs4ZnsMptfe7NN99sAqXyAiGlwZJWrmnOkG8USsvkddOpNC2jv/LKK00i9zXXXCM1Icp7bg07DF19XssCNQEsKSmJowIALqTVTjpKobkq2o8Hdf/3YvX7m6oxAADgWgRCAADAtQiEAACAaxEIAQAA1yIQAgAArkUgBABAJarSjRl1+/dBHyEAAMqhvW20R86RI0ekWbNm5nZNrH+F4LTjT15enll+Q38voa5mXxECIQAAyqFfttqr5ujRoyYYQt2QkJAg7dq1M7+fqiIQAgCgAjrqoF+6ugSFLi6K2qVrnOmCr9U1MkcgBABAJfRLV1dyP3c1dzgfydIAAMC1CIQAAIBrEQgBAADXIkeoghI93yq2AADAGXzf277v8coQCJXj1KlT5rJt27bV9bsBAAA1+D3eqFGjSveL8oYaMrmwa6X2jGjYsGG1Ns/SSFWDq0OHDklSUlK1vW6k47hx3Djf6jb+jXLc6sr5pmGNBkGtW7cOqc8QI0Ll0IPXpk0bCRf9xREIcdxqCucbx41zrW7j32j1HrdQRoJ8SJYGAACuRSAEAABci0CohsXHx8v06dPNJThunG91E/9OOWaca+75N0qyNAAAcC1GhAAAgGsRCAEAANciEAIAAK5FIAQAAFyLQChMZsyYIVdccYXpTN28eXO56aabZN++fQH7nD17Vh544AFp0qSJJCYmyvDhw+XYsWPiVqEcs4EDB5pO36W3++67T9zsD3/4g1x88cX+xmJ9+vSRv//97/7HOc/sHTfOtco9//zz5t/gxIkTOd+qeNw438p66qmnyvy9T0lJqfa/bQRCYbJx40bzC9qyZYusXbtW8vPz5frrr5fTp0/795k0aZK89957smLFCrO/Lulxyy23iFuFcszU2LFj5ejRo/7thRdeEDfTDuj6h3X79u3y2WefybXXXis33nij7N692zzOeWbvuCnOtfJt27ZN/vjHP5pgsjTON3vHjfMtuG7dugX8vf/444+r/1zTtcYQfhkZGbqmm3fjxo3mdmZmpjcuLs67YsUK/z579uwx+6SlpfErCXLM1IABA7wPPfQQx6cS5513nnf+/PmcZzaPG+daxU6dOuXt3Lmzd+3atQH/Jvm7Zu+4cb4FN336dO8ll1wS9LHqPNcYEaohWVlZ5jI5Odlc6v+F6ojHoEGD/PvokF+7du0kLS2tpj6Wo46Zz5IlS6Rp06Zy0UUXybRp0yQnJ6eWPmHdU1hYKMuWLTOjaDrVw3lm77j5cK4FpyO3w4YNC/j7pTjf7B03zrfyffPNN2bx1I4dO8qoUaPk4MGD1X6usehqDa1kr3PBV111lfnyVunp6eLxeKRx48YB+7Zo0cI85nbBjpm64447pH379uYfxldffSWPPPKIySNauXKluNnOnTvNF7jOmetc+V//+le58MILZceOHZxnNo6b4lwLTgPGzz//3EzxnIu/a/aOG+dbcL1795ZFixZJly5dzLTY008/LVdffbXs2rWrWs81AqEa+r8A/cWVntuEvWM2btw4//Xu3btLq1atJDU1Vf71r3/JT37yE9ceVv1DoUGPjqL95S9/kbvuusvMmcPecdNgiHOtrEOHDslDDz1kcvjq1avH6VWNx43zrawbbrjBf11zqjQw0v8Rfuutt6R+/fpSXZgaC7Px48fLqlWr5B//+IdJzvRp2bKl5OXlSWZmZsD+mvGuj7lZeccsGP2Hofbv3y9upv9n1KlTJ+nRo4epvrvkkktk9uzZnGc2j1swnGvF0xEZGRly+eWXS2xsrNk0cJwzZ465rv83zt8168dNp2Y53yqnoz8//elPzd/76vwOJRAKE6/Xa77Qdah9/fr10qFDh4DH9Q9vXFycrFu3zn+fTvHo/GfpHAU3qeyYBaP/N690ZAiBU4u5ubmcZzaPG+dacDr6qtOJ+u/Ot/Xs2dPkbviu83fN+nGLiYnhb1sIsrOzzei//r2v1u9QS6nVCNn999/vbdSokXfDhg3eo0eP+recnBz/Pvfdd5+3Xbt23vXr13s/++wzb58+fczmVpUds/3793ufeeYZc6wOHDjgfeedd7wdO3b09u/f3+tmU6dONZV1eky++uorczsqKsq7Zs0a8zjnmfXjxrkWunOrnzjfrB83zrfgpkyZYr4P9N/o5s2bvYMGDfI2bdrUVBRX57lGIBQmGmMG2xYuXOjf58yZM95f/epXpmQ3ISHBe/PNN5svfreq7JgdPHjQBD3Jycne+Ph4b6dOnby/+c1vvFlZWV43u+eee7zt27f3ejweb7Nmzbypqan+IEhxnlk/bpxr9gMhzjfrx43zLbgRI0Z4W7VqZf6Nnn/++ea2Bo3Vfa5F6X+sjSEBAABEBnKEAACAaxEIAQAA1yIQAgAArkUgBAAAXItACAAAuBaBEAAAcC0CIQAA4FoEQgAAwLUIhAAAgGsRCAEAANciEAIAAK5FIAQAAMSt/j94XCAYAMfF0AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "st = n_initial\n",
    "en = n_initial + n_bayes\n",
    "steps = np.arange(st, en)\n",
    "plt.plot(steps, best_fx_ard[st:en], label=\"ARD\")\n",
    "plt.plot(steps, best_fx_noard[st:en], label=\"no ARD\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing length scales\n",
    "\n",
    "You can obtain the learned length scale(s) using the Policy method `get_kernel_length_scale()`.\n",
    "\n",
    "- **ard=True**: Returns one value per dimension. A smaller value indicates that the dimension is more relevant (more sensitive to the objective).\n",
    "- **ard=False**: Returns a single scalar (isotropic kernel)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Kernel length scale ---\n",
      "ard=True:  per-dimension (smaller = more relevant)\n",
      "  dim 0 (relevant): 2.5764\n",
      "  dim 1 (relevant): 4.2119\n",
      "  dim 2 (irrelevant): 5.9329\n",
      "  dim 3 (irrelevant): 6.0144\n",
      "ard=False: single (isotropic): 2.200175177775015\n"
     ]
    }
   ],
   "source": [
    "ls_ard = policy_ard.get_kernel_length_scale()\n",
    "ls_noard = policy_noard.get_kernel_length_scale()\n",
    "\n",
    "print(\"--- Kernel length scale ---\")\n",
    "if ls_ard is not None:\n",
    "    print(\"ard=True:  per-dimension (smaller = more relevant)\")\n",
    "    for i in range(len(ls_ard)):\n",
    "        rel = \"relevant\" if weights[i]>0 else \"irrelevant\"\n",
    "        print(f\"  dim {i} ({rel}): {ls_ard[i]:.4f}\")\n",
    "else:\n",
    "    print(\"ard=True:  (not available)\")\n",
    "if ls_noard is not None:\n",
    "    print(\"ard=False: single (isotropic):\", ls_noard[0])\n",
    "else:\n",
    "    print(\"ard=False: (not available)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing permutation importance\n",
    "\n",
    "Permutation importance measures how much the prediction MSE increases when each dimension’s values are shuffled; that reflects the “importance” of that dimension.\n",
    "A larger value means the dimension has a stronger effect on the prediction.\n",
    "The pattern of importance depends on whether ARD is used or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Permutation importance ---\n",
      "ard=True:\n",
      "  dim 0 (relevant): mean = 83.1582, std = 22.8684\n",
      "  dim 1 (relevant): mean = 1.4901, std = 0.3866\n",
      "  dim 2 (irrelevant): mean = 0.1697, std = 0.2088\n",
      "  dim 3 (irrelevant): mean = 0.2912, std = 0.1730\n",
      "ard=False:\n",
      "  dim 0 (relevant): mean = 67.0282, std = 18.2652\n",
      "  dim 1 (relevant): mean = 16.3441, std = 4.8668\n",
      "  dim 2 (irrelevant): mean = 7.3731, std = 3.3190\n",
      "  dim 3 (irrelevant): mean = 6.2181, std = 2.3696\n",
      "(Higher mean => more important.)\n"
     ]
    }
   ],
   "source": [
    "imp_mean_ard, imp_std_ard = policy_ard.get_permutation_importance(n_perm=n_perm)\n",
    "imp_mean_noard, imp_std_noard = policy_noard.get_permutation_importance(n_perm=n_perm)\n",
    "\n",
    "print(\"--- Permutation importance ---\")\n",
    "print(\"ard=True:\")\n",
    "for i in range(D):\n",
    "    rel = \"relevant\" if weights[i]>0 else \"irrelevant\"\n",
    "    print(f\"  dim {i} ({rel}): mean = {imp_mean_ard[i]:.4f}, std = {imp_std_ard[i]:.4f}\")\n",
    "print(\"ard=False:\")\n",
    "for i in range(D):\n",
    "    rel = \"relevant\" if weights[i]>0 else \"irrelevant\"\n",
    "    print(f\"  dim {i} ({rel}): mean = {imp_mean_noard[i]:.4f}, std = {imp_std_noard[i]:.4f}\")\n",
    "print(\"(Higher mean => more important.)\")"
   ]
  }
 ],
 "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.11.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}