diff --git a/docs/sphinx/manual/en/source/notebook/index.rst b/docs/sphinx/manual/en/source/notebook/index.rst index cb0582b0..cb7cc36c 100644 --- a/docs/sphinx/manual/en/source/notebook/index.rst +++ b/docs/sphinx/manual/en/source/notebook/index.rst @@ -9,6 +9,7 @@ Here, the usage of PHYSBO is introduced through tutorials. tutorial_basic tutorial_Gaussian_process + tutorial_ard tutorial_interactive_mode tutorial_once_mode tutorial_multi_probe diff --git a/docs/sphinx/manual/en/source/notebook/tutorial_ard.ipynb b/docs/sphinx/manual/en/source/notebook/tutorial_ard.ipynb new file mode 100644 index 00000000..4ab56ce4 --- /dev/null +++ b/docs/sphinx/manual/en/source/notebook/tutorial_ard.ipynb @@ -0,0 +1,379 @@ +{ + "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": [ + "" + ] + }, + "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": [ + "
" + ] + }, + "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 +} diff --git a/docs/sphinx/manual/ja/source/notebook/index.rst b/docs/sphinx/manual/ja/source/notebook/index.rst index a26f6dd2..fa8340c9 100644 --- a/docs/sphinx/manual/ja/source/notebook/index.rst +++ b/docs/sphinx/manual/ja/source/notebook/index.rst @@ -9,6 +9,7 @@ tutorial_basic tutorial_Gaussian_process + tutorial_ard tutorial_interactive_mode tutorial_once_mode tutorial_multi_probe diff --git a/docs/sphinx/manual/ja/source/notebook/tutorial_ard.ipynb b/docs/sphinx/manual/ja/source/notebook/tutorial_ard.ipynb new file mode 100644 index 00000000..abbb17b2 --- /dev/null +++ b/docs/sphinx/manual/ja/source/notebook/tutorial_ard.ipynb @@ -0,0 +1,379 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# ARD (Automatic Relevance Determination) チュートリアル\n", + "\n", + "PHYSBOではサロゲートモデルにガウス過程 (GP) を用いており、カーネルとしてガウスカーネルを用います。\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", + "$\\lambda$ と $\\eta$ とはそれぞれカーネルの値と入力の長さスケールとを表すハイパーパラメータです。\n", + "\n", + "**ARD** カーネルでは、各入力次元 $i$ ごとに異なる長さスケール $\\eta_i$ を設定します。これらを別々に学習することで、目的関数に効く次元と効かない次元を自動的に識別し、効率的な最適化を行えます。\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", + "PHYSBOでは、 `policy.bayes_search` 関数でベイズ最適化を行うときに `ard=True` というキーワード引数を与えることでARDカーネルを利用できます。\n", + "\n", + "本チュートリアルでは、一部の次元だけが目的関数に効く問題を例に、ARD の有無による違いを比較します。" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## パッケージのインポート\n", + "\n", + "チュートリアルで使うパッケージをインポートします。" + ] + }, + { + "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": [ + "## このチュートリアルの問題設定\n", + "\n", + "入力の重みつき二乗和\n", + "\n", + "$$ f(\\vec{x}; \\vec{w}) = - \\sum_i w_i x_i^2 $$\n", + "\n", + "を目的関数とします。ここで重み $\\vec{w}$ を $ \\vec{w} = [5, 1, 0, 0]$ とすると、探索パラメータは $D=4$ 次元で、 $x_0$ の寄与が一番大きく、ついで$x_1$が寄与し、$x_2, x_3$ は全く目的関数に寄与しないパラメータとなります。\n", + "\n", + "探索候補点として、最大値を取る点 $\\vec{x} = \\vec{0}$ と、$D$ 次元正規分布から生成した $N-1 = 999$ 個の点を準備します。" + ] + }, + { + "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": [ + "## 共通パラメータの設定\n", + "\n", + "ARDの有無による効果を公平に比較するため、両方とも同じシード・同じステップ数で実行するようにパラメータを設定します。" + ] + }, + { + "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": [ + "## ard=True の場合\n", + "\n", + "まず、`ard=True` でベイズ最適化を実行します。\n", + "次元ごとに長さスケールを学習するため、関連する次元と無関係な次元を識別しやすくなります。" + ] + }, + { + "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": [ + "## ard=False の場合\n", + "\n", + "次に、`ard=False` で同じ設定のベイズ最適化を実行します。\n", + "等方カーネル（全次元で 1 つの長さスケール）を用いるため、無関係な次元の影響も同じ重みで入ります。" + ] + }, + { + "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": [ + "## 最良値の比較\n", + "\n", + "PHYSBOでは目的関数値が最大となる点を求める仕様です。\n", + "この問題では最大値が 0 であるため、見つかった最良値が 0 に近いほど、良い探索ができていることになります。" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "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": [ + "
" + ] + }, + "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": [ + "## 長さスケールの比較\n", + "\n", + "Policy の `get_kernel_length_scale()` メソッドを用いると、学習された長さスケールを取得できます。\n", + "\n", + "- **ard=True**: 次元ごとの値が返されます。小さいほどで重要であることを示します。\n", + "- **ard=False**: 1 つのスカラー（等方カーネル）が返されます。" + ] + }, + { + "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": [ + "## Permutation importance の比較\n", + "\n", + "各次元の値をシャッフルしたときに、予測の MSE がどれだけ増えるかで、その次元の「重要度」を測ります。\n", + "値が大きいほど、その次元が予測に効いていることを示します。\n", + "ARD の有無によって、重要度の出方も変わります。" + ] + }, + { + "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 +} diff --git a/examples/single_objective/ard_example.py b/examples/single_objective/ard_example.py new file mode 100644 index 00000000..31cabcb4 --- /dev/null +++ b/examples/single_objective/ard_example.py @@ -0,0 +1,118 @@ +# SPDX-License-Identifier: MPL-2.0 +# Copyright (C) 2020- The University of Tokyo +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. + +""" +Example: When ARD (Automatic Relevance Determination) is useful + +ARD is effective when the objective function depends on only a subset of +input dimensions. The GP kernel learns a separate length scale per dimension; +relevant dimensions get smaller length scales, irrelevant ones get larger +(and are effectively down-weighted). +""" + +import numpy as np +import physbo + +weights = np.array([5.0, 1.0, 0.0, 0.0]) +D = len(weights) +N = 1000 +M = 0 +np.random.seed(137) +test_X = np.random.randn(N, D) +test_X[0, :] = 0.0 + + +def simulator(actions: np.ndarray) -> np.ndarray: + """Objective function: f = -Σ_i w_i x_i^2.""" + X2 = test_X[actions, :] ** 2 + return -np.einsum("ai,i -> a", X2, weights) + + +# ---------- Run with ARD=True ---------- +n_initial = 20 +n_bayes = 30 +n_perm = 20 +score = "EI" +seed = 31415 + +print("=" * 60) +print("Running with ard=True") +print("=" * 60) +policy_ard = physbo.search.discrete.Policy(test_X) +policy_ard.set_seed(seed) +policy_ard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False) +policy_ard.bayes_search( + max_num_probes=n_bayes, + simulator=simulator, + score=score, + ard=True, + is_disp=False, + num_rand_basis=M, +) + +best_fx_ard, best_actions_ard = policy_ard.history.export_sequence_best_fx() +best_x_ard = test_X[best_actions_ard[-1], :] +print("\nBest value (ard=True):", best_fx_ard[-1]) +print("Best point:", best_x_ard) + +# ---------- Run with ARD=False ---------- +print("\n" + "=" * 60) +print("Running with ard=False") +print("=" * 60) +policy_noard = physbo.search.discrete.Policy(test_X) +policy_noard.set_seed(seed) +policy_noard.random_search(max_num_probes=n_initial, simulator=simulator, is_disp=False) +policy_noard.bayes_search( + max_num_probes=n_bayes, + simulator=simulator, + score=score, + ard=False, + is_disp=False, +) + +best_fx_noard, best_actions_noard = policy_noard.history.export_sequence_best_fx() +best_x_noard = test_X[best_actions_noard[-1], :] +print("\nBest value (ard=False):", best_fx_noard[-1]) +print("Best point:", best_x_noard) + +# ---------- Compare best values ---------- +print("\n" + "=" * 60) +print("Comparison: ard=True vs ard=False") +print("=" * 60) +print(f" Best value (ard=True): {best_fx_ard[-1]:.6f}") +print(f" Best value (ard=False): {best_fx_noard[-1]:.6f}") +print(f" (Higher is better; weights are {weights}).") + +# ---------- Compare kernel length scales ---------- +print("\n--- Kernel length scale ---") +ls_ard = policy_ard.get_kernel_length_scale() +ls_noard = policy_noard.get_kernel_length_scale() +if ls_ard is not None: + print(" ard=True: per-dimension length scale (smaller = more relevant)") + for i in range(len(ls_ard)): + rel = "relevant" if i < 2 else "irrelevant" + print(f" dim {i} ({rel}): {ls_ard[i]:.4f}") +else: + print(" ard=True: (not available)") +if ls_noard is not None: + print(" ard=False: single length scale (isotropic kernel):", ls_noard[0]) +else: + print(" ard=False: (not available)") + +# ---------- Permutation importance (both) ---------- +print("\n--- Permutation importance ---") +imp_mean_ard, imp_std_ard = policy_ard.get_permutation_importance(n_perm=n_perm) +imp_mean_noard, imp_std_noard = policy_noard.get_permutation_importance(n_perm=n_perm) +print(" ard=True:") +for i in range(D): + rel = "relevant" if i < 2 else "irrelevant" + print(f" dim {i} ({rel}): mean = {imp_mean_ard[i]:.4f}, std = {imp_std_ard[i]:.4f}") +print(" ard=False:") +for i in range(D): + rel = "relevant" if i < 2 else "irrelevant" + print(f" dim {i} ({rel}): mean = {imp_mean_noard[i]:.4f}, std = {imp_std_noard[i]:.4f}") +print(" (Higher mean => more important; both use the same GP predictor after optimization.)") diff --git a/src/physbo/gp/core/_model.py b/src/physbo/gp/core/_model.py index fcdacd7a..571118e2 100644 --- a/src/physbo/gp/core/_model.py +++ b/src/physbo/gp/core/_model.py @@ -8,7 +8,7 @@ import numpy as np from ... import blm -from .. import inf +from .. import cov, inf, lik, mean from . import learning from ._prior import Prior @@ -33,6 +33,35 @@ def __init__(self, lik, mean, cov, inf="exact"): self.params = self.cat_params(self.lik.params, self.prior.params) self.stats = () + @classmethod + def create_default(cls, ard=False, num_dim=None): + """ + Create a default GP model with Gaussian kernel, constant mean, and Gaussian likelihood. + + Parameters + ---------- + ard : bool + If True, use ARD (per-dimension length scale) for the Gaussian kernel. + num_dim : int or None + Input dimension. Required when ard is True; optional otherwise. + + Returns + ------- + Model + A new Model instance. + """ + if ard and num_dim is None: + raise ValueError( + "ard=True requires num_dim. " + "Provide the input dimension (e.g. from training data)." + ) + cov_kernel = cov.Gauss(num_dim=num_dim, ard=ard) + return cls( + lik=lik.Gauss(), + mean=mean.Const(), + cov=cov_kernel, + ) + def cat_params(self, lik_params, prior_params): """ Concatinate the likelihood and prior parameters diff --git a/src/physbo/search/discrete/_policy.py b/src/physbo/search/discrete/_policy.py index b812130e..81ab05af 100644 --- a/src/physbo/search/discrete/_policy.py +++ b/src/physbo/search/discrete/_policy.py @@ -15,6 +15,7 @@ from ._history import History from .. import utility from .. import score as search_score +from ... import gp from ...gp import Predictor as gp_predictor from ...blm import Predictor as blm_predictor from ...misc import SetConfig @@ -105,6 +106,7 @@ def __init__(self, test_X, config=None, initial_data=None, comm=None): else: self.config = config + self.ard = False if initial_data is not None: if len(initial_data) != 2: msg = "ERROR: initial_data should be 2-elements tuple or list (actions and objectives)" @@ -311,6 +313,7 @@ def bayes_search( num_rand_basis=0, optimizer=None, unify_method=None, + ard=False, ): """ Performing Bayesian optimization. @@ -327,7 +330,7 @@ def bayes_search( Base class is defined in physbo.predictor. If None, blm_predictor is defined. is_disp: bool - If true, process messages are outputted. + If true, process messages are outputted. simulator: callable Callable (function or object with ``__call__``) Here, action is an integer which represents the index of the candidate. @@ -343,6 +346,9 @@ def bayes_search( optimizer: optimizer object, optional This is for compatibility with the range-based Policies. This is not used. + ard: bool + If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel. + Default is False. Returns ------- @@ -360,6 +366,7 @@ def bayes_search( simulator = None is_rand_expans = num_rand_basis != 0 + self.ard = ard if training is not None: self.training = training @@ -450,7 +457,7 @@ def get_post_fmean(self, xs): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fmean()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm) predictor.prepare(self.training) return predictor.get_post_fmean(self.training, X) @@ -479,7 +486,7 @@ def get_post_fcov(self, xs, diag=True): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fcov()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm) predictor.prepare(self.training) return predictor.get_post_fcov(self.training, X, diag) @@ -487,6 +494,47 @@ def get_post_fcov(self, xs, diag=True): self._update_predictor() return self.predictor.get_post_fcov(self.training, X, diag) + def get_kernel_length_scale(self): + """ + Return the Gaussian kernel length scale(s) (width) of the predictor. + + With ARD, returns one length scale per input dimension; otherwise a single value. + + Returns + ------- + numpy.ndarray or None + Length scale(s). Shape (num_dim,) when ARD is used, (1,) otherwise. + None if the predictor is not set or not a GP with Gaussian kernel. + """ + if self.predictor is None: + return None + self._update_predictor() + try: + cov = self.predictor.model.prior.cov + except AttributeError: + return None + if not hasattr(cov, "width"): + return None + return np.atleast_1d(np.asarray(cov.width).flatten()) + + def get_num_dim(self): + """ + Return the input dimension (number of features) of the search space. + + Returns + ------- + int or None + The number of dimensions. None if neither training nor test data is set. + """ + if ( + self.training.X is not None + and self.training.X.shape[0] > 0 + ): + return self.training.X.shape[1] + if self.test.X is not None and self.test.X.shape[0] > 0: + return self.test.X.shape[1] + return None + def get_permutation_importance(self, n_perm: int, split_features_parallel=False): """ Calculating permutation importance of model @@ -508,7 +556,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False) if self.predictor is None: self._warn_no_predictor("get_post_fmean()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0) predictor.prepare(self.training) return predictor.get_permutation_importance( @@ -588,7 +636,7 @@ def get_score( if predictor is None: if self.predictor is None: self._warn_no_predictor("get_score()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm) predictor.prepare(training) else: @@ -899,9 +947,23 @@ def _init_predictor(self, is_rand_expans): If false, physbo.gp.Predictor is selected. """ if is_rand_expans: - self.predictor = blm_predictor(self.config) + self.predictor = self._make_blm_predictor() else: - self.predictor = gp_predictor(self.config) + self.predictor = self._make_gp_predictor() + + def _make_gp_predictor(self): + """Create a GP predictor, with ARD if self.ard is True.""" + ard = self.ard + num_dim = self.get_num_dim() + model = gp.core.Model.create_default(ard=ard, num_dim=num_dim) + return gp_predictor(self.config, model=model) + + def _make_blm_predictor(self): + """Create a BLM predictor, with ARD if self.ard is True.""" + ard = self.ard + num_dim = self.get_num_dim() + model = gp.core.Model.create_default(ard=ard, num_dim=num_dim) + return blm_predictor(self.config, model=model) def _learn_hyperparameter(self, num_rand_basis): self.predictor.fit(self.training, num_rand_basis, comm=self.mpicomm) diff --git a/src/physbo/search/discrete_multi/_policy.py b/src/physbo/search/discrete_multi/_policy.py index 22331355..f96821ec 100644 --- a/src/physbo/search/discrete_multi/_policy.py +++ b/src/physbo/search/discrete_multi/_policy.py @@ -14,8 +14,6 @@ from .. import discrete from .. import utility from .. import score_multi as search_score -from ...gp import Predictor as gp_predictor -from ...blm import Predictor as blm_predictor from ...misc import SetConfig from ..._variable import Variable @@ -42,6 +40,7 @@ def __init__( self.TS_candidate_num = None + self.ard = False if initial_data is not None: if len(initial_data) != 2: msg = "ERROR: initial_data should be 2-elements tuple or list (actions and objectives)" @@ -197,6 +196,7 @@ def bayes_search( num_rand_basis=0, optimizer=None, unify_method=None, + ard=False, ): """ Performing Bayesian optimization by using multi objective function @@ -230,6 +230,9 @@ def bayes_search( This is for compatibility with the range-based Policies. unify_method: callable, optional This is for compatibility with the unified-optimization Policies. + ard: bool, optional + If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel. + Default is False. Returns ------- @@ -246,6 +249,7 @@ def bayes_search( simulator = None is_rand_expans = False if num_rand_basis == 0 else True + self.ard = ard if training_list is not None: self.training = training_list @@ -253,11 +257,11 @@ def bayes_search( if predictor_list is None: if is_rand_expans: self.predictor_list = [ - blm_predictor(self.config) for i in range(self.num_objectives) + self._make_blm_predictor() for i in range(self.num_objectives) ] else: self.predictor_list = [ - gp_predictor(self.config) for i in range(self.num_objectives) + self._make_gp_predictor() for i in range(self.num_objectives) ] else: self.predictor_list = predictor_list @@ -359,7 +363,7 @@ def get_post_fmean(self, xs): self._warn_no_predictor("get_post_fmean()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) @@ -373,6 +377,35 @@ def get_post_fmean(self, xs): ] return np.array(fmean).T + def get_kernel_length_scale(self): + """ + Return the Gaussian kernel length scale(s) for each objective. + + Returns + ------- + list of numpy.ndarray or None + List of length num_objectives; each element is the length scale(s) + for that objective's GP, or None if not available. + """ + if self.predictor_list is None or all(p is None for p in self.predictor_list): + return None + self._update_predictor() + result = [] + for predictor in self.predictor_list: + if predictor is None: + result.append(None) + continue + try: + cov = predictor.model.prior.cov + except AttributeError: + result.append(None) + continue + if not hasattr(cov, "width"): + result.append(None) + continue + result.append(np.atleast_1d(np.asarray(cov.width).flatten())) + return result + def get_post_fcov(self, xs, diag=True): """ Calculate covariance of predictors (post distribution) @@ -395,7 +428,7 @@ def get_post_fcov(self, xs, diag=True): self._warn_no_predictor("get_post_fcov()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) @@ -441,7 +474,7 @@ def get_score( self._warn_no_predictor("get_score()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(training, objective_index=i) predictor_list.append(predictor) @@ -506,7 +539,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False) self._warn_no_predictor("get_post_fmean()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) diff --git a/src/physbo/search/discrete_unified/_policy.py b/src/physbo/search/discrete_unified/_policy.py index f9247f79..fcdc233a 100644 --- a/src/physbo/search/discrete_unified/_policy.py +++ b/src/physbo/search/discrete_unified/_policy.py @@ -14,8 +14,6 @@ from .. import discrete from .. import utility from .. import score as search_score -from ...gp import Predictor as gp_predictor -from ...blm import Predictor as blm_predictor from ...misc import SetConfig from ..._variable import Variable, normalize_t @@ -56,6 +54,8 @@ def __init__( self.config = SetConfig() else: self.config = config + + self.ard = False if initial_data is not None: if len(initial_data) != 2: @@ -203,6 +203,7 @@ def bayes_search( num_rand_basis=0, unify_method=None, optimizer=None, + ard=False, ): """ Performing Bayesian optimization by using unified objective function @@ -238,6 +239,9 @@ def bayes_search( optimizer: Optimizer object, optional This is for compatibility with the range-based Policies. This is not used. + ard: bool, optional + If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel. + Default is False. Returns ------- @@ -259,6 +263,7 @@ def bayes_search( if num_rand_basis < 0: raise ValueError("num_rand_basis must be non-negative") is_rand_expans = (num_rand_basis > 0) + self.ard = ard if training_list is not None: self.training = training_list @@ -362,7 +367,7 @@ def get_post_fmean(self, xs): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fmean()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training_unified, 0, comm=self.mpicomm) predictor.prepare(self.training_unified) return predictor.get_post_fmean(self.training_unified, X) @@ -391,7 +396,7 @@ def get_post_fcov(self, xs, diag=True): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fcov()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training_unified, 0, comm=self.mpicomm) predictor.prepare(self.training_unified) return predictor.get_post_fcov(self.training_unified, X, diag) @@ -462,7 +467,7 @@ def get_score( if predictor is None: if self.predictor is None: self._warn_no_predictor("get_score()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm) predictor.prepare(training) else: @@ -516,7 +521,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False) if self.predictor is None: self._warn_no_predictor("get_permutation_importance()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training_unified, 0) predictor.prepare(self.training_unified) return predictor.get_permutation_importance( @@ -676,9 +681,9 @@ def _learn_hyperparameter(self, num_rand_basis): def _initialize_predictor(self, is_rand_expans): if is_rand_expans: - self.predictor = blm_predictor(self.config) + self.predictor = self._make_blm_predictor() else: - self.predictor = gp_predictor(self.config) + self.predictor = self._make_gp_predictor() def _update_predictor(self): if self.new_data is not None: diff --git a/src/physbo/search/range/_policy.py b/src/physbo/search/range/_policy.py index 6161a571..31699034 100644 --- a/src/physbo/search/range/_policy.py +++ b/src/physbo/search/range/_policy.py @@ -14,6 +14,7 @@ from .. import utility from .. import score as search_score from ..optimize.random import Optimizer as RandomOptimizer +from ... import gp from ...gp import Predictor as gp_predictor from ...blm import Predictor as blm_predictor from ...misc import SetConfig @@ -64,6 +65,7 @@ def __init__( else: self.config = config + self.ard = False if initial_data is not None: if len(initial_data) != 2: msg = "ERROR: initial_data should be 2-elements tuple or list (X and objectives)" @@ -244,6 +246,7 @@ def bayes_search( interval=0, num_rand_basis=0, optimizer=None, + ard=False, ): """ Performing Bayesian optimization. @@ -276,6 +279,9 @@ def bayes_search( optimizer: Optimizer object Optimizer object for optimizing the acquisition function. If None, the default optimizer is used. + ard: bool + If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel. + Default is False. Returns ------- @@ -293,6 +299,7 @@ def bayes_search( simulator = None is_rand_expans = num_rand_basis != 0 + self.ard = ard if training is not None: self.training = training @@ -384,7 +391,7 @@ def get_post_fmean(self, xs): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fmean()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=0) predictor.prepare(self.training, objective_index=0) return predictor.get_post_fmean(self.training, X, objective_index=0) @@ -413,7 +420,7 @@ def get_post_fcov(self, xs, diag=True): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fcov()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=0) predictor.prepare(self.training, objective_index=0) return predictor.get_post_fcov(self.training, X, diag, objective_index=0) @@ -423,6 +430,40 @@ def get_post_fcov(self, xs, diag=True): self.training, X, diag, objective_index=0 ) + def get_kernel_length_scale(self): + """ + Return the Gaussian kernel length scale(s) (width) of the predictor. + + With ARD, returns one length scale per input dimension; otherwise a single value. + + Returns + ------- + numpy.ndarray or None + Length scale(s). Shape (num_dim,) when ARD is used, (1,) otherwise. + None if the predictor is not set or not a GP with Gaussian kernel. + """ + if self.predictor is None: + return None + self._update_predictor() + try: + cov = self.predictor.model.prior.cov + except AttributeError: + return None + if not hasattr(cov, "width"): + return None + return np.atleast_1d(np.asarray(cov.width).flatten()) + + def get_num_dim(self): + """ + Return the input dimension (number of features) of the search space. + + Returns + ------- + int + The number of dimensions (from min_X / max_X). + """ + return self.dim + def get_score( self, mode, *, xs=None, predictor=None, training=None, parallel=True, alpha=1 ): @@ -475,7 +516,7 @@ def get_score( if predictor is None: if self.predictor is None: self._warn_no_predictor("get_score()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm, objective_index=0) predictor.prepare(training, objective_index=0) else: @@ -709,9 +750,23 @@ def _init_predictor(self, is_rand_expans): If false, physbo.gp.Predictor is selected. """ if is_rand_expans: - self.predictor = blm_predictor(self.config) + self.predictor = self._make_blm_predictor() else: - self.predictor = gp_predictor(self.config) + self.predictor = self._make_gp_predictor() + + def _make_gp_predictor(self): + """Create a GP predictor, with ARD if self.ard is True.""" + ard = self.ard + num_dim = self.get_num_dim() + model = gp.core.Model.create_default(ard=ard, num_dim=num_dim) + return gp_predictor(self.config, model=model) + + def _make_blm_predictor(self): + """Create a BLM predictor, with ARD if self.ard is True.""" + ard = self.ard + num_dim = self.get_num_dim() + model = gp.core.Model.create_default(ard=ard, num_dim=num_dim) + return blm_predictor(self.config, model=model) def _learn_hyperparameter(self, num_rand_basis): self.predictor.fit( diff --git a/src/physbo/search/range_multi/_policy.py b/src/physbo/search/range_multi/_policy.py index 21e3ae1c..4ea65d98 100644 --- a/src/physbo/search/range_multi/_policy.py +++ b/src/physbo/search/range_multi/_policy.py @@ -15,8 +15,6 @@ from .. import utility from .. import score_multi as search_score from ..optimize.random import Optimizer as RandomOptimizer -from ...gp import Predictor as gp_predictor -from ...blm import Predictor as blm_predictor from ...misc import SetConfig from ..._variable import Variable @@ -55,6 +53,7 @@ def __init__( self.TS_candidate_num = None + self.ard = False if initial_data is not None: if len(initial_data) != 2: msg = "ERROR: initial_data should be 2-elements tuple or list (actions and objectives)" @@ -202,6 +201,7 @@ def bayes_search( num_rand_basis=0, optimizer=None, unify_method=None, + ard=False, ): """ Performing Bayesian optimization by using multi objective function @@ -237,6 +237,9 @@ def bayes_search( unify_method: callable, optional This is for compatibility with the unified-optimization Policies. This is not used. + ard: bool, optional + If True, use Automatic Relevance Determination (ARD) for the Gaussian kernel. + Default is False. Returns ------- @@ -254,6 +257,7 @@ def bayes_search( simulator = None is_rand_expans = False if num_rand_basis == 0 else True + self.ard = ard if training_list is not None: self.training = training_list @@ -261,11 +265,12 @@ def bayes_search( if predictor_list is None: if is_rand_expans: self.predictor_list = [ - blm_predictor(self.config) for i in range(self.num_objectives) + self._make_blm_predictor() + for i in range(self.num_objectives) ] else: self.predictor_list = [ - gp_predictor(self.config) for i in range(self.num_objectives) + self._make_gp_predictor() for i in range(self.num_objectives) ] else: self.predictor_list = predictor_list @@ -436,7 +441,7 @@ def get_post_fmean(self, xs): self._warn_no_predictor("get_post_fmean()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) @@ -450,6 +455,35 @@ def get_post_fmean(self, xs): ] return np.array(fmean).T + def get_kernel_length_scale(self): + """ + Return the Gaussian kernel length scale(s) for each objective. + + Returns + ------- + list of numpy.ndarray or None + List of length num_objectives; each element is the length scale(s) + for that objective's GP, or None if not available. + """ + if self.predictor_list is None or all(p is None for p in self.predictor_list): + return None + self._update_predictor() + result = [] + for predictor in self.predictor_list: + if predictor is None: + result.append(None) + continue + try: + cov = predictor.model.prior.cov + except AttributeError: + result.append(None) + continue + if not hasattr(cov, "width"): + result.append(None) + continue + result.append(np.atleast_1d(np.asarray(cov.width).flatten())) + return result + def get_post_fcov(self, xs, diag=True): """ Calculate covariance of predictors (post distribution) @@ -472,7 +506,7 @@ def get_post_fcov(self, xs, diag=True): self._warn_no_predictor("get_post_fcov()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) @@ -518,7 +552,7 @@ def get_score( self._warn_no_predictor("get_score()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm, objective_index=i) predictor.prepare(training, objective_index=i) predictor_list.append(predictor) @@ -574,7 +608,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False) self._warn_no_predictor("get_post_fmean()") predictor_list = [] for i in range(self.num_objectives): - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training, 0, objective_index=i) predictor.prepare(self.training, objective_index=i) predictor_list.append(predictor) diff --git a/src/physbo/search/range_unified/_policy.py b/src/physbo/search/range_unified/_policy.py index 57fd571c..df3bbe17 100644 --- a/src/physbo/search/range_unified/_policy.py +++ b/src/physbo/search/range_unified/_policy.py @@ -15,8 +15,6 @@ from .. import utility from .. import score as search_score from ..optimize.random import Optimizer as RandomOptimizer -from ...gp import Predictor as gp_predictor -from ...blm import Predictor as blm_predictor from ...misc import SetConfig from ..._variable import Variable @@ -54,6 +52,7 @@ def __init__( else: self.config = config + self.ard = False if initial_data is not None: if len(initial_data) != 2: msg = "ERROR: initial_data should be 2-elements tuple or list (actions and objectives)" @@ -192,6 +191,7 @@ def bayes_search( interval=0, num_rand_basis=0, optimizer=None, + ard=False, ): assert unify_method is not None, "unify_method must be provided" self.unify_method = unify_method @@ -207,15 +207,16 @@ def bayes_search( simulator = None is_rand_expans = False if num_rand_basis == 0 else True + self.ard = ard if training_list is not None: self.training = training_list if predictor is None: if is_rand_expans: - self.predictor = blm_predictor(self.config) + self.predictor = self._make_blm_predictor() else: - self.predictor = gp_predictor(self.config) + self.predictor = self._make_gp_predictor() else: self.predictor = predictor @@ -376,7 +377,7 @@ def get_post_fmean(self, xs): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fmean()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit( self.training_unified, 0, comm=self.mpicomm, objective_index=0 ) @@ -409,7 +410,7 @@ def get_post_fcov(self, xs, diag=True): X = self._make_variable_X(xs) if self.predictor is None: self._warn_no_predictor("get_post_fcov()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit( self.training_unified, 0, comm=self.mpicomm, objective_index=0 ) @@ -450,7 +451,7 @@ def get_score( if predictor is None: if self.predictor is None: self._warn_no_predictor("get_score()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(training, 0, comm=self.mpicomm, objective_index=0) predictor.prepare(training, objective_index=0) else: @@ -501,7 +502,7 @@ def get_permutation_importance(self, n_perm: int, split_features_parallel=False) if self.predictor is None: self._warn_no_predictor("get_permutation_importance()") - predictor = gp_predictor(self.config) + predictor = self._make_gp_predictor() predictor.fit(self.training_unified, 0) predictor.prepare(self.training_unified) return predictor.get_permutation_importance( diff --git a/tests/unit/test_policy.py b/tests/unit/test_policy.py index e72880e6..38e38d4d 100644 --- a/tests/unit/test_policy.py +++ b/tests/unit/test_policy.py @@ -97,6 +97,32 @@ def test_bayes_search(policy, mocker): assert get_actions_spy.call_count == N +def test_bayes_search_ard(policy, mocker): + """Test that ard=True creates a predictor with ARD kernel.""" + simulator = mocker.MagicMock(side_effect=lambda x: x) + N = 2 + policy.random_search(N, simulator=simulator) + policy.bayes_search(max_num_probes=N, simulator=simulator, score="TS", ard=True) + assert policy.predictor is not None + assert policy.predictor.model.prior.cov.ard is True + + +def test_get_kernel_length_scale(policy, mocker): + """Test get_kernel_length_scale returns None before fit and array after.""" + assert policy.get_kernel_length_scale() is None + + simulator = mocker.MagicMock(side_effect=lambda x: x) + policy.random_search(2, simulator=simulator) + policy.bayes_search(max_num_probes=2, simulator=simulator, score="TS", ard=True) + + length_scale = policy.get_kernel_length_scale() + assert length_scale is not None + assert length_scale.ndim == 1 + # X has 3 dimensions (see fixture) + assert len(length_scale) == 3 + assert np.all(length_scale > 0) + + def test_saveload(policy, X): simulator = lambda x: x