{ "cells": [ { "cell_type": "markdown", "id": "hearing-simon", "metadata": {}, "source": [ "# Warm-starting quantum optimization\n", "\n", "## Introduction\n", "\n", "Optimization problems with integer variables or constraints are often hard to solve. For example, the Quadratic Unconstrained Binary Optimization (QUBO) problem, i.e.\n", "\n", "\\begin{align}\n", "\\min_{x\\in\\{0,1\\}^n}x^T\\Sigma x + \\mu^Tx,\n", "\\end{align}\n", "\n", "is NP-Hard. Here, $\\Sigma$ is an $n\\times n$ matrix and $x$ is a vector of $n$ binary variables. Note that we could have added the linear term $\\mu$ to the diagonal as $x_i^2=x_i$ for $x_i\\in\\{0, 1\\}$. While QUBOs are hard to solve there exists many ways to relax them to problems that are easier to solve.\n", "For example, if $\\Sigma$ is semi-definite positive the QUBO can be relaxed and results in a convex Quadratic Program \n", "\n", "\\begin{align}\n", "\\min_{x\\in[0,1]^n}x^T\\Sigma x,\n", "\\end{align}\n", "\n", "which becomes easy to solve as $x$ now represents $n$ continuous variables bound to the range $[0, 1]$.\n", "Such relaxations can be leveraged to warm-start quantum optimization algorithms as shown in [1].\n", "\n", "## References\n", "\n", "[1] [D. J. Egger, J Marecek, S. Woerner, *Warm-starting quantum optimization*, arXiv:2009.10095](http://arxiv.org/abs/2009.10095)" ] }, { "cell_type": "code", "execution_count": 1, "id": "engaging-agreement", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import copy\n", "\n", "# Problem modelling imports\n", "from docplex.mp.model import Model\n", "\n", "# Qiskit imports\n", "from qiskit_algorithms import QAOA, NumPyMinimumEigensolver\n", "from qiskit_algorithms.optimizers import COBYLA\n", "from qiskit_algorithms.utils.algorithm_globals import algorithm_globals\n", "from qiskit.primitives import Sampler\n", "from qiskit_optimization.algorithms import MinimumEigenOptimizer, CplexOptimizer\n", "from qiskit_optimization import QuadraticProgram\n", "from qiskit_optimization.problems.variable import VarType\n", "from qiskit_optimization.converters.quadratic_program_to_qubo import QuadraticProgramToQubo\n", "from qiskit_optimization.translators import from_docplex_mp" ] }, { "cell_type": "markdown", "id": "regional-weekly", "metadata": {}, "source": [ "## Preliminaries: relaxing QUBOs\n", "\n", "First, we show how to relax a QUBO built with a semi-definite positive matrix to obtain an easy-to-solve QP." ] }, { "cell_type": "code", "execution_count": 2, "id": "southwest-stake", "metadata": {}, "outputs": [], "source": [ "def create_problem(mu: np.array, sigma: np.array, total: int = 3) -> QuadraticProgram:\n", " \"\"\"Solve the quadratic program using docplex.\"\"\"\n", "\n", " mdl = Model()\n", " x = [mdl.binary_var(\"x%s\" % i) for i in range(len(sigma))]\n", "\n", " objective = mdl.sum([mu[i] * x[i] for i in range(len(mu))])\n", " objective -= 2 * mdl.sum(\n", " [sigma[i, j] * x[i] * x[j] for i in range(len(mu)) for j in range(len(mu))]\n", " )\n", " mdl.maximize(objective)\n", " cost = mdl.sum(x)\n", " mdl.add_constraint(cost == total)\n", "\n", " qp = from_docplex_mp(mdl)\n", " return qp\n", "\n", "\n", "def relax_problem(problem) -> QuadraticProgram:\n", " \"\"\"Change all variables to continuous.\"\"\"\n", " relaxed_problem = copy.deepcopy(problem)\n", " for variable in relaxed_problem.variables:\n", " variable.vartype = VarType.CONTINUOUS\n", "\n", " return relaxed_problem" ] }, { "cell_type": "markdown", "id": "medieval-product", "metadata": {}, "source": [ "For this example, we use a positive semi-definite matrix $\\Sigma$ and a linear term $\\mu$ as defined below." ] }, { "cell_type": "code", "execution_count": 3, "id": "laden-number", "metadata": {}, "outputs": [], "source": [ "mu = np.array([3.418, 2.0913, 6.2415, 4.4436, 10.892, 3.4051])\n", "sigma = np.array(\n", " [\n", " [1.07978412, 0.00768914, 0.11227606, -0.06842969, -0.01016793, -0.00839765],\n", " [0.00768914, 0.10922887, -0.03043424, -0.0020045, 0.00670929, 0.0147937],\n", " [0.11227606, -0.03043424, 0.985353, 0.02307313, -0.05249785, 0.00904119],\n", " [-0.06842969, -0.0020045, 0.02307313, 0.6043817, 0.03740115, -0.00945322],\n", " [-0.01016793, 0.00670929, -0.05249785, 0.03740115, 0.79839634, 0.07616951],\n", " [-0.00839765, 0.0147937, 0.00904119, -0.00945322, 0.07616951, 1.08464544],\n", " ]\n", ")" ] }, { "cell_type": "markdown", "id": "aging-fortune", "metadata": {}, "source": [ "Using DOCPLEX we build a model with binary variables." ] }, { "cell_type": "code", "execution_count": 4, "id": "supreme-wallace", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem name: docplex_model1\n", "\n", "Maximize\n", " -2.15956824*x0^2 - 0.03075656*x0*x1 - 0.44910424*x0*x2 + 0.27371876*x0*x3\n", " + 0.04067172*x0*x4 + 0.0335906*x0*x5 - 0.21845774*x1^2 + 0.12173696*x1*x2\n", " + 0.008018*x1*x3 - 0.02683716*x1*x4 - 0.0591748*x1*x5 - 1.970706*x2^2\n", " - 0.09229252*x2*x3 + 0.2099914*x2*x4 - 0.03616476*x2*x5 - 1.2087634*x3^2\n", " - 0.1496046*x3*x4 + 0.03781288*x3*x5 - 1.59679268*x4^2 - 0.30467804*x4*x5\n", " - 2.16929088*x5^2 + 3.418*x0 + 2.0913*x1 + 6.2415*x2 + 4.4436*x3 + 10.892*x4\n", " + 3.4051*x5\n", "\n", "Subject to\n", " Linear constraints (1)\n", " x0 + x1 + x2 + x3 + x4 + x5 == 3 'c0'\n", "\n", " Binary variables (6)\n", " x0 x1 x2 x3 x4 x5\n", "\n" ] } ], "source": [ "qubo = create_problem(mu, sigma)\n", "print(qubo.prettyprint())" ] }, { "cell_type": "markdown", "id": "timely-sally", "metadata": {}, "source": [ "Such binary problems are hard to deal with but can be solved if the problem instance is small enough. Our example above has as solution" ] }, { "cell_type": "code", "execution_count": 5, "id": "contrary-bumper", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objective function value: 16.7689322\n", "variable values: x0=0.0, x1=0.0, x2=1.0, x3=1.0, x4=1.0, x5=0.0\n", "status: SUCCESS\n" ] } ], "source": [ "result = CplexOptimizer().solve(qubo)\n", "print(result.prettyprint())" ] }, { "cell_type": "markdown", "id": "established-stretch", "metadata": {}, "source": [ "We can create a relaxation of this problem in which the variables are no longer binary. Note that we use the `QuadraticProgramToQubo` converter to convert the constraint into a quadratic penalty term. We do this to remain consistent with the steps that the Qiskit optimization module applies internally. " ] }, { "cell_type": "code", "execution_count": 6, "id": "spectacular-african", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem name: docplex_model1\n", "\n", "Minimize\n", " 44.84880018*x0^2 + 85.40922044*x0*x1 + 85.82756812*x0*x2\n", " + 85.10474511999999*x0*x3 + 85.33779215999999*x0*x4 + 85.34487328*x0*x5\n", " + 42.90768968*x1^2 + 85.25672692*x1*x2 + 85.37044588*x1*x3 + 85.40530104*x1*x4\n", " + 85.43763867999999*x1*x5 + 44.65993794*x2^2 + 85.4707564*x2*x3\n", " + 85.16847247999999*x2*x4 + 85.41462863999999*x2*x5 + 43.89799534*x3^2\n", " + 85.52806848*x3*x4 + 85.34065100000001*x3*x5 + 44.28602462*x4^2\n", " + 85.68314192*x4*x5 + 44.85852282*x5^2 - 259.55339164*x0\n", " - 258.22669163999996*x1 - 262.37689163999994*x2 - 260.57899163999997*x3\n", " - 267.02739163999996*x4 - 259.54049163999997*x5 + 384.20308746\n", "\n", "Subject to\n", " No constraints\n", "\n", " Continuous variables (6)\n", " 0 <= x0 <= 1\n", " 0 <= x1 <= 1\n", " 0 <= x2 <= 1\n", " 0 <= x3 <= 1\n", " 0 <= x4 <= 1\n", " 0 <= x5 <= 1\n", "\n" ] } ], "source": [ "qp = relax_problem(QuadraticProgramToQubo().convert(qubo))\n", "print(qp.prettyprint())" ] }, { "cell_type": "markdown", "id": "absolute-creativity", "metadata": {}, "source": [ "The solution of this continuous relaxation is different from the solution to the binary problem but can be used to warm-start a solver when dealing with the binary problem." ] }, { "cell_type": "code", "execution_count": 7, "id": "moderate-photograph", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objective function value: -17.01205502568274\n", "variable values: x0=0.17524995761801201, x1=1.4803888163984595e-07, x2=0.9709053264087679, x3=0.7384168677494151, x4=0.9999999916475085, x5=0.14438904470168346\n", "status: SUCCESS\n" ] } ], "source": [ "sol = CplexOptimizer().solve(qp)\n", "print(sol.prettyprint())" ] }, { "cell_type": "code", "execution_count": 8, "id": "smoking-discretion", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.17524995761801201, 1.4803888163984595e-07, 0.9709053264087679, 0.7384168677494151, 0.9999999916475085, 0.14438904470168346]\n" ] } ], "source": [ "c_stars = sol.samples[0].x\n", "print(c_stars)" ] }, { "cell_type": "markdown", "id": "portuguese-profit", "metadata": {}, "source": [ "## QAOA\n", "\n", "Here, we illustrate how to warm-start the quantum approximate optimization algorithm (QAOA) by leveraging the relaxed problem shown above. \n", "\n", "### Standard QAOA\n", "\n", "First, we use standard QAOA to solve the QUBO. To do this, we convert the QUBO to `QuadraticProgram` class (note that the resulting problem is still a binary problem)." ] }, { "cell_type": "code", "execution_count": 9, "id": "recreational-packing", "metadata": {}, "outputs": [], "source": [ "algorithm_globals.random_seed = 12345\n", "qaoa_mes = QAOA(sampler=Sampler(), optimizer=COBYLA(), initial_point=[0.0, 1.0])\n", "exact_mes = NumPyMinimumEigensolver()" ] }, { "cell_type": "code", "execution_count": 10, "id": "pursuant-pendant", "metadata": {}, "outputs": [], "source": [ "qaoa = MinimumEigenOptimizer(qaoa_mes)" ] }, { "cell_type": "code", "execution_count": 11, "id": "painful-packing", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objective function value: 16.768932200000002\n", "variable values: x0=0.0, x1=0.0, x2=1.0, x3=1.0, x4=1.0, x5=0.0\n", "status: SUCCESS\n" ] } ], "source": [ "qaoa_result = qaoa.solve(qubo)\n", "print(qaoa_result.prettyprint())" ] }, { "cell_type": "markdown", "id": "arctic-wealth", "metadata": {}, "source": [ "### Warm-start QAOA\n", "\n", "Next, we compare this result to a warm-start QAOA in which we use the solution to the continuous relaxation of the problem. First, we create the initial state \n", "\n", "\\begin{align}\n", "|\\phi^*\\rangle=\\bigotimes_{i=0}^{n-1}R_y(\\theta_i)|0\\rangle_n .\n", "\\end{align}\n", "\n", "which is given by applying $R_y$ rotations with an angle $\\theta=2\\arcsin(\\sqrt{c^*_i})$ that depends on the solution to the relaxed problem.\n", "Here, $c^*_i$ the value of variable $i$ of the relaxed problem." ] }, { "cell_type": "code", "execution_count": 26, "id": "controversial-model", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAGwCAYAAAAJwO/qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAviElEQVR4nO3de1xVdb7/8ddGxM01RLCNolwCAlEgVBTHsUiaRMQ6FqkHOekx7SJp6khTjV1OYw5KdbxMo4426XFiKC0nJUctZCK6QYhykJ8KibqBnW0BE5QRZP/+YNrFAY3Lhu1afJ6Ph3/s7/qutT5r+34svuuy19KYTCYTQiiUjbULEKI7JMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0WytXYBon8kEzY3WrqLjbPqDRtP765UA36SaG+HwemtX0XHRi6GfXe+vV4YQQtEkwELRJMBC0STAQtEkwELRJMBC0eQ0moocLcvm15uiW7Vp7Rzx8ggkJiKJ+3/xJP36qeu/XF1bIwCIDp9NZNBUTJiouWTg0Nc72LR3GWfPl7D0wS3WLs+iJMAqFDA0gpjRc8yf4yc8wfw1Qez/aivzpqzC1cnDitVZloyB+wB7O0eCvMdjMpmovFBm7XIsSgLcR1T9K7guDm5WrsSyZAihQg2Nl7lYb8RkahkD7/18E6UVRwgaFomXR6C1y7OoPrEHNhqNpKSk4O/vj1arZdiwYSxZsoT6+nrmz5+PRqNh48aN1i7TYnYcfIEHX/Qg4aXBLHwtlL2fv8HEkTN4ae7frF2axal+D1xYWEhsbCwGgwFHR0dGjBhBZWUl69evp6ysjOrqagDCw8OtW6gFxY1byKTQBJqaGzldVURGdirGi3rs+mvNfVbtnEWzqZmVSe+Y276/XM2CtBAWTktjckSiNUrvNFXvgY1GI/Hx8RgMBpYvX05VVRUFBQUYDAZSU1PJzMwkLy8PjUZDaGiotcu1mKHuAUQExhAZFMvM6BRenreXE/o81u1+zNznyRlvUFyeS9aRdHPbhvcXEeI7UTHhBZUHePHixej1epKTk0lLS8PZ2dk8LSUlhbCwMJqamvDx8cHFxcWKlfasEJ8JxEQkkX00g+Lyz4CWg7nlCdvYuCcZ48VKPjm2i2Nl2Tw1Y5OVq+0c1Qa4pKSEjIwM3N3dWb16dbt9Ro8eDUBYWFir9tOnTzN9+nScnZ0ZOHAg//Ef/8GFCxd6vOaelBizEhubfmw/8Ly5bWzQFO4MfYjU9DlseO8JliVsxcVxkBWr7DzVBjg9PZ3m5mYSExNxcnJqt4+9vT3QOsCXLl0iOjoavV5Peno6W7ZsIScnh2nTptHc3NwrtfeEoe7+RIfN4kjpxxR9k2NuXxifRsWFUsYGxTIuOM6KFXaNagOclZUFQHR09HX76PV6oHWAt2zZQkVFBXv27GHatGkkJCTw9ttv88UXX/DBBx/0bNE9bPbk57DR2LD94I97YXs7Rzzd/PDVjbJiZV2n2rMQZ86cAcDb27vd6U1NTeTm5gKtA7xv3z4mTpzI8OHDzW1RUVH4+fmxd+9e7r///k7XMmbMGAwGQ6fmsbO1Z0vyqU7NE3bbXRxaa7rudO9bgzmw5lqnltlRAYEBXG260uX5dTod+fn5nZ5PtQGur68H4MqV9r/UjIwMjEYjzs7O+Pr6mtuPHz9OQkJCm/4hISEcP368S7UYDAYqKio6NY+2v0OX1mUtVZWVNDRe7vX1qjbAOp2OmpoaCgoKiIqKajWtqqqKFStWABAaGormJ78Hr6mpwdXVtc3y3NzcOHHiRJdr6Sw7W/surctaPIcM6fYeuCtUG+CYmBhKSkpITU3lnnvuITCw5RJqXl4eSUlJGI1GoHcuYHTlT+O1q733s/pXH8/u9jJOnTwlP6u3pJSUFAYNGsS5c+cICQlh1KhRBAQEEBkZiZ+fH3fffTfQ9hTawIEDqa2tbbO86upq3NzUdSOMGqg2wF5eXuTk5BAXF4dWq6W8vBw3Nzc2b95MZmYmJ0+eBNoGODg4uN2x7vHjxwkODu6V2kXHqXYIAS1h3LdvX5v2uro6ysvLsbGxYeTIka2mTZs2jWeffRa9Xo+XlxcAX375JWVlZaxdu7ZX6hYdp9o98I0UFxdjMpkICAjAwaH10f7ChQvx9PTkvvvuY9++fezatYvZs2cTGRnJfffdZ6WKxfX0yQAXFRUBbYcPAC4uLmRlZeHp6cmsWbN45JFHmDBhAvv27cPGpk9+XTc1VQ8hrudGAQa47bbb2h163Ez0351ibcbDXKw34qi9hRUz38JHF9KqT3NzM1v2/Zq8E3+nn40tLo6DWPrgnxjq7g/A+ZqzbHh/EXrjSWw0/YiPepz7Jz7Zahlr/jqXQ19v5/3/qsHJ3rW3Nq/DJMAKtW73o0wdt5B7x87lk2O7WJsxlz8syWvV5/PjH1BcnsvmZUex7defv3z0O97c/ywrk97BZDLx4vZ/Y2b0b7gzrOXCTc2lb1vNn1P0Hrb9+vfaNnVFn/ybmJWVhclkIi5OeTevANTUneekPp+YiJZfHv9y1AN8V3uOCmNpq34aNFxt+idXGxswmUxcbvgej1taDkyPnPqY/rYDzOEFGOh864/ruPQt6Vmv8Fj8a72wRV3XJ/fASvdd7TncXDzNDynRaDQMHjic87VnzcMDgPEj4iksO8zM/9JhP8AZ91uG8urj/wDgzPnj3OLowaqdszj33Ql0A314NP5VPAf5AfDargUsiFuDg9a5bQE3kT65B+4rTurzKTf8L+krK/jrykru8J9s/lXGtWtNFJZlkRizkk1LjzD69nt5eedDAHz45VYGuw7nDv+7rVl+h0iAFcjDdRjV31dx7VoTACaTifM1ZxnsOrxVv0Nf7yDc/26c7F2xsbHhnjEPU1h2GIDBA4fjP+QO84FfzOgkSisKaLrWyNGyw3xe/DfmvOLDnFd8AFj4WiilFUd6byM7SAKsQAOdBuM/NIKPCnYCkFO0G3dXr1bDBwBPNz8KS7NobLoKwBcl+/DRtVy4GRsUy3cX9Rgvttwl91XJhwwfHIxtv/488+9/4e3fnmPns+XsfLYcgC3LjuE/9I5e2sKOkzGwQj31wGbWZswlPesVHLQurHjozwC8+u4jRI2YzoSQ6Uz/xSLOni/h0dfDsLXpz0BnHU890PKbN3s7R5bM2MRz2+IAE47aW3gu8a9W3KKu0ZhMpuvfAS2spjfvRrMEecmLEF0gARaKJgEWiiYBFoomARaKJmchblLyruSOkQALRZMhhFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDR5tNRNQmm/gfs5vfUbOQnwTaK5UVmPkvo5vfWoKRlCCEWTAAtFkwALRZMAC0WTAAtFkwALRZMAC0WT88AqcrQsm19vim7VprVzxMsjkJiIJO7/xZPmd8uphbq2RgAQHT6byKCpmDBRc8nAoa93sGnvMs6eL2Hpg1usXZ5FSYBVKGBoBDGj55g/x094gvlrgtj/1VbmTVmFq5OHFauzLBkD9wH2do4EeY/HZDJReaHM2uVYlAS4j6j6V3BdHNysXIll9YkAG41GUlJS8Pf3R6vVMmzYMJYsWUJ9fT3z589Ho9GwceNGa5dpMQ2Nl7lYb6S27jtOVxWx/r1FlFYcIWhYJF4egdYuz6JUPwYuLCwkNjYWg8GAo6MjI0aMoLKykvXr11NWVkZ1dTUA4eHh1i3UgnYcfIEdB19o1TZx5Aye/Lc/WKminqPqPbDRaCQ+Ph6DwcDy5cupqqqioKAAg8FAamoqmZmZ5OXlodFoCA0NtXa5FhM3biGpCw6xav6HPDI1FWcHN4wX9dj115r7rNo5i5f/56FW831/uZqZ/+XJxwV/6e2Su0zVAV68eDF6vZ7k5GTS0tJwdnY2T0tJSSEsLIympiZ8fHxwcXGxYqWWNdQ9gIjAGCKDYpkZncLL8/ZyQp/Hut2Pmfs8OeMNistzyTqSbm7b8P4iQnwnMjki0Rpld4lqA1xSUkJGRgbu7u6sXr263T6jR48GICwszNz2Q+AjIyMZMGAAGmu8esfCQnwmEBORRPbRDIrLPwNaDuaWJ2xj455kjBcr+eTYLo6VZfPUjE1WrrZzVBvg9PR0mpubSUxMxMnJqd0+9vb2QOsAl5aWsnv3bnQ6HWPHju2VWntDYsxKbGz6sf3A8+a2sUFTuDP0IVLT57DhvSdYlrAVF8dBVqyy81Qb4KysLACio6Ov20ev1wOtAzxp0iSqqqr44IMPiImJ6dkie9FQd3+iw2ZxpPRjir7JMbcvjE+j4kIpY4NiGRccZ8UKu0a1AT5z5gwA3t7e7U5vamoiNzcXaB1gGxvVfiXMnvwcNhobth/8cS9sb+eIp5sfvrpRVqys61R7Gq2+vh6AK1eutDs9IyMDo9GIs7Mzvr6+PVrLmDFjMBgMN+xjZ2vPluRT3VpP2G13cWjt9d9b6X1rMAfWXOvWOjoqIDCAq03tf/ft0el05Ofnd3o9qg2wTqejpqaGgoICoqKiWk2rqqpixYoVAISGhvb4gZrBYKCiouKGfbT9HXq0ht5WVVlJQ+PlHl+PagMcExNDSUkJqamp3HPPPQQGtlyBysvLIykpCaPRCPTOBQydTvezfexs7Xu8jt7kOWRIp/fAXaHaAKekpPD2229z7tw5QkJCCAoKoqGhgdLSUmJjY/Hx8eHAgQOtxr89pSN/Gq9dtd5zIV59PNviyzx18pQ8F6I7vLy8yMnJIS4uDq1WS3l5OW5ubmzevJnMzExOnjwJ0CsBFj1HtXtggODgYPbt29emva6ujvLycmxsbBg5cqQVKhOWouoAX09xcTEmk4nAwEAcHNoePO3atQuA48ePt/rs4+PDmDFjeq9Q8bP6ZICLioqA6w8fEhIS2v388MMP89Zbb/VobaJzJMDtMJmufy61N+m/O8XajIe5WG/EUXsLK2a+hY8upE2//V9t46+Hf4+puZlw/7tZPOMNbPv17/K0v+f9mfdz1pmXb7yoZ5TfJF58+D3yThxga+bT5mm19edxc9bxx6cKevjbaJ9qD+Ju5OcCfLNYt/tRpo5byFtPn2Rm9NOszZjbpk9V9WneOrCS1x/PYftvSqmp+5bML7Z0a9qUsfPYvKzQ/M/NWcfkO1ruUBt7+72tpgUMjeDuO6x391qfDHBWVhYmk4m4uJv32n9N3XlO6vOJiWj5ceYvRz3Ad7XnqDCWtuqXc2wXUSOm4+aiQ6PRMG38YxwuTO/WtJ8qOfsltXXniQqZ3maa8WIlR059TMzoJEtvfof1yQArwXe153Bz8TQ/x0Gj0TB44HDO155t1e987VluHfjj/R46Nx9zn65O+6m/f7WNyaOTzMOOnzqY/xaRQVMZ6DS4G1vaPRJgcV1XrtaTXfhXYiPnt5lmMpk4kPcmU9qZ1pskwDcpD9dhVH9fxbVrTUBLYM7XnGWw6/BW/Qa7DufbmjPmz4bqcnOfrk77wSdH38X71hC8bx3Rpr5j3/yDq00NjLn93m5uafdIgG9SA50G4z80go8KdgKQU7Qbd1cvhrr7t+r3y1EP8PnxD6j+3oDJZGLfF5u4K3xWt6b94O952667h93/1TZ+NWYu/Wz6WXrTO6VPnkZTiqce2MzajLmkZ72Cg9aFFQ/9GYBX332EqBHTmRAyHc9Bfjz8q5d46g+/AFpuqZw2/lGALk8DOHf+BGWVhaz6zw/b1FV/5SK5Re+xZXlRz218B2lMN8tJzz7Omjfz9AR5yYsQHSABFoomARaKJgEWiiYBFoomARaKJqfRbhLysu+ukQALRZMhhFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQJsFA0CbBQNHm01E1KaT8x6q2fEP1fEuCbVHOjsh411VuPkvq/ZAghFE0CLBRNAiwUTQIsFE0CLBRNAiwUTQIsFE3OA6vI0bJsfr0pulWb1s4RL49AYiKSuP8XT5rfO6cW6toaAUB0+Gwig6ZiwkTNJQOHvt7Bpr3LOHu+hKUPbrF2eRYlAVahgKERxIyeY/4cP+EJ5q8JYv9XW5k3ZRWuTh5WrM6yZAzcB9jbORLkPR6TyUTlhTJrl2NREuA+oupfwXVxcLNyJZYlQwgVami8zMV6IyZTyxh47+ebKK04QtCwSLw8Aq1dnkX1iT2w0WgkJSUFf39/tFotw4YNY8mSJdTX1zN//nw0Gg0bN260dpkWs+PgCzz4ogcJLw1m4Wuh7P38DSaOnMFLc/9m7dIsTvV74MLCQmJjYzEYDDg6OjJixAgqKytZv349ZWVlVFdXAxAeHm7dQi0obtxCJoUm0NTcyOmqIjKyUzFe1GPXX2vus2rnLJpNzaxMesfc9v3lahakhbBwWhqTIxKtUXqnqXoPbDQaiY+Px2AwsHz5cqqqqigoKMBgMJCamkpmZiZ5eXloNBpCQ0OtXa7FDHUPICIwhsigWGZGp/DyvL2c0Oexbvdj5j5PzniD4vJcso6km9s2vL+IEN+JigkvqDzAixcvRq/Xk5ycTFpaGs7OzuZpKSkphIWF0dTUhI+PDy4uLlastGeF+EwgJiKJ7KMZFJd/BrQczC1P2MbGPckYL1byybFdHCvL5qkZm6xcbeeoNsAlJSVkZGTg7u7O6tWr2+0zevRoAMLCwsxtu3bt4oEHHsDb2xsHBweCgoJ47rnnqKur65W6e0pizEpsbPqx/cDz5raxQVO4M/QhUtPnsOG9J1iWsBUXx0FWrLLzVBvg9PR0mpubSUxMxMnJqd0+9vb2QOsAp6Wl0a9fP1555RX279/P448/zh//+EemTJlCc3Nzr9TeE4a6+xMdNosjpR9T9E2OuX1hfBoVF0oZGxTLuOA4K1bYNao9iMvKygIgOjr6un30ej3QOsB79+7Fw+PHK1V33nknHh4eJCYm8umnnzJp0qQeqrjnzZ78HIcL09l+8HnSHjsMtFzk8HTzw1c3ysrVdY1qA3zmzBkAvL29253e1NREbm4u0DrAPw3vD8aMGQNARUVFl2oZM2YMBoOhU/PY2dqzJflUp+YJu+0uDq29/mv/vG8N5sCaa51aZkcFBAZwtelKl+fX6XTk5+d3ej7VBri+vh6AK1fa/1IzMjIwGo04Ozvj6+t7w2UdPtyytwoODu5SLQaDodPh1/Z36NK6rKWqspKGxsu9vl7VBlin01FTU0NBQQFRUVGtplVVVbFixQoAQkND0dzggQYVFRWsXLmSKVOmdPlcsU6n6/Q8drb2XVqXtXgOGdLtPXBXqDbAMTExlJSUkJqayj333ENgYMsl1Ly8PJKSkjAajcCNL2DU1dVx3333YWdnx5tvvtnlWrryp/Ha1d57LsSrj2d3exmnTp6S50JYUkpKCoMGDeLcuXOEhIQwatQoAgICiIyMxM/Pj7vvvhtoPf79qStXrhAfH8/p06c5ePAgnp6evVm+6CDVBtjLy4ucnBzi4uLQarWUl5fj5ubG5s2byczM5OTJk0D7AW5sbOTBBx8kPz+f/fv3M2LEiN4uX3RQn3xbfV1dHS4uLmg0Gi5duoSDw48HTM3NzcyaNYsPPviADz/80Lyn7m29OYSwBGs9Wkq1Y+AbKS4uxmQyERgY2Cq8AIsWLeLdd9/lN7/5DQ4ODnzxxRfmabfddlu7p9mE9ah2CHEjRUVFQPvDh/379wPw+9//nqioqFb/MjMze7VO8fP65B74RgEuLy/v5Wos62pjA6v+Mosz3x5nQH97XJ0Gs3jGHxnq7t+m718Pp3Iofzu2/eyw669l0X3rCRoeaYWqu04CrEJTxy0kMigWjUbDntyNvPbuI21OlZVWFLL3szfY+uti7Ac48dHXO9m4J5mNi7+yTtFd1CcD/MN9Empk11/LuOCp5s/Bw8ez6x9pbfppNBqamhtpuFqP/QAn6hpqcb/FqzdLtYg+GeC+5P1P1xEVcl+b9tuGhPHAL5eStNoXZwc3+vcbwGtPfGKFCrunTx7E9RVvf/wKlcZS5se2vR+6qvo0nxa9x1tPl5L+Wz0zJi3ldztnWqHK7pEAq9S72Wl8+r/v8coj+9Hatb0x6NNju/H1HIX7LUMAuHfsPIrLc2lsutrbpXaLBFiFdv3jNQ4XppO64BBO9q7t9tEN8qO4PJcr/2z5pcmXx/fh5RFIf1srXI3oBhkDq8x3tXo271uOp5uf+UF/drYD2LD4S9468DyDXIYQH/UYE0f+GyfP5bFo3Rj62w5Aa+fIM//+tpWr77w+eSlZCeRScsfIEEIomgRYKJoEWCiaBFgomgRYKJqchbhJycu+O0YCLBRNhhBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELRJMBC0STAQtEkwELR5NFSNyn5TVzHSIBvUs2N8mipjpAhhFA0CbBQNAmwUDQJsFA0CbBQNAmwUDQ5jaYiR8uyza8V+IHWzhEvj0BiIpK4/xdP0q+fuv7L1bU1AoDo8NlEBk3FhImaSwYOfb2DTXuXcfZ8CUsf3GLt8ixKAqxCAUMjiBk9x/w5fsITzF8TxP6vtjJvyipcnTysWJ1lyRi4D7C3cyTIezwmk4nKC2XWLseiJMB9RNW/guvi4GblSixLhhAq1NB4mYv1RkymljHw3s83UVpxhKBhkXh5BFq7PIvqE3tgo9FISkoK/v7+aLVahg0bxpIlS6ivr2f+/PloNBo2btxo7TItZsfBF3jwRQ8SXhrMwtdC2fv5G0wcOYOX5v7N2qVZnOr3wIWFhcTGxmIwGHB0dGTEiBFUVlayfv16ysrKqK6uBiA8PNy6hVpQ3LiFTApNoKm5kdNVRWRkp2K8qMeuv9bcZ9XOWTSbmlmZ9I657fvL1SxIC2HhtDQmRyRao/ROU/Ue2Gg0Eh8fj8FgYPny5VRVVVFQUIDBYCA1NZXMzEzy8vLQaDSEhoZau1yLGeoeQERgDJFBscyMTuHleXs5oc9j3e7HzH2enPEGxeW5ZB1JN7dteH8RIb4TFRNeUHmAFy9ejF6vJzk5mbS0NJydnc3TUlJSCAsLo6mpCR8fH1xcXKxYac8K8ZlATEQS2UczKC7/DGg5mFuesI2Ne5IxXqzkk2O7OFaWzVMzNlm52s5RbYBLSkrIyMjA3d2d1atXt9tn9OjRAISFhZnbcnJyiImJwdPTkwEDBuDl5cXMmTMpKSnplbp7SmLMSmxs+rH9wPPmtrFBU7gz9CFS0+ew4b0nWJawFRfHQVassvNUG+D09HSam5tJTEzEycmp3T729vZA6wDX1NQwatQo1q9fz8GDB0lNTaW4uJioqCj0en2v1N4Thrr7Ex02iyOlH1P0TY65fWF8GhUXShkbFMu44DgrVtg1qj2Iy8rKAiA6Ovq6fX4I5E8DPH36dKZPn96q39ixY7n99tvZvXs3S5Ys6YFqe8fsyc9xuDCd7QefJ+2xw0DLRQ5PNz98daOsXF3XqDbAZ86cAcDb27vd6U1NTeTm5gKtA9yeQYNa/qza2t7cX1fYbXdxaO31X/vnfWswB9Zc68WKet7N/T/SDfX19QBcuXKl3ekZGRkYjUacnZ3x9fVtM/3atWs0Nzdz5swZnnnmGXQ6HQ899FCXahkzZgwGg6FT89jZ2rMl+VSX1mcNAYEBXG1q/7vuCJ1OR35+fqfnU22AdTodNTU1FBQUEBUV1WpaVVUVK1asACA0NBRNO78Hv/POO817aH9/f7KysvDw6NpNMAaDgYqKik7No+3v0KV1WUtVZSUNjZd7fb2qDXBMTAwlJSWkpqZyzz33EBjYcgk1Ly+PpKQkjEYjcP0LGNu2baO2tpbTp0+zdu1afvWrX5Gbm8vw4cM7XYtOp+v0PHa29p2ep6tefTy728vwHDKk23vgrlDtu5L1ej3h4eFcuHABW1tbgoKCaGhooLS0lNjYWJqbmzlw4ABbtmxhwYIFN1xWbW0tPj4+zJkzp9cuOV+7Ks+F6AjVnkbz8vIiJyeHuLg4tFot5eXluLm5sXnzZjIzMzl58iTw8wdwAK6urvj7+1NaWtrTZYtOUu0QAiA4OJh9+/a1aa+rq6O8vBwbGxtGjhz5s8s5f/48J06cYNy4cT1RpugGVQf4eoqLizGZTAQGBuLg0Ppgac6cOfj7+xMeHo6rqyunTp3i9ddfx9bWlqVLl1qpYnE9fTLARUVFQPvDh/Hjx7Njxw7WrVtHQ0MDw4YNIzo6mmefffa655SF9UiA/4/k5GSSk5N7u6RuudrYwKq/zOLMt8cZ0N8eV6fBLJ7xR4a6+7fp+8XxfWzZ92uuma7hqxvFiplv4ah1Ie/EAbZmPm3uV1t/HjdnHX98qqA3N6XTJMAqMXXcQiKDYtFoNOzJ3chr7z7S5vTYlX/W8eq783n18X8wfHAQG95P5i8fvczCaWsZe/u9jL39XnPf3745jbDbrn8Z/mah2rMQN5KVlYXJZCIuTnk3r7THrr+WccFTzRdkgoeP59ua8jb9vvp/+/EfcgfDBwcBMH3CExwuTG/Tz3ixkiOnPiZmdFKP1m0JfTLAavf+p+uICrmvTfv52rPcOvDHcfytA32o/r6Ka9eaWvU7mP8WkUFTGeg0uMdr7S4JsMq8/fErVBpLmR/b/j3QP8dkMnEg702mRM63cGU9QwKsIu9mp/Hp/77HK4/sR2vX9l6Kwa7D+bbmjPnztzXluLl4tnrc1LFv/sHVpgbG/GQ8fDOTAKvErn+8xuHCdFIXHMLJ3rXdPmNvn0JpRQFnz/8/AD747A3uCpvVqs/+r7bxqzFz6WfTr6dLtog+eRZCbb6r1bN533I83fzMD/ezsx3AhsVf8taB5xnkMoT4qMdw0DqzNGErL751P9eam/DRjSRl5nbzcuqvXCS36D22LC+y1qZ0mmpv5lE6uZmnY2QIIRRNAiwUTQIsFE0CLBRNAiwUTQIsFE1Oo92k5GXfHSMBFoomQwihaBJgoWgSYKFoEmChaBJgoWgSYKFoEmChaBJgoWgSYKFoEmChaBJgoWgSYKFoEmChaBJgoWgSYKFoEmChaBJgoWjyaKmblPykqGMkwDep5kZ5tFRHyBBCKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJqcB1aRo2XZ5ndk/EBr54iXRyAxEUnc/4snW72RSA3UtTUCgOjw2UQGTcWEiZpLBg59vYNNe5dx9nwJSx/cYu3yLEoCrEIBQyOIGT3H/Dl+whPMXxPE/q+2Mm/KKlydPKxYnWXJGLgPsLdzJMh7PCaTicoLZdYux6IkwH1E1b+C6+LgZuVKLEuGECrU0HiZi/VGTKaWMfDezzdRWnGEoGGReHkEWrs8i1L9HthoNJKSkoK/vz9arZZhw4axZMkS6uvrmT9/PhqNho0bN1q7TIvacfAFHnzRg4SXBrPwtVD2fv4GE0fO4KW5f7N2aRan6j1wYWEhsbGxGAwGHB0dGTFiBJWVlaxfv56ysjKqq6sBCA8Pt26hFhY3biGTQhNoam7kdFURGdmpGC/qseuvNfdZtXMWzaZmVia9Y277/nI1C9JCWDgtjckRidYovdNUuwc2Go3Ex8djMBhYvnw5VVVVFBQUYDAYSE1NJTMzk7y8PDQaDaGhodYu16KGugcQERhDZFAsM6NTeHneXk7o81i3+zFznydnvEFxeS5ZR9LNbRveX0SI70TFhBdUHODFixej1+tJTk4mLS0NZ2dn87SUlBTCwsJoamrCx8cHFxcXK1ba80J8JhATkUT20QyKyz8DWg7mlidsY+OeZIwXK/nk2C6OlWXz1IxNVq62c1QZ4JKSEjIyMnB3d2f16tXt9hk9ejQAYWFh111ObGwsGo2GF198sSfK7FWJMSuxsenH9gPPm9vGBk3hztCHSE2fw4b3nmBZwlZcHAdZscrOU2WA09PTaW5uJjExEScnp3b72NvbA9cP8DvvvENhYWFPldjrhrr7Ex02iyOlH1P0TY65fWF8GhUXShkbFMu44DgrVtg1qgxwVlYWANHR0dfto9frgfYD/P333/PUU0+RlpbWMwVayezJz2GjsWH7wR/3wvZ2jni6+eGrG2XFyrpOlWchzpw5A4C3t3e705uamsjNzQXaD/Bzzz1HYGAgiYmJzJkzp830zhozZgwGg6FT89jZ2rMl+VSn5gm77S4Orb3+a/+8bw3mwJprnVpmRwUEBnC16UqX59fpdOTn53d6PlUGuL6+HoArV9r/QjMyMjAajTg7O+Pr69tqWn5+Pn/605/4+uuvLVaPwWCgoqKiU/No+ztYbP29oaqykobGy72+XlUGWKfTUVNTQ0FBAVFRUa2mVVVVsWLFCgBCQ0PR/ORhBteuXePRRx8lOTmZkJAQi9bTWXa29hZbf2/wHDKk23vgrlBlgGNiYigpKSE1NZV77rmHwMCWy6d5eXkkJSVhNBqBthcwNm7cyLfffmvxsw5d+dN47WrvPRfi1cezu72MUydPyXMhLCUlJYVBgwZx7tw5QkJCGDVqFAEBAURGRuLn58fdd98NtB7/Go1GVq5cyfPPP09TUxO1tbXU1tYC0NDQQG1tLc3NzdbYHHEDqgywl5cXOTk5xMXFodVqKS8vx83Njc2bN5OZmcnJkyeB1gHW6/VcunSJRx99lIEDB5r/AaSmpjJw4EDOnj1rle0R19fn3lZfV1eHi4sLGo2GS5cu4eDgYG5v7099dHQ0Dz/8MHPnzmX8+PFotdo2fXpCbw4hLMFaj5ZS5Rj4RoqLizGZTAQGBprDC+Dk5MRdd93V7jw+Pj7XnSasS5VDiBspKioCbnwJWShHn9sDdzbAShlhPb3lV9RcMqDR2OCgdWbRfevxH3pHqz6G6nLWZsyltPIIuoG+bF5W2GY5JpOJlM2TOVVRwJ6Xa3un+G6QAKvEyqR3cLJ3BeDTovdZmzGXzcuOturjoHVh3pTfUd9wkTf3P9fucnZ/8jqeg27jVEVBT5dsEX1uCJGVlYXJZCIuTnk3rtzID+EFqG+4CLR92rSLgxsjfSeitXNsdxnlhmI+K97DrOjf9FCVltfn9sBqlpr+HxwtOwzAqvkfdmrepmuNvL5rAcsStmFj068nyusRfW4PrGZPz97B2789x9wpv+NPHz7dqXn/59BLTBw5A+9bg3uoup4hAVahX415mKOlh/m+/kKH5zn2zT/Yk7uBOa/4sPSNiVz+5/fMecWH2rrverDS7pMhhArUXaml4epl3G8ZAkDu/+7BxXEQzp14BsTrT/x4k7uhupzHXg9n57Plli7V4iTAKlDfcJGX/yeBfzZewUZjwy2OHrw8bx8ajYZX332EqBHTmRAynYarl5m3JpDGpn9S33CR2b/zIiYiiflT2//ZlRL0uUvJSiGXkjtGxsBC0STAQtEkwELRJMBC0STAQtHkLMRNSl723TESYKFoMoQQiiYBFoomARaKJgEWiiYBFoomARaKJgEWiiYBFoomARaKJgEWiiYBFoomARaKJgEWiiYBFoomARaKJgEWiiYBFoomT+a5SclPijpGAnyTam6UJ/N0hAwhhKJJgIWiSYCFokmAhaJJgIWiSYCFokmAhaLJeWAVOVqWza83Rbdq09o54uURSExEEvf/4kn69VPXf7m6tkYAEB0+m8igqZgwUXPJwKGvd7Bp7zLOni9h6YNbrF2eRUmAVShgaAQxo+eYP8dPeIL5a4LY/9VW5k1ZhauThxWrsywZA/cB9naOBHmPx2QyUXmhzNrlWJQEuI+o+ldwXTrx7jglkCGECjU0XuZivRGTqWUMvPfzTZRWHCFoWCReHoHWLs+i+sQe2Gg0kpKSgr+/P1qtlmHDhrFkyRLq6+uZP38+Go2GjRs3WrtMi9lx8AUefNGDhJcGs/C1UPZ+/gYTR87gpbl/s3ZpFqf6PXBhYSGxsbEYDAYcHR0ZMWIElZWVrF+/nrKyMqqrqwEIDw+3bqEWFDduIZNCE2hqbuR0VREZ2akYL+qx668191m1cxbNpmZWJr1jbvv+cjUL0kJYOC2NyRGJ1ii901S9BzYajcTHx2MwGFi+fDlVVVUUFBRgMBhITU0lMzOTvLw8NBoNoaGh1i7XYoa6BxARGENkUCwzo1N4ed5eTujzWLf7MXOfJ2e8QXF5LllH0s1tG95fRIjvRMWEF1Qe4MWLF6PX60lOTiYtLQ1nZ2fztJSUFMLCwmhqasLHxwcXFxcrVtqzQnwmEBORRPbRDIrLPwNaDuaWJ2xj455kjBcr+eTYLo6VZfPUjE1WrrZzVBvgkpISMjIycHd3Z/Xq9l9mPXr0aADCwsLMbdnZ2Wg0mjb/lD7ESIxZiY1NP7YfeN7cNjZoCneGPkRq+hw2vPcEyxK24uI4yIpVdp5qx8Dp6ek0NzeTmJiIk5NTu33s7e2B1gH+wR/+8AciIiLMnx0dHXum0F4y1N2f6LBZfHzkLxR9k8Mov18CsDA+jflrgxkbFMu44DgrV9l5qt0DZ2VlARAdHX3dPnq9Hmg/wCNGjGD8+PHmf6NGjeqZQnvR7MnPYaOxYfvBH/fC9naOeLr54atT5vapdg985swZALy9vdud3tTURG5uLtB+gC1pzJgxGAyGTs1jZ2vPluRTnZon7La7OLT2+q/98741mANrrnVqmR0VEBjA1aYrXZ5fp9ORn5/f6flUG+D6+noArlxp/0vNyMjAaDTi7OyMr69vm+kzZ87EaDQyaNAgpk+fzu9//3vc3d27VIvBYKCioqJT82j7O3RpXdZSVVlJQ+PlXl+vagOs0+moqamhoKCAqKioVtOqqqpYsWIFAKGhoWh+8kCDW265hRUrVjBp0iScnJz4/PPPWb16NV988QX5+flotVo6S6fTdXoeO1v7Ts9jTZ5DhnR7D9wVqn3V7OLFi9mwYQPDhg3jo48+IjCw5RJqXl4eSUlJfPPNNzQ2NrJo0aKfvQq3d+9epk+fzptvvsm8efN6o3yuXZXnQnSEag/iUlJSGDRoEOfOnSMkJIRRo0YREBBAZGQkfn5+3H333UDHxr/Tpk3D0dGxS2M00bNUG2AvLy9ycnKIi4tDq9VSXl6Om5sbmzdvJjMzk5MnTwKdO4DTWOPZSeKGVDsGBggODmbfvn1t2uvq6igvL8fGxoaRI0f+7HI++OAD6uvriYyM7IkyRTeoOsDXU1xcjMlkIjAwEAeH1kf7c+bMwc/Pj4iICPNB3Jo1awgPD2fWrFlWqlhcT58McFFREdD+8CEkJIS3336b//7v/+bKlSt4eXmxYMECXnjhBezsrHCUIm5IAvx/PPPMMzzzzDO9XVK36b87xdqMh7lYb8RRewsrZr6Fjy6kVZ+/5/2Z93PWmT8bL+oZ5TeJFx9+D4C/Hk7lUP52bPvZYddfy6L71hM0/OYeNkmAVWLd7keZOm4h946dyyfHdrE2Yy5/WJLXqs+UsfOYMvbH04AL0kYy+Y6WWydLKwrZ+9kbbP11MfYDnPjo651s3JPMxsVf9ep2dJZqz0LcSFZWFiaTibg45d280p6auvOc1OcTE9HyS+RfjnqA72rPUWEsve48JWe/pLbuPFEh04GWMyxNzY00XG25glnXUIv7LV49X3w39ck9sNp8V3sONxdP80NLNBoNgwcO53ztWYa6+7c7z9+/2sbk0UnY9usPwG1Dwnjgl0tJWu2Ls4Mb/fsN4LUnPum1beiqPrkH7uuuXK0nu/CvxEbON7dVVZ/m06L3eOvpUtJ/q2fGpKX8budMK1bZMRJgFfBwHUb191Vcu9YEgMlk4nzNWQa7Dm+3/ydH38X71hC8bx1hbvv02G58PUfhfssQAO4dO4/i8lwam672/AZ0gwRYBQY6DcZ/aAQfFewEIKdoN+6uXtcfPuRtY8pP9r4AukF+FJfncuWfdQB8eXwfXh6B9Le9uU8dyhhYJZ56YDNrM+aSnvUKDloXVjz0ZwBeffcRokZMZ8K/DtbOnT9BWWUhq/7zw1bzTxz5b5w8l8eidWPobzsArZ0jz/z7272+HZ2l2rvRlE7uRusYGUIIRZMAC0WTAAtFkwALRZMAC0WTsxA3KXnZd8dIgIWiyRBCKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKJoEWCiaBFgomgRYKNr/B9jzorFY4ulNAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit import QuantumCircuit\n", "\n", "thetas = [2 * np.arcsin(np.sqrt(c_star)) for c_star in c_stars]\n", "\n", "init_qc = QuantumCircuit(len(sigma))\n", "for idx, theta in enumerate(thetas):\n", " init_qc.ry(theta, idx)\n", "\n", "init_qc.draw(output=\"mpl\", style=\"clifford\")" ] }, { "cell_type": "markdown", "id": "tropical-armstrong", "metadata": {}, "source": [ "Next, we create the mixer operator for QAOA. When warm-starting QAOA we must ensure that the mixer operator has the initial state as ground state. We therefore chose the Hamiltonian\n", "\n", "\\begin{align}\n", "H_{M,i}^{(ws)}=\n", "\\begin{pmatrix}\n", "2c_i^*-1 & -2\\sqrt{c_i^*(1-c_i^*)} \\\\\n", "-2\\sqrt{c_i^*(1-c_i^*)} & 1-2c_i^*\n", "\\end{pmatrix}\n", "\\end{align}\n", "\n", "as mixer operator for qubit $i$. Once multiplied by $-i\\beta$ and exponentiated this mixer produces the following mixer circuit." ] }, { "cell_type": "code", "execution_count": 27, "id": "pacific-destiny", "metadata": { "tags": [ "nbsphinx-thumbnail" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAGwCAYAAAC0BO62AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNz0lEQVR4nO3deVxU9f7H8dewCaKEiIqKiqQEgmAuGKYWiSbuUaalpP5KWkRaVLrVNe2adlGqq1guV3P5VVxyvVfNJbciMsXcuMpPQkNlmXQENFACZH5/eJvkCsYyzPGc+TwfDx4P52x8vt/vOW8PZ86c0RmNRiNCCCFUy0bpAoQQQtSPBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicBLkQQqicndIFiLuP0QgVZUpXUTc29qDTKV1F/alxDMzd99IHNSdBLm5TUQb7FildRd2ExoCtg9JV1J8ax8DcfS99UHNyaUUIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVRObj8UZnP8zH6mLw2tNM3RwRnPFj6EdY9k1INTsbWVXa4hWfsYWGv7tdciobjQbk8R7DsEI0YKftHz1Q9rWbrlNc5fTOfVJ5YrXZ5VsPYxsLb2S5ALs+vctjthPcabXg/v8xLPzvdl+6EVTBo8F9cmLRSszjpY+xhYW/vlGrlocE4Ozvh2eACj0Uju5TNKl2OVrH0MtN5+CXJhEXn/OXhcGrspXIn1svYx0HL75dKKMLuSsmtcKTZgNN68PrnlwFIyc47i2y4YzxY+SpdnFax9DKyt/VYR5AaDgfnz57Nx40ays7Np0aIFERERzJs3j5iYGD755BMSEhKIjo5WulRNWLtrFmt3zao0rW9ABFMf+0ihiqyPtY+BtbVf80F+7NgxwsPD0ev1ODs706VLF3Jzc1m0aBFnzpwhPz8fgG7duilbqIYM7R1F/8DRlFeU8VNeGkn74zBcycbB3tG0zNxPx1JhrGBm5BemaVev5TM53p+oYfEM6D5OidI1oyZjkHY2mTdXht+2bvmNUioqbrBz/g1LlmxW1rYPavoaucFgYPjw4ej1eqZNm0ZeXh5HjhxBr9cTFxfHtm3bSE1NRafTERgYqHS5mtHWvTPdfcII9g1nTGgscyZt4XR2Kgs3vGBaZmrEx5zMSmHv0UTTtIRNU/Dv2FdVB9DdqiZj0NW7H1vmFlX6WRWbgYuzOxMenaNg9fVnbfugpoM8JiaG7OxsoqOjiY+Pp2nTpqZ5sbGxBAUFUV5ejpeXFy4uLgpWqm3+Xn0I6x7J/uNJnMz6Drj5htO00StZvDkaw5VcvjmxnhNn9vNKxFKFq9Wmqsbgv5WW/8o7ayMI8OrL0wPetHCFDUvr+6Bmgzw9PZ2kpCTc3d157733qlymR48eAAQFBVWa/tNPPzFixAiaNm1Ks2bNeOaZZ7h8+XKD16xl48JmYmNjy5qdb5um9fIdzEOBTxKXOJ6EjS/x2ugVuDg3V7BKbatqDG61cMMLlJaVMGPMassWZiFa3gc1G+SJiYlUVFQwbtw4mjRpUuUyTk5OQOUg/+WXXwgNDSU7O5vExESWL19OcnIyw4YNo6KiwiK1a1Fb906EBo3laOYe0s4mm6ZHDY8n53ImvXzD6e03VMEKta+6MQDY9O0iDqZv5Z2Jm3F0aKxQhQ1Ly/ugZoN87969AISGhla7THZ2NlA5yJcvX05OTg6bN29m2LBhjB49ms8//5zvv/+ef/3rXw1btMY9NeAtbHQ2rNn1+xmRk4Mzrd286ejRVcHKrEdVY3Ascx8rtr3OzMh1eLh5KVecBWh1H9TsXSvnzp0DoEOHDlXOLy8vJyUlBagc5Fu3bqVv3760b9/eNC0kJARvb2+2bNnCqFGjal1Lz5490ev1tV5PKQ52TiyP/rHW6wXd+zBfLTBWO79DK78GvxOis09nSsuvN+jvsARLjYE+P4t3P32SycMWEHTvw3Up1cTcfV+XPlB6H6xPH3h4eHD48OE6ravZIC8uLgbg+vWqOzUpKQmDwUDTpk3p2LGjafqpU6cYPXr0bcv7+/tz6tSpOtWi1+vJycmp07pKcLRX75/Webm5lJRdU7qMerPEGJSUXmPW6lGEdBnBqAfr/xkKc/e9GvdDpfY/zQa5h4cHBQUFHDlyhJCQkErz8vLymDFjBgCBgYHodDrTvIKCAlxdXW/bnpubG6dPn65zLWriYOekdAl11rpNG82ckTe05LQNnM07To4hg/3Hk26bv3L6KVo2a1/FmlUzd9+rcT+sTx/UJyc0G+RhYWGkp6cTFxfHwIED8fG5+bHc1NRUIiMjMRgMgGU+CFTXP5eUcqMU9i2y3O97/8X9ZtvWjxk/Yutgts0pxhJjMLBHJAN7RJpte+bue0vuh+baB5Xa/zT7ZmdsbCzNmzfnwoUL+Pv707VrVzp37kxwcDDe3t488sgjwO23HjZr1ozCwsLbtpefn4+bm/YetiOEUD/NBrmnpyfJyckMHToUR0dHsrKycHNzY9myZWzbto2MjAzg9iD38/Or8lr4qVOn8PPzs0jtQghRG5q9tAI3Q3nr1q23TS8qKiIrKwsbGxsCAgIqzRs2bBhvvvkm2dnZeHp6AnDw4EHOnDnDggULLFK3EELUhmbPyO/k5MmTGI1GOnfuTOPGld8Zj4qKonXr1owcOZKtW7eyfv16nnrqKYKDgxk5cqRCFQshRPWsMsjT0tKA2y+rALi4uLB3715at27N2LFjee655+jTpw9bt27FxsYqu0sIcZfT9KWV6twpyAHuvffeKi/JCKE2pWUlzP1sLOd+PkUjeydcm7QkJmIJbd07YTQa0el0rN01m0E9J9KqWQd0Oh0H079k+dbp2Nk60Mjeidixa+7KL2PIvvQjC5ImcKXYgLPjPcwYsxovD/9Ky1RUVLB863RST+/A1sYOF+fmvPrE32nr3gmAiwXnSdg0hWxDBjY6W4aHvMiovlMrbWP+Pyby1Q9r2PSXApo4uVqqebVilaeYfxTkom6yL/3Iy4v7MDHOhykLe5GlP1nlcgfTv+TFv3Xn+Q+6MTk+gF2H15jmlZb/SsKmaCbEdWby+1356+fjb1t/R+oqBs7QkfLvzQ3VFE0Z0juKVbGnWfbacUL8R/LBuueAm/eR/33b6xRdL+T0hUPEJUZytfgyizdN4aWRC1n22jG6dw5jy4ElCregags3PM+Q3lGsfj2DMaGvsyBp4m3LHDj1L05mpbDsteMsn3aC+zsN4JPtN5/saDQamb3mMcJ6PMOq2NOsnHGKh4KerLR+ctpG7GztLdGcerHKM/LfnsMizOu3A+vRXhP55sR6FiRN5KOXUystYzQaiUscT/wL+/FuE4g+P4v/WeBL34AIGjs2ZeWXf0Kn07E6NgOdTkf+1cqPNtDnZ7H94N/xa/+AJZumWg72jvT2G2J67df+AdZ/HQ9A/8AnaNWsA7HLBnA27zjznt2Og70jPxeew8ezJ9dLi8nI+YHg+27/8gmlFRRdJCP7MH+dvAuAfl0fZ/GmaHIMmaazbQAdOkrLf6W0rARbGzuulVylxT03b2I4+uMe7O0a8VDQ75/kbta01e+/45efSdw7j/jn97H90AoLtaxurPKMXJjfbwdWWPebZ9D9uj7OpcIL5Bgyb19Yp6OopBCAayVXcWncHHu7RlwvLWbHoZVMGjzX9GlbN5ffP+1WUVHBB+ueY8qoBOztGjV4m7Ro07cLCfG/+aZ9ctpGvjmxnkG9JjE85EU+XB/F1eLLGI1GTmal8OQ7rcjMPkLP+x5VuOrbXSq8gJtLa2xtb56L6nQ6WjZrz8XC85WWe6DLcILufZgxf/FgzF9aczRzDxMe/QsA5y6e4h7nFsz9dCwvfHg/s1c/Rt7ls6Z1P1g/mclD59PYsSl3O6s8Ixfmd6cDq9IZkk7Hn8cl8c6aCBwdnCm6XsCsZzZib+fAhdz/o2ljNxL3zuPIj7tpZO9E5MDZdO88AIAN33yAv9eD+Hj2UKSNd6OYhBByDFU/WGrJq0dp6drO9PrzPfPINWQy//k9APQNeIx+XSNYu2s297ULpn/gaNN/oAEd+/LPOVdJObmZVz/ux6rY03ft9eE7ycg+TJb+3yTOzKFxIxdWfvknFm54gT89/Sk3bpRz7MxeFkV/j5eHP1sOLGXOp0/y8cuH+fLgClq6tuf+To8o3YQakSAXNfJHgVFTN26U89med5k1YSOB3v05fSGVt1eNYPm0NG5UlPNzwTk6tOzCc0P+SmbOUV5fPpAV009SWHyJ5LQNfPDSN+ZqkiYsmnqgRsut2x/Pt//eyPyo3abnjf8W2s8Mml3lOjY2NvTrGsGqHW+RfSkD3/bBZqnZHFq4tiP/ah43bpRja2uH0WjkYsF5WrpWfjbMVz+spVunR0z/CQ3sOYE//X0QAC2btadTm/tNb5CG9YgkYdNLlN8o4/iZfaSd/YaD6b/f9BD1QSB/mfhPOrW93zKNrAUJclEjfxQY9naNanRgZeYe4/LVXAK9+wNwX7teuN/jSWbOUTq1vR8bnQ2P/Of7Eju1vR8Pt478lJdGjuFHfi7IYmJcZwDyf9Hzt/VR5F/NY3ifFxugxdqx/usP2Hcskbio3TU+qz6Yvo0B3cdxseA8Bb/oaXPLX1V3g2ZNWtKpbXd2H/mUR3tNJDltA+6unpX++gNo7ebNof/7ktEPTcfezoHv07fi5XHzQ4C9fMP5+7ZYDFdycL+nLYfSv6R9Sz/sbO154+nPKm1n4Awdy187cdf+VSJBLsyipgdWS9d25P+Sx7mf0+nQyo8cQyZ5l8/QrsV93OPsTrdOAzh8eie9/YaQl/8T+vyfaN/Kj+4+YZUCe9qSh4no9woPBoyycEvV5VJhNsu2TqO1mzfTl978khUHu0YkxBy843qbUxL4bPccbGxsmTIyAZfGd99zhl55fBkLkiaSuHcejR1dmPHkKgDeX/ccIV1G0Md/BCMenML5i+k8/2EQdjb2NGvqwSuP3/xOTicHZ16OWMpbK4cCRpwd7+Gtcf9QsEV1pzMajdU/hV1Ypbo+de7CxdMsSJrI1WuXTQdWx9Y3v3Xl1oNr79FEEvfOw0ZnQ4WxgqceeYNH7n8agLzLZ3l/3bNcKTZgo7NhfNjb9At8/LbfVV2Qh8YgTz+sp4EzdHW6Z9rcfa9kH9SVUvufBLm4jRoPoN9IkNefBHndKbX/yaUVIUQld/qqNHF3kvvIhRBC5STIhRBC5STIhRBC5eTNTnEboxEqypSuom5s7OGW79JWLTWOgbn7Xvqg5iTIhRBC5eTSihBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJwEuRBCqJyd0gUIy1PjdyGag5q+z1ONY1Tf/lVjm+/EkvubBLkVqiiDfYuUrsLyQmPA1kHpKmpGjWNU3/5VY5vvxJL7m1xaEUIIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZP7yIXZHD+zn+lLQytNc3RwxrOFD2HdIxn14FRsbWWXa0jWPgbW2n7ttUgoLrTbUwT7DsGIkYJf9Hz1w1qWbnmN8xfTefWJ5UqXZxWsfQysrf0S5MLsOrftTliP8abXw/u8xLPzfdl+aAWTBs/FtUkLBauzDtY+BtbWfrlGLhqck4Mzvh0ewGg0knv5jNLlWCVrHwOtt1+CXFhE3n8OHpfGbgpXYr2sfQy03H6rCHKDwUBsbCydOnXC0dGRdu3a8fLLL1NcXMyzzz6LTqdj8eLFSpepGSVl17hSbKCw6BI/5aWxaOMUMnOO4tsuGM8WPkqXZxWsfQysrf2av0Z+7NgxwsPD0ev1ODs706VLF3Jzc1m0aBFnzpwhPz8fgG7duilbqIas3TWLtbtmVZrWNyCCqY99pFBF1sfax8Da2q/pIDcYDAwfPhy9Xs+0adOYNWsWTZs2BWD+/Pm8/vrr2NnZodPpCAwMVLha7RjaO4r+gaMpryjjp7w0kvbHYbiSjYO9o2mZuZ+OpcJYwczIL0zTrl7LZ3K8P1HD4hnQfZwSpWtGTcYg7Wwyb64Mv23d8hulVFTcYOf8G5Ys2aysbR/U9KWVmJgYsrOziY6OJj4+3hTiALGxsQQFBVFeXo6XlxcuLi4KVqotbd07090njGDfcMaExjJn0hZOZ6eycMMLpmWmRnzMyawU9h5NNE1L2DQF/459VXUA3a1qMgZdvfuxZW5RpZ9VsRm4OLsz4dE5ClZff9a2D2o2yNPT00lKSsLd3Z333nuvymV69OgBQFBQkGnab8EfHBxMo0aN0KnlK2XuYv5efQjrHsn+40mczPoOuPmG07TRK1m8ORrDlVy+ObGeE2f280rEUoWr1aaqxuC/lZb/yjtrIwjw6svTA960cIUNS+v7oGaDPDExkYqKCsaNG0eTJk2qXMbJyQmoHOSZmZls2LABDw8PevXqZZFarcG4sJnY2NiyZufbpmm9fAfzUOCTxCWOJ2HjS7w2egUuzs0VrFLbqhqDWy3c8AKlZSXMGLPasoVZiJb3Qc0G+d69ewEIDQ2tdpns7GygcpD379+fvLw8/vWvfxEWFtawRVqRtu6dCA0ay9HMPaSdTTZNjxoeT87lTHr5htPbb6iCFWpfdWMAsOnbRRxM38o7Ezfj6NBYoQoblpb3Qc0G+blz5wDo0KFDlfPLy8tJSUkBKge5jY1mu0RxTw14CxudDWt2/X5G5OTgTGs3bzp6dFWwMutR1Rgcy9zHim2vMzNyHR5uXsoVZwFa3Qc1e9dKcXExANevX69yflJSEgaDgaZNm9KxY8cGraVnz57o9foG/R214WDnxPLoH82+3aB7H+arBcZq53do5afonRCdfTpTWl71/nC3qesY1XYM9PlZvPvpk0wetoCgex+uS6km9e1fc+yXd9M+WNv+8PDw4PDhw3X6XZoNcg8PDwoKCjhy5AghISGV5uXl5TFjxgwAAgMDG/wNTb1eT05OToP+jtpwtNfmn85/JC83l5Kya0qXUSOWGKOS0mvMWj2KkC4jGPVgdL23V9/+1dp+acn9TbNBHhYWRnp6OnFxcQwcOBAfn5uf5kpNTSUyMhKDwQBY5oNAHh4eDf47asPBzknpEhTRuk0bVZ2RN7TktA2czTtOjiGD/ceTbpu/cvopWjZrX+Pt1bd/tbZf1rY/6pMTOqPRWP3fISqWnZ1Nt27duHz5MnZ2dvj6+lJSUkJmZibh4eFUVFSwc+dOli9fzuTJk6vcxuzZs3nnnXfQWhfdKIV9i5SuwvJCY8DWQekqakaNY1Tf/lVjm+/EkvubZt/Z8/T0JDk5maFDh+Lo6EhWVhZubm4sW7aMbdu2kZGRAVR+o1MIIdRIs5dWAPz8/Ni6dett04uKisjKysLGxoaAgAAFKhNCCPPRdJBX5+TJkxiNRnx8fGjc+PY3WNavXw/AqVOnKr328vKiZ8+elitUCCFqwCqDPC0tDaj+ssro0aOrfD1hwgRWr17doLUJIURtSZBXQWtvbgrrVVpWwtzPxnLu51M0snfCtUlLYiKW0Na9E0ajEZ1Ox9pdsxnUcyKtmnVAp9NxMP1Llm+djp2tA43snYgdu0axZ3hnX/qRBUkTuFJswNnxHmaMWY2Xh/9ty20/tJJ/7PsrxooKunV6hJiIj7Gzta/zvB2pq9iUvNC0fcOVbLp692f2hI2knt7Jim2vm+YVFl/ErakHS1450sC9UT0JclErNT2w/mi5O82vbt7V4svMWDbAtI1fy66Rl3+WdW//XPX0WRc1+W0wtTWkdxTBvuHodDo2pyzmg3XP8f6L+0lO28DpC6mUlf/K6QuHWL3jz7w0ciGLN03hlSeW08NnIKt3zGTLgSW8OOJDRWpfuOF5hvSO4tFeE/nmxHoWJE3ko5dTKy2Tl/8Tq3fOZMnLR2jWtBVvrx7Jtu+XM/LBKXWeN7jXJAb3mmT6HZPjAxhw/80nIva671F63feoad6fPxlG0L3VPwrEEjR718qd7N27F6PRyNCh6nyugpJ+O7BWv57BmNDXWZA0sU7L3Wl+dfNcnJuz7LVjpp+hvaMIvi+8+ukS4jjYO9Lbb4jpQ29+7R/g54IsAPoHPkH/wCfYmfoJWw4sMT0w6ufCc/h49uR6aTEZOT/Q2s1bkdoLii6SkX2YsO43v0S5X9fHuVR4gRxDZqXlkk+sJ6TLCNxcPNDpdAx74AX2HUus17xbpZ8/SGHRRUL8R9w2z3All6M/7iGsR6S5m18rVhnkom5qemD90XJ3ml/T3wGwPXUlg4OfrfF0AZu+XUiI/0gAktM28s2J9QzqNYnhIS/y4foorhZfxmg0cjIrhSffaUVm9hF63nL2aUmXCi/g5tIaW9ubFw50Oh0tm7XnYuH5SstdLDxPq2a/P1PJw83LtExd591qx6GVDOgRabocc6tdh1cT7DuEZk1a1qOl9WeVl1ZE3dzpwGrr3qnGy91pvrPjPTX6HSezvqPoWgEP+A2rVGN107UqJiGEHEPVzydZ8upRWrq2M73+fM88cg2ZzH9+DwB9Ax6jX9cI1u6azX3tgukfONp05h7QsS//nHOVlJObefXjfqyKPU0TJ9cGb8/d5nppMfuP/YNFU7+/bZ7RaGRn6ie8NFL5TzFJkAuTPwqFu8mOQysZ2OMZU+D/0XStWjT1QI2WW7c/nm//vZH5UbtNj6n9LbSfGTS7ynVsbGzo1zWCVTveIvtSBr7tg81Sc021cG1H/tU8btwox9bWDqPRyMWC87R0rfzYgJau7cm9fMb0Wp+fZVqmrvN+883xdXRo5U+HVl1uq+/E2a8pLS9R7C+WW8mlFWGyaOoBNrxjqPKnpWu7SgcWUO2B9UfL3Wl+TX7H9V+L+PrEFzwa/D+Vfm91063d+q8/YN+xROImf1Xjs+qD6dsAuFhwnoJf9LS55a8hS2nWpCWd2nZn95FPgZvPhnF39az0lxncvPx24NS/yL+qx2g0svX7pTzcbWy95v1mxx0u020/tJJBPSdia2Nr7qbXmgS5qLGaHlh/tNyd5tfkd+w/noR36yDat/St9Hurm27NLhVms2zrNIquFzJ9aSjPf9CNqYt6/+F6m1MS+J/5vry5MpwpIxMUe+P4lceXse37ZUyM8+Ef+/7KjCdXAfD+uuf47uS/AGjd3JsJg97hlY8eZMJfO+Hq3IJhDzxfr3kAFy6e5kzuMR4OGnNbXcXXr5CStpHBve6OkwbNPjRLVK8+Dye6cPE0C5ImcvXaZRo7ujDjyVV0bH3zgfzvr3uOkC4j6OM/4o7L/dF2/mjdlxf3Ibz35Eq3h91p+m/koVk1M3CGjk1/Kaj1NXF5aFZlltzfJMitkNYOmJqSIK8ZCXLzsOT+Zh3vCAkhauxO37Aj7k5yjVwIIVROglwIIVROglwIIVROglwIIVRO7lqxQkYjVJQpXYXl2djDfz7MeNdT4xjVt3/V2OY7seT+JkEuhBAqJ5dWhBBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5eyULkDcndT4tVtq+iq3PyL9L31QGxLkokoVZbBvkdJV1E5oDNg6KF2FeUj/Sx/UhlxaEUIIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZMgF0IIlZP7yIXZHD+zn+lLQytNc3RwxrOFD2HdIxn14FRsbWWXa0jWPgbW2n7ttUgoLrTbUwT7DsGIkYJf9Hz1w1qWbnmN8xfTefWJ5UqXZxWsfQysrf0S5MLsOrftTliP8abXw/u8xLPzfdl+aAWTBs/FtUkLBauzDtY+BtbWfrlGLhqck4Mzvh0ewGg0knv5jNLlWCVrHwOtt1+CXFhE3n8OHpfGbgpXYr2sfQy03H65tCLMrqTsGleKDRiNN69PbjmwlMyco/i2C8azhY/S5VkFax8Da2u/VQS5wWBg/vz5bNy4kezsbFq0aEFERATz5s0jJiaGTz75hISEBKKjo5UuVRPW7prF2l2zKk3rGxDB1Mc+Uqgi62PtY2Bt7dd8kB87dozw8HD0ej3Ozs506dKF3NxcFi1axJkzZ8jPzwegW7duyhaqIUN7R9E/cDTlFWX8lJdG0v44DFeycbB3NC0z99OxVBgrmBn5hWna1Wv5TI73J2pYPAO6j1OidM2oyRiknU3mzZXht61bfqOUioob7Jx/w5Ilm5W17YOavkZuMBgYPnw4er2eadOmkZeXx5EjR9Dr9cTFxbFt2zZSU1PR6XQEBgYqXa5mtHXvTHefMIJ9wxkTGsucSVs4nZ3Kwg0vmJaZGvExJ7NS2Hs00TQtYdMU/Dv2VdUBdLeqyRh09e7HlrlFlX5WxWbg4uzOhEfnKFh9/VnbPqjpII+JiSE7O5vo6Gji4+Np2rSpaV5sbCxBQUGUl5fj5eWFi4uLgpVqm79XH8K6R7L/eBIns74Dbr7hNG30ShZvjsZwJZdvTqznxJn9vBKxVOFqtamqMfhvpeW/8s7aCAK8+vL0gDctXGHD0vo+qNkgT09PJykpCXd3d957770ql+nRowcAQUFBpmnr16/n8ccfp0OHDjRu3BhfX1/eeustioqKLFK3Vo0Lm4mNjS1rdr5tmtbLdzAPBT5JXOJ4Eja+xGujV+Di3FzBKrWtqjG41cINL1BaVsKMMastW5iFaHkf1GyQJyYmUlFRwbhx42jSpEmVyzg5OQGVgzw+Ph5bW1vmzZvH9u3befHFF1myZAmDBw+moqLCIrVrUVv3ToQGjeVo5h7SziabpkcNjyfncia9fMPp7TdUwQq1r7oxANj07SIOpm/lnYmbcXRorFCFDUvL+6Bmg3zv3r0AhIaGVrtMdnY2UDnIt2zZwhdffMG4ceN46KGHePnll1m8eDEpKSl8++23DVu0xj014C1sdDas2fX7GZGTgzOt3bzp6NFVwcqsR1VjcCxzHyu2vc7MyHV4uHkpV5wFaHUf1OxdK+fOnQOgQ4cOVc4vLy8nJSUFqBzkLVrc/tHdnj17ApCTk1OnWnr27Iler6/TukpxsHNiefSPtVon6N6H+WqBsdr5HVr5NeidEJ19OlNafr3Btm9Jdel/qP0Y6POzePfTJ5k8bAFB9z5cl1JNzN3/1rYPenh4cPjw4Tqtq9kgLy4uBuD69ao7NSkpCYPBQNOmTenYseMdt7Vv3z4A/Pz86lSLXq+v838CSnG0V9+f13m5uZSUXVO6DLOwRP+XlF5j1upRhHQZwagH6/8ZCnP3v+yDNafZIPfw8KCgoIAjR44QEhJSaV5eXh4zZswAIDAwEJ1OV+12cnJymDlzJoMHD67zveYeHh51Wk9JDnZOSpdQa63btNHUGXlDS07bwNm84+QYMth/POm2+Sunn6Jls/Y13p65+9/a9sH65ITOaDRW/3eIisXExJCQkEC7du3YvXs3Pj43P5abmppKZGQkZ8+epaysjClTprB48eIqt1FUVMTDDz+MXq8nNTWV1q1bW7IJirpRCvsWKV1F7YTGgK2D0lWYh/S/9EFtaPbNztjYWJo3b86FCxfw9/ena9eudO7cmeDgYLy9vXnkkUeAytfHb3X9+nWGDx/OTz/9xK5du6wqxIUQ6qLZIPf09CQ5OZmhQ4fi6OhIVlYWbm5uLFu2jG3btpGRkQFUHeRlZWU88cQTHD58mO3bt9OlSxdLly+EEDWm2WvkcPPNya1bt942vaioiKysLGxsbAgICKg077d7z/fs2cOXX35JcHCwpcoVQog60XSQV+fkyZMYjUZ8fHxo3LjyO+NTpkxh3bp1/OlPf6Jx48Z8//33pnn33ntvlbcnCiGEkjR7aeVO0tLSgKovq2zfvh2Av/71r4SEhFT62bZtm0XrFEKImrDKM/I7BXlWVpaFqxGi4ZSWlTD3s7Gc+/kUjeydcG3SkpiIJbR174TRaESn07F212wG9ZxIq2Yd0Ol0HEz/kuVbp2Nn60Ajeydix65R3Zcx3Knd/+0f++L46vAa7GwdcLB3ZMrIRfi2V9clVQly0eBqc1AdTP+S1Tv/TEVFBRUV5Yx+eAaDek5QoGrtGNI7imDfcHQ6HZtTFvPBuud4/8X9JKdt4PSFVMrKf+X0hUOs3vFnXhq5kMWbpvDKE8vp4TOQ1TtmsuXAEl4c8aHSzai16tp9q8ycY2z57mNWTD+JU6Mm7P7hUxZvjmZxzCFliq4jqwzy357DIiynJgeV0WgkLnE88S/sx7tNIPr8LP5ngS99AyJo7Ni06g2LO3Kwd6S33xDTa7/2D7D+63gA+gc+QatmHYhdNoCzeceZ9+x2HOwd+bnwHD6ePbleWkxGzg8E33f7l0/c7e7U7lvpdDrKK8ooKS3GqVETikoKcb/H05KlmoVVBrmwrJoeVADodBSVFAJwreQqLo2bY2/XyAJVWodN3y4kxH8kAMlpG/m/8wcZ1GsSAV4P8uH6KF4c8SFGo5GTWSnM/WwsTg5NeGnE35Qt2gxubfet7m0TxOP9XiXyvY40beyGvW0jPnjpGwUqrB8JcmFx1R1UOp2OP49L4p01ETg6OFN0vYBZz2zE3k4jH9dsADEJIeQYqn6w1JJXj9LStZ3p9ed75pFryGT+83sA6BvwGP26RrB212zuaxdM/8DRpsdVBHTsyz/nXCXl5GZe/bgfq2JP08TJtcHb0xD+u923ysv/iW/TNrL69Uzc72nD5pTFvPvpGP42RV1POpUgF/VWnzC51Y0b5Xy2511mTdhIoHd/Tl9I5e1VI1g+LY17nN0brH41WzT1QI2WW7c/nm//vZH5UbtNzxv/LbSfGTS7ynVsbGzo1zWCVTveIvtShureAISq232rb09soGPrrrjf0waAR3tN4qPNUykrL1XVCYQEuai3+oTJrTJzj3H5ai6B3v0BuK9dL9zv8SQz5yg9fAaatWZrsv7rD9h3LJG4qN01Pqs+mL6NAd3HcbHgPAW/6GlTxRvTd7uatNujuTc7D6/i+q9FODVqwsFTW/Fs4aOqEAcJcmEhNTmoWrq2I/+XPM79nE6HVn7kGDLJu3yGdi3us2yxGnKpMJtlW6fR2s2b6UtvfsmKg10jEmIO3nG9zSkJfLZ7DjY2tkwZmYBLYzdLlGs2d2r36p1v09ylDcNDXqBvwGNkXEhlysKe2Ns1wtHBmTee/lzh6mtPs08/FPVjzifPXSrM5um57Wjt5o1To5t3n9waJu+ve46QLiPo4z+CvUcTSdw7DxudDRXGCp565A0euf/pGv0eefqheQycoWPTXwpqfU1cnn6o3D4oZ+SiwbVw9bzjt7ZMG73C9O9H7n+KR+5/yhJlCaEZEuRCiEru9J+uuDtZ5bNWhBBCSyTIhRBC5STIhRBC5eSuFVEloxEqypSuonZs7OEO36OtKtL/0ge1IUEuhBAqJ5dWhBBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5STIhRBC5eyULkDcneT7EpUl/S99UBsS5KJKFWWwb5HSVdROaAzYOihdhXlI/0sf1IZcWhFCCJWTIBdCCJWTIBdCCJWTIBdCCJWTIBdCCJWTIBdCCJWT2w+F2Rw/s5/pS0MrTXN0cMazhQ9h3SMZ9eBUbG1ll2tI1j4G1tp+7bVIKC6021ME+w7BiJGCX/R89cNalm55jfMX03n1ieVKl2cVrH0MrK39EuTC7Dq37U5Yj/Gm18P7vMSz833ZfmgFkwbPxbVJCwWrsw7WPgbW1n65Ri4anJODM74dHsBoNJJ7+YzS5Vglax8DrbdfglxYRN5/Dh6Xxm4KV2K9rH0MtNx+ubQizK6k7BpXig0YjTevT245sJTMnKP4tgvGs4WP0uVZBWsfA2trv1UEucFgYP78+WzcuJHs7GxatGhBREQE8+bNIyYmhk8++YSEhASio6OVLlUT1u6axdpdsypN6xsQwdTHPlKoIutj7WNgbe3XfJAfO3aM8PBw9Ho9zs7OdOnShdzcXBYtWsSZM2fIz88HoFu3bsoWqiFDe0fRP3A05RVl/JSXRtL+OAxXsnGwdzQtM/fTsVQYK5gZ+YVp2tVr+UyO9ydqWDwDuo9TonTNqMkYpJ1N5s2V4betW36jlIqKG+ycf8OSJZuVte2Dmr5GbjAYGD58OHq9nmnTppGXl8eRI0fQ6/XExcWxbds2UlNT0el0BAYGKl2uZrR170x3nzCCfcMZExrLnElbOJ2dysINL5iWmRrxMSezUth7NNE0LWHTFPw79lXVAXS3qskYdPXux5a5RZV+VsVm4OLszoRH5yhYff1Z2z6o6SCPiYkhOzub6Oho4uPjadq0qWlebGwsQUFBlJeX4+XlhYuLi4KVapu/Vx/Cukey/3gSJ7O+A26+4TRt9EoWb47GcCWXb06s58SZ/bwSsVTharWpqjH4b6Xlv/LO2ggCvPry9IA3LVxhw9L6PqjZIE9PTycpKQl3d3fee++9Kpfp0aMHAEFBQaZpycnJhIWF0bp1axo1aoSnpydjxowhPT3dInVr1biwmdjY2LJm59umab18B/NQ4JPEJY4nYeNLvDZ6BS7OzRWsUtuqGoNbLdzwAqVlJcwYs9qyhVmIlvdBzQZ5YmIiFRUVjBs3jiZNmlS5jJOTE1A5yAsKCujatSuLFi1i165dxMXFcfLkSUJCQsjOzrZI7VrU1r0ToUFjOZq5h7SzyabpUcPjybmcSS/fcHr7DVWwQu2rbgwANn27iIPpW3ln4mYcHRorVGHD0vI+qNkg37t3LwChoaHVLvNbMN8a5CNGjODDDz9k9OjRPPTQQ4wbN46NGzdy5coVNmzY0LBFa9xTA97CRmfDml2/nxE5OTjT2s2bjh5dFazMelQ1Bscy97Fi2+vMjFyHh5uXcsVZgFb3Qc3etXLu3DkAOnToUOX88vJyUlJSgMpBXpXmzW/+qWVnp9nuMougex/mqwXGaud3aOWn6jsh1KC2Y6DPz+LdT59k8rAFBN37sAUqbFjWug9qNpmKi4sBuH79epXzk5KSMBgMNG3alI4dO942/8aNG1RUVHDu3DneeOMNPDw8ePLJJ+tUS8+ePdHr9XVaVykOdk4sj/5R6TJqpbNPZ0rLqx5vtbFE/5eUXmPW6lGEdBnBqAfr/xkKc/e/te2DHh4eHD58uE7rajbIPTw8KCgo4MiRI4SEhFSal5eXx4wZMwAIDAxEp9Pdtv5DDz1kOmPv1KkTe/fupUWLuj1oR6/Xk5OTU6d1leJor77rpHm5uZSUXVO6DLOwRP8np23gbN5xcgwZ7D+edNv8ldNP0bJZ+xpvz9z9L/tgzWk2yMPCwkhPTycuLo6BAwfi43PzY7mpqalERkZiMBiA6j8ItHLlSgoLC/npp59YsGABgwYNIiUlhfbta75j/8bDw6PO7VCKg52TxX7X+y/uN8t2Wrdpo6kz8oY2sEckA3tEmm175u5/a9sH65MTOqPRWP0FJRXLzs6mW7duXL58GTs7O3x9fSkpKSEzM5Pw8HAqKirYuXMny5cvZ/LkyXfcVmFhIV5eXowfP57FixdbqAXKulEK+xYpXUXthMaArYPSVZiH9L/0QW1o9q4VT09PkpOTGTp0KI6OjmRlZeHm5sayZcvYtm0bGRkZwB+/0Qng6upKp06dyMzMbOiyhRCi1jR7aQXAz8+PrVu33ja9qKiIrKwsbGxsCAgI+MPtXLx4kdOnT9O7d++GKFMIIepF00FenZMnT2I0GvHx8aFx48pvqIwfP55OnTrRrVs3XF1d+fHHH/nwww+xs7Pj1VdfVahiIYSonlUGeVpaGlD1ZZUHHniAtWvXsnDhQkpKSmjXrh2hoaG8+eab1d6TLoQQSpIg/y/R0dHyXHKhGaVlJcz9bCznfj5FI3snXJu0JCZiCW3dO2E0GtHpdKzdNZtBPSfSqlkHdDodB9O/ZPnW6djZOtDI3onYsWvu+i9juFM7/9v3p7ayfOt0bhhv0NGjKzPGrMbZ0YXU0ztZse1103KFxRdxa+rBkleOWLIpdSJBLsyuNgdV9qUfWZA0gSvFBpwd72HGmNV4efgDcLX4MjOWDTAt+2vZNfLyz7Ju1kVNfl1XQxnSO4pg33B0Oh2bUxbzwbrneP/F/SSnbeD0hVTKyn/l9IVDrN7xZ14auZDFm6bwyhPL6eEzkNU7ZrLlwBJeHPGh0s34Q9W181bXfy3i/XXP8v6LX9O+pS8Jm6L5bPccooYtoNd9j9LrvkdNy/75k2EE3Vv9Iz7uJpq9a+VO9u7di9FoZOhQdT4gRw2G9I5iVexplr12nBD/kXyw7rkql1u44XmG9I5i9esZjAl9nQVJE03zXJybs+y1Y6afob2jCL4vXEK8FhzsHentN8T0oTe/9g/wc0EWAP0Dn6B/4BPsTP2ELQeWmJ7893PhOXw8e3K9tJiMnB9o7eatYAtq5k7tvNWh/9tOpzb3076lLwAj+rzEvmOJty1nuJLL0R/3EGbG++wbklUGuWhYNT2oCooukpF9mLDu4wHo1/VxLhVeIMdQ9W2e21NXMjj42Qar2xps+nYhIf4jAUhO28g3J9YzqNckhoe8yIfro7hafBmj0cjJrBSefKcVmdlH6HnLWapa3NrOW10sPE+rZr+/19WqmRf5V/O4caO80nK7Dq8m2HcIzZq0bPBazcEqL60Iy6ruoLpUeAE3l9bY2t7cDXU6HS2btedi4fnbLsOczPqOomsFPOA3zCI1q0VMQgg5hqqfR7Lk1aO0dG1nev35nnnkGjKZ//weAPoGPEa/rhGs3TWb+9oF0z9wtOk/34COffnnnKuknNzMqx/3Y1XsaZo4uTZ4e8zhv9tZW0ajkZ2pn/DSSPV8GkmCXNRafcKjrnYcWsnAHs+YQl/ctGjqgRott25/PN/+eyPzo3abnjf+W2g/M2h2levY2NjQr2sEq3a8RfalDHzbB5ul5oZUVTtv1dK1PUcyvjK9/rkgq9LJBMCJs19TWl6iqr9E5KgQtVaf8LhVC9d2pj9rbW3tMBqNXCw4T0vXys+zuf5rEV+f+ILFMalmqd/arP/6A/YdSyQuaneNz6oPpm9jQPdxXCw4T8EvetpU8Ub13aYm7ex132AWb5rC+Yv/R/uWvvzru495OGhspWW2H1rJoJ4TsbWxtUDV5iFBLhpETQ6qZk1a0qltd3Yf+ZRHe00kOW0D7q6et11W2X88Ce/WQaY3qETNXSrMZtnWabR282b60pt3YDjYNSIh5uAd19ucksBnu+dgY2PLlJEJd/0bzHdq5+qdb9PcpQ3DQ16gsWNTXh29gtmrR3GjohwvjwBix6wxbaf4+hVS0jayfFqaUk2pE80+NEvUT30eWHSpMJun57ajtZs3To1ufuH1reHx/rrnCOkygj7+I7hw8TQLkiZy9dplGju6MOPJVXRsXfmbWl5e3Ifw3pMZ3GvSHX+vPDTLPAbO0LHpLwW1viYuD81Sbh+UM3Jhdi1cPe/4LS3TRq8w/btdy/v+8FLNwuiqv/VdCHGTBLkQopI7/Scs7k5yH7kQQqicBLkQQqicBLkQQqicBLkQQqic3H4oqmQ0QkWZ0lXUjo09/OfDiqon/S99UBsS5EIIoXJyaUUIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVROglwIIVTOTukCxN1JvmZLWdL/0ge1IUEuqlRRBvsWKV1F7YTGgK2D0lWYh/S/9EFtyKUVIYRQOQlyIYRQOQlyIYRQOQlyIYRQOQlyIYRQOQlyIYRQOQlyIYRQObmPXJjN8TP7mb40tNI0RwdnPFv4ENY9klEPTsXWVna5hmTtY2Ct7ddei4TiQrs9RbDvEIwYKfhFz1c/rGXpltc4fzGdV59YrnR5VsHax8Da2i9BLsyuc9vuhPUYb3o9vM9LPDvfl+2HVjBp8Fxcm7RQsDrrYO1jYG3tl2vkosE5OTjj2+EBjEYjuZfPKF2OVbL2MdB6+yXIhUXk/efgcWnspnAl1svax0DL7ZdLK8LsSsqucaXYgNF48/rklgNLycw5im+7YDxb+ChdnlWw9jGwtvZrPsgNBgPz589n48aNZGdn06JFCyIiIpg3bx4xMTF88sknJCQkEB0drXSpmrF21yzW7ppVaVrfgAimPvaRQhVZH2sfA2trv6aD/NixY4SHh6PX63F2dqZLly7k5uayaNEizpw5Q35+PgDdunVTtlCNGdo7iv6BoymvKOOnvDSS9sdhuJKNg72jaZm5n46lwljBzMgvTNOuXstncrw/UcPiGdB9nBKla0ZNxiDtbDJvrgy/bd3yG6VUVNxg5/wblizZrKxtH9TsNXKDwcDw4cPR6/VMmzaNvLw8jhw5gl6vJy4ujm3btpGamopOpyMwMFDpcjWlrXtnuvuEEewbzpjQWOZM2sLp7FQWbnjBtMzUiI85mZXC3qOJpmkJm6bg37Gvqg6gu1VNxqCrdz+2zC2q9LMqNgMXZ3cmPDpHwerrz9r2Qc0GeUxMDNnZ2URHRxMfH0/Tpk1N82JjYwkKCqK8vBwvLy9cXFwUrFT7/L36ENY9kv3HkziZ9R1w8w2naaNXsnhzNIYruXxzYj0nzuznlYilClerTVWNwX8rLf+Vd9ZGEODVl6cHvGnhChuW1vdBTQZ5eno6SUlJuLu7895771W5TI8ePQAICgqqdjvh4eHodDpmz57dEGValXFhM7GxsWXNzrdN03r5DuahwCeJSxxPwsaXeG30ClycmytYpbZVNQa3WrjhBUrLSpgxZrVlC7MQLe+DmgzyxMREKioqGDduHE2aNKlyGScnJ6D6IP/iiy84duxYQ5Voddq6dyI0aCxHM/eQdjbZND1qeDw5lzPp5RtOb7+hClaofdWNAcCmbxdxMH0r70zcjKNDY4UqbFha3gc1GeR79+4FIDQ0tNplsrOzgaqD/OrVq7zyyivEx8c3TIFW6qkBb2Gjs2HNrt/PiJwcnGnt5k1Hj64KVmY9qhqDY5n7WLHtdWZGrsPDzUu54ixAq/ugJu9aOXfuHAAdOnSocn55eTkpKSlA1UH+1ltv4ePjw7hx4xg/fvxt82urZ8+e6PX6em/HkhzsnFge/WOt1gm692G+WmCsdn6HVn4NeidEZ5/OlJZfb7DtW1Jd+h9qPwb6/Cze/fRJJg9bQNC9D9elVBNz97+17YMeHh4cPny4TutqMsiLi4sBuH696g5NSkrCYDDQtGlTOnbsWGne4cOH+fvf/84PP/xgtnr0ej05OTlm254lONqr78/rvNxcSsquKV2GWVii/0tKrzFr9ShCuoxg1IP1/xyFuftf9sGa02SQe3h4UFBQwJEjRwgJCak0Ly8vjxkzZgAQGBiITqczzbtx4wbPP/880dHR+Pv7m7UetXGwc1K6hFpr3aaNps7IG1py2gbO5h0nx5DB/uNJt81fOf0ULZu1r/H2zN3/1rYP1icnNBnkYWFhpKenExcXx8CBA/HxufmR3NTUVCIjIzEYDMDtHwRavHgxP//8s9nvUqnrn0tKulEK+xZZ5ne9/+J+s2znx4wfsXUwy6YUZ4n+H9gjkoE9Is22PXP3v+yDNafJNztjY2Np3rw5Fy5cwN/fn65du9K5c2eCg4Px9vbmkUceASpfHzcYDMycOZO3336b8vJyCgsLKSwsBKCkpITCwkIqKiqUaI4QQtyRJoPc09OT5ORkhg4diqOjI1lZWbi5ubFs2TK2bdtGRkYGUDnIs7Oz+eWXX3j++edp1qyZ6QcgLi6OZs2acf78eUXaI4QQd6LJSysAfn5+bN269bbpRUVFZGVlYWNjQ0BAgGl6p06d2Ldv323Lh4aGMmHCBCZOnKjKa91CCO3TbJBX5+TJkxiNRnx8fGjc+Pd3xZs0acLDDz9c5TpeXl7VzhNCCKVp8tLKnaSlpQF3/mi+EEKoidWdkdc2yI3G6j9cIMTdrrSshLmfjeXcz6doZO+Ea5OWxEQsoa17J4xGIzqdjrW7ZjOo50RaNeuATqfjYPqXLN86HTtbBxrZOxE7do0qvozh9eWDKPhFj05nQ2PHpkwZuYhObe+vtIw+P4sFSRPJzD2KR7OOLHvt2G3bMRqNxC4bwI85R9g8p9AyxdeTnJGLBvH68kFEvR/I8x9049WP+5GZc7TK5T7aHMP4eV4MnKEjM+dYtdvbkbqKgTN0pPx7c8MUrGFDekexKvY0y147Toj/SD5Y9xxw8z7yv297naLrhZy+cIi4xEiuFl9m8aYpvDRyIcteO0b3zmFsObBE4RbUzMzIL1g+7QTLXjvG4/1eY0HSxNuWaezowqTB7/Lm059Xu50N33xI6+b3NmCl5md1Qb53716MRiNDh6rz4ThqUZODCqBf4BN8+NK3tGpW9eMU4OZZ1PaDf8ev/QMNVK12Odg70ttviOmDb37tH+DngiwA+gc+Qf/AJ9iZ+glbDiwxPfnv58Jz+Hj25HppMRk5P9DazVvBFtRcEydX07+LS64AutuWcWnsRkDHvjg6OFe5jSz9Sb47uZmxoX9qoCobhtUFubCMmhxUAIHe/Wnh6lntdioqKvhg3XNMGZWAvV0jM1dpfTZ9u5AQ/5EAJKdt5JsT6xnUaxLDQ17kw/VRXC2+jNFo5GRWCk++04rM7CP0vO9RhauuubjEZ3j63Xas2TmTPz31v7Vat/xGGR+un8zLjy/Dxsa2gSpsGFZ3jVxYTlziMxw/c/OWzrnPflmnbWz45gP8vR7Ex7OHOUvTjJiEEHIMVT9YasmrR2np2s70+vM988g1ZDL/+T0A9A14jH5dI1i7azb3tQumf+Bo05l7QMe+/HPOVVJObubVj/uxKvZ0pf+c71avP7UWgF2H1/D3L19nXi32u//96h36BkTQoZUf+vysBqqwYUiQi1qraXjU56AC+En/b5LTNvDBS9/Ur2ANWzT1QI2WW7c/nm//vZH5UbtNzxv/LbSfGTS7ynVsbGzo1zWCVTveIvtSBr7tg81SsyUM6jmBhRte4Grx5Rp/UcSJs19zseA8//xuMTcqyrn261XGz/NicUwqrk1aNHDF9SNBLmqtpuHxm7ocVAD/PpvMzwVZTIzrDED+L3r+tj6K/Kt5DO/zYq1qsGbrv/6AfccSiYvaXeOz6oPp2xjQfRwXC85T8IueNu6dGrbIeiq6XkhJ6TXc72kDQMq/N+Pi3Jymjd1qvI0PX/r9yyb0+Vm88GE3Pn0zy9ylNggJcmF25jioAIb3ebFSYE9b8jAR/V7hwYBR5ixX0y4VZrNs6zRau3kzfenNL1pxsGtEQszBO663OSWBz3bPwcbGlikjE3Cp5dhZWnHJFeb872h+LbuOjc6Ge5xbMGfSVnQ6He+ve46QLiPo4z+CktJrTJrvQ1n5rxSXXOGpdz0J6x7Js0Oq/kpItZAgF2Z3p4MKqHRg/W398xz8v23k/6LnjRWP0rhRU9b8KVPhFmhHC1fPO37RQnXee26HKq6J/6ZVsw4sjjlU5bxpo1eY/u3o0JjEP2f/4fY83LxUcw85SJCLBnCngwoqH1ivPLGsxts116NGhdAaCXIhRCV1OYMXypL7yIUQQuUkyIUQQuUkyIUQQuV0Rnm8n6iC0QgVZUpXUTs29qCr+kkAqiP9L31QGxLkQgihcnJpRQghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVE6CXAghVM5O6QLE3Um+ZktZ0v/SB7UhQS6qVFEG+xYpXUXthMaArYPSVZiH9L/0QW3IpRUhhFA5CXIhhFA5CXIhhFA5CXIhhFA5CXIhhFA5CXIhhFA5CXIhhFA5uY9cmM3xM/uZvjS00jRHB2c8W/gQ1j2SUQ9OxdZWdrmGZO1jYK3t116LhOJCuz1FsO8QjBgp+EXPVz+sZemW1zh/MZ1Xn1iudHlWwdrHwNraL0EuzK5z2+6E9Rhvej28z0s8O9+X7YdWMGnwXFybtFCwOutg7WNgbe2Xa+SiwTk5OOPb4QGMRiO5l88oXY5VsvYx0Hr7JciFReT95+BxaeymcCXWy9rHQMvtl0srwuxKyq5xpdiA0Xjz+uSWA0vJzDmKb7tgPFv4KF2eVbD2MbC29ltFkBsMBubPn8/GjRvJzs6mRYsWREREMG/ePGJiYvjkk09ISEggOjpa6VI1Ye2uWazdNavStL4BEUx97COFKrI+1j4G1tZ+zQf5sWPHCA8PR6/X4+zsTJcuXcjNzWXRokWcOXOG/Px8ALp166ZsoRoytHcU/QNHU15Rxk95aSTtj8NwJRsHe0fTMnM/HUuFsYKZkV+Ypl29ls/keH+ihsUzoPs4JUrXjJqMQdrZZN5cGX7buuU3SqmouMHO+TcsWbJZWds+qOlr5AaDgeHDh6PX65k2bRp5eXkcOXIEvV5PXFwc27ZtIzU1FZ1OR2BgoNLlakZb98509wkj2DecMaGxzJm0hdPZqSzc8IJpmakRH3MyK4W9RxNN0xI2TcG/Y19VHUB3q5qMQVfvfmyZW1TpZ1VsBi7O7kx4dI6C1defte2Dmg7ymJgYsrOziY6OJj4+nqZNm5rmxcbGEhQURHl5OV5eXri4uChYqbb5e/UhrHsk+48ncTLrO+DmG07TRq9k8eZoDFdy+ebEek6c2c8rEUsVrlabqhqD/1Za/ivvrI0gwKsvTw9408IVNiyt74OaDfL09HSSkpJwd3fnvffeq3KZHj16ABAUFGSatn//fnQ63W0/cumlfsaFzcTGxpY1O982TevlO5iHAp8kLnE8CRtf4rXRK3Bxbq5gldpW1RjcauGGFygtK2HGmNWWLcxCtLwPajbIExMTqaioYNy4cTRp0qTKZZycnIDKQf6bjz76iAMHDph+/vd//7dB69W6tu6dCA0ay9HMPaSdTTZNjxoeT87lTHr5htPbb6iCFWpfdWMAsOnbRRxM38o7Ezfj6NBYoQoblpb3Qc0G+d69ewEIDQ2tdpns7Gyg6iDv0qULDzzwgOmna9euDVOoFXlqwFvY6GxYs+v3MyInB2dau3nT0UP61xKqGoNjmftYse11Zkauw8PNS7niLECr+6Bm71o5d+4cAB06dKhyfnl5OSkpKUDVQW5OPXv2RK/XN+jvMDcHOyeWR/9Yq3WC7n2YrxYYq53foZVfg94J0dmnM6Xl1xts+5ZUl/6H2o+BPj+Ldz99ksnDFhB078N1KdXE3P1vbfugh4cHhw8frtO6mg3y4uJiAK5fr7pTk5KSMBgMNG3alI4dO942f8yYMRgMBpo3b86IESP461//iru7e51q0ev15OTk1GldpTjaq+/P67zcXErKrildhllYov9LSq8xa/UoQrqMYNSD9f8Mhbn7X/bBmtNskHt4eFBQUMCRI0cICQmpNC8vL48ZM2YAEBgYiE6nM8275557mDFjBv3796dJkyYcOHCA9957j++//57Dhw/j6OhIbXl4eNSvMQpwsHNSuoRaa92mjabOyBtactoGzuYdJ8eQwf7jSbfNXzn9FC2bta/x9szd/9a2D9YnJ3RGo7H6v0NULCYmhoSEBNq1a8fu3bvx8bn5sdzU1FQiIyM5e/YsZWVlTJkyhcWLF99xW1u2bGHEiBF88sknTJo0yRLlK+5GKexbpHQVtRMaA7YOSldhHtL/0ge1odk3O2NjY2nevDkXLlzA39+frl270rlzZ4KDg/H29uaRRx4BanZ9fNiwYTg7O9f5+pUQQjQkzQa5p6cnycnJDB06FEdHR7KysnBzc2PZsmVs27aNjIwMoHZvdN56CUYIIe4Wmr1GDuDn58fWrVtvm15UVERWVhY2NjYEBAT84Xb+9a9/UVxcTHBwcEOUKYQQ9aLpIK/OyZMnMRqN+Pj40Lhx5XfGx48fj7e3N927dze92Tl//ny6devG2LFjFapYCCGqZ5VBnpaWBlR9WcXf35/PP/+cv/3tb1y/fh1PT08mT57MrFmzcHDQyDtpQghNkSD/L2+88QZvvPGGpUsSokGUlpUw97OxnPv5FI3snXBt0pKYiCW0de+E0WhEp9OxdtdsBvWcSKtmHdDpdBxM/5LlW6djZ+tAI3snYseuUcWXMWRf+pEFSRO4UmzA2fEeZoxZjZeHf6VldqSuYlPyQtNrw5Vsunr3Z/aEjQD8Y18cXx1eg52tAw72jkwZuQjf9nf/JVUJctEganJQXS2+zIxlA0yvfy27Rl7+WdbNuohLYzcOpn/J6p1/pqKigoqKckY/PINBPSdYuimqN6R3FMG+4eh0OjanLOaDdc/x/ov7SU7bwOkLqZSV/8rpC4dYvePPvDRyIYs3TeGVJ5bTw2cgq3fMZMuBJbw44kOlm/GHFm54niG9o3i010S+ObGeBUkT+ejl1ErLDO41icG9fr+FeHJ8AAPuv/nI2sycY2z57mNWTD+JU6Mm7P7hUxZvjmZxzCGLtqMurDLIf3sOi2g4NTmoXJybs+y1Y6bX6/bHc+Ls17g0dsNoNBKXOJ74F/bj3SYQfX4W/7PAl74BETR2bIqoGQd7R3r7DTG99mv/AOu/jgegf+ATtGrWgdhlAzibd5x5z27Hwd6RnwvP4ePZk+ulxWTk/EDwfbd/+cTdpqDoIhnZh/nr5F0A9Ov6OIs3RZNjyKSte6cq10k/f5DCoouE+I8Abt6VVl5RRklpMU6NmlBUUoj7PZ4Wa0N9aPb2Q6Gc3w6qsO7jgZsH1aXCC+QYMu+43vbUlQwOfvb3CTodRSWFAFwruYpL4+bY2zVqqLKtwqZvFxLiPxKA5LSNfHNiPYN6TWJ4yIt8uD6Kq8WXMRqNnMxK4cl3WpGZfYSe9z2qcNV/7FLhBdxcWmNre/PcVKfT0bJZey4Wnq92nR2HVjKgRyR2tvYA3NsmiMf7vUrkex156l1PNn7zIdGjEixSf31Z5Rm5aFh3OqiqOzs6mfUdRdcKeMBvmGmdP49L4p01ETg6OFN0vYBZz2zE3k7ecL5VTEIIOYaqHyy15NWjtHRtZ3r9+Z555Boymf/8HgD6BjxGv64RrN01m/vaBdM/cLTpsxIBHfvyzzlXSTm5mVc/7seq2NM0cXJt8PZYyvXSYvYf+weLpn5vmpaX/xPfpm1k9euZuN/Ths0pi3n30zH8bcq3ClZaMxLkotb+KDzqYsehlQzs8Ywp/G/cKOezPe8ya8JGAr37c/pCKm+vGsHyaWnc41y3h5dp0aKpB2q03Lr98Xz7743Mj9ptet74b6H9zKDZVa5jY2NDv64RrNrxFtmXMu7qN/1auLYj/2oeN26UY2trh9Fo5GLBeVq6Vv2smG+Or6NDK386tOpimvbtiQ10bN0V93vaAPBor0l8tHkqZeWld/0JhAS5qLU/Cg97u0a1Oqiu/1rE1ye+YHHM79fQM3OPcflqLoHe/QG4r10v3O/xJDPnKD18BpqvMVZg/dcfsO9YInFRu2t8Vn0wfRsDuo/jYsF5Cn7R06aav6TuFs2atKRT2+7sPvIpj/aaSHLaBtxdPav9C3DHf1/GAzyae7Pz8Cqu/1qEU6MmHDy1Fc8WPnd9iIMEuWgAtT2o9h9Pwrt1EO1b+pqmtXRtR/4veZz7OZ0OrfzIMWSSd/kM7VrcZ6lmaMKlwmyWbZ1Gazdvpi+9+SUrDnaNSIg5eMf1Nqck8NnuOdjY2DJlZAIujd0sUW69vPL4MhYkTSRx7zwaO7ow48lVALy/7jlCuoygz3/e1Lxw8TRnco8x93++rLR+34DHyLiQypSFPbG3a4SjgzNvPP25xdtRF5p9+qGon/o+ee7CxdMsSJrI1WuXTQdVx9Y3v4Hlvw+slxf3Ibz35Eq3hQHsPZpI4t552OhsqDBW8NQjb/DI/U9X+zvl6YfmMXCGjk1/Kaj1NXF5+qFy+6CckYsG0a7lfdVegpk2ekWl1wujq/5W90fuf4pH7n/K7LUJoTUS5EKISu70VWni7iT3kQshhMpJkAshhMpJkAshhMrJXSuiSkYjVJQpXUXt2NiDVr7ESfpf+qA2JMiFEELl5NKKEEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKonAS5EEKo3P8Dz3fhQkaSOoUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit.circuit import Parameter\n", "\n", "beta = Parameter(\"β\")\n", "\n", "ws_mixer = QuantumCircuit(len(sigma))\n", "for idx, theta in enumerate(thetas):\n", " ws_mixer.ry(-theta, idx)\n", " ws_mixer.rz(-2 * beta, idx)\n", " ws_mixer.ry(theta, idx)\n", "\n", "ws_mixer.draw(output=\"mpl\", style=\"clifford\")" ] }, { "cell_type": "markdown", "id": "requested-cooking", "metadata": {}, "source": [ "The initial state and mixer operator can then be passed to QAOA." ] }, { "cell_type": "code", "execution_count": 14, "id": "settled-mistress", "metadata": {}, "outputs": [], "source": [ "ws_qaoa_mes = QAOA(\n", " sampler=Sampler(),\n", " optimizer=COBYLA(),\n", " initial_state=init_qc,\n", " mixer=ws_mixer,\n", " initial_point=[0.0, 1.0],\n", ")" ] }, { "cell_type": "code", "execution_count": 15, "id": "wrapped-alberta", "metadata": {}, "outputs": [], "source": [ "ws_qaoa = MinimumEigenOptimizer(ws_qaoa_mes)" ] }, { "cell_type": "code", "execution_count": 16, "id": "aerial-parcel", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objective function value: 16.768932200000002\n", "variable values: x0=0.0, x1=0.0, x2=1.0, x3=1.0, x4=1.0, x5=0.0\n", "status: SUCCESS\n" ] } ], "source": [ "ws_qaoa_result = ws_qaoa.solve(qubo)\n", "print(ws_qaoa_result.prettyprint())" ] }, { "cell_type": "markdown", "id": "armed-chinese", "metadata": {}, "source": [ "### Analysis\n", "\n", "Both results appear to give the same result. However, when we look at the underlying probability distribution we observe that the warm-start QAOA has a much higher probability of sampling the optimal solution." ] }, { "cell_type": "code", "execution_count": 17, "id": "sharp-military", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['001110: value: 16.769, probability: 0.7%',\n", " '011010: value: 15.744, probability: 0.7%',\n", " '001011: value: 14.671, probability: 0.6%',\n", " '101010: value: 14.626, probability: 0.7%',\n", " '010110: value: 14.234, probability: 2.3%',\n", " '100110: value: 13.953, probability: 1.6%',\n", " '000111: value: 13.349, probability: 1.5%',\n", " '110010: value: 12.410, probability: 3.0%',\n", " '010011: value: 12.013, probability: 2.7%',\n", " '100011: value: 11.559, probability: 3.3%']" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def format_qaoa_samples(samples, max_len: int = 10):\n", " qaoa_res = []\n", " for s in samples:\n", " if sum(s.x) == 3:\n", " qaoa_res.append((\"\".join([str(int(_)) for _ in s.x]), s.fval, s.probability))\n", "\n", " res = sorted(qaoa_res, key=lambda x: -x[1])[0:max_len]\n", "\n", " return [(_[0] + f\": value: {_[1]:.3f}, probability: {1e2*_[2]:.1f}%\") for _ in res]\n", "\n", "\n", "format_qaoa_samples(qaoa_result.samples)" ] }, { "cell_type": "code", "execution_count": 18, "id": "political-dependence", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['001110: value: 16.769, probability: 72.8%',\n", " '001011: value: 14.671, probability: 1.2%',\n", " '101010: value: 14.626, probability: 2.6%',\n", " '100110: value: 13.953, probability: 0.3%',\n", " '000111: value: 13.349, probability: 0.1%',\n", " '100011: value: 11.559, probability: 0.0%']" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "format_qaoa_samples(ws_qaoa_result.samples)" ] }, { "cell_type": "markdown", "id": "animal-grave", "metadata": {}, "source": [ "## Warm-start QAOA\n", "\n", "The warm-start features above are available in the Qiskit optimization module as a single optimizer named `WarmStartQAOAOptimizer` which is illustrated below. This solver will solve a QUBO with a warm-start QAOA. It computes $c^*$ by relaxing the problem. This behavior is controlled by setting `relax_for_pre_solver` to `True`." ] }, { "cell_type": "code", "execution_count": 19, "id": "random-happiness", "metadata": {}, "outputs": [], "source": [ "from qiskit_optimization.algorithms import WarmStartQAOAOptimizer" ] }, { "cell_type": "code", "execution_count": 20, "id": "tracked-encoding", "metadata": {}, "outputs": [], "source": [ "qaoa_mes = QAOA(sampler=Sampler(), optimizer=COBYLA(), initial_point=[0.0, 1.0])\n", "ws_qaoa = WarmStartQAOAOptimizer(\n", " pre_solver=CplexOptimizer(), relax_for_pre_solver=True, qaoa=qaoa_mes, epsilon=0.0\n", ")" ] }, { "cell_type": "code", "execution_count": 21, "id": "insured-champagne", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objective function value: 16.768932200000002\n", "variable values: x0=0.0, x1=0.0, x2=1.0, x3=1.0, x4=1.0, x5=0.0\n", "status: SUCCESS\n" ] } ], "source": [ "ws_result = ws_qaoa.solve(qubo)\n", "print(ws_result.prettyprint())" ] }, { "cell_type": "code", "execution_count": 22, "id": "grave-initial", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['001110: value: 16.769, probability: 72.8%',\n", " '001011: value: 14.671, probability: 1.2%',\n", " '101010: value: 14.626, probability: 2.6%',\n", " '100110: value: 13.953, probability: 0.3%',\n", " '000111: value: 13.349, probability: 0.1%',\n", " '100011: value: 11.559, probability: 0.0%']" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "format_qaoa_samples(ws_result.samples)" ] }, { "cell_type": "code", "execution_count": 23, "id": "weird-dispatch", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

Version Information

SoftwareVersion
qiskit0.45.1
qiskit_algorithms0.2.1
qiskit_optimization0.7.0
System information
Python version3.10.13
Python compilerGCC 13.2.1 20230728 (Red Hat 13.2.1-1)
Python buildmain, Aug 28 2023 00:00:00
OSLinux
CPUs12
Memory (Gb)31.056411743164062
Wed Dec 06 10:26:59 2023 CET
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "

This code is a part of Qiskit

© Copyright IBM 2017, 2023.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import qiskit.tools.jupyter\n", "\n", "%qiskit_version_table\n", "%qiskit_copyright" ] }, { "cell_type": "code", "execution_count": null, "id": "ordinary-samuel", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.10.13" } }, "nbformat": 4, "nbformat_minor": 5 }