diff --git a/docs/guide_graph_based_optimization.ipynb b/docs/guide_graph_based_optimization.ipynb new file mode 100644 index 0000000..10e420d --- /dev/null +++ b/docs/guide_graph_based_optimization.ipynb @@ -0,0 +1,284 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "intro", + "metadata": {}, + "source": [ + "# Optimization with a graph surrogate\n", + "\n", + "Below is shown how an SBArchOpt problem can be optimized using a graph-based surrogate. The variables are processor count `n = 1 or 2` and common size `s in [0, 1]`. The objective `(n * s - 1.4)**2 + 0.1 * n` penalizes capacity error and processor count; its minimum is 0.2 for two processors of size 0.7. The code requires SBArchOpt's graph-surrogate support and the `arch_sbo` and notebook dependencies." + ] + }, + { + "cell_type": "code", + "id": "imports", + "metadata": {}, + "execution_count": 1, + "outputs": [], + "source": [ + "from operator import itemgetter\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import networkx as nx\n", + "import numpy as np\n", + "from pymoo.core.variable import Integer, Real\n", + "from pymoo.optimize import minimize\n", + "\n", + "from sb_arch_opt.problem import ArchOptProblemBase\n", + "from sb_arch_opt.algo.arch_sbo import get_sbo\n", + "from sb_arch_opt.algo.arch_sbo.graph import (\n", + " CompositeGraphKernel, GraphRepresentation, LdWloa, NodeLabeling,\n", + " SizingFeature, SizingKernel,\n", + ")\n", + "from sb_arch_opt.algo.arch_sbo.infill import get_default_infill\n", + "from sb_arch_opt.algo.arch_sbo.models import ModelFactory\n", + "\n", + "SEED = 42\n", + "N_INIT = 10\n", + "BATCH_SIZE = 2\n", + "MAX_EVALS = 24" + ] + }, + { + "cell_type": "markdown", + "id": "problem-text", + "metadata": {}, + "source": [ + "## Defining the Optimization Problem\n", + "\n", + "`ArchOptProblemBase` declares the design variables, and `_arch_evaluate` calculates objective values for a batch of designs. Both variables are always active. `_gen_all_discrete_x` lists the two processor-count choices; the sampler varies the continuous size when generating designs." + ] + }, + { + "cell_type": "code", + "id": "problem", + "metadata": {}, + "execution_count": 2, + "outputs": [], + "source": [ + "class ProcessorProblem(ArchOptProblemBase):\n", + " def __init__(self):\n", + " super().__init__([\n", + " Integer(bounds=(1, 2)),\n", + " Real(bounds=(0.0, 1.0)),\n", + " ])\n", + "\n", + " def _is_conditionally_active(self):\n", + " return [False, False]\n", + "\n", + " def _get_n_valid_discrete(self):\n", + " return 2\n", + "\n", + " def _gen_all_discrete_x(self):\n", + " x = np.array([[1, 0.5], [2, 0.5]])\n", + " return x, np.ones(x.shape, dtype=bool)\n", + "\n", + " def _arch_evaluate(self, x, is_active_out, f_out, g_out, h_out, *args, **kwargs):\n", + " n_processors, size = x.T\n", + " f_out[:, 0] = (n_processors * size - 1.4)**2 + 0.1 * n_processors\n", + "\n", + "\n", + "problem = ProcessorProblem()" + ] + }, + { + "cell_type": "markdown", + "id": "decoder-text", + "metadata": {}, + "source": [ + "## Representing Designs as Graphs\n", + "\n", + "The custom decoder creates a system node and one node per processor. Each node is identified by its role and instance index; constant and role labelings exclude the index. The `GraphDecoder` interface provides `labelings`, `sizing_features`, and `decode(x)`. Each decoded representation contains a graph, its node labels, and its sizing values. The size already lies in [0, 1] and needs no further scaling." + ] + }, + { + "cell_type": "code", + "id": "decoder", + "metadata": {}, + "execution_count": 3, + "outputs": [], + "source": [ + "class ProcessorDecoder:\n", + " def __init__(self):\n", + " self.constant_labeling = NodeLabeling(\"processor.constant\", self._constant_label)\n", + " self.role_labeling = NodeLabeling(\"processor.role\", itemgetter(0))\n", + " self.sizing_features = (SizingFeature(\"processor_size\", \"numeric\"),)\n", + "\n", + " @staticmethod\n", + " def _constant_label(node):\n", + " return \"node\"\n", + "\n", + " @property\n", + " def labelings(self):\n", + " return self.constant_labeling, self.role_labeling\n", + "\n", + " def decode(self, x):\n", + " representations = []\n", + " for n_processors, size in x:\n", + " graph = nx.DiGraph()\n", + " system = (\"system\", 0)\n", + " graph.add_node(system)\n", + " for index in range(int(n_processors)):\n", + " graph.add_edge(system, (\"processor\", index), kind=\"contains\")\n", + "\n", + " representations.append(GraphRepresentation(\n", + " graph=graph,\n", + " node_labels={\n", + " labeling: {node: labeling(node) for node in graph}\n", + " for labeling in self.labelings\n", + " },\n", + " sizing_values=np.array([size]),\n", + " ))\n", + " return representations\n", + "\n", + "\n", + "decoder = ProcessorDecoder()" + ] + }, + { + "cell_type": "markdown", + "id": "kernel-text", + "metadata": {}, + "source": [ + "## Constructing the Graph Surrogate\n", + "\n", + "`LdWloa` compares labeled graphs at depths 0, 1, and 2, while `SizingKernel` compares processor sizes. Their contributions are learned during training. The composite combines structure, sizing, and their product; `ModelFactory` creates the surrogate and input normalization." + ] + }, + { + "cell_type": "code", + "id": "surrogate", + "metadata": {}, + "execution_count": 4, + "outputs": [], + "source": [ + "# Define a hybrid graph / sizing kernel\n", + "kernel = CompositeGraphKernel(\n", + " branches=[\n", + " LdWloa(cutoff=2, labelings=decoder.labelings),\n", + " SizingKernel(decoder.sizing_features),\n", + " ],\n", + " composition=\"additive_interaction\",\n", + ")\n", + "\n", + "# Factory returns the graph surrogate that uses the specified kernels\n", + "surrogate, normalization = ModelFactory(problem).get_graph_kriging_model(\n", + " decoder=decoder,\n", + " kernel=kernel,\n", + " multi=False,\n", + " seed=SEED,\n", + " n_start=1,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "run-text", + "metadata": {}, + "source": [ + "## Running the Optimization\n", + "\n", + "`get_default_infill` selects the acquisition criterion, and `get_sbo` combines it with the surrogate. The run starts with ten evaluations and adds two designs per infill batch, stopping after 24 evaluations." + ] + }, + { + "cell_type": "code", + "id": "run", + "metadata": {}, + "execution_count": 5, + "outputs": [], + "source": [ + "# Get the infill criterion\n", + "infill, _ = get_default_infill(problem, n_parallel=BATCH_SIZE)\n", + "\n", + "# Define the SBO algorithm that uses the hybrid graph / sizing surrogate\n", + "algorithm = get_sbo(\n", + " surrogate_model=surrogate,\n", + " infill=infill,\n", + " init_size=N_INIT,\n", + " infill_size=BATCH_SIZE,\n", + " infill_pop_size=30,\n", + " infill_gens=10,\n", + " normalization=normalization,\n", + ")\n", + "result = minimize(\n", + " problem, algorithm,\n", + " termination=(\"n_evals\", MAX_EVALS),\n", + " seed=SEED,\n", + " verbose=False,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "result-text", + "metadata": {}, + "source": [ + "## Inspecting the Result\n", + "\n", + "The returned design vector gives the processor count and size. Compare it with the known optimum of two processors of size 0.7. The plot shows the best objective observed as evaluations accumulate; the dashed line marks the minimum of 0.2." + ] + }, + { + "cell_type": "code", + "id": "result", + "metadata": {}, + "execution_count": 6, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Processors: 2; size: 0.700\n", + "Best objective: 0.2000; known optimum: 0.2000\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAFzCAYAAADlp6mhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAASY1JREFUeJzt3Qd8VFX2wPGTDgESaggdCUIAUZAiHQQ0uoiwICooWLGsnQUVVxbUdQFF18Ii/m24rAguFgQBUZpYQEAEkY4ICCH09J75f84NM8wkGZhMJsxL8vt+PkPevLnz8vIyzJzce+65ATabzSYAAAAWEejvEwAAAHBGcAIAACyF4AQAAFgKwQkAALAUghMAAGApBCcAAMBSCE4AAIClBPv7BKwqLy9PDh8+LNWqVZOAgAB/nw4AAGWGllBLTk6W+vXrS2Bg8ftBCE7c0MCkUaNGJf39AABQYR08eFAaNmxY7OcRnLihPSb2CxsREVGy3w4AABVIUlKS+QPf/llaXAQnbtiHcjQwITgBAKD4vE2LICEWAABYCsEJAACwFIITAABgKeScAEAZK3OQlZXl79NABRcSEiJBQUGldnyCEwAoIzQo2bdvnwlQAH+rXr26REdHl0otMIITACgjRa3i4+PNX6s6RdObwlaAr16LaWlpcvToUXO/Xr164msEJwBQBuTk5JgPBK24GR4e7u/TQQVXuXJl81UDlKioKJ8P8RB6A0AZkJuba76Ghob6+1QAwx4kZ2dni68RnABAGcJaX6gIr0WCEwAAYCnknFwge44my7z1B8VqKocGy6iuTaR21TB/nwoAVFgPPvig/PHHH/LZZ5/5+1QsgeDkAjl4Ml3eWrNPrCg7N0+euCbW36cBoBx6+umnZe7cubJnzx7Hvk2bNsk111wjffv2lf/85z+mZkZFct9998nx48dl/vz5/j4VyyI4uUAa1QyXe3s1EyvZFp8ka3YflwMn0/x9KgAqiNWrV8v1118vI0eOlNdff50cmjOmT5/u31+M1dhQpMTERJteHv1aXi3afNjW5IlFtqEzvvP3qQA4j/T0dNu2bdvM17Lkb3/7my0mJsZsL1iwwFapUiXbxIkTXdrcdttttsGDB9seeOABW8OGDW01atSwjRw50paSkuJok5eXZ3v++edtjRo1soWGhtratGlj+9///ud4fMuWLbaAgADb0aNHzf3k5GRbcHCw7brrrnO0mTFjhi02Ntbj71nQ+c7BftxBgwbZRo8ebatfv76tWrVqtrvvvtvxe9ProZ8tzrevvvrKnIc+z/k41113nW3UqFG2Bg0a2KpWrWq74447bPHx8bZbb73VFhkZaatTp47t6aefdvn+AwYMsD3yyCMu+/7617/a4uLiSnzs4rwmS/oZSkJsBVY3Ij/P5EhShr9PBUAJqsa6u2ltFE/bFpwO6q6dt2bNmiU33HCDTJs2TSZNmlTocc21aNy4sWzdulW+/fZbWbVqlbz00kuOx7WXRZ/79ttvm9oaDzzwgNx0003y448/mscvueQSqVWrlumZUWvWrDEVTL/55hvHNGw9Zp8+fTz+ngWd7xzsFixYIFWrVpVffvnFHHP58uXy1FNPmcf+8Y9/yL333itDhw41xcz01r9//yK/36JFi6Rdu3aybds2+eqrr8zwWOvWraVLly5y4MAB+d///idTpkwxxy+u0jx2uRnWOXz4sPlFN2/e3PxCz0dLN+/du1eCg4PNC8td8Rdtt2vXLjMXW9vBVd2ISubr0aRM8x+EKYpA2TN58mS3j1188cUyYsQIx339YHVXk6JJkyZy++23O+6/+uqrpuhbQRMnTiz2OeqH3Z133ikvvvii+UAvSqdOneTxxx8325GRkTJ8+HATYNi98MILMm7cOLn66qvN/fvvv1+++OILs19zN/T9q1evXrJy5UoTBGlQcPPNN8unn34qP/30kzm+BiqvvPKKx9+zoPOdg13Dhg3Nz6qfTTVr1pR//vOfcscdd8jzzz/vKF7mic6dO8tjjz1mtjVo0J8vPT3dcQ179+4tV1xxhQms+vXr5/FxS/vYvuDXnpPMzEwTdcbExMiwYcOkbt268sYbb5zzOVOnTpUGDRrIgAEDzMVs1qyZLFmypFA7jVw1IImLizM3Tb46ceJEKf40ZU/UmZ6TrNw8OZXm+yI6AKDq1Kkj3bp1Mz0Pv/32m9tAylmNGjXk1KlTZjspKUkOHTpkPiydde3a1fzlb6e9IhqUKP165ZVXmg9Z3d6+fbscOXLE3Pfkexbk6Tmo9u3bu/zRrEFQRkaGWRepOJo3b+5yX3uCitrn7pz9dewy33Py3HPPmahMs7g14Pj4449NkNKxY0fzyyxIu+b0Qm3ZssW82PWv/b///e/mOdqTosGN0mNql9mMGTPknnvuMfu0a0rXpdBuP+QLCw6SmlVC5WRqliQkZZhtAGXL+PHj3T5WcP2dsWPHum1bsOf0kUceEV/R3oKlS5eaPyo1gNDeDf2j9Fzf3xMFe3z12A8//LD5TNEZQRqI6B+l2nuivfKxsbFmobqSfM/znYOvjuvuON4cu6iFIn117HLZc/LOO+/I3XffbQITpQGFjnm9++67RbbXSFTHwDQwsV9I7YJKTU013XZ2Op6pU9TsgYnSbikdk0TRQzvknQBlk5azd3fToW9P2xaczuuunbc0ONBebu2t0KBh9+7dHj83IiLCfE4UzO1Yt26dtGrVynFf3+Nr165t8jp0v/4xqgGLDtV8/fXXLvkmxeXpOSgNjOx5Lmr9+vUSFhYmF110kbmv17q0VpauUUTvT3GutVT04ETzTLSLrWAPiXaZ6S/WU/agxP5L14QtHVccOHCgCVo2b94sx44d8/HZl7+k2IREkmIBlC7N/9NETP0jVAOUnTt3evxczQ3RnBkNMnSI5c0335TFixebHBA7e97Jf//7XzOkozQY0nwSHeovSXDi6TmogwcPmn0nT540n1GaDKu1Tez5Jprfo0NB+riv9ezZUz7//HPZsGGDJCcnm04AzYspa/w2rGP/pWiykDONdD3NDdF2WlVPk5+0u05pIKIJX/qL1xelRtE6xml/wRb8fs75L3qz0xdeRRB9puckIenszw4ApUU/oPXDc8iQISZYWLFihUfPe+ihhyQlJcUk1uofti1atJB58+YVygHRoOSTTz5xCUR0+4MPPnDJN/GGp+cwaNAgExhoT45+vfHGG10SlzU5VmfING3a1Dyu275y5513mtSHa6+91gzrac6lJjprh0BZEqDzif3xjTVi1oBCX5j2CNf+y9f8kIIJRgVp8HDVVVeZ8T6NYrXLTemsH8090ch07dq1ZnxRAxZNxtIARaPIouhQ0DPPPFNof2JiouPY5dG/vtolry7fLcM7N5bJQ9r6+3QAuGFPqNRe4kqV8v+ogPVoIKBTuPWP4Yr8mkxKSjI9Vt5+hvptWKdRo0YmqisYzel9DSzORSNNnX2jY3rLli1z+cG1p0S7DjV/xZ74pDkqGrna57+7SyrTi2i/abdcxZpOzLAOAMAa/BacaAChc6sXLlzo2Kdz6rUXxHlOtQ7JaBdVwcBEc0u0K0ynOjnTgEd7VAoGPToFTAMXdzRZSYMc51tFEB1JITYAgLX4dSqxZlRrMZsJEyaYueKvvfaayQlxnmWjxWt0eEYr+GkuiU5F02li77//vksylU5Ls8/i0eGZHj16mKnK3bt3N9nUOt744Ycf+uXntLKoavacE3pOAKCktBIuynhworkmml+iCx7pDJu2bduaacTOvRYadGjVOnvPivaY6PhWwfLHmg2ti0mpyy67zNQ6efnll83xtRibftWcE7iKjswPTo6nZJnViUOCWNEAAFBBE2KtrqTJPGVFXp5NWk5YItm5Nvnuyb7SoLrnpZUBXDgkxMJqymVCLKwhMDCAoR2gDOHvSVSE1yLBCSjEBpQB9rVaSrIyMOBL9oUhC1YXLjerEsO/KGEPWJ+WotdZjlq3ST8MCq6bA1zIHhMNTLSumM6YdV7k0FcITuAITqgSC1iXlmavV6+eGePfv3+/v08HEA1MnBdS9CWCEzhm7DCdGLA2XXhPl+VgaAf+pr13pdFjYkdwgrM5J9Q6ASxPh3MoX4/yjkFLkHMCALAUghOczTlJpEosAMD/CE4g0WeCk9SsXEnJzOGKAAD8iuAEUiUsWKqF5acfHaH3BADgZwQnMKLOJMUeJSkWAOBnBCdwmU58hOAEAOBnBCcwqBILALAKghO4BCdHkzK5IgAAvyI4gcuMHRJiAQD+RnAC1yqxydQ6AQD4F8EJDAqxAQCsguAELrN1jiZnSl6ejasCAPAbghMYtauGSUCASE6eTU6kZnFVAAB+Q3ACIyQo0AQoitWJAQD+RHCCQkmxzNgBAPgTwQkKTSdmxg4AwJ8ITuDAjB0AgBUQnKBwcEKVWACAHxGcoHCVWBb/AwD4EcEJHKLsVWIJTgAAfkRwgkKF2AhOAAD+RHCCQsM6p9KyJSM7lysDAPALghM4RFYOkdDg/JfEseRMrgwAwC8ITuAQEBBAUiwAwO8ITuCCKrEAAH8jOIGbWicZXBkAgF8QnKDoEvYEJwAAPwkWi8jOzpaQkBCP2+fl5UlgYNGxVWZmpjmes6CgIKlcuXKJz7O8o0osAEAqes/Jc889J7Vq1ZJKlSpJq1at5Ouvvz5n+y+++EKuvPJKiYiIkKpVq0pcXJxs3brVpc0jjzwi1atXl+joaMetd+/epfyTlA91z9Q6oUosAKBCBif//ve/5cUXX5SPP/5YUlJSZMSIETJw4EDZu3dvke1zc3PljTfekEmTJsmxY8fk4MGDUrNmTbn66qslMTHRpe3gwYPNMe23H3/88QL9VGVb3WpUiQUAVODg5JVXXpG77rpL+vTpY4ZcJkyYIFFRUTJz5swi2+vQzKJFi0wviLavUaOGCW7i4+Nl3bp1hdrn5ORcgJ+i/FaJtdls/j4dAEAF5Lfg5MSJE7Jnzx7p1auXy34NPNauXevxcQ4dOmS+6tCQMw1iNIDR4Z8//elPsmPHDh+decXIOcnIzpOkdII7AEAFCk4SEhLM1zp16rjs1/tHjx716BhZWVkmv6RLly5y+eWXO/a3adNGlixZImlpafLrr79KaGioyVM5efKk22NpEm1SUpLLrSKqFBJkKsWqhGSmEwMAKmBCrM66KXhfK5Wej+af3HrrrXL48GGZN2+ey3MeeughE4zo7J9GjRrJ7NmzTU6KtnNn8uTJEhkZ6bjp8yr6dOIjiQQnAAALTyX+448/PD5ow4YNz9umfv365mvBXhK9X69evfMGJqNGjZLvv/9eVq9eLY0bNz5n+2rVqkmDBg3kt99+c9tm/PjxMmbMGMd97TmpqAFKVESY7ExIZsYOAMDawUlxPqg9SaTUqb46/LJ8+XK54YYbHL0mK1askLvvvttluEX322uU6PZtt91mgpJVq1ZJTEyMR/ktBw4cOOfPEBYWZm4423NylEJsAAArD+vs3r3bcZs+fbrp+dBpvRs2bDA33dYeD50e7Kknn3xS3nvvPTPcosHDww8/LOnp6XL//fc72jzwwAPSqVMnR9Bz5513yldffWUSXrV+iX2qsH1mjgYzOrVYgxzNa1m/fr0MGTLEJMzqMBA8n7FDrRMAgKV7Tpo3b+7YvvHGG2X+/PnStWtXx74OHTrIpZdeavI9/vKXv3h0TA0WNGn1mWeeMYFE27ZtTRE2+5CP0uJs4eHhZlsTWvX7qh49ergcS4Mi7VHR3o+//e1vMnXqVNm0aZOZbtyzZ0/58MMPTU0UnF+Uo4R9JpcLAHDBBdi8KGahwYImourQjLNTp06ZfJPU1FQp6zTnRBNjNZFWpyNXJF9tS5DR/9kglzaMlM8fdA0CAQAo7c9Qr2brNG3a1BRQK+jVV181j6FsqxuRn3vDbB0AQJlZ+E+DkEGDBsn//vc/6dixo8kF0byT/fv3y+eff+77s4RfEmKPp2RKTm6eBAf5fcY5AKAC8epT56qrrjLVXXWWjeaMaBLrsGHDzL5+/fr5/ixxQdWqGiZBgQGSZ9MAJYurDwCwfs+J0qRVTWRF+aOBSZ2qYWa2jq6xY5+9AwDAheB1f312drZ89913pvqqcz0RlA91mU4MAChLwYnWJNG1bHQIRyu12o0ePVoWLlzoy/ODn0SfSYrVnhMAACwfnDz22GOmxknBxfHGjh0rU6ZM8dW5wQKrExOcAADKRM6Jlo7fsWOHWe3XmRZR27hxo6/ODRYITo4kUogNAFAGek4yMjIkMDD/qc6rAcfHxzuquaJso+cEAFCmgpM+ffrIzJkzXYITrQo7btw46d+/v2/PEH6tdcKwDgCgTAzrTJs2TXr16iWLFy82Bdi0xsmaNWvMYzqDB2VfdOSZKrEkxAIAykLPSWxsrGzdulXi4uJk4MCBZpjn3nvvlc2bN0tMTIzvzxJ+W/wvOSNH0rLyV3wGAMCyPSe62m/79u1lwoQJvj8jWEK1sGAJDw2StKxcszrxRbW9rtcHAEDp95xojZM2bdrI5MmTTc0TlD+aS3R2xg61TgAAFg9Odu7cadbVeffdd80qxL1795a33npLTp8+7fszhN9XJz6aTHACALB4cNKiRQuzrs7u3bvlhx9+kMsuu0yefvppiY6OlqFDh/r+LOHXGTv0nAAAysTaOnZXXHGFvPbaa7JgwQKTKPvJJ5/45sxgoVonFGIDAJSR4GTPnj3y7LPPSsuWLaVbt24SEREhb775pu/ODn5FITYAgD94NQXj9ddflw8++EDWrVtnEmNvv/12ueWWW6Rx48a+P0P4jSMhllonAACrBydTp06V4cOHmyqx7dq18/1ZwVKF2KgSCwCwfHDyn//8R/r27ev7s4Ele06OJmWaSsDO6ygBAGCpnJOrr77afFihfIuqlh+cZOXmyam0bH+fDgCggvAqONHaJrt27fL92cBSQoMDpVaVULPNdGIAgKWDE61pokmwa9eulaSkJLO2jvMN5W+NHfJOAACWzjm54447zNeuXbsW+ThDPuVHdESYbI8nOAEAWDw4WbNmje/PBJYUHcl0YgBAGQhOevTo4fszgaWTYqkSCwCwfIXY7Oxs+e6772T27NmOfSdOnPDVecFiPSfknAAALB2cHDhwQC6//HLp16+fjBo1yrF/9OjRsnDhQl+eHyyyMjGzdQAAlg5OHnvsMZMMqzN1nI0dO1amTJniq3ODlQqxJTMLCwBg4ZyT1atXy44dOyQ0NL8Ghl3btm1l48aNvjo3WED0meDkeEqWZOXkmdonAACUJq8+abSWSWBg/lOdS5rHx8dLeHi4784OflcjPFRCgvJ/x8dSMv19OgCACsCr4KRPnz5m0T/n4CQ1NVXGjRsn/fv39+0Zwq8CAwMcM3bIOwEAWHZYZ9q0adKrVy9ZvHixKbg2bNgwR+0TncGD8pcUe+h0OjN2AADW7TmJjY2VrVu3SlxcnAwcONAM89x7772yefNmiYmJKdaxPvroI+nevbs0b95c/vznP8v27dvP2V7X9Ln//vvlsssukw4dOsgjjzwiR48eddv+5ZdfloYNG8qkSZOKdV44i+nEAADL95yoqKgomTBhguN+cnKyVK1atVjH+PTTT+WWW26R119/3cz+eemll6R3797y66+/Sp06dQq1z83NNQHMww8/bAKUtLQ0M0NIpzSvX79eKlXKH36w+/HHH82xQ0JC5PTp097+qBWefcbOkSRm7AAALNpzor0bTzzxhOO+9l5ERERI/fr15eeff/b4OM8995zcdtttct9995mekHfffdcME73xxhtFtg8KCpJffvnF9NJceuml0qVLF3nvvfdML44uQuhMpzmPGDFC3n77bYmMjPTmx0TB6cRJJMQCACwanIwZM0auuuoqxzDLO++8I4sWLZLhw4fL448/7tExNHjYtGmT4zgqODjY9IJ888037k/4zCwh50q1SntHnN1zzz0yePBgczz4ZjoxCbEAAMsO63z//fcyf/58s/3ll1/KoEGDZMCAASZ35KKLLvLoGIcOHTJfo6OjXfbXrVvX5K54QntZnnrqKZOv0qlTJ8d+7S3ROizvv/++xz9TZmamudkVLDBXkdl7TihhDwCwbM+J5nYcOXLEbC9ZskT69u1rtrOysgoVZnMnLy/P0VviTHtANLfEExqYrFixQubNm+f4vjrk9OSTT8oHH3wgYWH5pdc9MXnyZDP8Y781atTI4+dWlBL2BCcAAMsGJ9pLMmTIELnzzjvNFGKdsaOWL1/uMkxzLvaE1+PHj7vs1/tFJcMWpLNvpk+fbqYz6zo/dnoOKSkpZiaRztLR27Zt28zQk267C3zGjx8viYmJjtvBgwc9+jkqUs9JalauJGfkD6MBAGCp4ESDguuvv15ycnLMQn86c0dpoDJx4kSPjqHPadKkSaG6KN9++6107tz5nM995plnzMweDUx69Ojh8tgdd9whe/bsMQmy9tvFF18sN910k9nWpNqiaC+LJvU635CvSliwVAvL7+FKICkWAGDFnBMtUa8zbQqaMWNGsY7zwAMPmIUCdTqxrsujQY+ueKzJrHaaYPvTTz/J119/be4///zzpgicBiY9e/YsdMwqVaqYW8GhIp3mrD0n8E7dyEqSfDTFDO00jyrelHEAAC5InRMtVz9nzhxH0bTWrVubqbvFWVvnr3/9q1mPR6cEa4+G9lbMnTtXWrVq5Whz8uRJR36Lbj/99NPme+jMIGcvvvhioX3w7YydPUdTmLEDACh1ATad8lJMGzZsMHknmr+hPR5K649oD8UXX3zhkgPiCU2k1TyP2rVruywkqE6dOmUe11k8mkR7+PDhIo9Rs2ZNt4GRVpDVYZvi1DvR2TraXs+LIR6RMR/9LJ/8dEgev6al/KVP82L9fgEAFUtSCT9Dveo50aJpOn34tddec1Rl1RL2WrlVC6Rptdbi0Jk27pJga9So4VLjxJuhGXtODEpe6yQhkSqxAIDS5VVwouXlly5d6lIuXrf/+c9/SuPGjX15frAIStgDACw9W0cX97PngTjTfc2aNfPFecGyhdgoYQ8AsEhwosM29puWr9epuTpj5tixYyanQ7dvvPFGk+SK8oeViQEAlkuILZioei5e5NhaDgmxruIT06Xr5BUSFBggu/9xrQQGev56AABULEkXKiFWC6yh4qpTNUw0Ps3Ns8nx1EyJqnY23wgAAF/yODgpWIkVFUtwUKDUrhomx5IzJSGR4AQAYMEibOrEiROya9cuM4zTsmVLqVWrlu/ODJacTmyCk6QMaSue14wBAKDUZ+ukp6ebeiZaGK1bt27SvXt3s6379DGUT0wnBgBYNjgZN26crFixQubPn2/Kz+sUYt3WFYH1MZRPdSPCzNejSRRiAwBYbFhn3rx5smzZMmnfvr1j3+DBg80qw3FxcWYBP5TfKrFHCE4AAFbrOUlOTjaBSEG6T6cPobwP61CIDQBgseBEF/Z74YUXXOqZ6PaUKVOkQ4cOvjw/WEjdyPzghGEdAIDlhnWmTZsm11xzjXz88cfSqVMns08X+0tISDBr7qB8YlgHAGDZnhOdobN792655ZZbJDs7W3Jycsy27tPHUL4TYk+nZUtGdq6/TwcAUE55XedEpw5PmjTJt2cDS4usHCJhwYGSmZMnR5MypXGtcH+fEgCgHPKq5wQVk66vRK0TAEBpIziBV3knWiUWAIDSQHACr2bsEJwAAEoLwQmKpW61/KRYghMAgOWCE52l891338ns2bNdFgJE+RZ9pueEQmwAAEsFJwcOHDCF2Pr16yejRo1y7B89erQsXLjQl+cHi4my55wkknMCALBQcPLYY49J165dC5WqHzt2rKkSiwqQEJtMcAIAsFCdk9WrV8uOHTskNDTUZX/btm1l48aNvjo3WLlKbGKGWbJApxcDAOD3npOMjAwJDMx/qvOHU3x8vISHU5irPIs6UyVWC7Elpef4+3QAAOWQV8FJnz59ZObMmS7BSWpqqowbN0769+/v2zOEpVQKCZLq4SFm+wi1TgAAVlr4r1evXrJ48WLTtT9s2DBZs2aNeUxn8KD8D+3o+joanLSMrubv0wEAlDNe9ZzExsbK1q1bJS4uTgYOHGiGee69917ZvHmzxMTE+P4sYc0ZO/ScAACs0nOyadMmad++vUyYMMH3ZwTLiz6Td8J0YgCAZXpOtMZJmzZtZPLkyabmCSoW++J/TCcGAFgmONm5c6fccMMN8u6770rTpk2ld+/e8tZbb8np06d9f4awHMfKxImZ/j4VAEA55FVw0qJFC3nmmWdk9+7d8sMPP8hll10mTz/9tERHR8vQoUN9f5awFFYmBgBYeuG/K664Ql577TVZsGCBSZT95JNPfHNmsP6wDgmxAACrBSd79uyRZ599Vlq2bCndunWTiIgIefPNN313drCkupH5CbHHUzIlJzfP36cDAChnvJqt8/rrr8sHH3wg69atM4mxt99+u9xyyy3SuHFjr2b+aEG3hIQEU/5+zJgxUqNGDbfttdjbrFmz5Pvvv5fg4GDp0aOH+f4hIfmFwVReXp7Mnz9fvvzyS0lJSTHneM8995hhJ5RcrSphEhQYILl5NjmWkin1IitzWQEA/u05mTp1qvTs2dMEFlrvZPz48V4FJmvXrjULCGpgoYXcVq1aJd27d5e0tLQi22vQcckll5h1fQYMGGDOQWcMXXfddeYxu9tuu80UiNNKtoMHD5ZvvvnGzDA6cuSINz8uCtDAJKramenESSTFAgB8K8CmJV6LSQMB+9o6JdG3b1+JjIyUTz/91NxPTEyUBg0amIDjoYceKtReT/XkyZNSq1Ytx77169dL586dTaCj+S9KZw1Vr17d0SY9Pd0MOb399tsmcPGErris56bnpM+Fq8H//k5+PnhaZt7aQa65hB4pAIDvPkM9jjCOHz9ubkoDBPv9om6e0IBBezS0Z8NOfxBdm2fp0qVFPkfX8XEOTFRUVJT5qsM3ds6Bidq2bZvk5OSYWUbwjbpnCrEdTc7gkgIA/JNzUqdOHUfvhX3bHU86Yw4ePCi5ubnSqFEjl/0NGzaUlStXenpa8tJLL5mAxd5rYrdx40b529/+ZqI3TdydN2+eGUJyJzMz09zs9Hk4/3TiI4kEJwAAPwUnml9S1La3srKyzNfKlV2TKcPDwx2PnY8mxs6YMcMMC1WtWtXlsSZNmsijjz5qenLee+89+fvf/24WK3SXFKtDSVq7BcVbX4eViQEAfgtO2rVr59ieMmWKzJ07t8h2N998s9vHnNmHXnSIyNmJEyfOOVvHTr+HzsDRwEMXHyyodu3acs0115jtm266yQzp/Otf/zLJvEXRpF6dKeTcc1KwVweFe06OkhALAPAxr7JadYjE3XDORx995NExdPhGA4iCvTB6XyvOnot+D01s1ZL5I0eOPO/30tlAmmgbHx/vtk1YWJhJ2nG+wb3oSHpOAAAWCE50Fox9/Rz7tv2mPSCLFi2SevXqeXy8UaNGyTvvvONIotW6JBqcOM+oeeWVV+S+++5z3Nf6JRqQ/N///V+RM28yMjLMMZ3zXjTxVmf16Owg+DYhliqxAAC/FmFzHm4pauhFpxe7GzYpilaX/eWXX+Tiiy82N93+5z//aQqr2WkdFZ0mrHRK0ogRI8ysng8//NDc7B577DGJi4szvSQa4GjPTExMjAmcfv/9d3niiSc8nkYMz0vYJ2fkSFpWjoSHelXPDwCAktU52bBhg/naqVMn0xPhTIMCzdGoWbOmFJdO9dUKsa1bt5a6deu6PPbrr7+aAEOLs2VnZ8vy5cuLPIYWZ9OAxE4DGQ1sqlSpYgIf/Voc1Dk5N33ZtJn4paRl5cqKv/aWZnVcE5IBABVXUgnrnHhVhE17Ipo2bSrlGcHJ+fWdtkp+O54qH47uIl1jXOvPAAAqrqQLVYTNmfZuaHn4gnSfFldDxcDqxACA0uBVcKLTbnft2lVov+57+umnfXFeKANIigUAWCY4mTNnjtx6662F9uvKxM5Jqijf6jKdGABgleBEy7wXtXJwamqqyxo3KN/qVsufscN0YgCA34OT3r17m+EbnT1jp9u6lo2WiEfFKsSWQJVYAIAPeVWc4oUXXjC1SHSKbpcuXcy0Uq1For0pWvAMFSshlsX/AAB+7zmJjY01BdO0qJkO5egMndtvv122bNkirVq18ukJwvoJsUeTMzxaiRoAAE94XdZTy9Szim/FFnUm5yQ71yYnU7OkVtX8YAUAgAvec2LPMfnuu+9k9uzZLisKo+IIDQ6UWlVCzfaRpAx/nw4AoCIHJwcOHJDLL79c+vXrZxbvsxs9erQsXLjQl+eHMpJ3cpSkWACAP4MTXWSva9eupjyts7Fjx8qUKVN8dW4oQzN26DkBAPg152T16tWyY8cOCQ3N79K3a9u2rWzcuNFX54YygCqxAABL9JxkZGRIYGD+UwMCAhz74+PjJTw83HdnB8tjfR0AgCWCkz59+sjMmTNdghOdUjxu3Djp37+/b88QlkatEwCAJYZ1pk2bZirB6irEWt9i2LBhsmbNGvOYzuBBxRF9JiGWKrEAAL8XYdu6davExcXJwIEDzTDPvffeK5s3b5aYmBifnRysj2EdAIBlirBFRUXJhAkTfHs2KLMJsSdSsyQrJ8/UPgEA4IIEJ8ePHzdfa9eu7dh2p3LlylKlSpUSnRjKhppVQiU0KFCycvNMGfuGNUiIBgBcoOCkTp065qvmmNi3z6Vx48Ymafbaa68t2RnC0jQhOioiTP44lS4JSQQnAIALGJxs2rSpyG13pe2//vprefDBB2Xv3r0lO0OUibyT/OAk09+nAgCoSMFJu3btitx2R8vbf/XVV96fGcrcjJ0jiayvAwDwY0Ks1jWZM2eObN++3dxv3bq1jBgxwlGELSgoSFasWOGDU4TV6bCOSkgmOAEAlJxXUys2bNggzZo1k/Hjx5shHr09+eSTZhrxTz/95IPTQpmsdULPCQDAX8HJfffdJ4MGDZI//vhDVq5caW66rTVPtN4JKhYW/wMA+FKATaffFJNOFT548KCZVuxMpxjrLJ20tDQp63TF5cjISElMTJSIiAh/n46l/bD3hAx/a61UDQuWTk1r+Pt0UExNalWRCde1lqDAs+tkAYA/P0O9yjnR4ZsjR44UCk50nw73oGKJqVPFfLClZObIyp3H/H06KLZjEtcmWrrG1OLaAbAEj4MTLVFvN2bMGLnpppvkxRdflE6dOpnaJ5qHMnbsWLP4HyqWqIhKMv++rrLnaIq/TwXF9N91B2TzwdOy40gSwQmAsjesY1992BNejBRZDsM6qAheXrZTXluxR27q2Eim3nCpv08HQDmRdKGGdeyrDgMoP1pG579p7EhI9vepAEDxg5MePXp42hRAGRFbr5r5uutIsuTm2UiKBVC2i7CpEydOyK5du8wwTsuWLaVWLRLqgLKkaa0qEhYcKOnZuXLgZJpcVJsFOwGU0Ton6enppp5J3bp1pVu3btK9e3ezrfv0MQBlg86yalE3v/dk55Ekf58OAHgfnOiMHC1NP3/+fImPjzdTiHV7+fLlzNYBypjY6PzgZHs8eScAyvCwzrx582TZsmXSvn17x77BgwdLkyZNJC4uTqZPn+7xsRISEswaPfq1bdu2ZopycPC5T0uDoO+//96001yYnj17etUGgCbF2ntOCE4AlOGek+TkZBOIFKT7dPqQp3bv3m0CkqVLl0pISIhMnDhRrrnmGsnNzS2yfV5enlnteMqUKZKdnW1yXrRk/l/+8pditQFwVqt6Z2bsMKwDoCyXr9c8k169esnkyZMd9U/0ME888YR899135uaJIUOGmJL3q1atksDAQDlw4IA0b95c3nnnHRk5cmSh9vo9Nm/eLO3atXPs0x4c7a3ZsmWLCXQ8aeMJ6pygojiekikd//G16H/lX5+Jk/DQEuXJA4CU9DPUq56TadOmyYwZM6RFixYyYsQIc9PtmTNnmqqxntBejcWLF5vnamCidF2ePn36yGeffVbkczQQcg46lD3YOHTokMdtAJxVu2qYuemfKbsSqPILwP+8Ck6050SHZG655RYTZOTk5Jht3aePeUJ7STIzM806Pc70vh7HU7NmzTILEWoZ/ZK00XPRSM/5BlQUrc7UO9kRz+segP953X+rU4cnTZrk9Te2r1xcrVr+m6Kddv94uqrxypUrTZ7Kv/71L7c1Vjxpo3SI6plnninWzwCUFy3rVpM1u4/LDpJiAZTVnhNfsAclp0+fdtl/6tSpQgFLUXQmzvXXXy9PPvmkPPDAA163sRs/frwZG7PfDh48WKyfByjLYkmKBWAhfgtONL+katWqsmPHDpf9er9169bnfO4PP/xgZvU8/PDD8uyzz3rdxllYWJjptXG+ARWt1on2nJSHhTsBlG1+C040CXbo0KEmHyQjI8Ps09k0OtPnxhtvdKmpMnXqVMf9devWOYKO559/vshje9IGwFnNo6qaarGn07LlaHImlwZA2ZtKfPPNN8vcuXOL/VhBWlm2d+/eEhoaagq66eydAQMGyPvvv+9oc/fdd8vatWtl69atkpKSIg0bNjS9HMOHDy/0fbt06eJRG08wlRgVTf+XV8ueoyky645O0qdllL9PB0AZVtLPUK8rxBYVgGic89FHH3kcnERHR8vPP/8sS5YsMRVi77nnnkKrH2tAodOLzckGB7tNwrXnqXjSBkDRlWI1ONGhHYITAP5UrODEOXm1YCKrVmbVIZl69eoV6wR0iq8WY3Onf//+ju1KlSrJo48+es7jedIGQGGtoqvJF1viKWMPoGwFJzVq1Chy2zmPxDk/BEDZERud3/W6nVonAMpScLJ+/XrzVYuZ2bftdG2cRo0aSc2aNX17hgAu6AKAe4+lSHZunoQE+S1fHkAFV6zgpGPHjubrvn37pGnTpoUWA9SpwQDKpoY1KkvVsGBJycyR346lOoIVALjQvPrTKD093SzyZ/fII4+YbNz69eubBFcAZY+uS2UPSFihGECZC07GjBkjV111ldnetWuXWUV40aJFZuru448/7utzBOCHYmwA4C9eTSXWsvDz5883219++aUMGjTI1Cfp3r27XHTRRb4+RwAXuow9SbEAylrPiU7X1QJqSmuU9O3b12xnZWWZgmoAynbPyU56TgCUtZ4T7SXR2iQdOnSQNWvWmBL0avny5Y7hHgBljz3n5HBihiSmZUtkeIi/TwlABeRVz8n06dPNar85OTmycOFCiYrKL3WtgcrEiRN9fY4ALpCISiHSoHpls01SLIAy1XMSHh4uzz33XKH9M2bM8MU5AfDz0M6h0+myMyFZrmhWi98FgAvO6ypL2dnZplz97NmzHftOnDjhq/MC4Cex9fKHdrbHM2MHQBkKTg4cOCCXX3659OvXT0aNGuXYP3r0aDPMA6DsanmmjP3OI0n+PhUAFZRXwcljjz0mXbt2NUsiOxs7dqxMmTLFV+cGwE8LANpn7OTl2fgdACgbOSerV6+WHTt2FJo23LZtW9m4caOvzg2AHzStXUVCgwIlNStX/jiVLo1rhfN7AGD9npOMjAyzArG95LVdfHy8SZYFUHbpgn/No/LXyWLGDoAyE5z06dNHZs6c6RKcpKamyrhx46R///6+PUMAfkuKpYw9gDIzrDNt2jTp1auXLF68WGw2mwwbNszUOFE6gwdAeVljh6RYAGWk5yQ2Nla2bt0qcXFxMnDgQDPMc++998rmzZslJibG92cJ4IKKPTNjh54TAGWm50RpVdgJEyb49mwAWGpY5/fjqZKRnSuVQoL8fUoAKhCvi7A5y8zMlLS0NF8cCoAF1KkaJjWrhIrOJN6dkOLv0wFQwRQrOMnNzZUXX3xRbrnlFrPYn+abPPjgg1KlShWpVq2aXHvttVSJBcoBTXS3551sJ+8EgJWDEx3GefbZZ+XQoUPy0EMPmYqwy5Ytk1deeUVeffVV2bNnD0M9QDlboXgHZewBWDnnZM6cObJgwQLp27evrFixwpSv/+mnn6R9+/bmca0aO2TIkNI6VwAXUCt7GfsEZuwAsHDPifaY9OjRw2x3797dURXW7rLLLjNtAJSjWif0nACwcnCSk5PjKFkfFhZmvgYHn+180W3NSwFQ9l0cVU20xuKJ1Cw5lpzp79MBUIEUeyrxqlWrznkfQPlQOTRILqpVRX47nmqKsdWpVsffpwSggih2cHLllVee8z6A8pUUa4KT+GTpeTHBCQALBif79u0rvTMBYMlKsUu2HqFSLADrBidNmzYtvTMBYOEFAJmxA6CMVYgFUD7ZC7HtPpoiObl5/j4dABUEwQkAtxrVCJfw0CDJysmT30+kcqUAXBAEJwDcv0EEBkiLumfK2FPvBMAFQnAC4Jxanck72XkkmSsFwJpTiX0tOzvb1EpJSEgw1Wa1yuz5aBXaH3/80RR969ixo9SrV6/Idt9++63s2LFDBgwY4LYNgPPP2FEkxQKoEMHJiRMnzPo8KSkpJjB54IEHZOTIkTJ9+vQi2+sqyKNGjZI1a9bI5ZdfLmlpaWZ7ypQpZiFCuy+++EKeeOIJE7xs3rxZVq5cSXAClHABQIZ1AFSI4GT8+PGm50QDiCpVqsj69evliiuuMD0d1157bZHBSVxcnMyaNUuCgoLMPt2+6667zHOaNWtm9gUGBsrcuXOlevXq0qhRowv+cwHlccbOodPpkpSRLRGVQvx9SgDKOb/lnOTl5cm8efPkzjvvNIGJ6tSpk3Tp0kU+/PDDIp+jQcett97qCEzUn/70J3MsHb6x08DmkksuuQA/BVD+VQ8PleiISmZ7F3knAMpzcHLw4EFJSkqS1q1bu+xv06aNbN261ePjLF682AQtzqsjeyMzM9Ocj/MNQMFibCTFAijHwYn9w79GjRou+2vWrOlxYLBr1y4ZM2aMPProoyUevpk8ebJERkY6bgwHAWeRFAugQgQnlStXNl+Tk13/EtP79sfO5ffff5f+/fub2wsvvOCT/JfExETHTXt2ALjmnegCgABQbhNiGzduLCEhISbIKLi4YPPmzc/53P3790ufPn2kc+fOMmfOHJccFG+FhYWZGwD3wzpa60QT0wMCArhMAMpfz0loaKiZeaOzavTNTsXHx5tpvwMHDnSpVfLJJ5847h84cMAEJlrfRJ+r04UBlK5mtatKSFCAJGfmmFk7AFCa/PrJPnXqVOnWrZsMHTrUzNJ5//33pUOHDqbWiZ1OFV67dq0MGTJEMjIypG/fvqa+iQ7n6GN2PXv2lJYtW5rt3bt3y+rVq+XUqVOOuid79uwxtVH0BqB4QoMDJaZOVZMQq70nDWuEcwkBlM/gRGfqbNmyxQQZOlSjhdTuuOMOM9zjHHQ0aNDAbOfk5JheE7VhwwaXY8XGxjqCk2PHjpmARmkNFA1S9H5UVBTBCVCCvBMNTvTWr1VdriOAUhNgs4+pwIXOGNJZO5ocGxGRX74bqMjeWLVXpi7dIdddWk+mj6AHEkDpfYay8B+AYifFAkBpIjgB4JFWZxYA/O14qmTm5HLVAJQaghMAHqkbESaRlUMkN88me46mcNUAlBqCEwAe0domFGMDcCEQnADwmCM4OcLaUwBKD8EJAI/F1svPO2EBQAClieAEgBc9J8zYAVB6CE4AeKxF3fzg5FhyppxIyeTKASgVBCcAPFYlLFia1MovXU+9EwClheAEQLG0PNN7sp2hHQClhOAEgFdJsTuZsQOglBCcACiWViTFAihlBCcAiqXlmeBEc060WiwA+BrBCYBiaVKrilQKCZTMnDzZfyKVqwfA5whOABRLUGCAY0ox9U4AlAaCEwDFRjE2AKWJ4ARAscVGnyljH88aOwB8j+AEQLHRcwKgNBGcAPB6xs6Bk2mSmpnDFQTgUwQnAIqtVtUwiaoWZrZ3JrAIIADfIjgBUOJ6JwDgSwQnALzS6kwZe5JiAfgawQkAr7AAIIDSQnACwCux9c4O69hslLEH4DsEJwC80jyqqqkWm5ieLUeSMriKAHyG4ASAV8KCg6RZ7Spme0c8SbEAfIfgBIDXYu1JsczYAeBDBCcAfFApljL2AHyH4ARAiYMTap0A8CWCEwAlHtbZczRFsnLyuJIAfILgBIDX6kdWkmqVgiUnzyZ7j6VwJQH4BMEJAK8FBAQwtAPA5whOAPhkjZ3tJMUCKE/BycGDB2XDhg2SlORZxn9ubq5s375ddu/eLTk5OT47LoDii43OzzshKRZAuQhOMjIyZOjQodKyZUsZOXKkREdHy+uvv37O5zz//PPSoEED87yrr75amjZtKosWLSrxcQF4p9WZMvYUYgNQLoKTZ555Rn788UfZu3ev6QmZM2eOPPzww7Ju3Tq3PSbp6emybds2c9u3b5/cfffdctNNN8mRI0e8Pi4A77Womx+caAn702lZXEoAJRZg8+OKXdqjcf/998vEiRMd+9q2bSvdu3eXmTNnenSMhIQEc5wvvvhC/vSnP/nsuDoUFBkZKceOHZOIiPxua2eBgYESHBzsuJ+VlXXOpMGQkBCv2mZnZ7tdVK202qrQ0FCv2uowW15enk/a6vnqeZdmWw149eaLtvp60NeFVdrqNTjXsGdQUJC5+aJt/5dXy+HT6TLrjs7StXkdr46rrzF9rfmirfP/z9Jqe77/y7xHnMV7RMV7jzh58qTUqVNHEhMTi/wMPZ+z/8susMOHD5vAokOHDi77O3fuLJs2bfL4OOvXrzdfY2JiSnTczMxMc7Oz56m89NJLUqlSpULtL774YhkxYoTj/rRp09y+qTVp0kRuv/12x/1XX31V0tLSimxbv359GT16tOP+v//9b/PLLYr+4v/yl7847r/11lsmmCqKBlqPPvqo4/6sWbPMtSpKeHi4jBs3znH/gw8+kP3797v98H7qqacc9z/66COTC+SOc8D46aefmh4wd8aPH+94U9Ohu82bN7ttO3bsWKlSJX+dly+//NLkGrnzyCOPSPXq1c328uXL5YcffnDbVoPcqKgos71mzRpZvXq127bai6dDjmrt2rXy9ddfu2172223mSFJtXHjRlmyZInbtsOHD5cWLVqY7V9++UUWLFjgtu0NN9wgbdq0Mdvaazh//ny3bQcNGiTt2rUz23v27JEPP/zQbdtrr73W/B9SBw4ckPfff9/l8X76T2WR5XM3ycqGl0pgdEuz35Z6UnK3L3d73IB6rSWoQf752tITJffXZe7b1m0hQY0uy2+bmSq5vyx237ZOjAQ1uTy/bXam5G7+3H3bWk0k6KL8n82WmyO5mz5137ZGQwmK6eq4n7Phf+7bRkZL0MU9z7b96RORPDcfGlXrSHBsn7Ntf14gkuMm8AmvIcGt+59tu+ULkayi30+kUoQEXxJ3tu3WL0Uy3OTghYZL8KUDzrbd9rVI2qmi2waHSnC7QWfb7lglklL0e48EBknw5UMcd3N3rxFb4tme7kKH7jjsbNu9P4jt1B9u2wa1/7MEBOV/jOXu+1FsJ/a7b3vZ9RIQEpbfdv9PYju2133btn+SgLD895Pcg5vFlrDLfds2V0tA5cj8tod+FVu8+/e0oFb9JKBKTbOdd2Sn5P2xxW3bwBa9JTAi/70n7+geyTvg/vMrsHkPCaxeL7/t8d8l7/f17ts26yKBNRvltz15UPJ+W1sq7xFvvvmmlITfghONqlStWrVc9ut9+2Pnc/z4cXnooYfkxhtvNPklJTnu5MmTzXAQgJJZt++kbN29z2zXDkiVgYVje4efDpySn3/Lb1s9IF3+fI62W/5IlA2/57etGpApw87Rdlt8kqw9kN82TLJlRGX3bXclpMi3f+S3DZZcGXmOtr8dT5VVh/PbqjvO0fbAyXT5es3ZtrdWsklIfidcIfGJ6bLUqe3wSnlSyU3bYymZssip7Q1hOVLNzQD9qbQs+cyp7eCwLKnhpm1yRo7Md2p7XVim1HHTNiM7Tz50antNaLrUy/+juZDsXJu859S2f2i6NHLTVr3l1LZPaKpcdI62s77/XXIkv0GPkBS5+ByfaLPX7pdMye897hKSJK3O0Xbu+oOSYssPZDoGJ0rbs53OhczfeEhO2/I/W9oFn5L252j72abDctyW/8fmJcEnpdM52n7xS7wcyUs127FBJ6Tr2U7nQr789Yj8kZe/MnjzoOPS8xxtl28/Kr/n5feANA08KVfm/5iW47dhnZ07d0psbKysWLFCrrzySsd+DTZ036+//nrO52tvQr9+/cxf7suWLZNq1aqV6LhF9Zw0atSIYR2GdRjW8aB7NyktW+b8eEBSMnN0DFECAu3DL3n6BPf/kV3a2tz3LFiwrWmfm+OjtiIBgcG+b6vNg7xsqx9gttJom6sX2idttVfGPuxaNtoGSkBAoBdt8/Qf37d1+v95RbOa0vPiOi5NK+Swjn7w67jXoUOHXPbr/caNG5/zuRo46EwdvRBLly51BCYlOW5YWJi5FaTDCs7jpe540sabts55ImWhrfN4fFlo6/wfqry11f8Hnr7WStq2dmioPHx1K4+eD6DsCSyl9xO3xxA/0byGbt26yeefnx0LTk1NNeP0V111lWOfjnM554rYAxOlPSaaS+HNcQEAgDX5dbaOJhdqwPDXv/5VunbtamqR6PTgn3/+WapWrepIMtTkwq1bt5qE0z59+piA5d1333UJTDRBtW7duh4f19PZOt52SQEAUFEllfAz1G/DOqp3796ycuVKMyNF65LodN/Zs2e7BBAadNin6+kMFx3b032awOrsySeflOuuu87j4wIAAGvya8+JldFzAgCAfz5DLbG2DgAAgB3BCQAAsBSCEwAAYCkEJwAAwFIITgAAgKUQnAAAAEvxa50TK7PPsLavTgwAADxj/+z0tloJwYkbycnJjrV6AACAd5+lBZeZ8QRF2NzQVRUPHz5soj5dMPDgwYOUsb8A7KtBc70vDK73hcX15npXlNd3tWrVTGBSv359sxBgcdFz4oZezIYNGzq6prTCHWvsXDhc7wuL6831Ls94ffvnenvTY2JHQiwAALAUghMAAGApBCfnERYWJhMnTjRfUfq43hcW15vrXZ7x+i6715uEWAAAYCn0nAAAAEshOAEAAJZCcAIAACyFOifnceDAATl69Ki0aNGCOiel6Ndff5VTp0657KtZs6a0bt26NL9thZKeni6bNm2SBg0aSJMmTYpsc/r0admzZ49ER0ebOj/wXmJiomzdulWaNWsm9erVc3ksPj5e9u7dW+g5PXr04JJ7WTRTX7daNPOiiy6S0NDQIttpYc1Dhw5J8+bNpUaNGlxrL+Xm5srOnTvNdW7atKkEB7uGErt375aEhASXfVqU7bLLLvP8m9hQpLS0NNugQYNs4eHhttjYWFvlypVt//73v7lapSQuLs7WoEEDW/fu3R23xx9/nOvtAwkJCbZHH33UVq9ePVulSpVsTzzxRJHtXnzxRfN4q1atzOv9pptusmVmZvI7KKZ9+/bZRo8ebYuOjrYFBQXZXn/99UJtdJ9eY+fXu95QfC+//LJ572jRooUtJibGFhUVZZs7d65Lm+zsbNvIkSPN67t169a2sLAw2/PPP8/lLqbc3FzbpEmTzDXW94nGjRub29KlS13a3XbbbbY6deq4vLbvuuuuYn0vek7c+Pvf/27+yvztt9+kbt26Mn/+fLnxxhvliiuukA4dOngRa+J8br75Zpk2bRoXysf2799vSkrrX/G9evUqss3q1avl8ccfl6VLl8rVV19tegw7deokU6ZMMf8X4Dn9i7Jjx47y8ssvn7P3SZfF+Pbbb7m0JaQ9rhs3bjTv0+qVV16RkSNHmvdp7SFRL730kixZskS2bdtmelZWrFghV111lfk96esdnsnMzDS9U/oar169utl+6qmn5IYbbpDff/9datWq5Wir1/W///2veIuckyLoBX/vvffknnvucbzg9eLHxsaa/SgdKSkpsmHDBvPB6O1KlihMg4wxY8aYYTJ33n33XRN429+o9YNz1KhRZj+KJy4uzrx3VK1a9bxd47/88ots375dsrOzucxeevbZZx3v0+q+++4z1/PHH3907NPX8a233moCE9W3b1/p1q0br+9iqly5skyaNMkEJiogIMBcb33v3rx5c6FhZA0a9+3bZ4bdiovgpAi6aNGJEycK9ZDom7z2pqB0zJo1S+6++24zLtmmTRtZt24dl/oC0dd1wdd7586dTa9LwVwg+IbmnGhv7DXXXCO1a9eWmTNncml9QP/AUfZek9TUVNm1a1eRr2/ez0tu/fr15mtMTIzL/oULF8qdd95prrPmXX311VfFOi7BSRFOnjxpvjp3Udnv2x+Db912221y7Ngx+fnnn02y4OWXXy6DBw/mg/EC0dd1Ua93+2PwrUsuuUR27Nhhek00ANShiPvvv1+WLVvGpS4BXahVe62090o/FJU9uOb93PeOHDkijz76qOlldU6yHzhwoHkf194U/Tpo0CAZOnSo6RX3FMFJEUJCQszXjIyMQt1U7rLAUTLDhw832dyqUqVK5s1aX/irVq3i0l6g13xRr3fFa973+vTpY2YA2t1xxx1mWG3evHml8N0qhrS0NPOhqK/lOXPmOPbzfl469I8WDQK1V+SNN95weUwDEXswqDN5NJcwJydHvvjiC4+PT0JsEXS8XcfSdMqZM72vj6H0aX6EvqkU/B2gdOhfPUW93vV3oNOKUfo0b4LXu3c0kL7uuuvMB6YmuzrnV9WpU0fCw8N5P/ch7Y3q37+/+YNy8eLF5vqei76P6NTt4ry+6Tkpgl5w/Svm888/d+xLTk52ZHjD9xngGlU7W7lypUlq0+5vlD59XeuQgv4u7BYsWGD+wrf/5Qnf0TyIgsMRmmPF6937wESHhfU9WoMRZ4GBgSYB1vn9PCsry8ze4f3c+8BEAxKd3Vcw8VsTvQv2wmodK60xU5zXNwv/uaEvcu2y0umVOnb52muvmURZTaCqUqWKF79SuKNT0IYMGWKSYTWJTcfh//GPf5hsev2ARMlokGdPLtaxYb2ummEfEREhl156qeMNp127dubNQx/T4bTp06ebKcZdunThV1AMOnNBc6fUtddeK6NHjzavb0161Rl/Sj8sdVq3JtlrYKLTjnUYU39P9evX53oXg75P68wcnZHjHJjozBwtOKh++ukn6d69u8lH0Q/Wt99+2yTOak6E/l7gGQ069HWrPSDvvPOOS2DSsmVLc/319a/vGfp+3qpVK/P+/vzzz5tyBvq+4ukfOwQn57BmzRqZMWOGqRCrM0iefPJJiYqK8vDXiOLQioJ6rbUOgXZv6wwGzUPR4TWUjHZ1X3/99YX264yoN99803Ff33CmTp1q/srRiqYPPfSQ6UFE8WgNiLvuuqvQ/iuvvFKee+45R0+svt6/++47k9OjCeAPPvggVaiLSf9K7927d5GP6fXU2kl2Oq311VdfNa9zDRKfeOIJhumLSXun/vznPxf5mNZDspci+OOPP8wfNxqka+6J/o505k7BSrLnQnACAAAshZwTAABgKQQnAADAUghOAACApRCcAAAASyE4AQAAlkJwAgAALIXgBAAAWArBCQC/0uqd3377bbn7XgC8x8J/AIr0yy+/mGqxBWlFUy3H7itadlwrSvbo0cOnvwktT65rBWnZ8tL+XgB8i+AEQJHmzZtnyn0PGDDAZb+uLeXL4KS06Popx48fdwlOOnToYNZcAWBtBCcA3NJ1jubOnVvkY3v37jW9K4MHDy60TtLWrVvNGhy6VtKWLVvM/sjISGnbtq00bNjwvGta6YJi7du3d+zTNToSExMd66ic77h6Xnp+uoaN/fx1bRtdIystLa3Q+iy64F5CQoJcfPHFhVZO1cdsNptZJFEX/tSFzXTNoerVq7u0279/v/m5dS0RXStHe5gAeIfgBIDXK5RqAKKrSNtX21VPP/20WQlZH9NF8D777DPHAoS60J0uEKaLrrmjiw/q6tTOwcl///tf88FvD07Od1w9p3379plhHXs7DToKDutoQKKLTJ44ccIshLh27Vq56qqr5MMPP5SgoCDT5o033jCr1+rP26RJE7N6cHx8vFmx2f5z6yraL7zwgjmuBi+60vDHH38sMTExvLoALxCcAHBLP2gL9pzoisUaJOiHufZEfPDBBy6r7S5cuNAEE0oDFOdVTLXnQZdTHzZsmDRr1szrK3++4954442yYsUKM6zjrudH6UrjGoRoMKPDVRrQtGvXziwHf8899zjaaTCkq9rqEvDai9K/f3958cUXTTsNxPTnX7RokQls7L1H6enpXv98QEVHcALArdTUVEfPg51+eNt7MG699VbTs2APTj755BMJCwtzyVPRXgQNHo4ePWqGUCIiIswHfUmCE18cV4OMjz76SGbOnGkCE6X5KCNGjDABjXNw0qdPHxOYqICAAOnVq5d89dVX5n5gYKD5mXUoqV+/fua+Dg8B8B7BCQCvck7U8OHDzVDKDz/8IF27djW9KNp7oR/WSoc27rrrLjO80bhxY7Nfh1o0oCgJXxz32LFjJv+kYDCjx1y2bJnLvpo1a7rc1++nwzxKe15mz54tY8aMkcmTJ5vA7eabb5YbbrihRD8jUJERnADwWoMGDUyvggYl2uugQyl6s3vooYdMr4p+tatdu7bptXBHex7y8vJc9tkDgZIct6AaNWqYwEJzVpzpfT1WcQwaNMjcdDhnyZIlcvfdd8uBAwdMwAKg+CjCBqBEbrnlFjPtWHsPdMZMz549zX4datGcj5YtWzrarlq1yiSfni/g2bNnj+O+BirffPON476nx9UZPwWDGmchISHSuXNnMxTlfOxPP/20WHVQNLfk1KlTZluHcx5++GETqGhyLQDv0HMCoFgJsWro0KHmw13p8MUDDzwgzz77rOnJ0JwMpb0S1113nTzyyCPy2GOPmWEUrZsSHh5+ziuueSw6NDJ27FgTgMyfP9/MkKlfv36xjtuxY0dT62TGjBlmWEanEhf00ksvSd++fc0xNVDRHBRN6h0/frzHrwpt361bN7n++uvNlGadDaTDTjozCIB3CE4AFEnreuiQTcGEWKUfxPbgRBNRJ0yYYOqO3H777S7tdNaOJsx+//33pv7Hl19+afa1aNHCbWE0LZq2cuVKEyjo9GENQrRn4uDBg8U6rs7Y0V4NzYfR5FmdSlzwe2mejCbRzpo1y5S1159XpxE7D+toTZOCw0ytW7eWuLg4sx0VFWWq0dqPocNFmrOiAQsA7wTYijNICwAAUMrIOQEAAJZCcAIAACyF4AQAAFgKwQkAALAUghMAAGApBCcAAMBSCE4AAIClEJwAAABLITgBAACWQnACAAAsheAEAABYCsEJAAAQK/l/em7G3uM/hfgAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "best_x = result.X\n", + "print(f\"Processors: {int(best_x[0])}; size: {best_x[1]:.3f}\")\n", + "print(f\"Best objective: {result.F[0]:.4f}; known optimum: 0.2000\")\n", + "\n", + "values = result.pop.get(\"F\")[:, 0]\n", + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "ax.plot(np.arange(1, len(values) + 1), np.minimum.accumulate(values))\n", + "ax.axhline(0.2, color=\"grey\", linestyle=\"--\", label=\"Known optimum\")\n", + "ax.set(xlabel=\"Evaluations\", ylabel=\"Best objective observed\")\n", + "ax.legend()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/index.md b/docs/index.md index d101b7b..d9c9f5a 100644 --- a/docs/index.md +++ b/docs/index.md @@ -229,6 +229,16 @@ if __name__ == '__main__': DemoArchOptProblem().print_stats() ``` +A previously generated reference front can be stored as a NumPy array and loaded for use by experiment metrics: + +```python +from sb_arch_opt.pareto_front import load_pareto_front + +reference_pf = load_pareto_front('reference_front/pf.npy') +# Equivalently, pass the directory containing pf.npy: +reference_pf = load_pareto_front('reference_front') +``` + ### Architecture Optimization Problem with an Explicit Design Space Model If you are implementing an architecture optimization problem with design space hierarchy (i.e. some of your design diff --git a/mkdocs.yml b/mkdocs.yml index 1a08fc5..eeaae42 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -28,6 +28,7 @@ nav: - 'SMARTy': 'algo/smarty.md' - Tutorials: - SBArchOpt: tutorial.ipynb + - Graph-based Optimization: guide_graph_based_optimization.ipynb - 'Test Problems': 'test_problems.md' - API Reference: - Problem Definition: 'api/problem.md' diff --git a/sb_arch_opt/algo/arch_sbo/__init__.py b/sb_arch_opt/algo/arch_sbo/__init__.py index 0a0e47b..f051a9e 100644 --- a/sb_arch_opt/algo/arch_sbo/__init__.py +++ b/sb_arch_opt/algo/arch_sbo/__init__.py @@ -1 +1,2 @@ from .api import * +from .graph import * diff --git a/sb_arch_opt/algo/arch_sbo/graph/__init__.py b/sb_arch_opt/algo/arch_sbo/graph/__init__.py new file mode 100644 index 0000000..4f4d46c --- /dev/null +++ b/sb_arch_opt/algo/arch_sbo/graph/__init__.py @@ -0,0 +1,7 @@ +from .representation import * +from .kernels import * + +from sb_arch_opt.algo.arch_sbo.models import HAS_SMT + +if HAS_SMT: + from .surrogate import * diff --git a/sb_arch_opt/algo/arch_sbo/graph/kernels.py b/sb_arch_opt/algo/arch_sbo/graph/kernels.py new file mode 100644 index 0000000..7505121 --- /dev/null +++ b/sb_arch_opt/algo/arch_sbo/graph/kernels.py @@ -0,0 +1,736 @@ +""" +MIT License + +Copyright: (c) 2026, Deutsches Zentrum fuer Luft- und Raumfahrt e.V. +Contact: jasper.bussemaker@dlr.de + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" + +from abc import ABC, abstractmethod +from collections import Counter, OrderedDict, defaultdict +from dataclasses import dataclass +from itertools import combinations +import math +from typing import Any, Dict, Hashable, List, NamedTuple, Optional, Sequence, Tuple + +import numpy as np +from scipy import sparse + +from .representation import GraphRepresentation, NodeLabeling, SizingFeature + + +__all__ = [ + "ThetaParameter", + "GraphKernel", + "LdWloa", + "SizingKernel", + "EdgeMultiplicityKernel", + "CompositeGraphKernel", +] + + +class ThetaParameter(NamedTuple): + """Bounds, optimizer scale, and initial value of one kernel parameter.""" + + lower: float + upper: float + scale: str + initial: float + + +class GraphKernel(ABC): + """Compare decoded graphs for a graph surrogate. + + Values in ``theta`` follow the order from ``get_theta_parameters``. + ``fit_transform`` sets the training graphs that later ``transform`` calls + compare against. + """ + + def get_theta_parameters(self) -> Sequence[ThetaParameter]: + """Describe tunable parameters in the order used by ``theta``.""" + return () + + @abstractmethod + def fit_transform( + self, + graphs: Sequence[GraphRepresentation], + theta: Optional[np.ndarray] = None, + ) -> np.ndarray: + """Fit on training graphs and return their square kernel matrix. + + :param graphs: Training representations; rows and columns follow this order. + :param theta: Kernel parameters, or ``None`` for implementation defaults. + :return: Matrix with shape ``(n_train, n_train)``. + """ + pass + + @abstractmethod + def transform( + self, + graphs: Sequence[GraphRepresentation], + theta: Optional[np.ndarray] = None, + ) -> np.ndarray: + """Compare new graphs with the graphs from the last ``fit_transform``. + + :param graphs: New representations, in matrix-row order. + :param theta: Kernel parameters, or ``None`` for implementation defaults. + :return: Matrix with shape ``(n_new, n_train)``. + """ + pass + + +@dataclass(frozen=True) +class _WLFeatures: + counts: Tuple[Tuple[Counter, ...], ...] + + +@dataclass(frozen=True) +class _WLTrainingSnapshot: + features: Tuple[_WLFeatures, ...] + level_kernels: Tuple[Tuple[np.ndarray, ...], ...] + + +def _wl_features(structure_key, cutoff: int, intern) -> _WLFeatures: + """Compute WL counts using compact labels shared across kernel copies.""" + labelings, neighbors = structure_key + counts_by_labeling = [] + for labeling_index, raw_labels in enumerate(labelings): + current = tuple(intern(("base", labeling_index, label)) for label in raw_labels) + counts_by_level = [] + for level in range(cutoff + 1): + counts_by_level.append(Counter(current)) + if level < cutoff: + current = tuple( + intern(( + "wl", + current[index], + tuple(sorted(current[neighbor] for neighbor in adjacent)), + )) + for index, adjacent in enumerate(neighbors) + ) + counts_by_labeling.append(tuple(counts_by_level)) + return _WLFeatures(tuple(counts_by_labeling)) + + +class _WLCache: + """Shared by copies of one kernel; training data is never LRU-evicted.""" + + def __init__(self, cache_size: Optional[int]): + self.cache_size = cache_size + self._label_ids = {} + self.training_features = {} + self.candidate_features = OrderedDict() + self.training_snapshot: Optional[_WLTrainingSnapshot] = None + + def __deepcopy__(self, memo): + memo[id(self)] = self + return self + + def _intern(self, label) -> int: + if label not in self._label_ids: + self._label_ids[label] = len(self._label_ids) + return self._label_ids[label] + + def features(self, structure_key, cutoff: int, training: bool) -> _WLFeatures: + key = (cutoff, structure_key) + features = self.training_features.get(key) + if features is not None: + return features + + features = self.candidate_features.get(key) + if features is None: + features = _wl_features(structure_key, cutoff, self._intern) + if training: + self.training_features[key] = features + self.candidate_features.pop(key, None) + else: + self.candidate_features[key] = features + self.candidate_features.move_to_end(key) + if self.cache_size is not None and len(self.candidate_features) > self.cache_size: + self.candidate_features.popitem(last=False) + return features + + +def _same_objects(left, right): + return ( + left is not None + and right is not None + and len(left) == len(right) + and all(a is b for a, b in zip(left, right)) + ) + + +class LdWloa(GraphKernel): + """Learned-depth WLOA over one or more node labelings. + + Copies share theta-independent WL features and the latest training matrices. + ``cache_size`` limits transient candidate graphs; ``None`` disables eviction. + Training graph features are retained regardless of this limit. + """ + + def __init__( + self, + cutoff: int, + labelings: Sequence[NodeLabeling], + cache_size: Optional[int] = 8192, + ): + if cutoff < 0: + raise ValueError("cutoff must be non-negative") + if len(labelings) == 0: + raise ValueError("at least one node labeling is required") + if len({labeling.key for labeling in labelings}) != len(labelings): + raise ValueError("node labeling keys must be unique") + if cache_size is not None and cache_size < 0: + raise ValueError("cache_size must be non-negative or None") + + self.cutoff = cutoff + self.labelings = tuple(labelings) + self._cache = _WLCache(cache_size) + self._train_features: Optional[Tuple[_WLFeatures, ...]] = None + self._train_level_kernels = None + self._train_counts = None + self._transform_features = None + self._transform_level_kernels = None + self._transform_counts = None + + @property + def cache_size(self) -> Optional[int]: + return self._cache.cache_size + + def get_theta_parameters(self) -> Sequence[ThetaParameter]: + depth = [ThetaParameter(0.0, 1.0, "linear", 0.5)] * (self.cutoff + 1) + mixture = [ + ThetaParameter(0.0, 1.0, "linear", (remaining - 1) / remaining) + for remaining in range(len(self.labelings), 1, -1) + ] + return depth + mixture + + def fit_transform(self, graphs, theta=None): + features = tuple(self._build_features(graph, training=True) for graph in graphs) + counts = self._counts_by_labeling(features) + self._train_level_kernels = self._train_kernels(features, counts) + self._train_features = features + self._train_counts = counts + self._clear_transform_cache() + return self._combine_kernels( + self._train_level_kernels, + counts, + counts, + theta, + ) + + def transform(self, graphs, theta=None): + if self._train_features is None: + raise RuntimeError("fit_transform must be called before transform") + features = tuple(self._build_features(graph) for graph in graphs) + if _same_objects(features, self._transform_features): + counts = self._transform_counts + level_kernels = self._transform_level_kernels + else: + counts = self._counts_by_labeling(features) + level_kernels = tuple( + tuple( + self._kernel_matrix_for_level(left, right) + for left, right in zip(left_by_level, right_by_level) + ) + for left_by_level, right_by_level in zip(counts, self._train_counts) + ) + self._transform_features = features + self._transform_counts = counts + self._transform_level_kernels = level_kernels + return self._combine_kernels( + level_kernels, + counts, + self._train_counts, + theta, + ) + + def _build_features( + self, representation: GraphRepresentation, training: bool = False + ) -> _WLFeatures: + graph = representation.graph + nodes = tuple(graph.nodes) + if graph.is_directed(): + neighbors = { + node: tuple(set(graph.predecessors(node)) | set(graph.successors(node))) + for node in nodes + } + else: + neighbors = {node: tuple(graph.neighbors(node)) for node in nodes} + + node_index = {node: index for index, node in enumerate(nodes)} + structure_key = ( + tuple( + tuple(representation.node_labels[labeling][node] for node in nodes) + for labeling in self.labelings + ), + tuple( + tuple(sorted(node_index[neighbor] for neighbor in neighbors[node])) + for node in nodes + ), + ) + return self._cache.features(structure_key, self.cutoff, training) + + def _combine_kernels(self, level_kernels, left_counts, right_counts, theta): + depth_weights, labeling_weights = self._theta_parts(theta) + matrices = [ + self._labeling_kernel( + level_kernels[index], + left_counts[index], + right_counts[index], + depth_weights, + ) + for index in range(len(self.labelings)) + ] + shape = level_kernels[0][0].shape + result = np.zeros(shape, dtype=float) + for weight, matrix in zip(labeling_weights, matrices): + result += weight * matrix + return result + + @staticmethod + def _labeling_kernel(level_kernels, left_counts, right_counts, depth_weights): + matrix = np.zeros_like(level_kernels[0], dtype=float) + left_diag = np.zeros((len(left_counts[0]),), dtype=float) + right_diag = np.zeros((len(right_counts[0]),), dtype=float) + + for weight, level_kernel, left, right in zip( + depth_weights, + level_kernels, + left_counts, + right_counts, + ): + matrix += weight * level_kernel + left_diag += weight * np.fromiter( + (sum(count.values()) for count in left), + dtype=float, + count=len(left), + ) + right_diag += weight * np.fromiter( + (sum(count.values()) for count in right), + dtype=float, + count=len(right), + ) + + denominator = np.sqrt(np.outer(left_diag, right_diag)) + return np.divide(matrix, denominator, out=np.zeros_like(matrix), where=denominator > 0.0) + + def _counts_by_labeling(self, features): + return tuple( + tuple( + tuple(feature.counts[labeling_index][level] for feature in features) + for level in range(self.cutoff + 1) + ) + for labeling_index in range(len(self.labelings)) + ) + + def _train_kernels(self, features, counts): + if _same_objects(features, self._train_features): + return self._train_level_kernels + + snapshot = self._cache.training_snapshot + old_features = snapshot.features if snapshot is not None else None + old_kernels = snapshot.level_kernels if snapshot is not None else None + if _same_objects(features, old_features): + return old_kernels + + old_count = len(old_features) if old_features is not None else 0 + grows_previous_train = ( + old_count < len(features) + and _same_objects(features[:old_count], old_features) + and old_kernels is not None + ) + if not grows_previous_train: + kernels = tuple( + tuple(self._kernel_matrix_for_level(level, level) for level in by_level) + for by_level in counts + ) + else: + extended = [] + for labeling_index, by_level in enumerate(counts): + labeling_kernels = [] + for level, level_counts in enumerate(by_level): + old_kernel = old_kernels[labeling_index][level] + old_counts = level_counts[:old_count] + new_counts = level_counts[old_count:] + new_old = self._kernel_matrix_for_level(new_counts, old_counts) + new_new = self._kernel_matrix_for_level(new_counts, new_counts) + kernel = np.empty((len(features), len(features)), dtype=float) + kernel[:old_count, :old_count] = old_kernel + kernel[old_count:, :old_count] = new_old + kernel[:old_count, old_count:] = new_old.T + kernel[old_count:, old_count:] = new_new + labeling_kernels.append(kernel) + extended.append(tuple(labeling_kernels)) + kernels = tuple(extended) + + self._cache.training_snapshot = _WLTrainingSnapshot(features, kernels) + return kernels + + def _clear_transform_cache(self): + self._transform_features = None + self._transform_level_kernels = None + self._transform_counts = None + + @staticmethod + def _kernel_matrix_for_level(left_counts, right_counts): + labels = set() + for counts in left_counts: + labels.update(counts) + for counts in right_counts: + labels.update(counts) + label_index = {label: index for index, label in enumerate(labels)} + left = LdWloa._counts_to_csc(left_counts, label_index) + right = LdWloa._counts_to_csc(right_counts, label_index) + + matrix = np.zeros((len(left_counts), len(right_counts)), dtype=float) + for column in range(len(label_index)): + left_start, left_stop = left.indptr[column], left.indptr[column + 1] + right_start, right_stop = right.indptr[column], right.indptr[column + 1] + if left_start == left_stop or right_start == right_stop: + continue + left_rows = left.indices[left_start:left_stop] + right_rows = right.indices[right_start:right_stop] + left_values = left.data[left_start:left_stop] + right_values = right.data[right_start:right_stop] + matrix[np.ix_(left_rows, right_rows)] += np.minimum( + left_values[:, None], + right_values[None, :], + ) + return matrix + + @staticmethod + def _counts_to_csc(counts, label_index): + rows = [] + columns = [] + values = [] + for row, counter in enumerate(counts): + for label, value in counter.items(): + rows.append(row) + columns.append(label_index[label]) + values.append(float(value)) + return sparse.csc_matrix( + (values, (rows, columns)), + shape=(len(counts), len(label_index)), + ) + + def _theta_parts(self, theta): + if theta is None: + theta = [parameter.initial for parameter in self.get_theta_parameters()] + theta = np.clip(np.asarray(theta, dtype=float).ravel(), 0.0, 1.0) + n_depth = self.cutoff + 1 + if theta.size != n_depth + len(self.labelings) - 1: + raise ValueError("incorrect number of WLOA parameters") + + remaining = 1.0 + labeling_weights = [] + for fraction in theta[n_depth:]: + labeling_weights.append(remaining * (1.0 - fraction)) + remaining *= fraction + labeling_weights.append(remaining) + return theta[:n_depth], np.asarray(labeling_weights) + + +class SizingKernel(GraphKernel): + """Kernel over normalized numeric and categorical sizing features.""" + + def __init__( + self, + features: Sequence[SizingFeature], + power: float = 2.0, + ): + if len(features) == 0: + raise ValueError("at least one sizing feature is required") + if not 0.0 < power <= 2.0: + raise ValueError("sizing power must be in (0, 2]") + self.features = tuple(features) + self.power = power + self._numeric = np.array([feature.kind == "numeric" for feature in features]) + self._categorical = ~self._numeric + self._train_graphs = None + self._train_values: Optional[np.ndarray] = None + self._train_powered_distances = None + self._transform_graphs = None + self._transform_powered_distances = None + + def get_theta_parameters(self) -> Sequence[ThetaParameter]: + return [ThetaParameter(1e-6, 20.0, "log", 1.0)] * len(self.features) + + def fit_transform(self, graphs, theta=None): + graphs = tuple(graphs) + if not _same_objects(graphs, self._train_graphs): + values = self._arrays(graphs) + self._train_powered_distances = self._training_distances(graphs, values) + self._train_graphs = graphs + self._train_values = values + self._clear_transform_cache() + return self._kernel_from_distances(self._train_powered_distances, theta) + + def transform(self, graphs, theta=None): + if self._train_values is None: + raise RuntimeError("fit_transform must be called before transform") + graphs = tuple(graphs) + if not _same_objects(graphs, self._transform_graphs): + values = self._arrays(graphs) + self._transform_powered_distances = self._powered_distances( + values[:, None, :], + self._train_values[None, :, :], + ) + self._transform_graphs = graphs + return self._kernel_from_distances(self._transform_powered_distances, theta) + + def _arrays(self, graphs): + if len(graphs) == 0: + shape = (0, len(self.features)) + return np.empty(shape) + return np.asarray([graph.sizing_values for graph in graphs], dtype=float) + + def _training_distances(self, graphs, values): + old_graphs = self._train_graphs + old_count = len(old_graphs) if old_graphs is not None else 0 + grows_previous_train = ( + old_count < len(graphs) + and _same_objects(graphs[:old_count], old_graphs) + and self._train_powered_distances is not None + ) + if not grows_previous_train: + return self._powered_distances(values[:, None, :], values[None, :, :]) + + distances = np.empty((len(graphs), len(graphs), len(self.features)), dtype=float) + distances[:old_count, :old_count] = self._train_powered_distances + new_old = self._powered_distances( + values[old_count:, None, :], + values[None, :old_count, :], + ) + distances[old_count:, :old_count] = new_old + distances[:old_count, old_count:] = np.swapaxes(new_old, 0, 1) + distances[old_count:, old_count:] = self._powered_distances( + values[old_count:, None, :], + values[None, old_count:, :], + ) + return distances + + def _powered_distances(self, left_values, right_values): + distances = np.zeros(np.broadcast_shapes(left_values.shape, right_values.shape), dtype=float) + delta = np.abs(left_values - right_values) + distances[..., self._numeric] = delta[..., self._numeric] + distances[..., self._categorical] = ( + left_values[..., self._categorical] != right_values[..., self._categorical] + ) + + return distances**self.power + + def _kernel_from_distances(self, powered_distances, theta): + if theta is None: + theta = [parameter.initial for parameter in self.get_theta_parameters()] + theta = np.asarray(theta, dtype=float).ravel() + if theta.size != len(self.features): + raise ValueError(f"SizingKernel theta size must be {len(self.features)}, got {theta.size}") + return np.exp(-np.tensordot(powered_distances, theta, axes=([-1], [0]))) + + def _clear_transform_cache(self): + self._transform_graphs = None + self._transform_powered_distances = None + + +class EdgeMultiplicityKernel(GraphKernel): + """Kernel over multiplicities of a selected edge kind.""" + + def __init__(self, edge_kind: Hashable, endpoint_labeling: NodeLabeling, gamma0: float = 1.0): + self.edge_kind = edge_kind + self.endpoint_labeling = endpoint_labeling + self.gamma0 = gamma0 + self._train_features: Optional[Tuple[Dict[Tuple[Hashable, Hashable], float], ...]] = None + self._feature_cache = {} + self._train_distances = None + self._transform_features = None + self._transform_distances = None + + def get_theta_parameters(self) -> Sequence[ThetaParameter]: + return [ThetaParameter(1e-6, 20.0, "log", self.gamma0)] + + def fit_transform(self, graphs, theta=None): + features = tuple(self._features(graph) for graph in graphs) + if not _same_objects(features, self._train_features): + self._train_distances = self._training_distances(features) + self._train_features = features + self._clear_transform_cache() + return self._kernel(self._train_distances, theta) + + def transform(self, graphs, theta=None): + if self._train_features is None: + raise RuntimeError("fit_transform must be called before transform") + features = tuple(self._features(graph) for graph in graphs) + if not _same_objects(features, self._transform_features): + self._transform_distances = self._distances(features, self._train_features) + self._transform_features = features + return self._kernel(self._transform_distances, theta) + + def _features(self, representation): + cached = self._feature_cache.get(id(representation)) + if cached is not None and cached[0] is representation: + return cached[1] + + graph = representation.graph + counts = defaultdict(int) + edges = graph.edges(keys=True, data=True) if graph.is_multigraph() else graph.edges(data=True) + for edge in edges: + source, target, data = (edge[0], edge[1], edge[-1]) + if data.get("kind") == self.edge_kind: + counts[(source, target)] += 1 + + labels = representation.node_labels[self.endpoint_labeling] + values = defaultdict(float) + for (source, target), multiplicity in counts.items(): + values[(labels[source], labels[target])] += math.log1p(multiplicity) + features = dict(values) + self._feature_cache[id(representation)] = (representation, features) + return features + + def _training_distances(self, features): + old_features = self._train_features + old_count = len(old_features) if old_features is not None else 0 + grows_previous_train = ( + old_count < len(features) + and _same_objects(features[:old_count], old_features) + and self._train_distances is not None + ) + if not grows_previous_train: + return self._distances(features, features) + + distances = np.empty((len(features), len(features)), dtype=float) + distances[:old_count, :old_count] = self._train_distances + new_old = self._distances(features[old_count:], features[:old_count]) + distances[old_count:, :old_count] = new_old + distances[:old_count, old_count:] = new_old.T + distances[old_count:, old_count:] = self._distances( + features[old_count:], + features[old_count:], + ) + return distances + + @staticmethod + def _distances(left, right): + distances = np.empty((len(left), len(right)), dtype=float) + for i, values_left in enumerate(left): + for j, values_right in enumerate(right): + keys = set(values_left) | set(values_right) + distances[i, j] = sum( + abs(values_left.get(key, 0.0) - values_right.get(key, 0.0)) + for key in keys + ) + return distances + + def _kernel(self, distances, theta): + if theta is None: + gamma = self.gamma0 + else: + theta = np.asarray(theta, dtype=float).ravel() + if theta.size != 1: + raise ValueError(f"EdgeMultiplicityKernel theta size must be 1, got {theta.size}") + gamma = float(theta[0]) + return np.exp(-gamma * distances) + + def _clear_transform_cache(self): + self._transform_features = None + self._transform_distances = None + + +class CompositeGraphKernel(GraphKernel): + """Composition of structural, sizing, and edge graph-kernel branches.""" + + def __init__(self, branches: Sequence[GraphKernel], composition: str = "additive"): + if len(branches) == 0: + raise ValueError("at least one kernel branch is required") + if composition not in {"additive", "additive_interaction", "multiplicative"}: + raise ValueError(f"unsupported graph kernel composition: {composition!r}") + self.branches = tuple(branches) + self.composition = composition + + def get_theta_parameters(self) -> Sequence[ThetaParameter]: + parameters: List[ThetaParameter] = [] + if self._uses_weights(): + parameters.extend( + [ThetaParameter(0.01, 100.0, "linear", 1.0)] * (self._term_count() - 1) + ) + for branch in self.branches: + parameters.extend(branch.get_theta_parameters()) + return parameters + + def fit_transform(self, graphs, theta=None): + weights, branch_theta = self._theta_parts(theta) + matrices = [ + branch.fit_transform(graphs, values) + for branch, values in zip(self.branches, branch_theta) + ] + return self._compose(matrices, weights) + + def transform(self, graphs, theta=None): + weights, branch_theta = self._theta_parts(theta) + matrices = [ + branch.transform(graphs, values) + for branch, values in zip(self.branches, branch_theta) + ] + return self._compose(matrices, weights) + + def _uses_weights(self): + return self.composition in {"additive", "additive_interaction"} and self._term_count() > 1 + + def _term_count(self): + if self.composition == "additive_interaction": + return len(self.branches) + len(tuple(combinations(self.branches, 2))) + return len(self.branches) + + def _theta_parts(self, theta): + parameters = self.get_theta_parameters() + if theta is None: + theta = [parameter.initial for parameter in parameters] + theta = np.asarray(theta, dtype=float).ravel() + if theta.size != len(parameters): + raise ValueError( + f"CompositeGraphKernel theta size must be {len(parameters)}, got {theta.size}" + ) + + cursor = 0 + weights = None + if self._uses_weights(): + n_free = self._term_count() - 1 + free_weights = np.maximum(theta[:n_free], 0.0) + raw_weights = np.concatenate(([1.0], free_weights)) + weights = raw_weights / np.sum(raw_weights) + cursor = n_free + + branch_theta = [] + for branch in self.branches: + size = len(branch.get_theta_parameters()) + branch_theta.append(theta[cursor:cursor + size]) + cursor += size + return weights, branch_theta + + def _compose(self, matrices, weights): + if self.composition == "multiplicative": + return np.prod(matrices, axis=0) + + terms = list(matrices) + if self.composition == "additive_interaction": + terms.extend(left * right for left, right in combinations(matrices, 2)) + if weights is None: + return terms[0] + return sum(weight * term for weight, term in zip(weights, terms)) diff --git a/sb_arch_opt/algo/arch_sbo/graph/representation.py b/sb_arch_opt/algo/arch_sbo/graph/representation.py new file mode 100644 index 0000000..7028c17 --- /dev/null +++ b/sb_arch_opt/algo/arch_sbo/graph/representation.py @@ -0,0 +1,117 @@ +""" +MIT License + +Copyright: (c) 2026, Deutsches Zentrum fuer Luft- und Raumfahrt e.V. +Contact: jasper.bussemaker@dlr.de + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" + +from dataclasses import dataclass, field +from typing import Any, Callable, Hashable, Literal, Mapping, Protocol, Sequence, runtime_checkable + +import numpy as np + +try: + import networkx as nx + HAS_NETWORKX = True +except ImportError: + HAS_NETWORKX = False + + +__all__ = ["HAS_NETWORKX", "NodeLabeling", "SizingFeature", "GraphRepresentation", "GraphDecoder"] + + +@dataclass(frozen=True) +class NodeLabeling: + """A named rule for labeling nodes of a domain graph. + + :param key: Identifier used to match labels from a decoder with a kernel. + Use a different key when the meaning of the labels changes. + :param function: Maps a source-graph node to a hashable label. A decoder + can call it when building each graph representation. + """ + + key: str + function: Callable[[Any], Hashable] = field(compare=False, hash=False, repr=False) + + def __call__(self, node: Any) -> Hashable: + return self.function(node) + + +@dataclass(frozen=True) +class SizingFeature: + """Description of one sizing coordinate exported by a graph decoder. + + :param name: Human-readable name of the coordinate. + :param kind: ``numeric`` compares the difference between values; + ``categorical`` only checks whether values match. + """ + + name: str + kind: Literal["numeric", "categorical"] + + +def _empty_float_array() -> np.ndarray: + return np.empty((0,), dtype=float) + + +@dataclass(frozen=True) +class GraphRepresentation: + """Kernel inputs for one decoded design, independent of its source domain. + + :param graph: NetworkX graph whose nodes and edges describe the design. + :param node_labels: For each labeling declared by the decoder, a mapping + from every node in ``graph`` to its hashable label. + :param sizing_values: One-dimensional array with one value per decoder + sizing feature, in the order of ``GraphDecoder.sizing_features``. + + Do not change the graph, labels, or sizing values after decoding. Kernels + may reuse values computed from them. + """ + + graph: "nx.Graph" + node_labels: Mapping[NodeLabeling, Mapping[Hashable, Hashable]] + sizing_values: np.ndarray = field(default_factory=_empty_float_array) + + +@runtime_checkable +class GraphDecoder(Protocol): + """Domain adapter that converts design vectors into graph-kernel inputs. + + Implement this protocol for a new design space. Its labelings and sizing + features describe the data present in every representation it returns. + """ + + @property + def labelings(self) -> Sequence[NodeLabeling]: + """Labeling rules available in every decoded representation.""" + ... + + @property + def sizing_features(self) -> Sequence[SizingFeature]: + """Sizing coordinates, in the order used by ``sizing_values``.""" + ... + + def decode(self, x: np.ndarray) -> Sequence[GraphRepresentation]: + """Return one representation per row of ``x``, in the same order. + + :param x: Design vectors with shape ``(n_designs, n_variables)``. + """ + ... diff --git a/sb_arch_opt/algo/arch_sbo/graph/surrogate.py b/sb_arch_opt/algo/arch_sbo/graph/surrogate.py new file mode 100644 index 0000000..be76f3e --- /dev/null +++ b/sb_arch_opt/algo/arch_sbo/graph/surrogate.py @@ -0,0 +1,228 @@ +""" +MIT License + +Copyright: (c) 2026, Deutsches Zentrum fuer Luft- und Raumfahrt e.V. +Contact: jasper.bussemaker@dlr.de + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" + +from typing import Dict, Optional, Sequence, Tuple + +import numpy as np +from pymoo.util.normalization import Normalization +from smt.sampling_methods import LHS +from smt.surrogate_models.krg_based import KrgBased +from smt.surrogate_models.krg_based.hyperparam_optim import CobylaOptimizer +from smt.utils.kriging import cross_distances + +from .kernels import GraphKernel, ThetaParameter +from .representation import GraphDecoder, GraphRepresentation + + +__all__ = ["GraphKriging"] + + +class GraphKriging(KrgBased): + """SMT kriging model whose correlation is supplied by a graph kernel.""" + + name = "GraphKriging" + + def __init__( + self, + decoder: GraphDecoder, + kernel: GraphKernel, + normalization: Optional[Normalization] = None, + **kwargs, + ): + kwargs.setdefault("hyper_opt", "Cobyla") + super().__init__(**kwargs) + self.decoder = decoder + self.kernel = kernel + self.normalization = normalization + + self._representations_by_x: Dict[Tuple[float, ...], GraphRepresentation] = {} + self._train_kernel_matrix: Optional[np.ndarray] = None + self._train_kernel_theta: Optional[Tuple[float, ...]] = None + self._train_pairs: Optional[np.ndarray] = None + self._theta_parameters: Sequence[ThetaParameter] = () + + def _decode(self, x: np.ndarray): + x = np.asarray(x, dtype=float) + missing_keys = [] + missing_rows = [] + for row in x: + key = tuple(row) + if key not in self._representations_by_x: + missing_keys.append(key) + missing_rows.append(row) + + if missing_rows: + raw = np.asarray(missing_rows) + if self.normalization is not None: + raw = self.normalization.backward(raw) + for key, representation in zip(missing_keys, self.decoder.decode(raw)): + self._representations_by_x[key] = representation + + return tuple(self._representations_by_x[tuple(row)] for row in x) + + def _check_param(self): + requested_theta0 = np.asarray(self.options["theta0"], dtype=float).ravel() + self.options["theta0"] = np.array([float(requested_theta0[0])]) + super()._check_param() + + parameters = tuple(self.kernel.get_theta_parameters()) + if not parameters: + parameters = (ThetaParameter(0.01, 100.0, "linear", 1.0),) + theta0 = np.asarray([parameter.initial for parameter in parameters], dtype=float) + if requested_theta0.size > 1: + if requested_theta0.size != theta0.size: + raise ValueError( + f"GraphKriging expects theta0 with 1 or {theta0.size} values, " + f"got {requested_theta0.size}" + ) + theta0 = requested_theta0 + + self._theta_parameters = parameters + self.options["theta0"] = theta0 + self._theta0 = list(theta0) + self.n_param = theta0.size + + def _new_train(self): + self._train_kernel_matrix = None + self._train_kernel_theta = None + self._train_pairs = None + super()._new_train() + + def _to_optimizer_space(self, theta): + theta = np.asarray(theta, dtype=float).copy() + for index, parameter in enumerate(self._theta_parameters): + if parameter.scale == "log": + theta[index] = np.log10( + np.clip(theta[index], parameter.lower, parameter.upper) + ) + return theta + + def _from_optimizer_space(self, theta): + theta = np.asarray(theta, dtype=float).copy() + for index, parameter in enumerate(self._theta_parameters): + if parameter.scale == "log": + theta[index] = 10.0 ** theta[index] + return theta + + def _optimizer_bounds(self): + bounds = [] + for parameter in self._theta_parameters: + if parameter.scale == "log": + bounds.append([ + np.log10(parameter.lower), + np.log10(parameter.upper), + ]) + else: + bounds.append([parameter.lower, parameter.upper]) + return np.asarray(bounds) + + def _optimizer_constraints(self): + constraints = [] + for index, parameter in enumerate(self._theta_parameters): + lower = ( + np.log10(parameter.lower) + if parameter.scale == "log" + else parameter.lower + ) + upper = ( + np.log10(parameter.upper) + if parameter.scale == "log" + else parameter.upper + ) + constraints.append(lambda theta, i=index, bound=lower: theta[i] - bound) + constraints.append(lambda theta, i=index, bound=upper: bound - theta[i]) + return constraints + + def _optimize_hyperparam(self, D, use_multistart=True, limit=None): + del use_multistart, limit + self.D = D + + if self.options["hyper_opt"] not in {"Cobyla", "NoOp"}: + raise ValueError("GraphKriging supports hyper_opt='Cobyla' or 'NoOp'") + if self._eval_noise: + raise NotImplementedError("GraphKriging does not optimize the noise level") + no_optimization = self.options["hyper_opt"] == "NoOp" + self.noise0 = np.array(self.options["noise0"] if no_optimization else self._noise0) + theta0 = np.asarray(self.options["theta0"], dtype=float).ravel() + + if no_optimization: + optimal_theta = theta0 + else: + bounds = self._optimizer_bounds() + for index, parameter in enumerate(self._theta_parameters): + theta0[index] = np.clip(theta0[index], parameter.lower, parameter.upper) + theta0_optimizer = self._to_optimizer_space(theta0) + + def objective(theta): + native_theta = self._from_optimizer_space(theta) + return -self._reduced_likelihood_function(theta=native_theta)[0] + + random_start = bounds[:, 0] + self.rng.random(len(theta0)) * ( + bounds[:, 1] - bounds[:, 0] + ) + starts = [theta0_optimizer, random_start] + if self.options["n_start"] > 1: + sampler = LHS(xlimits=bounds, criterion="maximin", seed=self.rng) + starts.extend(sampler(self.options["n_start"])) + + result = CobylaOptimizer().optimize( + objective=objective, + theta_starts=np.vstack(starts), + constraints=self._optimizer_constraints(), + limit=max(12 * len(theta0), 50), + ) + optimal_theta = ( + theta0 + if result is None or "x" not in result + else self._from_optimizer_space(result["x"]) + ) + + optimal_value, optimal_parameters = self._reduced_likelihood_function( + theta=optimal_theta + ) + return optimal_value, optimal_parameters, optimal_theta + + def _matrix_data_corr(self, theta, *, x=None, kplsk_second_loop=False, **_): + if kplsk_second_loop: + raise NotImplementedError("GraphKriging does not support KPLSK") + + theta = np.asarray(theta, dtype=float).ravel() + if x is None: + x_train = self.training_points[None][0][0] + if self._train_pairs is None: + _, self._train_pairs = cross_distances(x_train) + theta_key = tuple(theta) + if self._train_kernel_theta != theta_key: + self._train_kernel_matrix = np.asarray( + self.kernel.fit_transform(self._decode(x_train), theta), + dtype=float, + ) + self._train_kernel_theta = theta_key + return self._train_kernel_matrix[ + self._train_pairs[:, 0], self._train_pairs[:, 1] + ].reshape(-1, 1) + + kernel_matrix = self.kernel.transform(self._decode(np.asarray(x)), theta) + return np.asarray(kernel_matrix, dtype=float).reshape(-1, 1) diff --git a/sb_arch_opt/algo/arch_sbo/models.py b/sb_arch_opt/algo/arch_sbo/models.py index 23c3fb8..43211cc 100644 --- a/sb_arch_opt/algo/arch_sbo/models.py +++ b/sb_arch_opt/algo/arch_sbo/models.py @@ -216,6 +216,41 @@ def get_md_kriging_model(self, kpls_n_comp: int = None, multi=True, ignore_hiera return surrogate, normalization + def get_graph_kriging_model(self, decoder, kernel, multi=True, ignore_hierarchy=False, + **kwargs_): + """Create a mixed-discrete kriging model using neutral graph representations.""" + check_dependencies() + from sb_arch_opt.algo.arch_sbo.graph.surrogate import GraphKriging + + normalization = self.get_md_normalization() + norm_ds_spec = self.create_smt_design_space_spec( + self.problem.design_space, + md_normalize=True, + ignore_hierarchy=ignore_hierarchy, + ) + kwargs = dict( + design_space=norm_ds_spec.design_space, + categorical_kernel=MixIntKernelType.GOWER, + hierarchical_kernel=MixHrcKernelType.ALG_KERNEL, + print_global=False, + print_training=False, + print_prediction=False, + print_problem=False, + print_solver=False, + ) + kwargs.update(kwargs_) + surrogate = GraphKriging( + decoder=decoder, + kernel=kernel, + normalization=normalization, + **kwargs, + ) + if ignore_hierarchy: + surrogate.supports['x_hierarchy'] = False + if multi: + surrogate = MultiSurrogateModel(surrogate) + return surrogate, normalization + @staticmethod def get_n_theta(problem: ArchOptProblemBase, surrogate: 'SurrogateModel') -> int: @@ -405,9 +440,6 @@ def train(self) -> None: model.train() - # rmse = np.linalg.norm(self.yt[:, i] - model.predict_values(self.xt)[:, 0], 2) - # print(f'TRAINED {i}: {rmse:.3g}') - if i == 0 and isinstance(model, KrgBased): try: theta0 = list(model.optimal_theta) diff --git a/sb_arch_opt/pareto_front.py b/sb_arch_opt/pareto_front.py index 9751443..32a7dbd 100644 --- a/sb_arch_opt/pareto_front.py +++ b/sb_arch_opt/pareto_front.py @@ -27,6 +27,7 @@ import pickle import hashlib import numpy as np +from pathlib import Path from typing import * import concurrent.futures import matplotlib.pyplot as plt @@ -48,7 +49,19 @@ # pymoo >= 0.6.1 from pymoo.operators.survival.rank_and_crowding.metrics import calc_crowding_distance -__all__ = ['CachedParetoFrontMixin'] +__all__ = ['CachedParetoFrontMixin', 'load_pareto_front'] + + +def load_pareto_front(path: Union[str, os.PathLike]) -> np.ndarray: + """Load a reference Pareto front from a NumPy file or a directory containing ``pf.npy``.""" + path = Path(path) + if path.is_dir(): + path = path / 'pf.npy' + + pareto_front = np.load(path, allow_pickle=False) + if pareto_front.ndim != 2: + raise ValueError(f'Expected a 2D Pareto front, got shape {pareto_front.shape}') + return pareto_front class CachedParetoFrontMixin(Problem): diff --git a/sb_arch_opt/tests/algo/test_arch_sbo_dependencies.py b/sb_arch_opt/tests/algo/test_arch_sbo_dependencies.py new file mode 100644 index 0000000..d524613 --- /dev/null +++ b/sb_arch_opt/tests/algo/test_arch_sbo_dependencies.py @@ -0,0 +1,38 @@ +import subprocess +import sys + +import pytest + + +@pytest.mark.parametrize("missing_modules", [("networkx",), ("smt",), ("networkx", "smt")]) +def test_arch_sbo_import_without_optional_dependencies(missing_modules): + result = subprocess.run( + [sys.executable, "-c", """ +import importlib.abc +import sys + +missing_modules = sys.argv[1:] + +class MissingDependency(importlib.abc.MetaPathFinder): + def find_spec(self, fullname, path=None, target=None): + if fullname.split('.')[0] in missing_modules: + raise ModuleNotFoundError(f'No module named {fullname!r}', name=fullname) + +sys.meta_path.insert(0, MissingDependency()) + +from sb_arch_opt.problem import ArchOptProblemBase +from sb_arch_opt.algo.arch_sbo import HAS_SMT, get_arch_sbo_gp +from sb_arch_opt.algo.arch_sbo.models import ModelFactory +from sb_arch_opt.algo.arch_sbo.graph import HAS_NETWORKX, GraphRepresentation, LdWloa, NodeLabeling + +for module in missing_modules: + assert module not in sys.modules +if 'smt' in missing_modules: + assert not HAS_SMT +if 'networkx' in missing_modules: + assert not HAS_NETWORKX +""", *missing_modules], + capture_output=True, + text=True, + ) + assert result.returncode == 0, result.stdout + result.stderr diff --git a/sb_arch_opt/tests/algo/test_arch_sbo_graph_kernels.py b/sb_arch_opt/tests/algo/test_arch_sbo_graph_kernels.py new file mode 100644 index 0000000..c567801 --- /dev/null +++ b/sb_arch_opt/tests/algo/test_arch_sbo_graph_kernels.py @@ -0,0 +1,358 @@ +import copy + +import numpy as np +import pytest + +from sb_arch_opt.algo.arch_sbo.graph import ( + HAS_NETWORKX, + CompositeGraphKernel, + EdgeMultiplicityKernel, + GraphRepresentation, + LdWloa, + NodeLabeling, + SizingFeature, + SizingKernel, +) + +if HAS_NETWORKX: + import networkx as nx + + +pytestmark = pytest.mark.skipif(not HAS_NETWORKX, reason="NetworkX dependency not installed") + + +TYPE = NodeLabeling("type", lambda node: node) +FAMILY = NodeLabeling("family", lambda node: node) + + +def make_graph(edges, type_labels, family_labels=None, sizing=()): + graph = nx.MultiDiGraph() + graph.add_nodes_from(type_labels) + for source, target, kind in edges: + graph.add_edge(source, target, kind=kind) + labels = {TYPE: type_labels} + if family_labels is not None: + labels[FAMILY] = family_labels + return GraphRepresentation( + graph=graph, + node_labels=labels, + sizing_values=np.asarray(sizing, dtype=float), + ) + + +def test_ld_wloa_single_labeling(): + graphs = [ + make_graph([(0, 1, "derives")], {0: "a", 1: "b"}), + make_graph([(0, 1, "derives")], {0: "a", 1: "b"}), + make_graph([], {0: "a", 1: "b"}), + ] + kernel = LdWloa(cutoff=1, labelings=[TYPE]) + + matrix = kernel.fit_transform(graphs) + + np.testing.assert_allclose(matrix, matrix.T) + np.testing.assert_allclose(np.diag(matrix), 1.0) + assert matrix[0, 1] == 1.0 + assert matrix[0, 2] < 1.0 + assert np.min(np.linalg.eigvalsh(matrix)) > -1e-10 + np.testing.assert_allclose(kernel.transform(graphs), matrix) + assert len(kernel.get_theta_parameters()) == 2 + + +def test_ld_wloa_multiple_labelings_share_depth_parameters(): + graphs = [ + make_graph([], {0: "same"}, {0: "family-a"}), + make_graph([], {0: "same"}, {0: "family-b"}), + ] + kernel = LdWloa(cutoff=2, labelings=[TYPE, FAMILY]) + + matrix = kernel.fit_transform(graphs) + + assert matrix[0, 1] < 1.0 + assert len(kernel.get_theta_parameters()) == 4 # 3 depth + 1 labeling-mixture + + +def test_ld_wloa_rejects_duplicate_labeling_keys(): + with pytest.raises(ValueError, match="labeling keys must be unique"): + LdWloa(1, [TYPE, NodeLabeling("type", lambda node: str(node))]) + + +def test_ld_wloa_reuses_level_kernels(monkeypatch): + graphs = [ + make_graph([(0, 1, "derives")], {0: "a", 1: "b"}), + make_graph([], {0: "a", 1: "b"}), + make_graph([(1, 0, "derives")], {0: "a", 1: "c"}), + ] + kernel = LdWloa(cutoff=1, labelings=[TYPE]) + calls = [] + kernel_matrix_for_level = kernel._kernel_matrix_for_level + + def track(left, right): + calls.append((len(left), len(right))) + return kernel_matrix_for_level(left, right) + + monkeypatch.setattr(kernel, "_kernel_matrix_for_level", track) + kernel.fit_transform(graphs[:2], theta=np.array([0.2, 0.8])) + first_call_count = len(calls) + + kernel.fit_transform(graphs[:2], theta=np.array([0.8, 0.2])) + assert len(calls) == first_call_count + + incremental = kernel.fit_transform(graphs, theta=np.array([0.8, 0.2])) + assert calls[first_call_count:] == [(1, 2), (1, 1)] * 2 + + expected = LdWloa(cutoff=1, labelings=[TYPE]).fit_transform( + graphs, + theta=np.array([0.8, 0.2]), + ) + np.testing.assert_allclose(incremental, expected) + + +def test_ld_wloa_reuses_transform_kernels(monkeypatch): + train = [make_graph([], {0: "a"}), make_graph([], {0: "b"})] + test = [make_graph([], {0: "c"})] + kernel = LdWloa(cutoff=1, labelings=[TYPE]) + kernel.fit_transform(train) + calls = [] + kernel_matrix_for_level = kernel._kernel_matrix_for_level + + def track(left, right): + calls.append((len(left), len(right))) + return kernel_matrix_for_level(left, right) + + monkeypatch.setattr(kernel, "_kernel_matrix_for_level", track) + kernel.transform(test, theta=np.array([0.2, 0.8])) + first_call_count = len(calls) + kernel.transform(test, theta=np.array([0.8, 0.2])) + + assert first_call_count == 2 + assert len(calls) == first_call_count + + +def test_ld_wloa_copies_share_graph_feature_cache(): + graph = make_graph([(0, 1, "derives")], {0: "a", 1: "b"}) + kernel = LdWloa(cutoff=1, labelings=[TYPE]) + copied = copy.deepcopy(kernel) + + original_features = kernel._build_features(graph) + copied_features = copied._build_features(graph) + + assert copied_features is original_features + + +def test_ld_wloa_kernel_copies_share_graph_features(): + first_graph = make_graph([(0, 1, "derives")], {0: "shared-a", 1: "shared-b"}) + same_graph = make_graph([(0, 1, "derives")], {0: "shared-a", 1: "shared-b"}) + template = LdWloa(cutoff=2, labelings=[TYPE]) + first = copy.deepcopy(template) + second = copy.deepcopy(template) + + assert first._build_features(first_graph) is second._build_features(same_graph) + + +def test_ld_wloa_kernel_copies_share_compact_labels_after_eviction(): + graph = make_graph([(0, 1, "derives")], {0: "a", 1: "b"}) + template = LdWloa(cutoff=2, labelings=[TYPE], cache_size=0) + first = copy.deepcopy(template) + second = copy.deepcopy(template) + + first_features = first._build_features(graph) + second_features = second._build_features(graph) + + assert first_features is not second_features + assert first_features.counts == second_features.counts + assert all( + isinstance(label, int) + for level in first_features.counts[0] + for label in level + ) + + +def test_ld_wloa_kernel_copies_reuse_and_extend_training_matrices(monkeypatch): + graphs = [ + make_graph([], {0: "matrix-a"}), + make_graph([], {0: "matrix-b"}), + make_graph([], {0: "matrix-c"}), + ] + template = LdWloa(cutoff=3, labelings=[TYPE], cache_size=1) + first = copy.deepcopy(template) + first.fit_transform(graphs[:2]) + + # Candidate searches may overflow the LRU without evicting training data. + first._build_features(make_graph([], {0: "candidate-a"})) + first._build_features(make_graph([], {0: "candidate-b"})) + assert len(template._cache.candidate_features) == 1 + + second = copy.deepcopy(template) + calls = [] + matrix_for_level = second._kernel_matrix_for_level + + def track(left, right): + calls.append((len(left), len(right))) + return matrix_for_level(left, right) + + monkeypatch.setattr(second, "_kernel_matrix_for_level", track) + second.fit_transform(graphs[:2]) + assert calls == [] + + incremental = second.fit_transform(graphs) + assert calls == [(1, 2), (1, 1)] * 4 + np.testing.assert_allclose(incremental, incremental.T) + np.testing.assert_allclose(np.diag(incremental), 1.0) + + +def test_ld_wloa_cache_size_validation(): + with pytest.raises(ValueError, match="cache_size must be non-negative"): + LdWloa(cutoff=1, labelings=[TYPE], cache_size=-1) + + unlimited = LdWloa(cutoff=0, labelings=[TYPE], cache_size=None) + unlimited._build_features(make_graph([], {0: "unlimited-a"})) + unlimited._build_features(make_graph([], {0: "unlimited-b"})) + assert unlimited.cache_size is None + assert len(unlimited._cache.candidate_features) == 2 + + +def test_sizing_kernel(): + features = [SizingFeature("size", "numeric"), SizingFeature("material", "categorical")] + graphs = [ + make_graph([], {}, sizing=[0.0, 0]), + make_graph([], {}, sizing=[1.0, 1]), + ] + kernel = SizingKernel(features) + + matrix = kernel.fit_transform(graphs, theta=np.ones(2)) + + np.testing.assert_allclose(np.diag(matrix), 1.0) + np.testing.assert_allclose(matrix[0, 1], np.exp(-2.0)) + np.testing.assert_allclose(kernel.transform(graphs, theta=np.ones(2)), matrix) + + +def test_sizing_kernel_rejects_invalid_parameters(): + graph = make_graph([], {}, sizing=[0.5]) + with pytest.raises(ValueError, match="power must be"): + SizingKernel([SizingFeature("size", "numeric")], power=2.5) + kernel = SizingKernel([SizingFeature("size", "numeric")]) + with pytest.raises(ValueError, match="theta size must be 1, got 2"): + kernel.fit_transform([graph], theta=np.array([1.0, 2.0])) + + +def test_sizing_kernel_reuses_and_extends_distances(): + features = [SizingFeature("size", "numeric")] + graphs = [ + make_graph([], {}, sizing=[0.0]), + make_graph([], {}, sizing=[0.5]), + make_graph([], {}, sizing=[1.0]), + ] + kernel = SizingKernel(features, power=1.9) + kernel.fit_transform(graphs[:2], theta=np.array([0.5])) + cached_distances = kernel._train_powered_distances + + kernel.fit_transform(graphs[:2], theta=np.array([2.0])) + assert kernel._train_powered_distances is cached_distances + + incremental = kernel.fit_transform(graphs, theta=np.array([2.0])) + expected = SizingKernel(features, power=1.9).fit_transform( + graphs, + theta=np.array([2.0]), + ) + np.testing.assert_allclose(incremental, expected) + + +def test_edge_multiplicity_kernel(): + graphs = [ + make_graph([(0, 1, "connects")], {0: "provided", 1: "needed"}), + make_graph( + [(0, 1, "connects"), (0, 1, "connects")], + {0: "provided", 1: "needed"}, + ), + ] + kernel = EdgeMultiplicityKernel("connects", TYPE) + + matrix = kernel.fit_transform(graphs, theta=np.ones(1)) + + np.testing.assert_allclose(np.diag(matrix), 1.0) + np.testing.assert_allclose(matrix[0, 1], np.exp(-abs(np.log(2.0) - np.log(3.0)))) + np.testing.assert_allclose(kernel.transform(graphs, theta=np.ones(1)), matrix) + + +def test_edge_multiplicity_kernel_rejects_extra_parameter(): + graph = make_graph([], {0: "provided"}) + with pytest.raises(ValueError, match="theta size must be 1, got 2"): + EdgeMultiplicityKernel("connects", TYPE).fit_transform( + [graph], theta=np.array([1.0, 2.0]) + ) + + +def test_edge_multiplicity_kernel_reuses_and_extends_distances(): + graphs = [ + make_graph([(0, 1, "connects")], {0: "provided", 1: "needed"}), + make_graph( + [(0, 1, "connects"), (0, 1, "connects")], + {0: "provided", 1: "needed"}, + ), + make_graph([], {0: "provided", 1: "needed"}), + ] + kernel = EdgeMultiplicityKernel("connects", TYPE) + kernel.fit_transform(graphs[:2], theta=np.array([0.5])) + cached_distances = kernel._train_distances + + kernel.fit_transform(graphs[:2], theta=np.array([2.0])) + assert kernel._train_distances is cached_distances + + incremental = kernel.fit_transform(graphs, theta=np.array([2.0])) + expected = EdgeMultiplicityKernel("connects", TYPE).fit_transform( + graphs, + theta=np.array([2.0]), + ) + np.testing.assert_allclose(incremental, expected) + + +def test_composite_graph_kernel(): + features = [SizingFeature("size", "numeric")] + graphs = [ + make_graph([], {0: "a"}, sizing=[0.0]), + make_graph([], {0: "b"}, sizing=[1.0]), + ] + + structure_expected = LdWloa(0, [TYPE]).fit_transform(graphs, theta=np.array([0.5])) + sizing_expected = SizingKernel(features).fit_transform(graphs, theta=np.array([1.0])) + kernel = CompositeGraphKernel( + [LdWloa(0, [TYPE]), SizingKernel(features)], + composition="additive", + ) + + matrix = kernel.fit_transform(graphs, theta=np.array([1.0, 0.5, 1.0])) + + np.testing.assert_allclose(matrix, 0.5 * structure_expected + 0.5 * sizing_expected) + np.testing.assert_allclose(kernel.transform(graphs, theta=np.array([1.0, 0.5, 1.0])), matrix) + + +def test_composite_graph_kernel_rejects_invalid_configuration(): + graph = make_graph([], {0: "a"}) + with pytest.raises(ValueError, match="unsupported graph kernel composition"): + CompositeGraphKernel([LdWloa(0, [TYPE])], composition="unknown") + + kernel = CompositeGraphKernel([LdWloa(0, [TYPE])]) + with pytest.raises(ValueError, match="theta size must be 1, got 2"): + kernel.fit_transform([graph], theta=np.array([0.5, 0.5])) + + +def test_composite_graph_kernel_clamps_negative_weights(): + features = [SizingFeature("size", "numeric")] + graphs = [ + make_graph([], {0: "a"}, sizing=[0.0]), + make_graph([], {0: "b"}, sizing=[1.0]), + ] + structure_expected = LdWloa(0, [TYPE]).fit_transform( + graphs, + theta=np.array([0.5]), + ) + kernel = CompositeGraphKernel( + [LdWloa(0, [TYPE]), SizingKernel(features)], + composition="additive", + ) + + matrix = kernel.fit_transform(graphs, theta=np.array([-1.0, 0.5, 1.0])) + + np.testing.assert_allclose(matrix, structure_expected) + assert np.min(np.linalg.eigvalsh(matrix)) > -1e-10 diff --git a/sb_arch_opt/tests/algo/test_arch_sbo_graph_representation.py b/sb_arch_opt/tests/algo/test_arch_sbo_graph_representation.py new file mode 100644 index 0000000..8a3119a --- /dev/null +++ b/sb_arch_opt/tests/algo/test_arch_sbo_graph_representation.py @@ -0,0 +1,70 @@ +import numpy as np +import pytest + +from sb_arch_opt.algo.arch_sbo.graph import HAS_NETWORKX, GraphDecoder, GraphRepresentation, NodeLabeling, SizingFeature + +if HAS_NETWORKX: + import networkx as nx + + +pytestmark = pytest.mark.skipif(not HAS_NETWORKX, reason="NetworkX dependency not installed") + + +def node_kind(node): + return node["kind"] + + +def test_node_labeling(): + labeling = NodeLabeling("example.kind", node_kind) + + assert labeling({"kind": "component"}) == "component" + assert labeling == NodeLabeling("example.kind", lambda _node: "other") + + +def test_graph_representation(): + labeling = NodeLabeling("example.kind", node_kind) + graph = nx.MultiDiGraph() + graph.add_nodes_from([0, 1]) + graph.add_edge(0, 1, kind="connects") + graph.add_edge(0, 1, kind="connects") + labels = {labeling: {0: "source", 1: "target"}} + + representation = GraphRepresentation( + graph=graph, + node_labels=labels, + sizing_values=np.array([0.5, 1.0]), + ) + + assert representation.graph is graph + assert representation.node_labels is labels + assert representation.graph.number_of_edges(0, 1) == 2 + np.testing.assert_array_equal(representation.sizing_values, [0.5, 1.0]) + + +def test_graph_representation_without_sizing(): + representation = GraphRepresentation(nx.Graph(), {}) + + assert representation.sizing_values.shape == (0,) + + +def test_graph_decoder_protocol(): + labeling = NodeLabeling("example.kind", node_kind) + + class Decoder: + labelings = (labeling,) + sizing_features = (SizingFeature("size", "numeric"),) + + def decode(self, x): + return [ + GraphRepresentation( + graph=nx.Graph(), + node_labels={labeling: {}}, + sizing_values=np.asarray(row[:1], dtype=float), + ) + for row in x + ] + + decoder = Decoder() + + assert isinstance(decoder, GraphDecoder) + assert len(decoder.decode(np.array([[0.25], [0.75]]))) == 2 diff --git a/sb_arch_opt/tests/algo/test_arch_sbo_graph_surrogate.py b/sb_arch_opt/tests/algo/test_arch_sbo_graph_surrogate.py new file mode 100644 index 0000000..627f85e --- /dev/null +++ b/sb_arch_opt/tests/algo/test_arch_sbo_graph_surrogate.py @@ -0,0 +1,141 @@ +import numpy as np +import pytest + +from sb_arch_opt.algo.arch_sbo.graph import ( + HAS_NETWORKX, + CompositeGraphKernel, + GraphRepresentation, + LdWloa, + NodeLabeling, + SizingFeature, + SizingKernel, +) +from sb_arch_opt.algo.arch_sbo.models import HAS_SMT, MultiSurrogateModel + +if HAS_NETWORKX: + import networkx as nx + +if HAS_SMT: + from sb_arch_opt.algo.arch_sbo.graph import GraphKriging + try: + from smt.design_space import DesignSpace, IntegerVariable + except ImportError: + from smt.utils.design_space import DesignSpace, IntegerVariable + + +pytestmark = pytest.mark.skipif( + not HAS_NETWORKX or not HAS_SMT, + reason="GraphKriging dependencies not installed", +) + + +LABELING = NodeLabeling("kind", lambda node: node) + + +class ExampleDecoder: + labelings = (LABELING,) + sizing_features = (SizingFeature("size", "numeric"),) + + def decode(self, x): + representations = [] + for row in x: + graph = nx.Graph() + graph.add_node(0) + labels = {0: "root"} + if row[0] >= 0.5: + graph.add_edge(0, 1) + labels[1] = "leaf" + representations.append( + GraphRepresentation( + graph, + {LABELING: labels}, + np.asarray(row[:1]), + ) + ) + return representations + + +def test_graph_kriging_defaults_to_cobyla(): + model = GraphKriging(ExampleDecoder(), LdWloa(0, [LABELING])) + + assert model.options["hyper_opt"] == "Cobyla" + + +def test_graph_kriging_train_and_predict(): + decoder = ExampleDecoder() + kernel = CompositeGraphKernel( + [ + LdWloa(1, [LABELING]), + SizingKernel(decoder.sizing_features), + ], + composition="additive_interaction", + ) + model = GraphKriging( + decoder, + kernel, + hyper_opt="NoOp", + poly="constant", + corr="squar_exp", + print_global=False, + print_training=False, + print_prediction=False, + print_problem=False, + print_solver=False, + design_space=DesignSpace([IntegerVariable(0, 3)]), + ) + x = np.array([[0.0], [1.0], [2.0], [3.0]]) + y = np.array([[0.0], [0.5], [0.5], [0.0]]) + + model.set_training_values(x, y) + model.train() + prediction = model.predict_values(np.array([[0.0], [3.0]])) + variance = model.predict_variances(np.array([[0.0], [3.0]])) + + assert prediction.shape == (2, 1) + assert variance.shape == (2, 1) + assert np.all(np.isfinite(prediction)) + assert np.all(np.isfinite(variance)) + assert len(model.optimal_theta) == len(kernel.get_theta_parameters()) + np.testing.assert_array_equal( + model.optimal_theta, + [parameter.initial for parameter in kernel.get_theta_parameters()], + ) + + warm_start = np.linspace(0.2, 0.6, len(kernel.get_theta_parameters())) + model.options["theta0"] = warm_start + model.train() + np.testing.assert_array_equal(model.optimal_theta, warm_start) + + +def test_multi_graph_kriging_reuses_wl_matrices_across_retraining(): + decoder = ExampleDecoder() + model = GraphKriging( + decoder, + CompositeGraphKernel( + [LdWloa(1, [LABELING]), SizingKernel(decoder.sizing_features)], + composition="additive", + ), + hyper_opt="NoOp", + poly="constant", + print_global=False, + print_training=False, + print_prediction=False, + print_problem=False, + print_solver=False, + design_space=DesignSpace([IntegerVariable(0, 3)]), + ) + multi = MultiSurrogateModel(model) + x = np.arange(4.0).reshape(-1, 1) + y = np.column_stack((x[:, 0], x[:, 0] ** 2)) + + multi.set_training_values(x[:3], y[:3]) + multi.train() + old_first = multi._models[0].kernel.branches[0]._train_level_kernels + assert multi._models[1].kernel.branches[0]._train_level_kernels is old_first + + multi.set_training_values(x, y) + multi.train() + new_first = multi._models[0].kernel.branches[0]._train_level_kernels + assert multi._models[1].kernel.branches[0]._train_level_kernels is new_first + for old_level, new_level in zip(old_first[0], new_first[0]): + np.testing.assert_array_equal(new_level[:3, :3], old_level) diff --git a/sb_arch_opt/tests/test_problem.py b/sb_arch_opt/tests/test_problem.py index 9aaa3b7..92df07d 100644 --- a/sb_arch_opt/tests/test_problem.py +++ b/sb_arch_opt/tests/test_problem.py @@ -443,6 +443,17 @@ def test_cached_pareto_front_mixin(problem: ArchOptTestProblemBase, discrete_pro assert os.path.exists(discrete_problem._pf_cache_path()) +def test_load_pareto_front(tmp_path): + from sb_arch_opt.pareto_front import load_pareto_front + + pf = np.array([[0., 1.], [.5, .5], [1., 0.]]) + pf_path = tmp_path / 'pf.npy' + np.save(pf_path, pf) + + assert np.array_equal(load_pareto_front(pf_path), pf) + assert np.array_equal(load_pareto_front(tmp_path), pf) + + def test_failing_evaluations(failing_problem: ArchOptTestProblemBase): out = failing_problem.evaluate(np.random.random((4, 5)), return_as_dictionary=True) is_failed = failing_problem.get_failed_points(out) diff --git a/setup.py b/setup.py index 63caeb8..4e654e7 100644 --- a/setup.py +++ b/setup.py @@ -66,6 +66,7 @@ def _get_readme(): 'smt~=2.2,!=2.4,!=2.10.0', 'numba', 'scikit-learn', + 'networkx', ], # 'ota': [ # pip install -r requirements-ota.txt # 'open_turb_arch @ git+https://github.com/jbussemaker/OpenTurbofanArchitecting@pymoo_optional#egg=open_turb_arch',