diff --git a/.pylintrc b/.pylintrc index b26aeee4f..23aacadaa 100644 --- a/.pylintrc +++ b/.pylintrc @@ -129,6 +129,7 @@ disable=R, wrong-import-order, xrange-builtin, zip-builtin-not-iterating, + invalid-name, [REPORTS] diff --git a/docs/api/linprog.rst b/docs/api/linprog.rst new file mode 100644 index 000000000..927233fee --- /dev/null +++ b/docs/api/linprog.rst @@ -0,0 +1,12 @@ +Linear programming +================== + +.. currentmodule:: optax.linprog + +.. autosummary:: + rhpdhg + + +Restarted Halpern primal-dual hybrid gradient method +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +.. autofunction:: rhpdhg diff --git a/docs/gallery.rst b/docs/gallery.rst index 5f7a3b134..5d96431b2 100644 --- a/docs/gallery.rst +++ b/docs/gallery.rst @@ -209,7 +209,7 @@ .. only:: html .. image:: /images/examples/linear_assignment_problem.png - :alt: + :alt: Linear assignment problem. :doc:`_collections/examples/linear_assignment_problem` @@ -219,6 +219,23 @@ +.. raw:: html + +
+ +.. only:: html + + .. image:: /images/examples/linear_programming.png + :alt: Linear programming. + + :doc:`_collections/examples/linear_programming` + +.. raw:: html + +
Linear programming.
+
+ + .. raw:: html diff --git a/docs/images/examples/linear_programming.png b/docs/images/examples/linear_programming.png new file mode 100644 index 000000000..95fdd200c Binary files /dev/null and b/docs/images/examples/linear_programming.png differ diff --git a/docs/index.rst b/docs/index.rst index 694dadc97..6bacaeef8 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -54,6 +54,7 @@ for instructions on installing JAX. :caption: 📖 Reference :maxdepth: 2 + api/linprog api/assignment api/optimizers api/transformations diff --git a/examples/linear_programming.ipynb b/examples/linear_programming.ipynb new file mode 100644 index 000000000..980e33b02 --- /dev/null +++ b/examples/linear_programming.ipynb @@ -0,0 +1,229 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "205fbe7e-4a73-4ee0-b785-311b870357cf", + "metadata": {}, + "source": [ + "# Linear programming\n", + "\n", + "[Linear programming](https://en.wikipedia.org/wiki/Linear_programming) is one of the most important problems in optimization.\n", + "\n", + "A linear program is an optimization problem of the following form:\n", + "\n", + "$$\n", + "\\begin{align*}\n", + " \\text{minimize} \\quad & c \\cdot x \\\\\n", + " \\text{subject to} \\quad\n", + " & A x = b \\\\\n", + " & G x \\leq h\n", + "\\end{align*}\n", + "$$\n", + "\n", + "where:\n", + "- $c \\in \\mathbb{R}^d$ is a cost vector\n", + "- $A \\in \\mathbb{R}^{n \\times d}$ is an equality constraint matrix\n", + "- $b \\in \\mathbb{R}^n$ is an equality constraint vector\n", + "- $G \\in \\mathbb{R}^{m \\times d}$ is an inequality constraint matrix\n", + "- $h \\in \\mathbb{R}^m$ is an inequality constraint vector\n", + "\n", + "A linear program solver returns a solution $x \\in \\mathbb{R}^d$ to this problem, if one exists.\n", + "\n", + "Optax has a solver based on the [restarted Halpern primal-dual hybrid gradient (RHPDHG) method](https://arxiv.org/abs/2407.16144), which is a [matrix-free](https://en.wikipedia.org/wiki/Matrix-free_methods) [primal-dual](https://en.wikipedia.org/wiki/Duality_(optimization)) algorithm." + ] + }, + { + "cell_type": "markdown", + "id": "05d9b905-46aa-477d-ae96-632e797faa9a", + "metadata": {}, + "source": [ + "## Example\n", + "\n", + "Consider the following problem:\n", + "\n", + "$$\n", + "\\begin{align*}\n", + "\\text{maximize} \\quad 2 x + y & \\\\\n", + "\\text{subject to} \\quad\n", + "3 x + y &\\leq 21 \\\\\n", + "x + y &\\leq 9 \\\\\n", + "x + 4 y &\\leq 24\n", + "\\end{align*}\n", + "$$\n", + "\n", + "Note that this is a maximization problem.\n", + "\n", + "First, let's put this problem into the matrix form we described in the introduction." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "e527600d-b231-4c8c-b02a-c73971042fda", + "metadata": {}, + "outputs": [], + "source": [ + "from jax import numpy as jnp\n", + "\n", + "# We are trying to maximize rather than minimize, so we use a minus sign here.\n", + "c = -jnp.array([2, 1])\n", + "\n", + "# Our problem has no equality constraints, so we use a zero-size A and zero-size b.\n", + "A = jnp.zeros([0, 2])\n", + "b = jnp.zeros(0)\n", + "\n", + "G = jnp.array([[3, 1], [1, 1], [1, 4]])\n", + "h = jnp.array([21, 9, 24])" + ] + }, + { + "cell_type": "markdown", + "id": "535c5d85-3fa3-4f8b-9b34-35542d81107b", + "metadata": {}, + "source": [ + "Next, let's import optax and solve it." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d49955a5-573e-4be9-b759-104c18d49bbe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[5.999964 2.999992] -14.99992\n" + ] + } + ], + "source": [ + "import optax\n", + "\n", + "x = optax.linprog.rhpdhg(c, A, b, G, h, 1_000_000)['primal']\n", + "print(x, c @ x)" + ] + }, + { + "cell_type": "markdown", + "id": "469ff08a-6a17-486f-a5f8-9c4ef79dd37c", + "metadata": {}, + "source": [ + "Up to a small numerical error, the solution is $(6, 3)$, with a profit of $15$.\n", + "\n", + "Finally, let's plot the solution:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9688e718-b05a-4869-8a6f-309b6ca1b4de", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABK8UlEQVR4nO3deVQUV74H8G8DssoiCAIKgWeUoCIYcX9vJGjihlExosYYcM2ojGbQ0ZAJilvQRIwmZnRMjsgYNeqMGrOMBonEjSiKuIyIxOA2AbcorRIQ6X5/dLpjI0s3vVRV9/dzTh/p6urqW4D15f7qVl2ZUqlUgoiIyMxshG4AERFZJwYQEREJggFERESCYAAREZEgGEBERCQIBhAREQmCAURERIJgABERkSDshG5AbQqFAj///DNcXV0hk8mEbg4REelJqVTi/v378Pf3h41N/f0c0QXQzz//jICAAKGbQUREBrp27RratGlT7+uiCyBXV1cAqoa7ubkJ3BordWUHcGwy4Po/wMB8gD1RItKDXC5HQECA5nheH9EFkLrs5ubmxgASynNxwPlEoOYnQHEZaBEudIuISIIaO43CQQj0tGaugN8g1ddXtwvbFiKyWAwgqlvgKNW/V3cAvGE6EZkAA4jq1joGsHUE7hcD984I3RoiskCiOwdEIqEuw13fpSrD8TxQnZRKJR4/foyamhqhm0JkNra2trCzszP4UhkGENUvcNRvAbQD6LyEo+FqefToEUpLS1FRUSF0U4jMztnZGX5+frC3t2/yNhhAVL/aZTj2gjQUCgVKSkpga2sLf39/2Nvb88JpsgpKpRKPHj3CrVu3UFJSgnbt2jV4sWlDGEBUP5bh6vXo0SMoFAoEBATA2dlZ6OYQmZWTkxOaNWuGK1eu4NGjR3B0dGzSdjgIgRrG0XANaupffkRSZ4zfff7voYZxNBwRmQgDiBrGi1KJyEQYQNQ4luEsRlRUFN58802hm9EkGzduhIeHh+Z5amoqIiIiBGuPVCQkJGD48OFCN6NODCBqHMtwFmPnzp1YvHix0M0wijlz5iA7O1vzXMwHWn0YO1hXr16NjRs36vUemUyG3bt3G60N9WEAUeNYhrMYnp6ejd6hWCqaN28OLy8voZshmOrqap3Wc3d31+o5igkDiHTDMlzDlErg8UNhHnr8PGqX4IKCgvDuu+9i4sSJcHV1RWBgINavX6/1nmvXriEuLg4eHh7w9PTEsGHDcPnyZc3rNTU1SEpKgoeHB7y8vDB37lzEx8dr9UaCgoKwatUqre1GREQgNTVV83zlypUICwuDi4sLAgICMH36dDx48KDefXmyp5CamorMzEx88cUXkMlkkMlkyMnJQXR0NBITE7Xed+vWLdjb22v1nmr78ssv0a1bNzg6OqJly5YYMWKE5rW7d+/i9ddfR4sWLeDs7IxBgwahuLhY87q6VLhv3z6EhoaiefPmGDhwIEpLSzXr5OTkoHv37nBxcYGHhwf69OmDK1euYOPGjVi4cCFOnz6t2Q9170Umk2Ht2rV4+eWX4eLigqVLl6KmpgaTJk1CcHAwnJycEBISgtWrV2vtS+2eYVRUFGbOnIm5c+fC09MTvr6+Wj+HoKAgAMCIESMgk8k0z02B1wGRbnhRasNqKoDtzYX57LgHgJ1Lk9+enp6OxYsX4+2338Y///lPTJs2DX379kVISAiqq6sxYMAA9OrVC4cOHYKdnR2WLFmCgQMH4syZM7C3t0d6ejo2btyIDRs2IDQ0FOnp6di1axeio6P1aoeNjQ0+/PBDBAcH46effsL06dMxd+5c/O1vf2v0vXPmzEFhYSHkcjkyMjIAqHp7kydPRmJiItLT0+Hg4AAA+Oyzz9C6det62/f1119jxIgR+Otf/4p//OMfePToEb755hvN6wkJCSguLsaePXvg5uaGefPmYfDgwTh//jyaNWsGAKioqMCKFSuwadMm2NjY4LXXXsOcOXOwefNmPH78GMOHD8eUKVOwdetWPHr0CMePH4dMJsPo0aNx7tw57N27F/v37weg6sGopaamYtmyZVi1ahXs7OygUCjQpk0b7NixA15eXjh69CimTp0KPz8/xMXF1fv9yszMRFJSEo4dO4bc3FwkJCSgT58+ePHFF5GXlwcfHx9kZGRg4MCBsLW1bfT731QMININL0q1WIMHD8b06dMBAPPmzcMHH3yAAwcOICQkBNu2bYNCocCnn36qudNDRkYGPDw8kJOTg5deegmrVq1CcnIyYmNjAQDr1q3Dvn379G5H7Z7ZkiVL8Mc//lGnAGrevDmcnJxQVVUFX19fzfLY2FgkJibiiy++0ByQN27ciISEhHrvXLF06VKMGTMGCxcu1CwLD1f9vquD58iRI+jduzcAYPPmzQgICMDu3bsxapSqUlBdXY1169ahbdu2AIDExEQsWrQIgGqytvLycsTExGheDw0N1doXOzs7rf1Qe/XVVzFhwgStZU+2Mzg4GLm5udi+fXuDAdS5c2csWLAAANCuXTusWbMG2dnZePHFF+Ht7Q0A8PDwqLMNxsQAIt3x3nD1s3VW9USE+mwDdO7cWfO1TCaDr68vbt68CQA4ffo0fvzxx6fOG1VWVuLSpUsoLy9HaWkpevTooXnNzs4OkZGRUOpZqt2/fz/S0tJw4cIFyOVyPH78GJWVlaioqGjy3SYcHR0xfvx4bNiwAXFxccjPz8e5c+ewZ8+eet9TUFCAKVOm1PlaYWEh7OzstPbXy8sLISEhKCws1CxzdnbWhAsA+Pn5ab6nnp6eSEhIwIABA/Diiy+if//+iIuLg5+fX6P7ExkZ+dSyjz/+GBs2bMDVq1fx66+/4tGjR40OYnjyZ167febEc0CkO46Gq59MpiqDCfEw8A8Bddno912RQaFQAAAePHiArl27oqCgQOtx8eJFvPrqqzp/ho2NzVOB9ORJ9MuXLyMmJgadO3fGv/71L5w8eRIff/wxANVtjwwxefJkZGVl4fr168jIyEB0dDSeeeaZetd3cnIy6POAur+nT+5/RkYGcnNz0bt3b2zbtg3t27fHDz/80Oh2XVy0S62ff/455syZg0mTJuHbb79FQUEBJkyY0Oj3rKGfuTkxgEh3HA1ndZ5//nkUFxfDx8cHzz77rNbD3d0d7u7u8PPzw7FjxzTvefz4MU6ePKm1HW9vb62T8HK5HCUlJZrnJ0+ehEKhQHp6Onr27In27dvj559/1qut9vb2dU6LERYWhsjISHzyySfYsmULJk6c2OB2OnfuXO8AhdDQUDx+/Fhrf+/cuYOioiJ06NBBr/Z26dIFycnJOHr0KDp16oQtW7Y0uB91UZcCp0+fji5duuDZZ5/FpUuX9GpHXZo1a2aWKUb0DqCDBw9i6NCh8Pf3f2qseHV1NebNm6cZyeLv74/XX39d718kEjGOhrMq48aNQ8uWLTFs2DAcOnQIJSUlyMnJwcyZM3H9+nUAwKxZs7Bs2TLs3r0bFy5cwPTp03Hv3j2t7URHR2PTpk04dOgQzp49i/j4eK2T288++yyqq6vx0Ucf4aeffsKmTZuwbt06vdoaFBSEM2fOoKioCLdv39bqYU2ePBnLli2DUqnUGtFWlwULFmDr1q1YsGABCgsLcfbsWSxfvhyA6nzJsGHDMGXKFBw+fBinT5/Ga6+9htatW2PYsGE6tbOkpATJycnIzc3FlStX8O2336K4uFhzHigoKAglJSUoKCjA7du3UVVVVe+22rVrhxMnTmDfvn24ePEiUlJSkJeXp1M7GhIUFITs7GyUlZXh7t27Bm+vPnoH0MOHDxEeHq7pHj+poqIC+fn5SElJQX5+Pnbu3ImioiK8/PLLRmksiQDLcFbF2dkZBw8eRGBgIGJjYxEaGopJkyahsrISbm5uAIDZs2dj/PjxiI+PR69eveDq6vrUQT45ORl9+/ZFTEwMhgwZguHDh2udIwkPD8fKlSuxfPlydOrUCZs3b0ZaWppebZ0yZQpCQkIQGRkJb29vHDlyRPPa2LFjYWdnh7FjxzZ65+aoqCjs2LEDe/bsQUREBKKjo3H8+HHN6xkZGejatStiYmLQq1cvKJVKfPPNN0+Vterj7OyMCxcuYOTIkWjfvj2mTp2KGTNm4I033gAAjBw5EgMHDsQLL7wAb29vbN26td5tvfHGG4iNjcXo0aPRo0cP3LlzRzOgxBDp6enIyspCQEAAunTpYvD26iNT6num8Mk3y2TYtWtXg1cf5+XloXv37rhy5QoCAwMb3aZcLoe7uzvKy8s1v+AkMgdjVYMROr4NhC8VujWCqKysRElJCYKDg5t8K3pLlpCQgHv37pnlanpdXL58GW3btkVeXh6ef/55oZtjERr6P6Drcdzk54DKy8shk8nqvRK3qqoKcrlc60EixzIcSUR1dTXKysrwzjvvoGfPngwfkTFpAFVWVmLevHkYO3ZsvSmYlpamOZnp7u6OgIAAUzaJjIFlOJKII0eOwM/PD3l5eXqfUyLTM9l1QNXV1YiLi4NSqcTatWvrXS85ORlJSUma53K5nCEkdrwolRqh780vTSUqKkrv65HIfEzSA1KHz5UrV5CVldVgDdDBwQFubm5aD5IAluGIyEBGDyB1+BQXF2P//v1Wfbdai8YyHBEZSO8S3IMHD/Djjz9qnqvHq3t6esLPzw+vvPIK8vPz8dVXX6GmpgZlZWUAVLefsLe3N17LSVgswxGRgfTuAZ04cQJdunTRjA1PSkpCly5dMH/+fPz3v//Fnj17cP36dURERMDPz0/zOHr0qNEbTwJjGY6IDKB3D6ixk3o84WdFOEUDERmA94KjpuO94SSn9oR0UqKe6E3N2FNXWyoxT1XOACLDsAwnKTt37sTixYuFboZRzJkzR+umoWI+0OrD2MG6evVqvYfF177Pp6lwPiAyDMtwkuLp6Sl0E4ymefPmaN5coFloRaC6ulqn+889OaOq2LAHRIZhGU5SapfggoKC8O6772LixIlwdXVFYGAg1q9fr/Wea9euIS4uDh4eHvD09MSwYcNw+fJlzes1NTVISkqCh4cHvLy8MHfuXMTHx2v1RoKCgrBq1Sqt7UZERCA1NVXzfOXKlZo76QcEBGD69Ol48KD+Sf6e7CmkpqYiMzMTX3zxBWQyGWQyGXJychAdHY3ExESt9926dQv29vb1TrkAAF9++SW6desGR0dHtGzZUuvmqnfv3sXrr7+OFi1awNnZGYMGDUJxcbHmdXWpcN++fQgNDUXz5s0xcOBArekocnJy0L17d7i4uMDDwwN9+vTBlStXsHHjRixcuBCnT5/W7Ie69yKTybB27Vq8/PLLcHFxwdKlS1FTU4NJkyYhODgYTk5OCAkJwerVq7X2pXbPMCoqCjNnzsTcuXPh6ekJX19frZ9DUFAQAGDEiBGQyWSa56bAACLDsQyn2u/HD4V5GPg9T09PR2RkJE6dOoXp06dj2rRpKCoqAqD6K3vAgAFwdXXFoUOHcOTIEc0BVT3pWXp6OjZu3IgNGzbg8OHD+OWXX7Br1y6922FjY4MPP/wQ//nPf5CZmYnvvvsOc+fO1em9c+bMQVxcnOZAX1pait69e2Py5MnYsmWL1pQGn332GVq3bo3o6Og6t/X1119jxIgRGDx4ME6dOoXs7Gx0795d83pCQgJOnDiBPXv2IDc3F0qlEoMHD9aa/qGiogIrVqzApk2bcPDgQVy9ehVz5swBoJovafjw4ejbty/OnDmD3NxcTJ06FTKZDKNHj8bs2bPRsWNHzX6MHj1as93U1FSMGDECZ8+excSJE6FQKNCmTRvs2LED58+fx/z58/H2229j+/aG/xjMzMyEi4sLjh07hvfeew+LFi1CVlYWAGimc8jIyEBpaalRpneoD0twZDiW4YCaCmC7QOWguAeqmVGbaPDgwZpb+M+bNw8ffPABDhw4gJCQEGzbtg0KhQKffvopZL/NvJqRkQEPDw/k5OTgpZdewqpVq5CcnIzY2FgAwLp167Bv3z6921G7Z7ZkyRL88Y9/xN/+9rdG39u8eXM4OTmhqqoKvr6+muWxsbFITEzEF198gbi4OACqHkpCQoJmf2pbunQpxowZg4ULF2qWhYerfqeLi4uxZ88ezURwALB582YEBARg9+7dGDVK9cdYdXU11q1bp5lyIjExEYsWLQKgut1YeXk5YmJiNK+r5wJS74udnZ3Wfqi9+uqrmDBhgtayJ9sZHByM3NxcbN++XbO/dencuTMWLFgAQDWn0Jo1a5CdnY0XX3wR3t7eAAAPD48622BM7AGR4ViGk7TOnTtrvpbJZPD19cXNmzcBAKdPn8aPP/4IV1dXzTkXT09PVFZW4tKlSygvL0dpaSl69Oih2YadnR0iIyP1bsf+/fvRr18/tG7dGq6urhg/fjzu3LmDioqKJu+bo6Mjxo8fjw0bNgAA8vPzce7cOSQkJNT7noKCAvTr16/O1woLC2FnZ6e1v15eXggJCUFhYaFmmbOzs9Z8R35+fprvqaenJxISEjBgwAAMHToUq1ev1irPNaSu7+vHH3+Mrl27wtvbG82bN8f69etx9erVBrfz5M+8dvvMiT0gMo7AUb/dFWEH0HkJUM9flxbL1lnVExHqsw1Q+0S2TCaDQqEAoLrzSdeuXbF58+an3qf+S1kXNjY2T10j+GTJ6vLly4iJicG0adOwdOlSeHp64vDhw5g0aRIePXoEZ+em7+PkyZMRERGB69evIyMjA9HR0XjmmWfqXd/JyanJn6VW1/f0yf3PyMjAzJkzsXfvXmzbtg3vvPMOsrKy0LNnzwa36+Ki3dP9/PPPMWfOHKSnp2smA3z//fe1pgzXtX3qn7k5MYDIOKy9DCeTGVQGE6vnn38e27Ztg4+PT703Cvbz88OxY8fwhz/8AYDqHMfJkye15t7x9vbW+itfLpejpKRE8/zkyZNQKBRIT0+HjY2qMNPYeYza7O3tUVNT89TysLAwREZG4pNPPsGWLVuwZs2aBrfTuXNnZGdnP1XqAlSlssePH+PYsWOaEtydO3dQVFSEDh066NVe9R1lkpOT0atXL2zZsgU9e/asdz/qoi4FPjkL6qVLl/RqR12aNWumcxsMwRIcGQfLcBZp3LhxaNmyJYYNG4ZDhw6hpKQEOTk5mDlzJq5fvw4AmDVrFpYtW4bdu3fjwoULmD59Ou7du6e1nejoaGzatAmHDh3C2bNnER8fD1tbW83rzz77LKqrq/HRRx/hp59+wqZNm/SevycoKAhnzpxBUVERbt++rdXDmjx5MpYtWwalUvnUdOG1LViwAFu3bsWCBQtQWFiIs2fPYvny5QBU50uGDRuGKVOm4PDhwzh9+jRee+01tG7dGsOGDdOpnSUlJUhOTkZubi6uXLmCb7/9FsXFxZrzQEFBQZp7bN6+fVtrAEVt7dq1w4kTJ7Bv3z5cvHgRKSkpRhk0EBQUhOzsbJSVleHu3bsGb68+DCAyHo6GszjOzs44ePAgAgMDERsbi9DQUEyaNAmVlZWaHtHs2bMxfvx4xMfHa8pAtQ/yycnJ6Nu3L2JiYjBkyBAMHz5c6xxJeHg4Vq5cieXLl6NTp07YvHkz0tLS9GrrlClTEBISgsjISHh7e+PIkSOa18aOHQs7OzuMHTu20SnUo6KisGPHDuzZswcRERGIjo7G8ePHNa9nZGSga9euiImJQa9evaBUKvHNN9/odE0OoPqeXrhwASNHjkT79u0xdepUzJgxA2+88QYAYOTIkRg4cCBeeOEFeHt7Y+vWrfVu64033kBsbCxGjx6NHj164M6dO1q9oaZKT09HVlYWAgICNPf9NAWZUmQ3b9N1LnESoer7wE4foKYSGFRg0WW4yspKlJSUIDg4uNEDmjVKSEjAvXv3zHI1vS4uX76Mtm3bIi8vj9NyG0lD/wd0PY6zB0TGwzIciUx1dTXKysrwzjvvoGfPngwfkWEAkXGxDEcicuTIEfj5+SEvL0/vc0pkehwFR8Zl7aPhCAD0vvmlqTQ2fQwJiz0gMi6W4YhIRwwgMj6W4YhIBwwgMr7aZTgiojowgMj4WIYjIh0wgMg0WIYjokYwgMg0WIYjokYwgMg0WIYTpdozokqJeqZRtSdnRKX61Z4RVUwYQGQ6LMOJzs6dO7F48WKhm2EUc+bM0ZpWW8wHWn0YO1hXr16t93VZMpnMLLdR4oWoZDq8KFV0PD09hW6C0agnyLNW1dXVOt0A1d3d3QytaRr2gMh0WIYTndoluKCgILz77ruYOHEiXF1dERgYiPXr12u959q1a4iLi4OHhwc8PT0xbNgwXL58WfN6TU0NkpKS4OHhAS8vL8ydOxfx8fFavZGgoCCsWrVKa7sRERFITU3VPF+5ciXCwsLg4uKCgIAATJ8+HQ8e1D/J35M9hdTUVGRmZuKLL76ATCaDTCZDTk4OoqOjkZiYqPW+W7duwd7eXqv3VNuXX36Jbt26wdHRES1bttS6u/fdu3fx+uuvo0WLFnB2dsagQYNQXFyseV1dKty3bx9CQ0PRvHlzDBw4UGs+pJycHHTv3h0uLi7w8PBAnz59cOXKFWzcuBELFy7E6dOnNfuh7r3IZDKsXbsWL7/8MlxcXLB06VLU1NRg0qRJCA4OhpOTE0JCQrB69WqtfandM4yKisLMmTMxd+5ceHp6wtfXV+vnEBQUBAAYMWIEZDKZ5rkpMIDItKylDKdUAo8fCvMw8Puanp6OyMhInDp1CtOnT8e0adNQVFQEQPVX9oABA+Dq6opDhw7hyJEjmgPqo0ePNO/fuHEjNmzYgMOHD+OXX37Brl279G6HjY0NPvzwQ/znP/9BZmYmvvvuO8ydO1en986ZMwdxcXGaA31paSl69+6NyZMnY8uWLVpz6nz22Wdo3bo1oqOj69zW119/jREjRmDw4ME4deoUsrOz0b17d83rCQkJOHHiBPbs2YPc3FwolUoMHjxYa/6hiooKrFixAps2bcLBgwdx9epVzJkzB4Bqwr7hw4ejb9++OHPmDHJzczF16lTIZDKMHj0as2fPRseOHTX7MXr0aM12U1NTMWLECJw9exYTJ06EQqFAmzZtsGPHDpw/fx7z58/H22+/3ehkfpmZmXBxccGxY8fw3nvvYdGiRcjKygIAzXxCGRkZKC0tNcr8QvVhCY5My1rKcDUVwHaBykFxDwyajXXw4MGaOWTmzZuHDz74AAcOHEBISAi2bdsGhUKBTz/9FLLfplnPyMiAh4cHcnJy8NJLL2HVqlVITk5GbGwsAGDdunXYt2+f3u2o3TNbsmQJ/vjHP+Jvf/tbo+9t3rw5nJycUFVVBV9fX83y2NhYJCYm4osvvkBcXBwAVQ8lISFBsz+1LV26FGPGjMHChQs1y8LDVb+3xcXF2LNnj2YmUgDYvHkzAgICsHv3bowapfqDq7q6GuvWrdPMeZSYmIhFixYBUE1VUF5ejpiYGM3r6sno1PtiZ2entR9qr7766lMztT7ZzuDgYOTm5mL79u2a/a1L586dsWDBAgCqSe3WrFmD7OxsvPjii5qp1j08POpsgzGxB0SmxTKc6HXu3FnztUwmg6+vL27evAkAOH36NH788Ue4urpqzrl4enqisrISly5dQnl5OUpLS9GjRw/NNuzs7BAZGal3O/bv349+/fqhdevWcHV1xfjx43Hnzh1UVFQ0ed8cHR0xfvx4bNiwAQCQn5+Pc+fOISEhod73FBQUoF+/fnW+VlhYCDs7O6399fLyQkhICAoLCzXLnJ2dtSbc8/Pz03xPPT09kZCQgAEDBmDo0KFYvXq1VnmuIXV9Xz/++GN07doV3t7eaN68OdavX4+rV682uJ0nf+a122dO7AGR6QWOAq7vUpXhOi8B6vnLU9JsnVU9EaE+2wC1T2TLZDIoFAoAwIMHD9C1a1ds3rz5qfep/1LWhY2NzVN3pX6yZHX58mXExMRg2rRpWLp0KTw9PXH48GFMmjQJjx49grNz0/dx8uTJiIiIwPXr15GRkYHo6Gg888wz9a7v5OTU5M9Sq+t7+uT+Z2RkYObMmdi7dy+2bduGd955B1lZWejZs2eD23Vx0e7pfv7555gzZw7S09M1s9G+//77OHbsmN7tU//MzYkBRKZnDWU4mcygMphYPf/889i2bRt8fHzqndnSz88Px44dwx/+8AcAqnMcJ0+e1Jr8zdvbW+uvfLlcjpKSEs3zkydPQqFQID09HTY2qsJMY+cxarO3t0dNTc1Ty8PCwhAZGYlPPvkEW7ZswZo1axrcTufOnZGdnf1UqQtQlcoeP36MY8eOaUpwd+7cQVFRETp06KBXe7t06YIuXbogOTkZvXr1wpYtW9CzZ89696Mu6lLgk9NwX7p0Sa921KVZs2Y6t8EQLMGR6bEMJ1njxo1Dy5YtMWzYMBw6dAglJSXIycnBzJkzcf36dQDArFmzsGzZMuzevRsXLlzA9OnTce/ePa3tREdHY9OmTTh06BDOnj2L+Ph42Nraal5/9tlnUV1djY8++gg//fQTNm3apPcEckFBQThz5gyKiopw+/ZtrR7W5MmTsWzZMiiVSq0RbXVZsGABtm7digULFqCwsBBnz57F8uXLAajOlwwbNgxTpkzB4cOHcfr0abz22mto3bo1hg0bplM7S0pKkJycjNzcXFy5cgXffvstiouLNeeBgoKCUFJSgoKCAty+fVtrAEVt7dq1w4kTJ7Bv3z5cvHgRKSkpRhk0EBQUhOzsbJSVleHu3bsGb68+DCAyD2sZDWdhnJ2dcfDgQQQGBiI2NhahoaGYNGkSKisrNT2i2bNnY/z48YiPj9eUgWof5JOTk9G3b1/ExMRgyJAhGD58uNY5kvDwcKxcuRLLly9Hp06dsHnzZqSlpenV1ilTpiAkJASRkZHw9vbGkSNHNK+NHTsWdnZ2GDt2LBwdHRvcTlRUFHbs2IE9e/YgIiIC0dHROH78uOb1jIwMdO3aFTExMejVqxeUSiW++eYbna7JAVTf0wsXLmDkyJFo3749pk6dihkzZuCNN94AAIwcORIDBw7ECy+8AG9vb2zdurXebb3xxhuIjY3F6NGj0aNHD9y5c0erN9RU6enpyMrKQkBAALp06WLw9uojU4psukC5XA53d3eUl5fX2+UnCaq+D+z0AWoqgUEFki/DVVZWoqSkBMHBwY0e0KxRQkIC7t27Z5ar6XVx+fJltG3bFnl5eVqlQWq6hv4P6HocZw+IzINlOBJAdXU1ysrK8M4776Bnz54MH5FhAJH5sAxHZnbkyBH4+fkhLy9P73NKZHocBUfmYw2j4QgA9L75palERUU9NfybxIM9IDIfluGI6AkMIDIvluGI6DcMIDIvzpRKRL9hAJF5sQxHRL9hAJH5sQxHRGAAkRBYhiMiNCGADh48iKFDh8Lf37/OecOVSiXmz58PPz8/ODk5oX///lqzBRKxDGd5KioqMHLkSLi5uUEmk+HevXt1zoJK9CS9rwN6+PAhwsPDMXHiRM0EVE9677338OGHHyIzMxPBwcFISUnBgAEDcP78ed6yhH5noVM0yBaadz+UC8RRwszMzMShQ4dw9OhRtGzZEu7u7sjLy9OaPkAmk2HXrl1a00OTddM7gAYNGoRBgwbV+ZpSqcSqVavwzjvvaO4M+49//AOtWrXC7t27MWbMGMNaS5aDF6VKwqNHj2Bvb9/oepcuXUJoaCg6deqkWabPfEFknYx6DqikpARlZWXo37+/Zpm7uzt69OiB3NzcOt9TVVUFuVyu9QAAFDU8ZwdJHMtwgoiKikJiYiISExPh7u6Oli1bIiUlRXO3gKCgICxevBivv/463NzcMHXqVADAv/71L3Ts2BEODg4ICgpCenq61jbT09Nx8OBByGQyREVFabalLsEFBQUBAEaMGAGZTKZ5TtbNqAFUVlYGAGjVqpXW8latWmleqy0tLQ3u7u6aR0BAgOqF038FClcas3kkNhwNJ4jMzEzY2dnh+PHjWL16NVauXIlPP/1U8/qKFSsQHh6OU6dOISUlBSdPnkRcXBzGjBmDs2fPIjU1FSkpKZrb7ezcuRNTpkxBr169UFpaip07dz71meo5ajIyMlBaWmqUOWtI+gS/F1xycjKSkpI0z+Vy+e8hdGq26t/QpDreSZLHMpwgAgIC8MEHH0AmkyEkJARnz57FBx98gClTpgBQTR43e/Zszfrjxo1Dv379kJKSAgBo3749zp8/j/fffx8JCQnw9PSEs7Mz7O3t4evrW+dnqstxHh4e9a5D1seoPSD1L9aNGze0lt+4caPeXzoHBwe4ublpPQAAHeap/j01mz0hS8UynCB69uwJ2RODPnr16oXi4mLNFMyRkZFa6xcWFqJPnz5ay/r06aP1HqKmMGoABQcHw9fXF9nZ2Zplcrkcx44dQ69evfTbWMdkoNN81dcMIcvFMpzoPDlyjciU9C7BPXjwAD/++KPmuXruck9PTwQGBuLNN9/EkiVL0K5dO80wbH9/f/2HXspkQFiq6utzi1iOs1Qsw5ndsWPHtJ7/8MMPaNeuHWxtbetcPzQ0VGt6a0A1z0779u3rfU9dmjVrxh4TadG7B3TixAl06dJFM094UlISunTpgvnzVb2VuXPn4k9/+hOmTp2Kbt264cGDB9i7d2/TrgFShxB7QpaLZTizu3r1KpKSklBUVIStW7fio48+wqxZs+pdf/bs2cjOzsbixYtx8eJFZGZmYs2aNZgzZ45enxsUFITs7GyUlZXh7t27hu4GWQC9A0g9wVPth3pEjEwmw6JFi1BWVobKykrs378f7du3b3oLGUKWj2U4s3r99dfx66+/onv37pgxYwZmzZqlGW5dl+effx7bt2/H559/jk6dOmH+/PlYtGgREhIS9Prc9PR0ZGVlISAgQPMHLFk3mVJk0wXK5XK4u7ujvLz89wEJgOrAdDZVVY4DgC7pLMdZiur7wE4foKYSGFQgiTJcZWUlSkpKEBwcLKk7fERFRSEiIoK3yCGDNfR/oN7jeC3SuRkpe0KWi2U4IqsknQACGEKWjGU4Iqsj+IWoeuPoOMvE0XBmkZOTI3QTiDSk1QNSY0/I8rAMR2R1pBlAAEPIErEMR2RVpBtAAEPI0nCmVCKrIu0AAhhCloRlOCKrIv0AAhhCloRlOCKrYRkBBDCELAXLcERWw3ICCGAIWQKW4UTr8uXLkMlkKCgoMHhbMpkMu3fvNng7JG2WFUAAQ8gSWFEZrkahRO6lO/ii4L/IvXQHNQrL2t/U1FREREQ8tby0tBSDBg0yf4NIVKR3IaoueLGqtFnJRal7z5Vi4ZfnUVpeqVnm5+6IBUM7YGAnPwFbZnqcFZUAS+wBqbEnJF1WUIbbe64U0z7L1wofACgrr8S0z/Kx91ypyT77n//8J8LCwuDk5AQvLy/0798fDx8+hEKhwKJFi9CmTRs4ODggIiICe/furXc7GzduhIeHh9ay3bt3a2Zb3bhxIxYuXIjTp09DJpNBJpNp3TX/yRLc2bNnER0drWnT1KlT8eDBA83rCQkJGD58OFasWAE/Pz94eXlhxowZqK6uNtr3hczPcgMIYAhJmQWX4WoUSiz88jzq2iv1soVfnjdJOa60tBRjx47FxIkTUVhYiJycHMTGxkKpVGL16tVIT0/HihUrcObMGQwYMAAvv/wyiouLm/RZo0ePxuzZs9GxY0eUlpaitLQUo0ePfmq9hw8fYsCAAWjRogXy8vKwY8cO7N+/H4mJiVrrHThwAJcuXcKBAweQmZmJjRs3agKNpMmyAwhgCEmVBY+GO17yy1M9nycpAZSWV+J4yS9G/+zS0lI8fvwYsbGxCAoKQlhYGKZPn47mzZtjxYoVmDdvHsaMGYOQkBAsX77coKkbnJyc0Lx5c9jZ2cHX1xe+vr5wcnJ6ar0tW7agsrIS//jHP9CpUydER0djzZo12LRpE27cuKFZr0WLFlizZg2ee+45xMTEYMiQIcjOzm7qt4JEwPIDCGAISZEFl+Fu3q8/fJqynj7Cw8PRr18/hIWFYdSoUfjkk09w9+5dyOVy/Pzzz+jTp4/W+n369EFhYaHR2/GkwsJChIeHw8XFRetzFQoFioqKNMs6duyoNQW4n58fbt68adK2kWlZRwABDCEpstAynI+rbhPY6bqePmxtbZGVlYV///vf6NChAz766COEhISgpKRE723Z2Nig9nyWpjwn06xZM63nMpkMCoXCZJ9Hpmc9AQQwhKTGQstw3YM94efuCFk9r8ugGg3XPdjTJJ8vk8nQp08fLFy4EKdOnYK9vT2ys7Ph7++PI0eOaK175MgRdOjQoc7teHt74/79+3j48KFmWe1rhOzt7VFTU9Nge0JDQ3H69Gmt7Rw5cgQ2NjYICQnRc+9ISqwrgACGkJRYaBnO1kaGBUNVB/XaIaR+vmBoB9ja1BdRTXfs2DG8++67OHHiBK5evYqdO3fi1q1bCA0NxV/+8hcsX74c27ZtQ1FREd566y0UFBRg1qxZdW6rR48ecHZ2xttvv41Lly5hy5YtTw0KCAoKQklJCQoKCnD79m1UVVU9tZ1x48bB0dER8fHxOHfuHA4cOIA//elPGD9+PFq1amX07wGJh/UFEMAQkhILLcMN7OSHta89D1937TKbr7sj1r72vMmuA3Jzc8PBgwcxePBgtG/fHu+88w7S09MxaNAgzJw5E0lJSZg9ezbCwsKwd+9e7NmzB+3atatzW56envjss8/wzTffICwsDFu3bkVqaqrWOiNHjsTAgQPxwgsvwNvbG1u3bn1qO87Ozti3bx9++eUXdOvWDa+88gr69euHNWvWmOJbQCIiU9Yu4gpMLpfD3d0d5eXlcHNzM+2HKZXA2VTVxaoA0CWdF6uKTfV9YKcPUFMJDCoQzUWplZWVKCkpQXBwMBwdm36upkahxPGSX3DzfiV8XFVlN1P0fIiMraH/A7oexy3zTgi64h0TxE9dhru+S1WGE0kAGYutjQy92noJ3QwiQVhnCe5JLMeJn4WW4YisHQMIYAiJnYWOhiOydgwgNYaQeFnoaDgia8cAehJDSLxYhiOyOAyg2hhC4iTSMpzIBpESmY0xfvcZQHVhCImPyMpw6tvCVFRUCNwSImGof/dr3yJJH9Y9DLshHKItPoGjfhuOvQPovET1MxKIra0tPDw8NDfDdHZ21syDQ2TJlEolKioqcPPmTXh4eGjdIFZfDKCGMITERWQzpapn9eQdmckaeXh4GDyzLQOoMQwh8RDZRakymQx+fn7w8fHhzJxkVZo1a2ZQz0eNAaQLhpB4iKgMp2Zra2uU/4xE1oaDEHTFgQniINLRcESkPwaQPhhCwhPZaDgiajoGkL4YQsLjRalEFoEB1BQMIWGxDEdkERhATcUQEg7LcEQWgQFkCIaQcFiGI5I8BpChGELCYBmOSPIYQMbAEDI/luGIJI8BZCwMIfNjGY5I0oweQDU1NUhJSUFwcDCcnJzQtm1bLF682DpuW88QMi+W4Ygkzei34lm+fDnWrl2LzMxMdOzYESdOnMCECRPg7u6OmTNnGvvjxIe37TEfkd0bjoj0Y/Qe0NGjRzFs2DAMGTIEQUFBeOWVV/DSSy/h+PHjxv4o8WJPyHxYhiOSLKMHUO/evZGdnY2LFy8CAE6fPo3Dhw9j0KBBda5fVVUFuVyu9bAIDCHzYBmOSLKMXoJ76623IJfL8dxzz8HW1hY1NTVYunQpxo0bV+f6aWlpWLhwobGbIQ4sx5key3BEkmX0HtD27duxefNmbNmyBfn5+cjMzMSKFSuQmZlZ5/rJyckoLy/XPK5du2bsJgmLPSHTYxmOSJKM3gP6y1/+grfeegtjxowBAISFheHKlStIS0tDfHz8U+s7ODjAwcHB2M0QF/aETEtkM6USkW6M3gOqqKiAjY32Zm1tbaFQKIz9UdLCnpDp8KJUIkkyegANHToUS5cuxddff43Lly9j165dWLlyJUaMGGHsj5IehpDpsAxHJDkypZGvEL1//z5SUlKwa9cu3Lx5E/7+/hg7dizmz58Pe3v7Rt8vl8vh7u6O8vJyuLm5GbNp4qFUAmdTVeU4AOiSznKcoarvAzt9gJpKYFABy3BEAtL1OG70ADKUVQQQwBAyhYOxqtFwHd8GwpcK3Roiq6XrcZz3ghMKy3HGxzIckaQwgITEEDIuXpRKJCkMIKExhIyHo+GIJIUBJAYMIeNhGY5IMhhAYsEQMg6W4YgkgwEkJgwhw7EMRyQZDCCxYQgZjmU4IklgAIkRQ8gwLMMRSQIDSKwYQk3HMhyRJDCAxIwh1HQswxGJHgNI7BhCTcMyHJHoMYCkgCGkP5bhiESPASQVDCH9sQxHJGoMIClhCOmHZTgiUWMASQ1DSHcswxGJGgNIihhCumMZjki0GEBSxRDSDctwRKLFAJIyhlDjWIYjEi0GkNQxhBrHMhyRKDGALAFDqGEswxGJEgPIUjCE6scyHJEoMYAsCUOofizDEYkOA8jSMITqxjIckegwgCwRQ+hpLMMRiQ4DyFIxhJ7GMhyRqDCALBlDSBvLcESiwgCydAyh37EMRyQqDCBrwBD6HctwRKLBALIWDCEVluGIRIMBZE0YQizDEYkIA8jaMIRYhiMSCQaQNbL2EGIZjkgUGEDWyppDiGU4IlFgAFkzaw4hluGIBMcAsnbWGkIswxEJjgFE1hlCLMMRCY4BRCrWGEIswxEJSrwBxAOC+VlbCLEMRyQo8QbQl+2Aw2OA4nVA+QUGkrlYUwixDEckKPEGUOUt4Oo2IG8a8HUosMuPgWQu1hRCLMMRCcZO6AbU64V/AxV5wM0c4PZRoPKGKpCublO97tgK8IkCWkWp/nULUR04yTjUIQQA5xapQggAQpMEa5JJ1C7DtQgXukVEVsMkPaD//ve/eO211+Dl5QUnJyeEhYXhxIkT+m3EuzcQlgL0ywZeuQf0PwiELQJaRasOGOpAYg/JdKyhJ8QyHJFgZEqlcY/Ud+/eRZcuXfDCCy9g2rRp8Pb2RnFxMdq2bYu2bds2+n65XA53d3eUl5fDzc2t7pVqqoA7x4EbOb/3kGoqtddhD8l4lErgbKqqJwQAXdItqyd0eStw9FXAtR0QU8TfEyID6XQchwkC6K233sKRI0dw6NChJr1f14ZrYSCZniWHUPV9YKeP6ndmUAHLcEQGEiyAOnTogAEDBuD69ev4/vvv0bp1a0yfPh1Tpkypc/2qqipUVVVpNTwgIEC/AKqNgWQalhxCB2OB67uAjm8D4UuFbg2RpAkWQI6OjgCApKQkjBo1Cnl5eZg1axbWrVuH+Pj4p9ZPTU3FwoULn1puUADVxkAyHksNIZbhiIxGsACyt7dHZGQkjh49qlk2c+ZM5OXlITc396n1TdIDagwDyTCWGEIswxEZja4BZPRh2H5+fujQoYPWstDQUPzrX/+qc30HBwc4ODgYuxkNs3UAfP5P9UBK3YHEYd/1s8Qh2urRcNd3qUbDMYCITM7oAdSnTx8UFRVpLbt48SKeeeYZY3+U8TCQ9GeJIRQ46rcA2gF0XmJdP08iARi9BJeXl4fevXtj4cKFiIuLw/HjxzFlyhSsX78e48aNa/T9TRoFZ2os2dXPkspxLMMRGYVg54AA4KuvvkJycjKKi4sRHByMpKSkekfB1SbKAKqNgaTNkkKIo+GIDCZoABlCEgFUGwPJckKIo+GIDMYAEpK1BpIlhBDLcEQGYwCJiTUFkiWEEMtwRAZhAImZpQeS1EOIZTgigzCApMQSA0nKIcQyHJFBGEBSZimBJOUQYhmOqMkYQJZEyoEk1RBiGY6oyRhAlkxqgSTFEGIZjqjJGEDWRAqBJMUQYhmOqEkYQNZMrIEktRBiGY6oSRhA9DsxBZKUQohlOKImYQBR/YQOJCmFEMtwRHpjAJHuhAgkqYQQy3BEemMAUdOZK5CkEEIswxHpjQFExmPKQJJCCLEMR6QXBhCZjrEDSewhxDIckV4YQGQ+xggkMYcQy3BEemEAkXCaGkiAeEOIZTginTGASDz0CSSfvkD5OaD4b6rlYgkhluGIdMYAIvHSJZDsXIDHD1Vfh84FIpYJe9BnGY5IZwwgkg6dAskV8B8s7M1VWYYj0gkDiKRLE0gHgJ8ygIeXn15HiHvZsQxHpBMGEFkGpRI4/Q5w/l3V8+bPAr9eF+ZedizDEelE1+O4nRnbRKQ/mQwIXwLY2KlGxz34EQhfDnj30i7ZVd4Arm5TPQDTBFIzV8BvkKoMd3U7A4jIQOwBkTQ0dJ2QOe9lxzIcUaNYgiPLo+vFqqYMJJbhiBrFACLL1JQ7Jhg7kDgajqhBDCCyXIbetsfQQGIZjqhBDCCybMa8d5y+geTZDcjqAyiqWIYjqgMDiCyfqW5gqksg2TioAshvIPD8B8JcGEskUgwgsg7muIu20FOYE0kMA4ish7mncqipUoXR9zGA8vHvvaEnMZDIivFCVLIeMhkQlqr6+twi4NRs1demCiFbB8B/ANB6qGo0XMgsoHWM+S+MJZI49oDIcpi7J1TfaDiW7MjKsQRH1smcIaTrRakMJLIyDCCyXuYMoaZclMpAIgvHACLrZq4QMsZFqQwksjAMICJzhJAp7g3HQCKJYwARAeYJIVPfG46BRBLDACJSM3UImfvecAwkEjkGENGTTBlCQk/RwEAikWEAEdVmyhAS0xQNIg8kpVIJGcPPojGAiOpiqhAS8xQNAgTSr9W/4tLdSyi+U4wff/kRxb+o/v318a/4bMRnaOvZtun7Q6InmgBatmwZkpOTMWvWLKxatarR9RlAZHKmCCGhy3D6MHIgHSg5gBM/n9CETPEvxbguv/7Uev6u/siJz0E7r3ZG3iESG1HcCy4vLw9///vf0blzZ1N+DJF+THHvuGaugN8gVRnu6nZxB5CtA+Dzf6oHUuoOJD3uZXeq7BTm7p/b4EcyfKguNqba8IMHDzBu3Dh88sknaNGihak+hqhp1CHUab7q+anZQOFKw7YZOEr179Udql6WVKgDKSwF6JcNvHIP6H8QCFsEtIoGbB1/D6S8acDXocAuP+DwGKB4Hf7Y7gX4uPjUu3mGD9XHZD2gGTNmYMiQIejfvz+WLFlS73pVVVWoqvr9VvZyudxUTSLSZuyeUOsY1cH6fjFw74y4e0EN0aOHJC/ZhjXlQNWvdZfnGD7UEJME0Oeff478/Hzk5eU1um5aWhoWLlxoimYQNc6YISSlMpw+6ggkeekBrMldgfSS7/HL48cAnu7x+Ts4IafvG2hnV6PqEYppYAaJgtFLcNeuXcOsWbOwefNmODo6Nrp+cnIyysvLNY9r164Zu0lEDTNmOU6qZTgdyavkePdoOoK3jMNfz2f/Fj6Ap4MbBrX6vZfjbwvk+P6KdkULnirZofyCRX5vSH9GHwW3e/dujBgxAra2tpplNTU1kMlksLGxQVVVldZrtXEUHAnGGKPjpDQaTg/yKjnWHF+D9Nx0/PLrL5rlnk6emNNrDhK7J8LWxhbBq4NhZ2OHnJj30a7qkiivQyLTE2wUXL9+/XD27FmtZRMmTMBzzz2HefPmNRg+RIIyRjnOwspwugSPq4OrZvmKF1egZ5ueT5zzMXyUHVkus1yIGhUVhYiICF4HRNJgaE9IzBel6kjf4NGLyO/UQIYTxXVARJJkaE9IwqPhTBo8aka+Domki7fiIaqPIT0hMd0bTgdmCR5dsYckeaK5FY++GEAkKk0NIYmU4UQVPPVhIEkOA4jIWJoSQiIfDSeJ4KkPA0n0GEBExtSUEBJhGU7SwVMfBpLoMICIjE3fEBJRGc4ig6c+DCTBMYCITEGfEBJBGc6qgqc+DCSzYwARmYo+ISRQGY7B0wAGkskxgIhMSdcQMnMZjsHTBAwko2MAEZmaLiFkpjIcg8eIGEgGYwARmYMuIWTCMlxDwTO712wkdk+EmwP/HxmEgaQ3BhCRuTQWQiYowzF4BMRAahQDiMicGgohI5bhGDwixEB6CgOIyNwaCiEDy3AMHglhIDGAiARRXwg1sQzH4LEAVhhIDCAiodQVQs9O0asMx+CxYFYQSAwgIiHVFUK3DjdahmPwWCELDCQGEJHQaodQ0DjUlGzFceVLuNnxU/i4OaJ7sCdsbWQMHvqdBQQSA4hIDJ4Iob3lvbDw5zdQWt1S83IrN3t0fvYidl6ez+ChukkwkBhARGKhVGLvN6sx7dCzUP1n+/3AoIQCgAy37N/Fr7a5DB5qnAQCiQFEJBI1CiX+d/l3KC3/FU+Gj5oSCihldzF10HXM7MHgIT2JMJB0PY7bmawFRAQAOF7yC0rLK1FX+ACADDaQKb3Qr/UQhg/pz9YB8Pk/1QMpdQdS5Q3g6jbVAxC8h6TGACIysZv3KxtfSY/1iBokoUBiABGZmI+ro1HXI9KLiAOJAURkYt2DPeHn7oiy8krUdcJVBgV8nWvQPdjT7G0jKySiQOIgBCIz2HuuFNM+ywcArRCS/fZs7TPvYmDUK/XPrEpkLkYY1MBRcEQis/dcKRZ+ef63AQkqfvZ3saDHLQws/y14Gprem0gITQgkOfzg7uHBACISkxqFEsdLfsHNe/fgc2o0ujudgu2gk6pb9DQ2vTeRGOgQSHKFN9zH3+IwbCIxsbWRoVdbLwBewMNA4PpJ4NoOoPMS1QrnFgGnZqu+ZgiRGOlyDqnilk6bsjFlO4moAYGjVP9e3aH6NywV6DRf9fWp2UDhSkGaRaQXdSCFpQD9soFX7gEv/FuntzKAiITSOgawdQTuFwP3zqhO4jKESOpsHQDv3jqtygAiEkozV8BvkOrrq9tV/zKEyIowgIiE9GQZTj0eiCFEVoIBRCSk2mU4NYYQWQEGEJGQ6irDqTGEyMIxgIiEVlcZTo0hRBaMAUQktPrKcGoMIbJQDCAioTVUhlNjCJEFYgARiUFDZTg1hhBZGAYQkRg0VoZTYwiRBWEAEYmBLmU4NYYQWQgGEJFY6FKGU2MIkQVgABGJha5lODWGEEkcA4hILPQpw6kxhEjCjB5AaWlp6NatG1xdXeHj44Phw4ejqKjI2B9DZJn0KcOpMYRIooweQN9//z1mzJiBH374AVlZWaiursZLL72Ehw8fGvujiCyPvmU4NYYQSZDRZ0Tdu3ev1vONGzfCx8cHJ0+exB/+8AdjfxyRZVGX4a7vUpXhWoTr/l51CAGcWZUkweTngMrLywEAnp6edb5eVVUFuVyu9SCyak0pw6mxJ0QSYtIAUigUePPNN9GnTx906tSpznXS0tLg7u6ueQQEBJiySUTi19QynBpDiCTCpAE0Y8YMnDt3Dp9//nm96yQnJ6O8vFzzuHbtmimbRCR+TRkNVxtDiCTAZAGUmJiIr776CgcOHECbNm3qXc/BwQFubm5aDyKrZ0gZTo0hRCJn9ABSKpVITEzErl278N133yE4ONjYH0Fk+Qwtw6kxhEjEjB5AM2bMwGeffYYtW7bA1dUVZWVlKCsrw6+//mrsjyKyXMYow6kxhEikjB5Aa9euRXl5OaKiouDn56d5bNu2zdgfRWTZjFGGU2MIkQgZ/TogpaH/UYhIpXYZTp9rgurC64RIZHgvOCKxMmYZTo09IRIRBhCRmBmzDKfGECKRYAARiZmxRsPVxhAiEWAAEYmZKcpwagwhEhgDiEjsTFGGU2MIkYAYQERiZ6oynBpDiATCACISO1OW4dQYQiQABhCRFJiyDKfGECIzYwARSYGpy3BqDCEyIwYQkRSYowynxhAiM2EAEUmFOcpwagwhMgMGEJFUmKsMp8YQIhNjABFJhTnLcGoMITIhBhCRlJizDKfGECITYQARSYm5y3BqDCEyAQYQkZQIUYZTYwiRkTGAiKRGiDKcGkOIjIgBRCQ1QpXh1BhCZCQMICKpEbIMp8YQIiNgABFJkZBlODWGEBmIAUQkRUKX4dQYQmQABhCRFImhDKfGEKImYgARSZUYynBqDCFqAgYQkVSJpQynxhAiPTGAiKRKTGU4NYYQ6YEBRCRlYirDqTGESEcMICIpE1sZTo0hRDpgABFJmRjLcGoMIWoEA4hI6sRYhlNjCFEDGEBEUifWMpwaQ4jqwQAikjoxl+HUGEJUBwYQkSUQcxlOjSFEtTCAiCyB2MtwagwhegIDiMgSSKEMp8YQot8wgIgshRTKcGoMIQIDiMhySKUMp8YQsnoMICJLIaUynBpDyKoxgIgsiZTKcGoMIavFACKyJFIrw6kxhKwSA4jIkkixDKfGELI6DCAiSyPFMpwaQ8iqMICILI1Uy3BqDCGrYbIA+vjjjxEUFARHR0f06NEDx48fN9VHEdGTpFyGU2MIWQWTBNC2bduQlJSEBQsWID8/H+Hh4RgwYABu3rxpio8jotqkXIZTYwhZPJME0MqVKzFlyhRMmDABHTp0wLp16+Ds7IwNGzaY4uOIqDapl+HUGEIWzegB9OjRI5w8eRL9+/f//UNsbNC/f3/k5uY+tX5VVRXkcrnWg4gMZAllODWGkMUyegDdvn0bNTU1aNWqldbyVq1aoays7Kn109LS4O7urnkEBAQYu0lE1skSynBqDCGLJPgouOTkZJSXl2se165dE7pJRJbBUspwagwhi2P0AGrZsiVsbW1x48YNreU3btyAr6/vU+s7ODjAzc1N60FERmBJZTg1hpBFMXoA2dvbo2vXrsjOztYsUygUyM7ORq9evYz9cUTUEEsqw6kxhCyGnSk2mpSUhPj4eERGRqJ79+5YtWoVHj58iAkTJpji44ioPrXLcC3ChW6RcahDCADOLVKFEACEJgnWJNKfSQJo9OjRuHXrFubPn4+ysjJERERg7969Tw1MICITU5fhru9SleEsJYAAhpAFkCmV4uqXy+VyuLu7o7y8nOeDiIzh8lbg6KuAazsgpkh14LYkSiVwNlUVQgDQJZ0hJDBdj+OCj4IjIhOztNFwtfGckGQxgIgsnSWOhquNISRJDCAia2CJo+FqYwhJDgOIyBpYehlOjSEkKQwgImtgDWU4NYaQZDCAiKyFNZTh1BhCksAAIrIW1lKGU2MIiR4DiMhaWFMZTo0hJGoMICJrYk1lODWGkGgxgIisibWV4dQYQqLEACKyJtZYhlNjCIkOA4jI2lhjGU6NISQqDCAia2OtZTg1hpBoMICIrI01l+HUGEKiwAAiskbWXIZTYwgJjgFEZI2svQynxhASFAOIyBqxDPc7hpBgGEBE1opluN8xhARhJ3QDalPPEC6XywVuCZGFc/0DUOUAVBQDV48CLcKEbpHwnkkCHlQB55cDR2YDDyqBkEShWyU56uO3spE/bGTKxtYws+vXryMgIEDoZhARkYGuXbuGNm3a1Pu66AJIoVDg559/hqurK2QymUk/Sy6XIyAgANeuXYObm5tJP8tcLG2fLG1/AO6TVHCfmk6pVOL+/fvw9/eHjU39Z3pEV4KzsbFpMDFNwc3NzWJ+wdQsbZ8sbX8A7pNUcJ+axt3dvdF1OAiBiIgEwQAiIiJBWHUAOTg4YMGCBXBwcBC6KUZjaftkafsDcJ+kgvtkeqIbhEBERNbBqntAREQkHAYQEREJggFERESCYAAREZEgrDaAPv74YwQFBcHR0RE9evTA8ePHhW5Sk6WlpaFbt25wdXWFj48Phg8fjqKiIqGbZVTLli2DTCbDm2++KXRTDPLf//4Xr732Gry8vODk5ISwsDCcOHFC6GY1WU1NDVJSUhAcHAwnJye0bdsWixcvbvQeYGJx8OBBDB06FP7+/pDJZNi9e7fW60qlEvPnz4efnx+cnJzQv39/FBcXC9NYHTW0T9XV1Zg3bx7CwsLg4uICf39/vP766/j5558FaatVBtC2bduQlJSEBQsWID8/H+Hh4RgwYABu3rwpdNOa5Pvvv8eMGTPwww8/ICsrC9XV1XjppZfw8OFDoZtmFHl5efj73/+Ozp07C90Ug9y9exd9+vRBs2bN8O9//xvnz59Heno6WrRoIXTTmmz58uVYu3Yt1qxZg8LCQixfvhzvvfcePvroI6GbppOHDx8iPDwcH3/8cZ2vv/fee/jwww+xbt06HDt2DC4uLhgwYAAqKyvN3FLdNbRPFRUVyM/PR0pKCvLz87Fz504UFRXh5ZdfFqClAJRWqHv37soZM2ZontfU1Cj9/f2VaWlpArbKeG7evKkEoPz++++FborB7t+/r2zXrp0yKytL2bdvX+WsWbOEblKTzZs3T/m///u/QjfDqIYMGaKcOHGi1rLY2FjluHHjBGpR0wFQ7tq1S/NcoVAofX19le+//75m2b1795QODg7KrVu3CtBC/dXep7ocP35cCUB55coV8zTqCVbXA3r06BFOnjyJ/v37a5bZ2Nigf//+yM3NFbBlxlNeXg4A8PT0FLglhpsxYwaGDBmi9fOSqj179iAyMhKjRo2Cj48PunTpgk8++UToZhmkd+/eyM7OxsWLFwEAp0+fxuHDhzFo0CCBW2a4kpISlJWVaf3uubu7o0ePHhZzrABUxwuZTAYPDw+zf7bobkZqardv30ZNTQ1atWqltbxVq1a4cOGCQK0yHoVCgTfffBN9+vRBp06dhG6OQT7//HPk5+cjLy9P6KYYxU8//YS1a9ciKSkJb7/9NvLy8jBz5kzY29sjPj5e6OY1yVtvvQW5XI7nnnsOtra2qKmpwdKlSzFu3Dihm2awsrIyAKjzWKF+TeoqKysxb948jB07VpAbrlpdAFm6GTNm4Ny5czh8+LDQTTHItWvXMGvWLGRlZcHR0VHo5hiFQqFAZGQk3n33XQBAly5dcO7cOaxbt06yAbR9+3Zs3rwZW7ZsQceOHVFQUIA333wT/v7+kt0na1FdXY24uDgolUqsXbtWkDZYXQmuZcuWsLW1xY0bN7SW37hxA76+vgK1yjgSExPx1Vdf4cCBA2af0sLYTp48iZs3b+L555+HnZ0d7Ozs8P333+PDDz+EnZ0dampqhG6i3vz8/NChQwetZaGhobh69apALTLcX/7yF7z11lsYM2YMwsLCMH78ePz5z39GWlqa0E0zmPp4YInHCnX4XLlyBVlZWYJNN2F1AWRvb4+uXbsiOztbs0yhUCA7Oxu9evUSsGVNp1QqkZiYiF27duG7775DcHCw0E0yWL9+/XD27FkUFBRoHpGRkRg3bhwKCgpga2srdBP11qdPn6eGx1+8eBHPPPOMQC0yXEVFxVMTjtna2kKhUAjUIuMJDg6Gr6+v1rFCLpfj2LFjkj1WAL+HT3FxMfbv3w8vLy/B2mKVJbikpCTEx8cjMjIS3bt3x6pVq/Dw4UNMmDBB6KY1yYwZM7BlyxZ88cUXcHV11dSn3d3d4eTkJHDrmsbV1fWpc1guLi7w8vKS7LmtP//5z+jduzfeffddxMXF4fjx41i/fj3Wr18vdNOabOjQoVi6dCkCAwPRsWNHnDp1CitXrsTEiROFbppOHjx4gB9//FHzvKSkBAUFBfD09ERgYCDefPNNLFmyBO3atUNwcDBSUlLg7++P4cOHC9foRjS0T35+fnjllVeQn5+Pr776CjU1NZrjhaenJ+zt7c3bWLOPuxOJjz76SBkYGKi0t7dXdu/eXfnDDz8I3aQmA1DnIyMjQ+imGZXUh2ErlUrll19+qezUqZPSwcFB+dxzzynXr18vdJMMIpfLlbNmzVIGBgYqHR0dlf/zP/+j/Otf/6qsqqoSumk6OXDgQJ3/d+Lj45VKpWoodkpKirJVq1ZKBwcHZb9+/ZRFRUXCNroRDe1TSUlJvceLAwcOmL2tnI6BiIgEYXXngIiISBwYQEREJAgGEBERCYIBREREgmAAERGRIBhAREQkCAYQEREJggFERESCYAAREZEgGEBERCQIBhAREQmCAURERIL4f2YpA7hPuhYMAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from matplotlib import pyplot as plt\n", + "\n", + "def plot_lp(c, A, b, G, h, x):\n", + " fig, ax = plt.subplots()\n", + "\n", + " for Ai, bi in zip(A, b):\n", + " ax.axline((0, bi / Ai[1]), (bi / Ai[0], 0), label=\"equality constraint\", c=\"purple\")\n", + "\n", + " for Gi, hi in zip(G, h):\n", + " ax.axline((0, hi / Gi[1]), (hi / Gi[0], 0), label=\"inequality constraint\", c=\"orange\")\n", + "\n", + " plt.arrow(x[0], x[1], -c[0].item(), -c[1].item(), width=0.1, ec=\"none\", fc=\"green\", label=\"profit\", zorder=2)\n", + "\n", + " ax.plot(*x, \"o\", label=\"solution\")\n", + "\n", + " ax.legend()\n", + " ax.set(aspect=\"equal\", xlim=(-1, 13), ylim=(-1, 13))\n", + " plt.show()\n", + "\n", + "plot_lp(c, A, b, G, h, x)" + ] + }, + { + "cell_type": "markdown", + "id": "5090f3a1-92b7-4835-b953-6bb571ec2828", + "metadata": {}, + "source": [ + "As you can see, the solution goes \"as far as it can go\" in the direction of increasing profit.\n", + "\n", + "If we change the profit vector, we end up in a different place, because we want to maximize in a different direction:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "05284d21-3d97-4a11-a84a-0e88830ce00b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[4.0000005 5. ] -19.0\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMJklEQVR4nO3deVgUV9o28LsBWWURRAEDgSRKUBFQ4sb7TgyauGFUjKgxRtxHZdRBR0MmKK5oIkajGX1NPpExatSJW5bRIJG4hCiKuIyIqLgl4BYFlQGxu78/OtWhka3ppaq779919SVdXV11CrBvzlOn6siUSqUSRERERmYldgOIiMgyMYCIiEgUDCAiIhIFA4iIiETBACIiIlEwgIiISBQMICIiEgUDiIiIRGEjdgOqUygU+PXXX+Hs7AyZTCZ2c4iISEtKpRIPHz6Ej48PrKxq7+dILoB+/fVX+Pr6it0MIiLS0Y0bN/Dcc8/V+rrkAsjZ2RmAquEuLi4it8ZCXdsBHBsPOL8A9MkB2BMlIi2UlpbC19dX/XleG8kFkFB2c3FxYQCJ5eUY4HwcIL8CKK4CzULEbhERmaD6TqNwEAI9q4kz4N1X9fX17eK2hYjMFgOIauY3VPXv9R0Ab5hORAbAAKKatYoCrO2BhwXAgzNit4aIzJDkzgGRRAhluJu7VGU4ngeqkVKpxNOnTyGXy8VuCpHRWFtbw8bGRudLZRhAVDu/ob8H0A6gwyKOhqvmyZMnKCoqQllZmdhNITI6R0dHeHt7w9bWttHbYABR7aqX4dgLUlMoFCgsLIS1tTV8fHxga2vLC6fJIiiVSjx58gR37txBYWEhWrduXefFpnVhAFHtWIar1ZMnT6BQKODr6wtHR0exm0NkVA4ODmjSpAmuXbuGJ0+ewN7evlHb4SAEqhtHw9WpsX/5EZk6ffzu838P1Y2j4YjIQBhAVDdelEpEBsIAovqxDGc2evTogRkzZojdjEbZuHEj3Nzc1M+TkpIQGhoqWntMRWxsLAYNGiR2M2rEAKL6sQxnNnbu3ImFCxeK3Qy9mDVrFjIyMtTPpfxBqw19B+uqVauwceNGrd4jk8mwe/duvbWhNgwgqh/LcGbD3d293jsUm4qmTZvCw8ND7GaIprKyskHrubq6avQcpYQBRA3DMlzdlErg6WNxHlr8PKqX4Pz9/bFkyRKMHTsWzs7O8PPzw/r16zXec+PGDcTExMDNzQ3u7u4YOHAgrl69qn5dLpcjPj4ebm5u8PDwwOzZszF69GiN3oi/vz9Wrlypsd3Q0FAkJSWpn69YsQLBwcFwcnKCr68vpkyZgkePHtV6LFV7CklJSUhLS8OePXsgk8kgk8mQmZmJyMhIxMXFabzvzp07sLW11eg9Vff111/jlVdegb29PZo3b47BgwerX7t//z7effddNGvWDI6Ojujbty8KCgrUrwulwv379yMoKAhNmzZFnz59UFRUpF4nMzMTnTt3hpOTE9zc3BAREYFr165h48aNmD9/Pk6fPq0+DqH3IpPJsHbtWrz55ptwcnLC4sWLIZfLMW7cOAQEBMDBwQGBgYFYtWqVxrFU7xn26NED06ZNw+zZs+Hu7g4vLy+Nn4O/vz8AYPDgwZDJZOrnhsDrgKhheFFq3eRlwPam4uw75hFg49Tot6ekpGDhwoV4//338a9//QuTJ0/Gq6++isDAQFRWVqJ3797o1q0bDh8+DBsbGyxatAh9+vTBmTNnYGtri5SUFGzcuBEbNmxAUFAQUlJSsGvXLkRGRmrVDisrK3zyyScICAjAlStXMGXKFMyePRv/+Mc/6n3vrFmzkJeXh9LSUqSmpgJQ9fbGjx+PuLg4pKSkwM7ODgDwxRdfoFWrVrW279tvv8XgwYPx97//Hf/85z/x5MkTfPfdd+rXY2NjUVBQgL1798LFxQVz5sxBv379cP78eTRp0gQAUFZWhuXLl2PTpk2wsrLCO++8g1mzZmHz5s14+vQpBg0ahAkTJmDr1q148uQJjh8/DplMhmHDhuHcuXPYt28fDhw4AEDVgxEkJSVh6dKlWLlyJWxsbKBQKPDcc89hx44d8PDwwE8//YSJEyfC29sbMTExtX6/0tLSEB8fj2PHjiErKwuxsbGIiIjA66+/juzsbLRo0QKpqano06cPrK2t6/3+NxYDiBqGF6WarX79+mHKlCkAgDlz5uDjjz/GwYMHERgYiG3btkGhUODzzz9X3+khNTUVbm5uyMzMxBtvvIGVK1ciISEB0dHRAIB169Zh//79Wrejes9s0aJF+POf/9ygAGratCkcHBxQUVEBLy8v9fLo6GjExcVhz5496g/kjRs3IjY2ttY7VyxevBjDhw/H/Pnz1ctCQlS/70LwHD16FN27dwcAbN68Gb6+vti9ezeGDlVVCiorK7Fu3Tq8+OKLAIC4uDgsWLAAgGqytpKSEkRFRalfDwoK0jgWGxsbjeMQvP322xgzZozGsqrtDAgIQFZWFrZv315nAHXo0AHz5s0DALRu3Rpr1qxBRkYGXn/9dXh6egIA3NzcamyDPjGAqOF4b7jaWTuqeiJi7VsHHTp0UH8tk8ng5eWF27dvAwBOnz6NS5cuPXPeqLy8HJcvX0ZJSQmKiorQpUsX9Ws2NjYIDw+HUstS7YEDB5CcnIwLFy6gtLQUT58+RXl5OcrKyhp9twl7e3uMGjUKGzZsQExMDHJycnDu3Dns3bu31vfk5uZiwoQJNb6Wl5cHGxsbjeP18PBAYGAg8vLy1MscHR3V4QIA3t7e6u+pu7s7YmNj0bt3b7z++uvo1asXYmJi4O3tXe/xhIeHP7Ps008/xYYNG3D9+nX897//xZMnT+odxFD1Z169fcbEc0DUcBwNVzuZTFUGE+Oh4x8CQtnoj0ORQaFQAAAePXqETp06ITc3V+Nx8eJFvP322w3eh5WV1TOBVPUk+tWrVxEVFYUOHTrgq6++wsmTJ/Hpp58CUN32SBfjx49Heno6bt68idTUVERGRuL555+vdX0HBwed9gfU/D2tevypqanIyspC9+7dsW3bNrRp0wY///xzvdt1ctIstX755ZeYNWsWxo0bh++//x65ubkYM2ZMvd+zun7mxsQAoobjaDiL07FjRxQUFKBFixZ46aWXNB6urq5wdXWFt7c3jh07pn7P06dPcfLkSY3teHp6apyELy0tRWFhofr5yZMnoVAokJKSgq5du6JNmzb49ddftWqrra1tjdNiBAcHIzw8HJ999hm2bNmCsWPH1rmdDh061DpAISgoCE+fPtU43nv37iE/Px9t27bVqr1hYWFISEjATz/9hPbt22PLli11HkdNhFLglClTEBYWhpdeegmXL1/Wqh01adKkiVGmGNE6gA4dOoQBAwbAx8fnmbHilZWVmDNnjnoki4+PD959912tf5FIwjgazqKMHDkSzZs3x8CBA3H48GEUFhYiMzMT06ZNw82bNwEA06dPx9KlS7F7925cuHABU6ZMwYMHDzS2ExkZiU2bNuHw4cM4e/YsRo8erXFy+6WXXkJlZSVWr16NK1euYNOmTVi3bp1WbfX398eZM2eQn5+Pu3fvavSwxo8fj6VLl0KpVGqMaKvJvHnzsHXrVsybNw95eXk4e/Ysli1bBkB1vmTgwIGYMGECjhw5gtOnT+Odd95Bq1atMHDgwAa1s7CwEAkJCcjKysK1a9fw/fffo6CgQH0eyN/fH4WFhcjNzcXdu3dRUVFR67Zat26NEydOYP/+/bh48SISExORnZ3doHbUxd/fHxkZGSguLsb9+/d13l5ttA6gx48fIyQkRN09rqqsrAw5OTlITExETk4Odu7cifz8fLz55pt6aSxJAMtwFsXR0RGHDh2Cn58foqOjERQUhHHjxqG8vBwuLi4AgJkzZ2LUqFEYPXo0unXrBmdn52c+5BMSEvDqq68iKioK/fv3x6BBgzTOkYSEhGDFihVYtmwZ2rdvj82bNyM5OVmrtk6YMAGBgYEIDw+Hp6cnjh49qn5txIgRsLGxwYgRI+q9c3OPHj2wY8cO7N27F6GhoYiMjMTx48fVr6empqJTp06IiopCt27doFQq8d133z1T1qqNo6MjLly4gCFDhqBNmzaYOHEipk6dikmTJgEAhgwZgj59+uC1116Dp6cntm7dWuu2Jk2ahOjoaAwbNgxdunTBvXv31ANKdJGSkoL09HT4+voiLCxM5+3VRqbU9kxh1TfLZNi1a1edVx9nZ2ejc+fOuHbtGvz8/OrdZmlpKVxdXVFSUqL+BSeJORStGozQ7n0gZLHYrRFFeXk5CgsLERAQ0Ohb0Zuz2NhYPHjwwChX0zfE1atX8eKLLyI7OxsdO3YUuzlmoa7/Aw39HDf4OaCSkhLIZLJar8StqKhAaWmpxoMkjmU4MhGVlZUoLi7GBx98gK5duzJ8JMagAVReXo45c+ZgxIgRtaZgcnKy+mSmq6srfH19Ddkk0geW4chEHD16FN7e3sjOztb6nBIZnsGuA6qsrERMTAyUSiXWrl1b63oJCQmIj49XPy8tLWUISR0vSqV6aHvzS0Pp0aOH1tcjkfEYpAckhM+1a9eQnp5eZw3Qzs4OLi4uGg8yASzDEZGO9B5AQvgUFBTgwIEDFn23WrPGMhwR6UjrEtyjR49w6dIl9XNhvLq7uzu8vb3x1ltvIScnB9988w3kcjmKi4sBqG4/YWtrq7+Wk7hYhiMiHWndAzpx4gTCwsLUY8Pj4+MRFhaGuXPn4pdffsHevXtx8+ZNhIaGwtvbW/346aef9N54EhnLcESkA617QPWd1OMJPwvCKRqISAe8Fxw1Hu8NZ3KqT0hnSoSJ3gT6nrraXEl5qnIGEOmGZTiTsnPnTixcuFDsZujFrFmzNG4aKuUPWm3oO1hXrVql9bD46vf5NBTOB0S6YRnOpLi7u4vdBL1p2rQpmjYVaRZaCaisrGzQ/eeqzqgqNewBkW5YhjMp1Utw/v7+WLJkCcaOHQtnZ2f4+flh/fr1Gu+5ceMGYmJi4ObmBnd3dwwcOBBXr15Vvy6XyxEfHw83Nzd4eHhg9uzZGD16tEZvxN/fHytXrtTYbmhoKJKSktTPV6xYob6Tvq+vL6ZMmYJHj2qf5K9qTyEpKQlpaWnYs2cPZDIZZDIZMjMzERkZibi4OI333blzB7a2trVOuQAAX3/9NV555RXY29ujefPmGjdXvX//Pt599100a9YMjo6O6Nu3LwoKCtSvC6XC/fv3IygoCE2bNkWfPn00pqPIzMxE586d4eTkBDc3N0RERODatWvYuHEj5s+fj9OnT6uPQ+i9yGQyrF27Fm+++SacnJywePFiyOVyjBs3DgEBAXBwcEBgYCBWrVqlcSzVe4Y9evTAtGnTMHv2bLi7u8PLy0vj5+Dv7w8AGDx4MGQymfq5ITCASHcsw6mO++ljcR46fs9TUlIQHh6OU6dOYcqUKZg8eTLy8/MBqP7K7t27N5ydnXH48GEcPXpU/YEqTHqWkpKCjRs3YsOGDThy5Ah+++037Nq1S+t2WFlZ4ZNPPsF//vMfpKWl4YcffsDs2bMb9N5Zs2YhJiZG/UFfVFSE7t27Y/z48diyZYvGlAZffPEFWrVqhcjIyBq39e2332Lw4MHo168fTp06hYyMDHTu3Fn9emxsLE6cOIG9e/ciKysLSqUS/fr105j+oaysDMuXL8emTZtw6NAhXL9+HbNmzQKgmi9p0KBBePXVV3HmzBlkZWVh4sSJkMlkGDZsGGbOnIl27dqpj2PYsGHq7SYlJWHw4ME4e/Ysxo4dC4VCgeeeew47duzA+fPnMXfuXLz//vvYvr3uPwbT0tLg5OSEY8eO4cMPP8SCBQuQnp4OAOrpHFJTU1FUVKSX6R1qwxIc6Y5lOEBeBmwXqRwU80g1M2oj9evXT30L/zlz5uDjjz/GwYMHERgYiG3btkGhUODzzz+H7PeZV1NTU+Hm5obMzEy88cYbWLlyJRISEhAdHQ0AWLduHfbv3691O6r3zBYtWoQ///nP+Mc//lHve5s2bQoHBwdUVFTAy8tLvTw6OhpxcXHYs2cPYmJiAKh6KLGxserjqW7x4sUYPnw45s+fr14WEqL6nS4oKMDevXvVE8EBwObNm+Hr64vdu3dj6FDVH2OVlZVYt26desqJuLg4LFiwAIDqdmMlJSWIiopSvy7MBSQci42NjcZxCN5++22MGTNGY1nVdgYEBCArKwvbt29XH29NOnTogHnz5gFQzSm0Zs0aZGRk4PXXX4enpycAwM3NrcY26BN7QKQ7luFMWocOHdRfy2QyeHl54fbt2wCA06dP49KlS3B2dlafc3F3d0d5eTkuX76MkpISFBUVoUuXLupt2NjYIDw8XOt2HDhwAD179kSrVq3g7OyMUaNG4d69eygrK2v0sdnb22PUqFHYsGEDACAnJwfnzp1DbGxsre/Jzc1Fz549a3wtLy8PNjY2Gsfr4eGBwMBA5OXlqZc5OjpqzHfk7e2t/p66u7sjNjYWvXv3xoABA7Bq1SqN8lxdavq+fvrpp+jUqRM8PT3RtGlTrF+/HtevX69zO1V/5tXbZ0zsAZF++A39/a4IO4AOi4Ba/ro0W9aOqp6IWPvWQfUT2TKZDAqFAoDqziedOnXC5s2bn3mf8JdyQ1hZWT1zjWDVktXVq1cRFRWFyZMnY/HixXB3d8eRI0cwbtw4PHnyBI6OjT/G8ePHIzQ0FDdv3kRqaioiIyPx/PPP17q+g4NDo/clqOl7WvX4U1NTMW3aNOzbtw/btm3DBx98gPT0dHTt2rXO7To5afZ0v/zyS8yaNQspKSnqyQA/+ugjjSnDG9o+4WduTAwg0g9LL8PJZDqVwaSqY8eO2LZtG1q0aFHrjYK9vb1x7Ngx/OlPfwKgOsdx8uRJjbl3PD09Nf7KLy0tRWFhofr5yZMnoVAokJKSAisrVWGmvvMY1dna2kIulz+zPDg4GOHh4fjss8+wZcsWrFmzps7tdOjQARkZGc+UugBVqezp06c4duyYugR379495Ofno23btlq1V7ijTEJCArp164YtW7aga9eutR5HTYRSYNVZUC9fvqxVO2rSpEmTBrdBFyzBkX6wDGeWRo4ciebNm2PgwIE4fPgwCgsLkZmZiWnTpuHmzZsAgOnTp2Pp0qXYvXs3Lly4gClTpuDBgwca24mMjMSmTZtw+PBhnD17FqNHj4a1tbX69ZdeegmVlZVYvXo1rly5gk2bNmk9f4+/vz/OnDmD/Px83L17V6OHNX78eCxduhRKpfKZ6cKrmzdvHrZu3Yp58+YhLy8PZ8+exbJlywCozpcMHDgQEyZMwJEjR3D69Gm88847aNWqFQYOHNigdhYWFiIhIQFZWVm4du0avv/+exQUFKjPA/n7+6vvsXn37l2NARTVtW7dGidOnMD+/ftx8eJFJCYm6mXQgL+/PzIyMlBcXIz79+/rvL3aMIBIfzgazuw4Ojri0KFD8PPzQ3R0NIKCgjBu3DiUl5ere0QzZ87EqFGjMHr0aHUZqPqHfEJCAl599VVERUWhf//+GDRokMY5kpCQEKxYsQLLli1D+/btsXnzZiQnJ2vV1gkTJiAwMBDh4eHw9PTE0aNH1a+NGDECNjY2GDFiRL1TqPfo0QM7duzA3r17ERoaisjISBw/flz9empqKjp16oSoqCh069YNSqUS3333XYOuyQFU39MLFy5gyJAhaNOmDSZOnIipU6di0qRJAIAhQ4agT58+eO211+Dp6YmtW7fWuq1JkyYhOjoaw4YNQ5cuXXDv3j2N3lBjpaSkID09Hb6+vur7fhqCTCmxm7c1dC5xkqDKh8DOFoC8HOiba9ZluPLychQWFiIgIKDeDzRLFBsbiwcPHhjlavqGuHr1Kl588UVkZ2dzWm49qev/QEM/x9kDIv1hGY4kprKyEsXFxfjggw/QtWtXho/EMIBIv1iGIwk5evQovL29kZ2drfU5JTI8joIj/bL00XAEAFrf/NJQ6ps+hsTFHhDpF8twRNRADCDSP5bhiKgBGECkf9XLcERENWAAkf6xDEdEDcAAIsNgGY6I6sEAIsNgGY6I6sEAIsNgGU6Sqs+IakqEmUYFVWdEpdpVnxFVShhAZDgsw0nOzp07sXDhQrGboRezZs3SmFZbyh+02tB3sK5atUrr67JkMplRbqPEC1HJcHhRquS4u7uL3QS9ESbIs1SVlZUNugGqq6urEVrTOOwBkeGwDCc51Utw/v7+WLJkCcaOHQtnZ2f4+flh/fr1Gu+5ceMGYmJi4ObmBnd3dwwcOBBXr15Vvy6XyxEfHw83Nzd4eHhg9uzZGD16tEZvxN/fHytXrtTYbmhoKJKSktTPV6xYgeDgYDg5OcHX1xdTpkzBo0e1T/JXtaeQlJSEtLQ07NmzBzKZDDKZDJmZmYiMjERcXJzG++7cuQNbW1uN3lN1X3/9NV555RXY29ujefPmGnf3vn//Pt599100a9YMjo6O6Nu3LwoKCtSvC6XC/fv3IygoCE2bNkWfPn005kPKzMxE586d4eTkBDc3N0RERODatWvYuHEj5s+fj9OnT6uPQ+i9yGQyrF27Fm+++SacnJywePFiyOVyjBs3DgEBAXBwcEBgYCBWrVqlcSzVe4Y9evTAtGnTMHv2bLi7u8PLy0vj5+Dv7w8AGDx4MGQymfq5ITCAyLAspQynVAJPH4vz0PH7mpKSgvDwcJw6dQpTpkzB5MmTkZ+fD0D1V3bv3r3h7OyMw4cP4+jRo+oP1CdPnqjfv3HjRmzYsAFHjhzBb7/9hl27dmndDisrK3zyySf4z3/+g7S0NPzwww+YPXt2g947a9YsxMTEqD/oi4qK0L17d4wfPx5btmzRmFPniy++QKtWrRAZGVnjtr799lsMHjwY/fr1w6lTp5CRkYHOnTurX4+NjcWJEyewd+9eZGVlQalUol+/fhrzD5WVlWH58uXYtGkTDh06hOvXr2PWrFkAVBP2DRo0CK+++irOnDmDrKwsTJw4ETKZDMOGDcPMmTPRrl079XEMGzZMvd2kpCQMHjwYZ8+exdixY6FQKPDcc89hx44dOH/+PObOnYv333+/3sn80tLS4OTkhGPHjuHDDz/EggULkJ6eDgDq+YRSU1NRVFSkl/mFasMSHBmWpZTh5GXAdpHKQTGPdJqNtV+/fuo5ZObMmYOPP/4YBw8eRGBgILZt2waFQoHPP/8cst+nWU9NTYWbmxsyMzPxxhtvYOXKlUhISEB0dDQAYN26ddi/f7/W7ajeM1u0aBH+/Oc/4x//+Ee9723atCkcHBxQUVEBLy8v9fLo6GjExcVhz549iImJAaDqocTGxqqPp7rFixdj+PDhmD9/vnpZSIjq97agoAB79+5Vz0QKAJs3b4avry92796NoUNVf3BVVlZi3bp16jmP4uLisGDBAgCqqQpKSkoQFRWlfl2YjE44FhsbG43jELz99tvPzNRatZ0BAQHIysrC9u3b1cdbkw4dOmDevHkAVJParVmzBhkZGXj99dfVU627ubnV2AZ9Yg+IDItlOMnr0KGD+muZTAYvLy/cvn0bAHD69GlcunQJzs7O6nMu7u7uKC8vx+XLl1FSUoKioiJ06dJFvQ0bGxuEh4dr3Y4DBw6gZ8+eaNWqFZydnTFq1Cjcu3cPZWVljT42e3t7jBo1Chs2bAAA5OTk4Ny5c4iNja31Pbm5uejZs2eNr+Xl5cHGxkbjeD08PBAYGIi8vDz1MkdHR40J97y9vdXfU3d3d8TGxqJ3794YMGAAVq1apVGeq0tN39dPP/0UnTp1gqenJ5o2bYr169fj+vXrdW6n6s+8evuMiT0gMjy/ocDNXaoyXIdFQC1/eZo0a0dVT0Ssfeug+olsmUwGhUIBAHj06BE6deqEzZs3P/M+4S/lhrCysnrmrtRVS1ZXr15FVFQUJk+ejMWLF8Pd3R1HjhzBuHHj8OTJEzg6Nv4Yx48fj9DQUNy8eROpqamIjIzE888/X+v6Dg4Ojd6XoKbvadXjT01NxbRp07Bv3z5s27YNH3zwAdLT09G1a9c6t+vkpNnT/fLLLzFr1iykpKSoZ6P96KOPcOzYMa3bJ/zMjYkBRIZnCWU4mUynMphUdezYEdu2bUOLFi1qndnS29sbx44dw5/+9CcAqnMcJ0+e1Jj8zdPTU+Ov/NLSUhQWFqqfnzx5EgqFAikpKbCyUhVm6juPUZ2trS3kcvkzy4ODgxEeHo7PPvsMW7ZswZo1a+rcTocOHZCRkfFMqQtQlcqePn2KY8eOqUtw9+7dQ35+Ptq2batVe8PCwhAWFoaEhAR069YNW7ZsQdeuXWs9jpoIpcCq03BfvnxZq3bUpEmTJg1ugy5YgiPDYxnOZI0cORLNmzfHwIEDcfjwYRQWFiIzMxPTpk3DzZs3AQDTp0/H0qVLsXv3bly4cAFTpkzBgwcPNLYTGRmJTZs24fDhwzh79ixGjx4Na2tr9esvvfQSKisrsXr1aly5cgWbNm3SegI5f39/nDlzBvn5+bh7965GD2v8+PFYunQplEqlxoi2msybNw9bt27FvHnzkJeXh7Nnz2LZsmUAVOdLBg4ciAkTJuDIkSM4ffo03nnnHbRq1QoDBw5sUDsLCwuRkJCArKwsXLt2Dd9//z0KCgrU54H8/f1RWFiI3Nxc3L17V2MARXWtW7fGiRMnsH//fly8eBGJiYl6GTTg7++PjIwMFBcX4/79+zpvrzYMIDIOSxkNZ2YcHR1x6NAh+Pn5ITo6GkFBQRg3bhzKy8vVPaKZM2di1KhRGD16tLoMVP1DPiEhAa+++iqioqLQv39/DBo0SOMcSUhICFasWIFly5ahffv22Lx5M5KTk7Vq64QJExAYGIjw8HB4enri6NGj6tdGjBgBGxsbjBgxAvb29nVup0ePHtixYwf27t2L0NBQREZG4vjx4+rXU1NT0alTJ0RFRaFbt25QKpX47rvvGnRNDqD6nl64cAFDhgxBmzZtMHHiREydOhWTJk0CAAwZMgR9+vTBa6+9Bk9PT2zdurXWbU2aNAnR0dEYNmwYunTpgnv37mn0hhorJSUF6enp8PX1RVhYmM7bq41MKbHpAktLS+Hq6oqSkpJau/xkgiofAjtbAPJyoG+uyZfhysvLUVhYiICAgHo/0CxRbGwsHjx4YJSr6Rvi6tWrePHFF5Gdna1RGqTGq+v/QEM/x9kDIuNgGY5EUFlZieLiYnzwwQfo2rUrw0diGEBkPCzDkZEdPXoU3t7eyM7O1vqcEhkeR8GR8VjCaDgCAK1vfmkoPXr0eGb4N0kHe0BkPCzDEVEVDCAyLpbhiOh3DCAyLs6USkS/YwCRcbEMR0S/YwCR8bEMR0RgAJEYWIYjIjQigA4dOoQBAwbAx8enxnnDlUol5s6dC29vbzg4OKBXr14aswUSsQxnfsrKyjBkyBC4uLhAJpPhwYMHNc6CSlSV1tcBPX78GCEhIRg7dqx6AqqqPvzwQ3zyySdIS0tDQEAAEhMT0bt3b5w/f563LKE/mOkUDbL5xj0O5TxplDDT0tJw+PBh/PTTT2jevDlcXV2RnZ2tMX2ATCbDrl27NKaHJsumdQD17dsXffv2rfE1pVKJlStX4oMPPlDfGfaf//wnWrZsid27d2P48OG6tZbMBy9KNQlPnjyBra1tvetdvnwZQUFBaN++vXqZNvMFkWXS6zmgwsJCFBcXo1evXuplrq6u6NKlC7Kysmp8T0VFBUpLSzUeAID8uufsIBPHMpwoevTogbi4OMTFxcHV1RXNmzdHYmKi+m4B/v7+WLhwId599124uLhg4sSJAICvvvoK7dq1g52dHfz9/ZGSkqKxzZSUFBw6dAgymQw9evRQb0sowfn7+wMABg8eDJlMpn5Olk2vAVRcXAwAaNmypcbyli1bql+rLjk5Ga6uruqHr6+v6oXTfwfyVuizeSQ1HA0nirS0NNjY2OD48eNYtWoVVqxYgc8//1z9+vLlyxESEoJTp04hMTERJ0+eRExMDIYPH46zZ88iKSkJiYmJ6tvt7Ny5ExMmTEC3bt1QVFSEnTt3PrNPYY6a1NRUFBUV6WXOGjJ9ot8LLiEhAfHx8ernpaWlf4TQqZmqf4Pia3gnmTyW4UTh6+uLjz/+GDKZDIGBgTh79iw+/vhjTJgwAYBq8riZM2eq1x85ciR69uyJxMREAECbNm1w/vx5fPTRR4iNjYW7uzscHR1ha2sLLy+vGvcplOPc3NxqXYcsj157QMIv1q1btzSW37p1q9ZfOjs7O7i4uGg8AABt56j+PTWTPSFzxTKcKLp27QpZlUEf3bp1Q0FBgXoK5vDwcI318/LyEBERobEsIiJC4z1EjaHXAAoICICXlxcyMjLUy0pLS3Hs2DF069ZNu421SwDaz1V9zRAyXyzDSU7VkWtEhqR1Ce7Ro0e4dOmS+rkwd7m7uzv8/PwwY8YMLFq0CK1bt1YPw/bx8dF+6KVMBgQnqb4+t4DlOHPFMpzRHTt2TOP5zz//jNatW8Pa2rrG9YOCgjSmtwZU8+y0adOm1vfUpEmTJuwxkQate0AnTpxAWFiYep7w+Ph4hIWFYe5cVW9l9uzZ+Mtf/oKJEyfilVdewaNHj7Bv377GXQMkhBB7QuaLZTiju379OuLj45Gfn4+tW7di9erVmD59eq3rz5w5ExkZGVi4cCEuXryItLQ0rFmzBrNmzdJqv/7+/sjIyEBxcTHu37+v62GQGdA6gIQJnqo/hBExMpkMCxYsQHFxMcrLy3HgwAG0adOm8S1kCJk/luGM6t1338V///tfdO7cGVOnTsX06dPVw61r0rFjR2zfvh1ffvkl2rdvj7lz52LBggWIjY3Var8pKSlIT0+Hr6+v+g9YsmwypcSmCywtLYWrqytKSkr+GJAAqD6YziapynEAEJbCcpy5qHwI7GwByMuBvrkmUYYrLy9HYWEhAgICTOoOHz169EBoaChvkUM6q+v/QK2f49WYzs1I2RMyXyzDEVkk0wkggCFkzliGI7I4ol+IqjWOjjNPHA1nFJmZmWI3gUjNtHpAAvaEzA/LcEQWxzQDCGAImSOW4YgsiukGEMAQMjecKZXIoph2AAEMIXPCMhyRRTH9AAIYQuaEZTgii2EeAQQwhMwFy3BEFsN8AghgCJkDluEk6+rVq5DJZMjNzdV5WzKZDLt379Z5O2TazCuAAIaQObCgMpxcoUTW5XvYk/sLsi7fg1xhXseblJSE0NDQZ5YXFRWhb9++xm8QSYrpXYjaELxY1bRZyEWp+84VYf7X51FUUq5e5u1qj3kD2qJPe28RW2Z4nBWVAHPsAQnYEzJdFlCG23euCJO/yNEIHwAoLinH5C9ysO9ckcH2/a9//QvBwcFwcHCAh4cHevXqhcePH0OhUGDBggV47rnnYGdnh9DQUOzbt6/W7WzcuBFubm4ay3bv3q2ebXXjxo2YP38+Tp8+DZlMBplMpnHX/KoluLNnzyIyMlLdpokTJ+LRo0fq12NjYzFo0CAsX74c3t7e8PDwwNSpU1FZWam37wsZn/kGEMAQMmVmXIaTK5SY//V51HRUwrL5X583SDmuqKgII0aMwNixY5GXl4fMzExER0dDqVRi1apVSElJwfLly3HmzBn07t0bb775JgoKChq1r2HDhmHmzJlo164dioqKUFRUhGHDhj2z3uPHj9G7d280a9YM2dnZ2LFjBw4cOIC4uDiN9Q4ePIjLly/j4MGDSEtLw8aNG9WBRqbJvAMIYAiZKjMeDXe88Ldnej5VKQEUlZTjeOFvet93UVERnj59iujoaPj7+yM4OBhTpkxB06ZNsXz5csyZMwfDhw9HYGAgli1bptPUDQ4ODmjatClsbGzg5eUFLy8vODg4PLPeli1bUF5ejn/+859o3749IiMjsWbNGmzatAm3bt1Sr9esWTOsWbMGL7/8MqKiotC/f39kZGQ09ltBEmD+AQQwhEyRGZfhbj+sPXwas542QkJC0LNnTwQHB2Po0KH47LPPcP/+fZSWluLXX39FRESExvoRERHIy8vTezuqysvLQ0hICJycnDT2q1AokJ+fr17Wrl07jSnAvb29cfv2bYO2jQzLMgIIYAiZIjMtw7VwbtgEdg1dTxvW1tZIT0/Hv//9b7Rt2xarV69GYGAgCgsLtd6WlZUVqs9nachzMk2aNNF4LpPJoFAoDLY/MjzLCSCAIWRqzLQM1znAHd6u9pDV8roMqtFwnQPcDbJ/mUyGiIgIzJ8/H6dOnYKtrS0yMjLg4+ODo0ePaqx79OhRtG3btsbteHp64uHDh3j8+LF6WfVrhGxtbSGXy+tsT1BQEE6fPq2xnaNHj8LKygqBgYFaHh2ZEssKIIAhZErMtAxnbSXDvAGqD/XqISQ8nzegLaytaouoxjt27BiWLFmCEydO4Pr169i5cyfu3LmDoKAg/O1vf8OyZcuwbds25Ofn47333kNubi6mT59e47a6dOkCR0dHvP/++7h8+TK2bNnyzKAAf39/FBYWIjc3F3fv3kVFRcUz2xk5ciTs7e0xevRonDt3DgcPHsRf/vIXjBo1Ci1bttT794Ckw/ICCGAImRIzLcP1ae+Nte90hJerZpnNy9Uea9/paLDrgFxcXHDo0CH069cPbdq0wQcffICUlBT07dsX06ZNQ3x8PGbOnIng4GDs27cPe/fuRevWrWvclru7O7744gt89913CA4OxtatW5GUlKSxzpAhQ9CnTx+89tpr8PT0xNatW5/ZjqOjI/bv34/ffvsNr7zyCt566y307NkTa9asMcS3gCREpqxexBVZaWkpXF1dUVJSAhcXF8PuTKkEziapLlYFgLAUXqwqNZUPgZ0tAHk50DdXMhellpeXo7CwEAEBAbC3b/y5GrlCieOFv+H2w3K0cFaV3QzR8yHSt7r+DzT0c9w874TQULxjgvQJZbibu1RlOIkEkL5YW8nQ7UUPsZtBJArLLMFVxXKc9JlpGY7I0jGAAIaQ1JnpaDgiS8cAEjCEpMtMR8MRWToGUFUMIeliGY7I7DCAqmMISZNEy3ASG0RKZDT6+N1nANWEISQ9EivDCbeFKSsrE7klROIQfver3yJJG5Y9DLsuHKItPX5Dfx+OvQPosEj1MxKJtbU13Nzc1DfDdHR0VM+DQ2TOlEolysrKcPv2bbi5uWncIFZbDKC6MISkRWIzpQqzevKOzGSJ3NzcdJ7ZlgFUH4aQdEjsolSZTAZvb2+0aNGCM3OSRWnSpIlOPR8BA6ghGELSIaEynMDa2lov/xmJLA0HITQUByZIg0RHwxGR9hhA2mAIic+Io+Gu3L9i0O0TWToGkLYYQuIzwkWpO/N2Ii03zSDbJiIVBlBjMITEZeAy3JX7VzB2z1i9b5eINDGAGoshJB4DluEqnlZg2L+GoaSiRK/bJaJnMYB0wRASj4HKcLPTZ+PEryf0tj0iqh0DSFcMIXEYoAy3M28nPjn+iV62RUT1YwDpA0PI+PRchuN5HyLjYwDpC0PI+PRUhuN5HyJx6D2A5HI5EhMTERAQAAcHB7z44otYuHChZdy2niFkXHoqw/G8D5E49H4rnmXLlmHt2rVIS0tDu3btcOLECYwZMwaurq6YNm2avncnPbxtj/Ho4d5wPO9DJB6994B++uknDBw4EP3794e/vz/eeustvPHGGzh+/Li+dyVd7AkZjw5lOKVSieO/HEdzx+YGaBgR1UfvAdS9e3dkZGTg4sWLAIDTp0/jyJEj6Nu3b43rV1RUoLS0VONhFhhCxqFDGU4mk2Fpr6UomlmEgYEDDdRAIqqN3gPovffew/Dhw/Hyyy+jSZMmCAsLw4wZMzBy5Mga109OToarq6v64evrq+8miYchZHh6GA0ngwxZN7MAAC80ewEHRh3ApE6T2DMiMjC9B9D27duxefNmbNmyBTk5OUhLS8Py5cuRllbzfbUSEhJQUlKifty4cUPfTRIXQ8jwdBwNd+jaIdx+rJpUbmjboej5Qk+si1qHv3T5iz5bSUTV6H0Qwt/+9jd1LwgAgoODce3aNSQnJ2P06NHPrG9nZwc7Ozt9N0NaODDBsHScKXXH+R3qr2Paxei7dURUC733gMrKymBlpblZa2trKBQKfe/KtLAnZDg6lOHkCjm+yvsKgKr8FuYVpu/WEVEt9B5AAwYMwOLFi/Htt9/i6tWr2LVrF1asWIHBgwfre1emhyFkOI0sw1Uvv8kkMMMqkaXQewlu9erVSExMxJQpU3D79m34+Phg0qRJmDt3rr53ZZpYjjOMRpbhWH4jEo9MKbFbFJSWlsLV1RUlJSVwcXERuzmGo1QCZ5NUIQQAYSkMIV0dilZdlNrufSBkcb2ryxVy+Kzwwe3Ht/FCsxdw6S+X2AMi0oOGfo7zXnBiYTlO/7Qsw7H8RiQuBpCYGEL6peVFqSy/EYmLASQ2hpD+aDEajqPfiMTHAJIChpD+NLAMx/IbkfgYQFLBENKPBpbhWH4jEh8DSEoYQrprQBmO5TciaWAASQ1DSHf1lOFYfiOSBgaQFDGEdFNPGY7lNyJpYABJFUOo8eoow7H8RiQdDCApYwg1Xi1lOJbfiKSDASR1DKHGqaUMx/IbkXQwgEwBQ0h7NZThWH4jkhYGkKlgCGmvWhmO5TciaWEAmRKGkHaqleFYfiOSFgaQqWEINVyVMpz82pcsvxFJDAPIFDGEGu73Mtyh85tYfiOSGAaQqWIINczvZbgdt39RL2L5jUgaGECmjCFUvybOkHv1wVePVE9ZfiOSDgaQqWMI1euQbTvclqu+Hhr0FstvRBLBADIHDKE67bhbrP46xo+9HyKpYACZC4ZQjeQKOb7K/xoA8EITIOy/9U/VTUTGwQAyJwyhZ2hcfNoUkN34V50zpRKR8TCAzA1DSIPGxaeutvXOlEpExsMAMkcMIQA13PvNv5/qhVpmSiUi42IAmSuG0LP3fnv+9+t/apkplYiMiwFkziw8hJ6591s9M6USkXExgMydhYZQjVMv1DFTKhEZHwPIElhgCNU69UItM6USkfExgCyFhYVQrVMvsAxHJBkMIEtiISFU58ynLMMRSQYDyNJYQAjVO/Mpy3BEksAAskRmHkL1znzKMhyRJDCALJWZhlCd5TcBy3BEksAAsmRmGEL1lt8ELMMRiY4BZOnMLIS2/+ePHk2dM5+yDEckOgYQmU0IyRVy7LywE0ADZj5lGY5IdAwgUjGDEGpw+U3AMhyRqKQbQPxAMD4TD6EGl98ELMMRiUq6AfR1a+DIcKBgHVBygYFkLCYaQlqV3wQswxGJSroBVH4HuL4NyJ4MfBsE7PJmIBmLCYaQ1uU3ActwRKKxEbsBtXrt30BZNnA7E7j7E1B+SxVI17epXrdvCbToAbTsofrXJVD1wUn6IYQQAJxboAohAAiKF61JddG6/CaoXoZrFmKA1hFRTQzSA/rll1/wzjvvwMPDAw4ODggODsaJEye024hndyA4EeiZAbz1AOh1CAheALSMVH1gCIHEHpLhmEhPqFHlNwHLcESi0XsP6P79+4iIiMBrr72Gf//73/D09ERBQQGaNWvW+I1a2wEt/lf1QCIgrwDuHQduZbKHZGgm0BNqdPlN4DcUuLlLVYbrsIi/J0RGovcAWrZsGXx9fZGamqpeFhAQoN+dMJCMS+Ih1Ojym4BlOCJRyJRK/daq2rZti969e+PmzZv48ccf0apVK0yZMgUTJkyocf2KigpUVFSon5eWlsLX1xclJSVwcXFpXCNqCiR5ueY6DCTtKZXA2SRVCAFAWIroISRXyOGzwge3H9/GC81ewKW/XNK+BwQAh6JVvaB27wMhi/XfUCILUlpaCldX13o/x/V+DujKlStYu3YtWrdujf3792Py5MmYNm0a0tLSalw/OTkZrq6u6oevr6/ujRB6SDyHpF8SPCekc/lNwNFwREan9x6Qra0twsPD8dNPP6mXTZs2DdnZ2cjKynpmfYP0gOrDHpJuJNQTmvzNZKw7uQ4AcHLiSXT07ti4DVU+BHa2UP0e9M1lGY5IBw3tAen9HJC3tzfatm2rsSwoKAhfffVVjevb2dnBzs5O382oG88h6UYi54R0Gv1WnTAa7uYu1Wg4BhCRwek9gCIiIpCfn6+x7OLFi3j++ef1vSv9YSBpTwIhpLfym4Cj4YiMSu8B9Ne//hXdu3fHkiVLEBMTg+PHj2P9+vVYv369vndlOAykhhE5hHQe/VYdR8MRGZXezwEBwDfffIOEhAQUFBQgICAA8fHxtY6Cq66htUNR8RySJhHOCelt9Ft1HA1HpLOGfo4bJIB0YRIBVB0DyeghdLDwICL/GQkAmBMxB0t7LdXPhq9uBX56G3BuDUTlm9fPiMhIRBuEYJFYsjN6OU7v5TcBy3BERsMAMgRLDSQjhZBeR79Vx9FwREbDADIGSwokI4SQ3ke/VcfRcERGwQASg7kHkoFDyGDlNwHLcERGwQCSAnMMJAOFkEHLbwKW4YiMggEkReYSSAYIIYOX3wQswxEZHAPIFJhyIOk5hAxefhOwDEdkcAwgU2RqgaSnEDJK+U3AMhyRwTGAzIEpBJIeQsho5TcBy3BEBsUAMkdSDSQdQ8ho5TcBy3BEBsUAsgRSCqRGhpBRy28CluGIDIoBZInEDqRGhJDRy28CluGIDIYBROIEkpYhZPTym4BlOCKDYQDRs4wVSA0MIVHKbwKW4YgMhgFE9TNkIDUghEQrvwlYhiMyCAYQaU/fgVRLCMkD/4rjhb9h5aEjsJMHo8LqP8YtvwlYhiMyCAYQ6U5fgVQlhPZl/gvzd7yAorImADrCCx0Bq/u4ddcb8Dby8bEMR2QQnBGVDE+bGWNbvIp95+9g8s8dofrFrFruUkIGGda+0xF92hs5hThTKlGDcUpukq46AkmutML/XPh/KKpsDs3wUZEB8HK1x5E5kbC2MmIIVD4EdrZQtbNvLntBRHXglNwkXXWU7I5fuISiSs9a36oEUFRSjuOF99DtxeZGazLLcET6ZyV2A4jUgRSciNutlzToLbd/GAscGQ4UrANKLgDG6Mj7DVX9e32HcfZHZObYAyJJaeFs37D1lFeB62eNey87joYj0isGEElK5wB3eLvao7ikHDX1MWRQwMvxKTq/uQa486Nx72XHMhyRXnEQAknOvnNFmPxFDgBohJDs92drn1+CPj3e+uNiVW1G2ekaSBwNR1QvjoIjk7bvXBHmf30eRSV/BIm3qz3mdTiLPiV/VS0IS6n53nGGDCSOhiOqFwOITJ5cocTxwt9w+2E5Wjjbo3OAO6xlAM4mqe6YANQeQhob0nMgHYpWleHavQ+ELNbhCInMEwOIzJdSqX0IVaVrILEMR1QnBhCZN11DqCptA8n9FSA9AlBUsAxHVAMGEJk/fYZQVQ0JJCs7VQB59wE6fmyYYd9EJooBRJbBUCFUlTFH2RGZAQYQWQ5jhFBV8gpVGP0YBSif/tEbqoqBRBaM94Ijy6Hl9N46s7YDfHoDrQaoRsMFTlfdJcGQU5gTmSH2gMh8GLsnVNtoOJbsyMKxBEeWyZgh1NCLUhlIZGEYQGS5jBlCjbkolYFEZo4BRJbNWCGkj4tSGUhkZhhARMYIIUPcG46BRCaOAUQEGCeEDH1vOAYSmRgGEJHA0CFk7HvDMZBI4hhARFUZMoTEnqKBgUQSwwAiqs6QISSlKRoYSCQyBhBRTQwVQlKeooGBREbW0M9xK0M3ZOnSpZDJZJgxY4ahd0VUP+G2Pe3nqp6fmgnkrdB9u62iAGt74GEB8OCM7tvTJ2s7oMX/AsGJQM8M4K0HQK9DQPACoGWkqt3CrYOyJwPfBgG7vIEjw4GCdUDJBVVwE+mZQe8Fl52djf/7v/9Dhw4dDLkbIu0Y4t5xTZwB776qMtz17dKeI0gIpBb/CyCx5h4S72VHRmCwHtCjR48wcuRIfPbZZ2jWrJmhdkPUOIboCfkNVf17fYdp9RjYQyKRGKwHNHXqVPTv3x+9evXCokWLal2voqICFRV/3Mq+tLTUUE0i0qTvnlD1MpyUe0F1YQ+JjMQgAfTll18iJycH2dnZ9a6bnJyM+fPnG6IZRPXTZwiZUhlOGwwkMhC9j4K7ceMGwsPDkZ6erj7306NHD4SGhmLlypXPrF9TD8jX15ej4Mi49DU6Tsqj4QyFo+yoGtGGYe/evRuDBw+GtbW1eplcLodMJoOVlRUqKio0Xmtsw4n0Th8hJPZFqVLAQLJ4ogXQw4cPce3aNY1lY8aMwcsvv4w5c+agffv2db6fAUSi0kcISemiVClgIFkc0abkdnZ2fiZknJyc4OHhUW/4EIlOH+eE/Ib+fh5oB9BhET9IeQ6JamHQ64CITJKuIWQuo+EMhYFEv+OteIhqo0s5jmW4xmPJzuTxXnBE+tDYELLE0XCGwkAyOQwgIn1pTAhxNJzhMJAkjwFEpE+NCSGW4YyDgSQ5DCAifdM2hFiGEwcDSXQMICJD0CaEWIaTBgaS0TGAiAxFmxBiGU56GEgGxwAiMqSGhhDLcNLHQNI7BhCRoTUkhFiGMz0MJJ0xgIiMoSEhxDKcaWMgaY0BRGQs9YUQy3DmhYFULwYQkTHVFUIsw5k3BtIzGEBExlZXCLEMZzkYSAwgIlHUFkIsw1kuCwwkBhCRWGoKoZcmsAxHKhYQSAwgIjHVFEJ3jrAMR88yw0BiABGJrXoI+Y8Erm5mGY7qZgaBxAAikoLqISSzAZRPWYajhjPBQGIAEUlF9RACWIajxjOBQGIAEUlJ9RCyaw5E35ZU2YRMlAQDiQFEJDVKJXA6ATi/TPX85ZlAx+XitonMjwQCiQFEJEVKJfBNW+DhBdXzhsysSqQLEQKJAUQkVYVbgKyRfzxnCJExGSGQGEBEUlX5EPjKE1BU/LGMIURiMUAgMYCIpEy4N5zn/6guUAUYQiQNeggkBhCRlAn3hmv6EuA/Aji3ULWcIURS04hAKoU3XN3cGEBEklR1ioY+p1S9ofqm9yaSggYEUqnCE66j7tT7OW5j2JYSUY2aOAPefVXBc2MH0GGRavm5BcCpmaqvGUIkRdZ2QIv/VT2QWHMgld1p0KasDNlOIqqD31DVv9d3qP4NTgLaz1V9fWomkLdClGYRaUUIpOBEoGcG8NYD4LV/N+itDCAisbSKAqztgYcFwIMzqpO4DCEyddZ2gGf3Bq3KACISi1CGA4Dr21X/MoTIgjCAiMRUtQwnjAdiCJGFYAARial6GU7AECILwAAiElNNZTgBQ4jMHAOISGw1leEEDCEyYwwgIrHVVoYTMITITDGAiMRWVxlOwBAiM8QAIpKCuspwAoYQmRkGEJEU1FeGEzCEyIwwgIikoCFlOAFDiMwEA4hIKhpShhMwhMgMMICIpKKhZTgBQ4hMHAOISCq0KcMJGEJkwvQeQMnJyXjllVfg7OyMFi1aYNCgQcjPz9f3bojMkzZlOAFDiEyU3gPoxx9/xNSpU/Hzzz8jPT0dlZWVeOONN/D48WN974rI/GhbhhMwhMgE6X1G1H379mk837hxI1q0aIGTJ0/iT3/6k753R2Reqs6Uen070Cyk4e8VQgjgzKpkEgx+DqikpAQA4O7uXuPrFRUVKC0t1XgQWbTGlOEE7AmRCTFoACkUCsyYMQMRERFo3759jeskJyfD1dVV/fD19TVkk4ikr7FlOAFDiEyEQQNo6tSpOHfuHL788sta10lISEBJSYn6cePGDUM2iUj6GjMarjqGEJkAgwVQXFwcvvnmGxw8eBDPPfdcrevZ2dnBxcVF40Fk8XQpwwkYQiRxeg8gpVKJuLg47Nq1Cz/88AMCAgL0vQsi86drGU7AECIJ03sATZ06FV988QW2bNkCZ2dnFBcXo7i4GP/973/1vSsi86WPMpyAIUQSpfcAWrt2LUpKStCjRw94e3urH9u2bdP3rojMmz7KcAKGEEmQ3q8DUur6H4WIVKqX4bS5JqgmvE6IJIb3giOSKn2W4QTsCZGEMICIpEyfZTgBQ4gkggFEJGX6Gg1XHUOIJIABRCRlhijDCRhCJDIGEJHUGaIMJ2AIkYgYQERSZ6gynIAhRCJhABFJnSHLcAKGEImAAURkCgxZhhMwhMjIGEBEpsDQZTgBQ4iMiAFEZAqMUYYTMITISBhARKbCGGU4AUOIjIABRGQqjFWGEzCEyMAYQESmwphlOAFDiAyIAURkSoxZhhMwhMhAGEBEpsTYZTgBQ4gMgAFEZErEKMMJGEKkZwwgIlMjRhlOwBAiPWIAEZkascpwAoYQ6QkDiMjUiFmGEzCESA8YQESmSMwynIAhRDpiABGZIrHLcAKGEOmAAURkiqRQhhMwhKiRGEBEpkoKZTgBQ4gagQFEZKqkUoYTMIRISwwgIlMlpTKcgCFEWmAAEZkyKZXhBAwhaiAGEJEpk1oZTsAQogZgABGZMimW4QQMIaoHA4jI1EmxDCdgCFEdGEBEpk6qZTgBQ4hqwQAiMnVSLsMJGEJUAwYQkTmQchlOwBCiahhAROZA6mU4AUOIqmAAEZkDUyjDCRhC9DsGEJG5MIUynIAhRGAAEZkPUynDCRhCFo8BRGQuTKkMJ2AIWTQGEJE5MaUynIAhZLEYQETmxNTKcAKGkEViABGZE1MswwkYQhaHAURkbkyxDCdgCFkUBhCRuTHVMpyAIWQxDBZAn376Kfz9/WFvb48uXbrg+PHjhtoVEVVlymU4AUPIIhgkgLZt24b4+HjMmzcPOTk5CAkJQe/evXH79m1D7I6IqjPlMpyAIWT2DBJAK1aswIQJEzBmzBi0bdsW69atg6OjIzZs2GCI3RFRdaZehhMwhMya3gPoyZMnOHnyJHr16vXHTqys0KtXL2RlZT2zfkVFBUpLSzUeRKQjcyjDCRhCZkvvAXT37l3I5XK0bNlSY3nLli1RXFz8zPrJyclwdXVVP3x9ffXdJCLLZA5lOAFDyCyJPgouISEBJSUl6seNGzfEbhKReTCXMpyAIWR29B5AzZs3h7W1NW7duqWx/NatW/Dy8npmfTs7O7i4uGg8iEgPzKkMJ2AImRW9B5CtrS06deqEjIwM9TKFQoGMjAx069ZN37sjorqYUxlOwBAyGzaG2Gh8fDxGjx6N8PBwdO7cGStXrsTjx48xZswYQ+yOiGpTvQzXLETsFumHEEIAcG6BKoQAIChetCaR9gwSQMOGDcOdO3cwd+5cFBcXIzQ0FPv27XtmYAIRGZhQhru5S1WGM5cAAhhCZkCmVEqrX15aWgpXV1eUlJTwfBCRPlzdCvz0NuDcGojKV31wmxOlEjibpAohAAhLYQiJrKGf46KPgiMiAzO30XDV8ZyQyWIAEZk7cxwNVx1DyCQxgIgsgTmOhquOIWRyGEBElsDcy3AChpBJYQARWQJLKMMJGEImgwFEZCksoQwnYAiZBAYQkaWwlDKcgCEkeQwgIkthSWU4AUNI0hhARJbEkspwAoaQZDGAiCyJpZXhBAwhSWIAEVkSSyzDCRhCksMAIrI0lliGEzCEJIUBRGRpLLUMJ2AISQYDiMjSWHIZTsAQkgQGEJElsuQynIAhJDoGEJElsvQynIAhJCoGEJElYhnuDwwh0TCAiCwVy3B/YAiJwkbsBlQnzBBeWloqckuIzJzzn4AKO6CsALj+E9AsWOwWie/5eOBRBXB+GXB0JvCoHAiME7tVJkf4/FbW84eNTFnfGkZ28+ZN+Pr6it0MIiLS0Y0bN/Dcc8/V+rrkAkihUODXX3+Fs7MzZDKZQfdVWloKX19f3LhxAy4uLgbdl7GY2zGZ2/EAPCZTwWNqPKVSiYcPH8LHxwdWVrWf6ZFcCc7KyqrOxDQEFxcXs/kFE5jbMZnb8QA8JlPBY2ocV1fXetfhIAQiIhIFA4iIiERh0QFkZ2eHefPmwc7OTuym6I25HZO5HQ/AYzIVPCbDk9wgBCIisgwW3QMiIiLxMICIiEgUDCAiIhIFA4iIiERhsQH06aefwt/fH/b29ujSpQuOHz8udpMaLTk5Ga+88gqcnZ3RokULDBo0CPn5+WI3S6+WLl0KmUyGGTNmiN0Unfzyyy9455134OHhAQcHBwQHB+PEiRNiN6vR5HI5EhMTERAQAAcHB7z44otYuHBhvfcAk4pDhw5hwIAB8PHxgUwmw+7duzVeVyqVmDt3Lry9veHg4IBevXqhoKBAnMY2UF3HVFlZiTlz5iA4OBhOTk7w8fHBu+++i19//VWUtlpkAG3btg3x8fGYN28ecnJyEBISgt69e+P27dtiN61RfvzxR0ydOhU///wz0tPTUVlZiTfeeAOPHz8Wu2l6kZ2djf/7v/9Dhw4dxG6KTu7fv4+IiAg0adIE//73v3H+/HmkpKSgWbNmYjet0ZYtW4a1a9dizZo1yMvLw7Jly/Dhhx9i9erVYjetQR4/foyQkBB8+umnNb7+4Ycf4pNPPsG6detw7NgxODk5oXfv3igvLzdySxuurmMqKytDTk4OEhMTkZOTg507dyI/Px9vvvmmCC0FoLRAnTt3Vk6dOlX9XC6XK318fJTJyckitkp/bt++rQSg/PHHH8Vuis4ePnyobN26tTI9PV356quvKqdPny52kxptzpw5yv/5n/8Ruxl61b9/f+XYsWM1lkVHRytHjhwpUosaD4By165d6ucKhULp5eWl/Oijj9TLHjx4oLSzs1Nu3bpVhBZqr/ox1eT48eNKAMpr164Zp1FVWFwP6MmTJzh58iR69eqlXmZlZYVevXohKytLxJbpT0lJCQDA3d1d5JboburUqejfv7/Gz8tU7d27F+Hh4Rg6dChatGiBsLAwfPbZZ2I3Syfdu3dHRkYGLl68CAA4ffo0jhw5gr59+4rcMt0VFhaiuLhY43fP1dUVXbp0MZvPCkD1eSGTyeDm5mb0fUvuZqSGdvfuXcjlcrRs2VJjecuWLXHhwgWRWqU/CoUCM2bMQEREBNq3by92c3Ty5ZdfIicnB9nZ2WI3RS+uXLmCtWvXIj4+Hu+//z6ys7Mxbdo02NraYvTo0WI3r1Hee+89lJaW4uWXX4a1tTXkcjkWL16MkSNHit00nRUXFwNAjZ8Vwmumrry8HHPmzMGIESNEueGqxQWQuZs6dSrOnTuHI0eOiN0Undy4cQPTp09Heno67O3txW6OXigUCoSHh2PJkiUAgLCwMJw7dw7r1q0z2QDavn07Nm/ejC1btqBdu3bIzc3FjBkz4OPjY7LHZCkqKysRExMDpVKJtWvXitIGiyvBNW/eHNbW1rh165bG8lu3bsHLy0ukVulHXFwcvvnmGxw8eNDoU1ro28mTJ3H79m107NgRNjY2sLGxwY8//ohPPvkENjY2kMvlYjdRa97e3mjbtq3GsqCgIFy/fl2kFunub3/7G9577z0MHz4cwcHBGDVqFP76178iOTlZ7KbpTPg8MMfPCiF8rl27hvT0dNGmm7C4ALK1tUWnTp2QkZGhXqZQKJCRkYFu3bqJ2LLGUyqViIuLw65du/DDDz8gICBA7CbprGfPnjh79ixyc3PVj/DwcIwcORK5ubmwtrYWu4lai4iIeGZ4/MWLF/H888+L1CLdlZWVPTPhmLW1NRQKhUgt0p+AgAB4eXlpfFaUlpbi2LFjJvtZAfwRPgUFBThw4AA8PDxEa4tFluDi4+MxevRohIeHo3Pnzli5ciUeP36MMWPGiN20Rpk6dSq2bNmCPXv2wNnZWV2fdnV1hYODg8itaxxnZ+dnzmE5OTnBw8PDZM9t/fWvf0X37t2xZMkSxMTE4Pjx41i/fj3Wr18vdtMabcCAAVi8eDH8/PzQrl07nDp1CitWrMDYsWPFblqDPHr0CJcuXVI/LywsRG5uLtzd3eHn54cZM2Zg0aJFaN26NQICApCYmAgfHx8MGjRIvEbXo65j8vb2xltvvYWcnBx88803kMvl6s8Ld3d32NraGrexRh93JxGrV69W+vn5KW1tbZWdO3dW/vzzz2I3qdEA1PhITU0Vu2l6ZerDsJVKpfLrr79Wtm/fXmlnZ6d8+eWXlevXrxe7STopLS1VTp8+Xenn56e0t7dXvvDCC8q///3vyoqKCrGb1iAHDx6s8f/O6NGjlUqlaih2YmKismXLlko7Oztlz549lfn5+eI2uh51HVNhYWGtnxcHDx40els5HQMREYnC4s4BERGRNDCAiIhIFAwgIiISBQOIiIhEwQAiIiJRMICIiEgUDCAiIhIFA4iIiETBACIiIlEwgIiISBQMICIiEgUDiIiIRPH/ATrwGMrebnApAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "c = -jnp.array([1, 3])\n", + "x = optax.linprog.rhpdhg(c, A, b, G, h, 1_000_000)['primal']\n", + "print(x, c @ x)\n", + "plot_lp(c, A, b, G, h, x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16ad881e-bee1-4981-8b71-ee2f1fcdab62", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/optax/__init__.py b/optax/__init__.py index 0840b1d5a..f623db676 100644 --- a/optax/__init__.py +++ b/optax/__init__.py @@ -19,6 +19,7 @@ from optax import assignment from optax import contrib +from optax import linprog from optax import losses from optax import monte_carlo from optax import perturbations @@ -364,6 +365,7 @@ "lion", "linear_onecycle_schedule", "linear_schedule", + "linprog", "log_cosh", "lookahead", "LookaheadParams", diff --git a/optax/_src/alias.py b/optax/_src/alias.py index b2c2d1865..d251e2341 100644 --- a/optax/_src/alias.py +++ b/optax/_src/alias.py @@ -2482,7 +2482,7 @@ def lbfgs( ... ) ... params = optax.apply_updates(params, updates) ... print('Objective function: ', f(params)) - Objective function: 7.5166864 + Objective function: 7.516686... Objective function: 7.460699e-14 Objective function: 2.6505726e-28 Objective function: 0.0 diff --git a/optax/linprog/__init__.py b/optax/linprog/__init__.py new file mode 100644 index 000000000..6e27365a2 --- /dev/null +++ b/optax/linprog/__init__.py @@ -0,0 +1,19 @@ +# Copyright 2024 DeepMind Technologies Limited. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== +"""The linear programming sub-package.""" + +# pylint:disable=g-importing-member + +from optax.linprog._rhpdhg import solve_general as rhpdhg diff --git a/optax/linprog/_rhpdhg.py b/optax/linprog/_rhpdhg.py new file mode 100644 index 000000000..01a98b155 --- /dev/null +++ b/optax/linprog/_rhpdhg.py @@ -0,0 +1,215 @@ +# Copyright 2024 DeepMind Technologies Limited. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== +"""The restarted Halpern primal-dual hybrid gradient method.""" + +from jax import lax, numpy as jnp +from optax import tree_utils as otu + + +def solve_canonical( + c, A, b, iters, reflect=True, restarts=True, tau=None, sigma=None +): + r"""Solves a linear program using the restarted Halpern primal-dual hybrid + gradient (RHPDHG) method. + + Minimizes :math:`c \cdot x` subject to :math:`A x = b` and :math:`x \geq 0`. + + See also `MPAX `_. + + Args: + c: Cost vector. + A: Equality constraint matrix. + b: Equality constraint vector. + iters: Number of iterations to run the solver for. + reflect: Use reflection. See paper for details. + restarts: Use restarts. See paper for details. + tau: Primal step size. See paper for details. + sigma: Dual step size. See paper for details. + + Returns: + A dictionary whose entries are as follows: + - primal: The final primal solution. + - dual: The final dual solution. + - primal_iterates: The primal iterates. + - dual_iterates: The dual iterates. + + Examples: + >>> from jax import numpy as jnp + >>> import optax + >>> c = -jnp.array([2, 1]) + >>> A = jnp.zeros([0, 2]) + >>> b = jnp.zeros(0) + >>> G = jnp.array([[3, 1], [1, 1], [1, 4]]) + >>> h = jnp.array([21, 9, 24]) + >>> x = optax.linprog.rhpdhg(c, A, b, G, h, 1_000_000)['primal'] + >>> print(x[0]) + 5.99... + >>> print(x[1]) + 2.99... + + References: + Haihao Lu, Jinwen Yang, `Restarted Halpern PDHG for Linear Programming + `_, 2024 + Haihao Lu, Zedong Peng, Jinwen Yang, `MPAX: Mathematical Programming in JAX + `_, 2024 + """ + + if tau is None or sigma is None: + A_norm = jnp.linalg.norm(A, axis=(0, 1), ord=2) + if tau is None: + tau = 1 / (2 * A_norm) + if sigma is None: + sigma = 1 / (2 * A_norm) + + def T(z): + # primal dual hybrid gradient (PDHG) + x, y = z + xn = x + tau * (y @ A - c) + xn = xn.clip(min=0) + yn = y + sigma * (b - A @ (2 * xn - x)) + return xn, yn + + def H(z, k, z0): + # Halpern PDHG + Tz = T(z) + if reflect: + zc = otu.tree_sub(otu.tree_scalar_mul(2, Tz), z) + else: + zc = Tz + kp2 = k + 2 + zn = otu.tree_add( + otu.tree_scalar_mul((k + 1) / kp2, zc), + otu.tree_scalar_mul(1 / kp2, z0), + ) + return zn, Tz + + def update(carry, _): + z, k, z0, d0 = carry + zn, Tz = H(z, k, z0) + + if restarts: + d = otu.tree_l2_norm(otu.tree_sub(z, Tz), squared=True) + restart = d <= d0 * jnp.exp(-2) + new_carry = otu.tree_where( + restart, + (zn, 0, zn, d), + (zn, k + 1, z0, d0), + ) + else: + new_carry = zn, k + 1, z0, d0 + + return new_carry, z + + def run(): + m, n = A.shape + x = jnp.zeros(n) + y = jnp.zeros(m) + z0 = x, y + d0 = otu.tree_l2_norm(otu.tree_sub(z0, T(z0)), squared=True) + (z, _, _, _), zs = lax.scan(update, (z0, 0, z0, d0), length=iters) + x, y = z + xs, ys = zs + return { + "primal": x, + "dual": y, + "primal_iterates": xs, + "dual_iterates": ys, + } + + return run() + + +def general_to_canonical(c, A, b, G, h): + """Converts a linear program from general form to canonical form. + + The solution to the new linear program will consist of the concatenation of + - the positive part of x + - the negative part of x + - slacks + + That is, we go from + + Minimize c · x subject to + A x = b + G x ≤ h + + to + + Minimize c · (x⁺ - x⁻) subject to + A (x⁺ - x⁻) = b + G (x⁺ - x⁻) + s = h + x⁺, x⁻, s ≥ 0 + + Args: + c: Cost vector. + A: Equality constraint matrix. + b: Equality constraint vector. + G: Inequality constraint matrix. + h: Inequality constraint vector. + + Returns: + A triple (c', A', b') representing the corresponding canonical form. + """ + c_can = jnp.concatenate([c, -c, jnp.zeros(h.size)]) + G_ = jnp.concatenate([G, -G, jnp.eye(h.size)], 1) + A_ = jnp.concatenate([A, -A, jnp.zeros([b.size, h.size])], 1) + A_can = jnp.concatenate([A_, G_], 0) + b_can = jnp.concatenate([b, h]) + return c_can, A_can, b_can + + +def solve_general( + c, A, b, G, h, iters, reflect=True, restarts=True, tau=None, sigma=None +): + r"""Solves a linear program using the restarted Halpern primal-dual hybrid + gradient (RHPDHG) method. + + Minimizes :math:`c \cdot x` subject to :math:`A x = b` and :math:`G x \leq h`. + + See also `MPAX `_. + + Args: + c: Cost vector. + A: Equality constraint matrix. + b: Equality constraint vector. + G: Inequality constraint matrix. + h: Inequality constraint vector. + iters: Number of iterations to run the solver for. + reflect: Use reflection. See paper for details. + restarts: Use restarts. See paper for details. + tau: Primal step size. See paper for details. + sigma: Dual step size. See paper for details. + + Returns: + A dictionary whose entries are as follows: + - primal: The final primal solution. + - slacks: The final primal slack values. + - canonical_result: The result for the canonical program that was used + internally to find this solution. See paper for details. + + References: + Haihao Lu, Jinwen Yang, `Restarted Halpern PDHG for Linear Programming + `_, 2024 + Haihao Lu, Zedong Peng, Jinwen Yang, `MPAX: Mathematical Programming in JAX + `_, 2024 + """ + canonical = general_to_canonical(c, A, b, G, h) + result = solve_canonical(*canonical, iters, reflect, restarts, tau, sigma) + x_pos, x_neg, slacks = jnp.split(result["primal"], [c.size, c.size * 2]) + return { + "primal": x_pos - x_neg, + "slacks": slacks, + "canonical_result": result, + } diff --git a/optax/linprog/_rhpdhg_test.py b/optax/linprog/_rhpdhg_test.py new file mode 100644 index 000000000..d2d2832b9 --- /dev/null +++ b/optax/linprog/_rhpdhg_test.py @@ -0,0 +1,95 @@ +# Copyright 2024 DeepMind Technologies Limited. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== +"""Tests for the restarted Halpern primal-dual hybrid gradient method.""" + +from functools import partial + +from absl.testing import absltest +from absl.testing import parameterized +import jax +from jax import numpy as jnp +import numpy as np +import cvxpy as cp + +from optax.linprog import rhpdhg + + +def solve_cvxpy(c, A, b, G, h): + x = cp.Variable(c.size) + constraints = [] + if A.shape[0] > 0: + constraints.append(A @ x == b) + if G.shape[0] > 0: + constraints.append(G @ x <= h) + objective = cp.Minimize(c @ x) + problem = cp.Problem(objective, constraints) + problem.solve(solver='GLPK') + return x.value, problem.status + + +class RHPDHGTest(parameterized.TestCase): + + def setUp(self): + super().setUp() + self.f = jax.jit(partial(rhpdhg, iters=100_000)) + + @parameterized.parameters( + dict(n_vars=n_vars, n_eq=n_eq, n_ineq=n_ineq) + for n_vars in range(8) + for n_eq in range(n_vars) + for n_ineq in range(8) + if n_eq + n_ineq >= n_vars + # Make sure set of solvable LPs with these shapes is not null in measure. + ) + def test_hungarian_algorithm(self, n_vars, n_eq, n_ineq): + # Find a solvable LP. + while True: + + c = np.random.normal(size=(n_vars,)) + A = np.random.normal(size=(n_eq, n_vars)) + b = np.random.normal(size=(n_eq,)) + G = np.random.normal(size=(n_ineq, n_vars)) + h = np.random.normal(size=(n_ineq,)) + + # For numerical testing purposes, constrain x to [-limit, limit]. + limit = 5 + G = jnp.concatenate([G, jnp.eye(n_vars), -jnp.eye(n_vars)]) + h = jnp.concatenate([h, jnp.full(n_vars * 2, limit)]) + + r, status = solve_cvxpy(c, A, b, G, h) + + if status == 'optimal': + break + + result = self.f(c, A, b, G, h) + x = result['primal'] + + rtol = 1e-2 + atol = 1e-2 + + with self.subTest('x approximately satisfies equality constraints'): + np.testing.assert_allclose(A @ x, b, rtol=rtol, atol=atol) + + with self.subTest('x approximately satisfies inequality constraints'): + np.testing.assert_allclose((G @ x).clip(min=h), h, rtol=rtol, atol=atol) + + with self.subTest('x is approximately as good as the reference solution'): + cx = c @ x + cr = c @ r + np.testing.assert_allclose(cx.clip(min=cr), cr, rtol=rtol, atol=atol) + + +if __name__ == '__main__': + absltest.main() diff --git a/pyproject.toml b/pyproject.toml index 7bccb6db3..1d21ce9ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,7 +49,8 @@ test = [ "dm-tree>=0.1.7", "flax>=0.5.3", "scipy>=1.7.1", - "scikit-learn" + "scikit-learn", + "cvxpy[GLPK]", ] examples = [