docs/source/tutorials/qaoa_portfolio_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6e228419-e835-4860-ac03-946eb02d08bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorcircuit as tc\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ddb8a88-779a-43f7-ae14-115463bd87f5",
   "metadata": {},
   "source": [
    "### Solving QUBO problems using QAOA\n",
    "\n",
    "Here we show how to solve a quadratic unconstrained binary optimization (QUBO) problem using QAOA.  Later on below we will extend this to show how to solve binary Markowitz portfolio optimization problems.\n",
    "\n",
    "Consider minimizing the following 2x2 QUBO objective function:\n",
    "\n",
    "$\\begin{pmatrix}x_1 & x_2\\end{pmatrix}\\begin{pmatrix}-5& -2 \\\\-2 & 6\\end{pmatrix}\\begin{pmatrix}x_1\\\\x_2\\end{pmatrix} = -5x_1^2 -4x_1x_2 +6x_2^2$\n",
    "\n",
    "Clearly this is minimized at $(x_1,x_2) = (1,0)$, with corresponding objective function value of $-5$\n",
    "\n",
    "We first convert this to an Ising Hamiltonian by mapping $x_i\\rightarrow \\frac{I-Z_i}{2}$\n",
    "\n",
    "This gives\n",
    "\n",
    "$$-\\frac{5}{4}(I-Z_1)^2 -\\frac{4}{4}(I-Z_1)(I-Z_2) + \\frac{6}{4}(I-Z_2)^2 $$\n",
    "\n",
    "which simplifies to\n",
    "\n",
    "$$-\\frac{1}{2}I +\\frac{7}{2}Z_1   -2Z_2 -Z_1Z_2$$ \n",
    "\n",
    "The $-I/2$ term is simply a constant offset, so we can solve the problem by finding the minimum of \n",
    "$$\\langle \\psi | \\frac{7}{2}Z_1 -2Z_2 -Z_1Z_2 |\\psi\\rangle$$ \n",
    "Note that the minimum should correspond to the computational basis state $|10\\rangle$, and the corresponding true objective function value should be $-4.5$ (ignoring the offset value of $-1/2$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4006848a-1a2f-407a-9f80-63e75ea0d3a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we first manually encode the terms (-7/2) Z_1 - 2 Z_2 - Z_1Z_2 as:\n",
    "pauli_terms = [\n",
    "    [1, 0],\n",
    "    [0, 1],\n",
    "    [1, 1],\n",
    "]  # see the TensorCircuit whitepaper for 'pauli structures'\n",
    "weights = [7 / 2, -2, -1]\n",
    "\n",
    "# see below for a function to generate the pauli terms and weights from the QUBO matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d197cf4a-1bad-4470-a846-998bfe68ba3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we define the QAOA ansatz of depth nlayers\n",
    "def QAOA_from_Ising(params, nlayers, pauli_terms, weights):\n",
    "    nqubits = len(pauli_terms[0])\n",
    "    c = tc.Circuit(nqubits)\n",
    "    for i in range(nqubits):\n",
    "        c.h(i)\n",
    "    for j in range(nlayers):\n",
    "        for k in range(len(pauli_terms)):\n",
    "            term = pauli_terms[k]\n",
    "            index_of_ones = []\n",
    "            for l in range(len(term)):\n",
    "                if term[l] == 1:\n",
    "                    index_of_ones.append(l)\n",
    "            if len(index_of_ones) == 1:\n",
    "                c.rz(index_of_ones[0], theta=2 * weights[k] * params[2 * j])\n",
    "            elif len(index_of_ones) == 2:\n",
    "                c.exp1(\n",
    "                    index_of_ones[0],\n",
    "                    index_of_ones[1],\n",
    "                    unitary=tc.gates._zz_matrix,\n",
    "                    theta=weights[k] * params[2 * j],\n",
    "                )\n",
    "            else:\n",
    "                raise ValueError(\"Invalid number of Z terms\")\n",
    "\n",
    "        for i in range(nqubits):\n",
    "            c.rx(i, theta=params[2 * j + 1])  # mixing terms\n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb38c120-500a-44cc-96ec-fb5ceb11032d",
   "metadata": {},
   "source": [
    "For a general state that is the output of a quantum circuit c, we first define the corresponding loss with respect to the Ising Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "88cec9cf-3ab6-4a4c-b743-ed95ee8c3817",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ising_loss(c, pauli_terms, weights):\n",
    "    loss = 0.0\n",
    "    for k in range(len(pauli_terms)):\n",
    "        term = pauli_terms[k]\n",
    "        index_of_ones = []\n",
    "\n",
    "        for l in range(len(term)):\n",
    "            if term[l] == 1:\n",
    "                index_of_ones.append(l)\n",
    "\n",
    "        if len(index_of_ones) == 1:\n",
    "            delta_loss = weights[k] * c.expectation_ps(z=[index_of_ones[0]])\n",
    "\n",
    "        else:\n",
    "            delta_loss = weights[k] * c.expectation_ps(\n",
    "                z=[index_of_ones[0], index_of_ones[1]]\n",
    "            )\n",
    "\n",
    "        loss += delta_loss\n",
    "\n",
    "    return K.real(loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30a3aa96-7823-4337-9b2f-5170502bb893",
   "metadata": {},
   "source": [
    "For the particular case of a circuit corresponding to a QAOA ansatz this is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "26e4bec2-ce5b-4d0d-9e06-c80fda20619f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QAOA_loss(nlayers, pauli_terms, weights, params):\n",
    "    c = QAOA_from_Ising(params, nlayers, pauli_terms, weights)\n",
    "    return Ising_loss(c, pauli_terms, weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ca8b7a76-5c4f-4ad8-9173-90126b8bb93b",
   "metadata": {},
   "outputs": [],
   "source": [
    "K = tc.set_backend(\"tensorflow\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d5b1897f-cd77-4cd3-bd3b-ece64e6004fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QAOA_solve(pauli_terms, weights, nlayers, iterations):\n",
    "    print_every = 100\n",
    "    learning_rate = 1e-2\n",
    "\n",
    "    loss_val_grad = K.value_and_grad(partial(QAOA_loss, nlayers, pauli_terms, weights))\n",
    "    loss_val_grad_jit = K.jit(loss_val_grad)\n",
    "\n",
    "    opt = K.optimizer(tf.keras.optimizers.Adam(learning_rate))\n",
    "\n",
    "    params = K.implicit_randn(shape=[2 * nlayers], stddev=0.5)\n",
    "    print(f\"initial params: {params}\")\n",
    "    for i in range(iterations):\n",
    "        loss, grads = loss_val_grad_jit(params)\n",
    "        if i % print_every == 0:\n",
    "            print(K.numpy(loss))\n",
    "        params = opt.update(grads, params)\n",
    "\n",
    "    return params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2bc533da-b4f2-4ffb-b486-c65880a30a6d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial params: [ 0.39931756 -0.49578992 -0.22545011 -0.40585193]\n",
      "-2.1728685\n",
      "-4.177884\n",
      "-4.2291102\n",
      "-4.2291365\n",
      "-4.229136\n"
     ]
    }
   ],
   "source": [
    "iterations = 500\n",
    "nlayers = 2\n",
    "final_params = QAOA_solve(pauli_terms, weights, nlayers, iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e93d8cbe-2884-4f0a-80f8-3b600194927b",
   "metadata": {},
   "source": [
    "We note that for nlayers=2 and 500 iterations, the objective function does not in this case (although it depends on the initial parameters)converge to the true value of $-4.5$.  However, the we see below that the final wavefunction does have large overlap with the desired state $|10\\rangle$, so measuring the output of the QAOA algorithm will, with high probability, output the correct answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "294ea9ce-5064-4176-94d0-8dbb7d1707f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_output(c):\n",
    "    n = c._nqubits\n",
    "    N = 2**n\n",
    "    x_label = r\"$\\left|{0:0\" + str(n) + r\"b}\\right>$\"\n",
    "    labels = [x_label.format(i) for i in range(N)]\n",
    "    plt.bar(range(N), c.probability())\n",
    "    plt.xticks(range(N), labels, rotation=70);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fc1353ab-7a7a-4cdc-931c-3b90417c4961",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdrklEQVR4nO3df3DX9X3A8VcSIAEx0UpJhKamFVv1kGChpqCt/SMVN0bLRk/OUqGs02urzhrdlVgkqB3RbuXoTVqunuy6cx70x3Q/YNQay7UrdMzgDp2KqxaJ2oSgZxLiSJZ8v/ujIy4jKF8E33y/PB53+YNPPp/knXvlQ575fD/fb4qy2Ww2AAASKU69AADg1CZGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqVGpF3A0MplMvPLKK3H66adHUVFR6uUAAEchm81GT09PTJo0KYqLj3z9Iy9i5JVXXonq6urUywAAjkFbW1u8733vO+L78yJGTj/99Ij43RdTXl6eeDUAwNHo7u6O6urqoZ/jR5IXMXLooZny8nIxAgB55u1usXADKwCQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqVGpFwCQq5plm1Iv4ZS15+65qZdAAXJlBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFLHFCNr166NmpqaKCsri7q6utixY8db7r9mzZr48Ic/HGPHjo3q6uq4+eab4+DBg8e0YACgsOQcIxs3boyGhoZoamqKnTt3Rm1tbcyZMyf27ds34v4PPvhgLFu2LJqamuKZZ56J+++/PzZu3Bi33XbbO148AJD/co6R1atXx7XXXhtLly6NCy+8MNatWxfjxo2L9evXj7j/tm3b4tJLL43Pfe5zUVNTE1dccUVcffXVb3s1BQA4NeQUI/39/dHa2hr19fVvfoDi4qivr4/t27ePeMzs2bOjtbV1KD5eeOGF2Lx5c/z+7//+ET9PX19fdHd3D3sDAArTqFx23r9/fwwODkZlZeWw7ZWVlfHss8+OeMznPve52L9/f1x22WWRzWZjYGAgvvSlL73lwzTNzc1xxx135LI0ACBPnfBn02zdujVWrVoV3/nOd2Lnzp3xd3/3d7Fp06a46667jnhMY2NjdHV1Db21tbWd6GUCAInkdGVkwoQJUVJSEh0dHcO2d3R0RFVV1YjH3H777XHNNdfEn/zJn0RExEUXXRS9vb1x3XXXxde//vUoLj68h0pLS6O0tDSXpQEAeSqnKyNjxoyJGTNmREtLy9C2TCYTLS0tMWvWrBGPeeONNw4LjpKSkoiIyGazua4XACgwOV0ZiYhoaGiIJUuWxMyZM+OSSy6JNWvWRG9vbyxdujQiIhYvXhyTJ0+O5ubmiIiYN29erF69Oi6++OKoq6uLX//613H77bfHvHnzhqIEADh15RwjCxcujM7OzlixYkW0t7fH9OnTY8uWLUM3te7du3fYlZDly5dHUVFRLF++PF5++eV473vfG/PmzYs///M/P35fBQCQt4qyefBYSXd3d1RUVERXV1eUl5enXg6QWM2yTamXcMrac/fc1Esgjxztz29/mwYASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJDUMcXI2rVro6amJsrKyqKuri527Njxlvu//vrrcf3118fZZ58dpaWl8aEPfSg2b958TAsGAArLqFwP2LhxYzQ0NMS6deuirq4u1qxZE3PmzIndu3fHxIkTD9u/v78/PvWpT8XEiRPjRz/6UUyePDlefPHFOOOMM47H+gGAPJdzjKxevTquvfbaWLp0aURErFu3LjZt2hTr16+PZcuWHbb/+vXr47XXXott27bF6NGjIyKipqbmna0aACgYOT1M09/fH62trVFfX//mBygujvr6+ti+ffuIx/zDP/xDzJo1K66//vqorKyMqVOnxqpVq2JwcPCIn6evry+6u7uHvQEAhSmnGNm/f38MDg5GZWXlsO2VlZXR3t4+4jEvvPBC/OhHP4rBwcHYvHlz3H777fGtb30rvvGNbxzx8zQ3N0dFRcXQW3V1dS7LBADyyAl/Nk0mk4mJEyfG9773vZgxY0YsXLgwvv71r8e6deuOeExjY2N0dXUNvbW1tZ3oZQIAieR0z8iECROipKQkOjo6hm3v6OiIqqqqEY85++yzY/To0VFSUjK07YILLoj29vbo7++PMWPGHHZMaWlplJaW5rI0ACBP5XRlZMyYMTFjxoxoaWkZ2pbJZKKlpSVmzZo14jGXXnpp/PrXv45MJjO07bnnnouzzz57xBABAE4tOT9M09DQEPfdd198//vfj2eeeSa+/OUvR29v79CzaxYvXhyNjY1D+3/5y1+O1157LW666aZ47rnnYtOmTbFq1aq4/vrrj99XAQDkrZyf2rtw4cLo7OyMFStWRHt7e0yfPj22bNkydFPr3r17o7j4zcaprq6On/zkJ3HzzTfHtGnTYvLkyXHTTTfF1772teP3VQAAeasom81mUy/i7XR3d0dFRUV0dXVFeXl56uUAidUs25R6CaesPXfPTb0E8sjR/vz2t2kAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASOqYYmTt2rVRU1MTZWVlUVdXFzt27Diq4zZs2BBFRUUxf/78Y/m0AEAByjlGNm7cGA0NDdHU1BQ7d+6M2tramDNnTuzbt+8tj9uzZ0/ceuut8fGPf/yYFwsAFJ6cY2T16tVx7bXXxtKlS+PCCy+MdevWxbhx42L9+vVHPGZwcDAWLVoUd9xxR3zwgx98RwsGAApLTjHS398fra2tUV9f/+YHKC6O+vr62L59+xGPu/POO2PixInxxS9+8dhXCgAUpFG57Lx///4YHByMysrKYdsrKyvj2WefHfGYf/mXf4n7778//v3f//2oP09fX1/09fUN/bu7uzuXZQIAeeSEPpump6cnrrnmmrjvvvtiwoQJR31cc3NzVFRUDL1VV1efwFUCACnldGVkwoQJUVJSEh0dHcO2d3R0RFVV1WH7P//887Fnz56YN2/e0LZMJvO7TzxqVOzevTvOPffcw45rbGyMhoaGoX93d3cLEgAoUDnFyJgxY2LGjBnR0tIy9PTcTCYTLS0tccMNNxy2//nnnx9PPvnksG3Lly+Pnp6e+Pa3v33EwCgtLY3S0tJclgYA5KmcYiQioqGhIZYsWRIzZ86MSy65JNasWRO9vb2xdOnSiIhYvHhxTJ48OZqbm6OsrCymTp067PgzzjgjIuKw7QDAqSnnGFm4cGF0dnbGihUror29PaZPnx5btmwZuql17969UVzshV0BgKNTlM1ms6kX8Xa6u7ujoqIiurq6ory8PPVygMRqlm1KvYRT1p6756ZeAnnkaH9+u4QBACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABI6phiZO3atVFTUxNlZWVRV1cXO3bsOOK+9913X3z84x+PM888M84888yor69/y/0BgFNLzjGycePGaGhoiKampti5c2fU1tbGnDlzYt++fSPuv3Xr1rj66qvjZz/7WWzfvj2qq6vjiiuuiJdffvkdLx4AyH9F2Ww2m8sBdXV18dGPfjTuvffeiIjIZDJRXV0dN954Yyxbtuxtjx8cHIwzzzwz7r333li8ePFRfc7u7u6oqKiIrq6uKC8vz2W5QAGqWbYp9RJOWXvunpt6CeSRo/35ndOVkf7+/mhtbY36+vo3P0BxcdTX18f27duP6mO88cYb8d///d/xnve854j79PX1RXd397A3AKAw5RQj+/fvj8HBwaisrBy2vbKyMtrb24/qY3zta1+LSZMmDQua/6+5uTkqKiqG3qqrq3NZJgCQR97VZ9PcfffdsWHDhnjooYeirKzsiPs1NjZGV1fX0FtbW9u7uEoA4N00KpedJ0yYECUlJdHR0TFse0dHR1RVVb3lsX/5l38Zd999dzz66KMxbdq0t9y3tLQ0SktLc1kaAJCncroyMmbMmJgxY0a0tLQMbctkMtHS0hKzZs064nHf/OY346677ootW7bEzJkzj321AEDByenKSEREQ0NDLFmyJGbOnBmXXHJJrFmzJnp7e2Pp0qUREbF48eKYPHlyNDc3R0TEPffcEytWrIgHH3wwampqhu4tGT9+fIwfP/44fikAQD7KOUYWLlwYnZ2dsWLFimhvb4/p06fHli1bhm5q3bt3bxQXv3nB5bvf/W709/fHZz/72WEfp6mpKVauXPnOVg8A5L2cX2ckBa8zAvxfXmckHa8zQi5OyOuMAAAcb2IEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJIalXoBABARUbNsU+olnLL23D036ed3ZQQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApI4pRtauXRs1NTVRVlYWdXV1sWPHjrfc/4c//GGcf/75UVZWFhdddFFs3rz5mBYLABSenGNk48aN0dDQEE1NTbFz586ora2NOXPmxL59+0bcf9u2bXH11VfHF7/4xXjiiSdi/vz5MX/+/Hjqqafe8eIBgPyXc4ysXr06rr322li6dGlceOGFsW7duhg3blysX79+xP2//e1vx5VXXhl/9md/FhdccEHcdddd8ZGPfCTuvffed7x4ACD/jcpl5/7+/mhtbY3GxsahbcXFxVFfXx/bt28f8Zjt27dHQ0PDsG1z5syJhx9++Iifp6+vL/r6+ob+3dXVFRER3d3duSz3qExt+slx/5gcnafumJN6CeSpTN8bqZdwyjoR/w8fYq7pnKi5Hvq42Wz2LffLKUb2798fg4ODUVlZOWx7ZWVlPPvssyMe097ePuL+7e3tR/w8zc3Ncccddxy2vbq6OpflcpKrWJN6BUCunLeF6UTPtaenJyoqKo74/pxi5N3S2Ng47GpKJpOJ1157Lc4666woKipKuLKTS3d3d1RXV0dbW1uUl5enXg7HibkWLrMtXGY7smw2Gz09PTFp0qS33C+nGJkwYUKUlJRER0fHsO0dHR1RVVU14jFVVVU57R8RUVpaGqWlpcO2nXHGGbks9ZRSXl7um78AmWvhMtvCZbaHe6srIofkdAPrmDFjYsaMGdHS0jK0LZPJREtLS8yaNWvEY2bNmjVs/4iIn/70p0fcHwA4teT8ME1DQ0MsWbIkZs6cGZdcckmsWbMment7Y+nSpRERsXjx4pg8eXI0NzdHRMRNN90Ul19+eXzrW9+KuXPnxoYNG+Lxxx+P733ve8f3KwEA8lLOMbJw4cLo7OyMFStWRHt7e0yfPj22bNkydJPq3r17o7j4zQsus2fPjgcffDCWL18et912W5x33nnx8MMPx9SpU4/fV3GKKi0tjaampsMe0iK/mWvhMtvCZbbvTFH27Z5vAwBwAvnbNABAUmIEAEhKjAAASYkRACApMQIAJCVG8lQmk0m9BE4Qsy1cZluYzPWd89TeApLNZv3tngJltoXLbAuTuebmpPxDeYysq6srWltb48knn4yampr4wAc+EC+99FKcddZZMW3atBg7dmzqJXKMzLZwmW1hMtfjy5WRPNLa2hr33XdfdHZ2xq5duyKbzcZFF10Uu3fvjueffz6++tWvxi233BITJ06MTCYz7JVwObmZbeEy28JkrseXGMkzPT090d/fH2eddVbs27cv2tvbY8qUKfHzn/88vvvd78aUKVPitttui7POOiv1UsmR2RYusy1M5nr8iJE8cqTHIA9t37p1a9x5553R1tYWq1atij/6oz+KkpKSBCslV/93tiPN2Wzzl/O2MDlnjy8xkod6enqitLQ0xowZM+L777nnnnj00UdjwYIF8aUvfeldXh3H6mhueDPb/OW8LTzO2ePHDax5or+/P7Zu3Rrf/OY3Y8yYMTF27NioqamJz3zmM/GJT3wiIt48MW655ZYYHByMvXv3DtvOye2tZnToMeebb77ZbPOI87awOWePH1dGTnKHvqH/4i/+ItavXx8XXHBBVFVVRX9/f7zwwgvR2dkZtbW1ceONN0ZdXd2wYwcGBmLUqFG++U9Sh2a7devWePXVV+NjH/tYTJw4MUaPHn3YvoODg8Mu8Zrtyc15W5icsyeOGMkTkydPjqamprjuuusiIqK3tzeee+65+MUvfhEbNmyIbDYbDzzwQJx77rlxaKS+4fPDpEmTor29PS688ML49Kc/HXPnzo0Pf/jDceaZZw79Z/bDH/4wXnzxxbj11lv9Z5ZHnLeFyTl7/HmuUR546aWXory8PKZNmza07bTTTouLL744/vRP/zQeeuiheP311+M73/lORPzuPzPf+Pnh6aefjnHjxsUPfvCD+IM/+IP467/+67j88stj/vz58Vd/9Vexa9eu6OzsjKamphgYGIiICL8/5AfnbWFyzp4gWU56b7zxRnbu3LnZ2bNnZ9va2rKZTOawfe6///7sRRddlH3jjTcSrJBj9cgjj2Q//elPZ5944omhbb/61a+yixYtyp522mnZ8ePHZy+//PJsUVFRtrOzM5vNZkecPycf521hcs6eGK6M5IGxY8fG8uXL4+DBg/GVr3wlfvrTn0ZXV9dQdUdE7NmzJ8aNGxdjx471dxLySG1tbdxwww3x/ve/PyJ+95h0XV1dPPDAA3HgwIF44IEHYteuXTFnzpyYMGFCDA4O+u05TzhvC5Nz9sRwz0geOHTT1KOPPhorV66Mbdu2xXnnnRdXXHFFTJo0KbZs2RKvvvpq3HHHHbFgwYLDbpwiP2T/93HlTCYTmUwmRo363ZPdqqqq4p577oklS5aYbR762c9+FitXroxf/OIXztsC5Zx958RIHvqP//iP+Nu//dvYvHlzlJeXxznnnBOLFy+OT33qU6mXxnHwf//zevLJJ2PBggWxe/duv13luaeffjo2bNgQ//RP/xTjx4933uap7NvcjOqcPTZi5CT3+uuvxxNPPBG//e1vo7+/P2bPnh0f+tCHht7f09MT48eP902fhw7Ntr29PTKZTHzsYx+Lc889d9g+fX190dnZGe973/v8fYs8NTAwECUlJcPO0Z6enhg3bpzfmAvQf/3Xf8X+/fujurraOZsDMXISOvQN/Nhjj8Xq1atj8+bN8f73vz+qqqriwIEDcf7558fChQtj7ty5MW7cON/weeStZtvb2xvnn39+LFq0KH7v934vSktLUy+XHBz6jbmzszMGBwdj4sSJw87LbDYbg4ODUVxc7HzNI283V44PMXISOvQD6yMf+UhMnz49Vq5cGaWlpdHa2hq7du2Kf/3Xf42XX345vvCFL8RXvvKV1MslB2ZbuA7N9g//8A+jq6srFi5cGJdddlmcc845MX78+GH7Pv300zEwMDDsab+cnHKZ6zPPPBP9/f1RW1ubaLV5LMEzeDgKr7zySvaMM87I/uY3vznsfc8991y2sbExW1RUlH3ggQfe/cXxjpht4XrppZeyo0ePzs6ePTtbWlqaraioyF511VXZH/zgB9nnn38+e/DgwWwmk8nOnz8/+41vfCObzXraZz4w1xPPtaaT1MGDB+ODH/xg/PjHPz7sfeedd16sWrUqrrvuunjkkUeiv78/wQo5VmZbeLL/e4F527ZtMWvWrHj44Yfj4MGDce+998ZLL70UCxcujE9+8pNx6623xt/8zd/E3//938dnPvOZxKvm7Zjru0eMnIQymUx84AMfiCuvvDK+//3vx49//OM4cODAYftNmTIldu3adcS/AsrJx2wL06GbU2tqauLKK68cmunnP//5+OUvfxn79++PG2+8MR599NFYunRp1NbWxtSpU71M+EnOXN897hk5ibW1tcUtt9wSDz30UMyePTuuueaa+OhHPxqjR4+O3bt3x5133hlXXXVVNDY2Dv0RJvKD2RamgYGB2LdvX0yaNGno30VFRcOeNTNlypS44YYb4qtf/arZ5glzfRekfIyIw23cuDH79NNPD9v2y1/+MrtgwYLsaaedlj3ttNOyF198cfb000/P3nTTTdkDBw5ks1mPT+YDsy1cI812cHAwOzAwkM1mfzfDgYGB7L/9279li4qKsvv27RvazsnLXN89roycRLLZbCxYsCAef/zxOOecc2LBggWxaNGieO973zu0z2OPPRa//e1vY+bMmTFlyhSvU5AnzLZw/f/ZXnXVVbFo0aJ4z3veM7TPoZcEb2tri61bt3plzjxgru8uMXKSee211+I3v/lNPPLII7Fly5Z48cUXY9q0afH5z38+PvvZzx72/Hbf+PnDbAvXSLOtra2Na665JhYsWDDs/oH+/n73AuUJc333iJGT1MDAQLz66qvx1FNPxebNm6OlpSUOHDgQn/jEJ+KP//iP47LLLhva99AI3TCVH8y2cP3/2T722GNx4MCB+OQnPxlf+MIX4tJLLx3a12zzh7meeGIkD/T19UV7e3s8/vjj8Y//+I+xffv2KCsri3nz5sXixYuHvTw8+cVsC5fZFiZzPTHESJ7p7e2NvXv3xs9//vP453/+5/jVr34VCxYsiJUrVw67/4D8Y7aFy2wLk7keP2IkT2Uymeju7o7//M//jK6urrjsssuirKws9bI4Dsy2cJltYTLXd06MAABJeQVWACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApP4H157oSnD133MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = QAOA_from_Ising(final_params, nlayers, pauli_terms, weights)\n",
    "\n",
    "print_output(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d155c5e5-5843-4bba-9edc-595a18cb0c9a",
   "metadata": {},
   "source": [
    "### General Case\n",
    "For the general QUBO case, we wish to minimize\n",
    "$$ x^T Q x$$\n",
    "where $x\\in\\{0,1\\}^n$ and $Q\\in\\mathbb{R}^{n\\times n}$ is a real symmetric matrix.\n",
    "\n",
    "This maps to an Ising Hamiltonian \n",
    "$$\\frac{1}{2}\\left(\\sum_{i=1}^n C_{ii} + \\sum_{i<j}C_{ij}\\right) I - \\frac{1}{2}\\sum_{i=1}^n \\left(\\sum_{j=1}^n C_{ij} \\right)Z_i +\\frac{1}{2}\\sum_{i<j}C_{ij}Z_iZ_j$$\n",
    "\n",
    "Below is a simple function which can perform this mapping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a01ab8a2-d8cd-440d-ad70-7444e8b6a009",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QUBO_to_Ising(Q):\n",
    "    # input is nxn symmetric numpy array corresponding to QUBO matrix Q\n",
    "\n",
    "    n = Q.shape[0]\n",
    "\n",
    "    offset = np.triu(Q, 0).sum() / 2\n",
    "    pauli_terms = []\n",
    "    weights = []\n",
    "\n",
    "    weights = -np.sum(Q, axis=1) / 2\n",
    "\n",
    "    for i in range(n):\n",
    "        term = np.zeros(n)\n",
    "        term[i] = 1\n",
    "        pauli_terms.append(term)\n",
    "\n",
    "    for i in range(n - 1):\n",
    "        for j in range(i + 1, n):\n",
    "            term = np.zeros(n)\n",
    "            term[i] = 1\n",
    "            term[j] = 1\n",
    "            pauli_terms.append(term)\n",
    "\n",
    "            weight = Q[i][j] / 2\n",
    "            weights = np.concatenate((weights, weight), axis=None)\n",
    "\n",
    "    return pauli_terms, weights, offset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c49427-8555-4fb8-aad7-3e374d75659c",
   "metadata": {},
   "source": [
    "### Solving portfolio optimization problems with QAOA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1744c60e-3914-4c92-a653-735052ce8453",
   "metadata": {},
   "source": [
    "In a simple boolean Markowitz portfolio optimization problem, we wish to solve \n",
    "\n",
    "$$\\min_{x\\in\\{0,1\\}^n}\\quad q x^T \\Sigma x - \\mu^T x$$\n",
    "subject to \n",
    "$$ 1^T x = B$$\n",
    "where \n",
    "* $n$: number of assets under consideration\n",
    "* $q > 0 $: risk-appetite\n",
    "* $\\Sigma \\in \\mathbb{R}^{n\\times n}$: covariance matrix of the assets\n",
    "* $\\mu\\in\\mathbb{R}^n$: mean return of the assets\n",
    "* $B$: budget (i.e., total number of assets out of $n$ that can be selected)\n",
    "\n",
    "Our first step is to convert this constrained quadratic programming problem into a QUBO.  We do this by adding a penalty factor $t$ and consider the alternative problem:\n",
    "\n",
    "$$ \\min_{x\\in\\{0,1\\}^n}\\quad q x^T \\Sigma x - \\mu^Tx  + t(1^Tx-B)^2$$\n",
    "\n",
    "The variables in the linear terms $\\mu^Tx = \\mu_1 x_1 + \\mu_2 x_2+\\ldots$ can all be squared (since they are boolean variables), i.e. we can consider\n",
    "\n",
    "$$\\min_{x\\in\\{0,1\\}^n}\\quad q x^T \\Sigma x - \\sum_{i=1}^n\\mu_i x_i^2  + t(1^Tx-B)^2$$\n",
    "which is a QUBO defined by the matrix $Q$ \n",
    "\n",
    "$$ Q = q\\Sigma -\\mu\\begin{pmatrix}1 & \\\\ & 1\\\\ & & \\ddots\\end{pmatrix} + t\\begin{pmatrix}1 -2B & 1 & \\ldots & 1 \\\\\n",
    "1 & 1-2B & 1 & \\ldots \\\\1 & 1 & 1-2B \\\\\n",
    "\\vdots\\end{pmatrix}$$\n",
    "i.e., we wish to mimimize\n",
    "$$ x^T Q X + tB$$\n",
    "and we ignore the constant term $t B$.\n",
    "We can now solve this by QAOA as above.\n",
    "\n",
    "Let us first define a function to convert portfolio data into a QUBO matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3080b901-fb6c-4bda-8348-c96540cbc39a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QUBO_from_portfolio(cov, mean, q, B, t):\n",
    "    # cov: n x n covariance numpy array\n",
    "    # mean: numpy array of means\n",
    "    n = cov.shape[0]\n",
    "    R = np.diag(mean)\n",
    "    S = np.ones((n, n)) - 2 * B * np.diag(np.ones(n))\n",
    "\n",
    "    Q = q * cov - R + t * S\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4cdcb0e-15a2-461c-b487-084488486c67",
   "metadata": {},
   "source": [
    "We can test this using the qiskit_finance package to generate some stock covariance and mean data:\n",
    "\n",
    "*Note that this was tested with qiskit version 0.39.3 and qiskit-finance version 0.3.4.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2168c69e-73ce-4306-8a39-4ddc475acc8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "from qiskit_finance.data_providers import RandomDataProvider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6a7bc671-a496-4cd8-b954-50a280b5dd85",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_assets = 4\n",
    "seed = 123\n",
    "\n",
    "# Generate expected return and covariance matrix from (random) time-series\n",
    "stocks = [(\"TICKER%s\" % i) for i in range(num_assets)]\n",
    "data = RandomDataProvider(\n",
    "    tickers=stocks,\n",
    "    start=datetime.datetime(2016, 1, 1),\n",
    "    end=datetime.datetime(2016, 1, 30),\n",
    "    seed=seed,\n",
    ")\n",
    "data.run()\n",
    "\n",
    "mu = data.get_period_return_mean_vector()\n",
    "sigma = data.get_period_return_covariance_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6dc53d4-7ed0-436d-aa1f-8674c56e756e",
   "metadata": {},
   "source": [
    "Using this mean and covariance data, we can now define our portfolio optimization problem, convert it to a QUBO matrix, and then extract the pauli terms and weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "3f6edcd5-3c10-49fc-86ea-160fc6d3187e",
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 0.5\n",
    "budget = 3  # Note that in this example, there are 4 assets, but a budget of only 3\n",
    "penalty = 3\n",
    "\n",
    "Q = QUBO_from_portfolio(sigma, mu, q, budget, penalty)\n",
    "portfolio_pauli_terms, portfolio_weights, portfolio_offset = QUBO_to_Ising(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "809b90fa-7047-4c88-b862-355da4f58a50",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0000: 21.006979417669037\n",
      "0001: 6.006208358124514\n",
      "0010: 6.006857249462996\n",
      "0011: -2.994037697463167\n",
      "0100: 6.007889613170697\n",
      "0101: -2.992836964752989\n",
      "0110: -2.992179512275861\n",
      "0111: -5.9930299775811875\n",
      "1000: 5.992965725313347\n",
      "1001: -3.007905195444355\n",
      "1010: -3.0070278423618397\n",
      "1011: -6.008022650501182\n",
      "1100: -3.0060506769683\n",
      "1101: -6.006877116105166\n",
      "1110: -6.005991201884008\n",
      "1111: -3.006941528402507\n"
     ]
    }
   ],
   "source": [
    "# Brute force search over classical results for comparison before we run QAOA\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        for k in range(2):\n",
    "            for l in range(2):\n",
    "                x = np.array([i, j, k, l])\n",
    "                print(f\"{i}{j}{k}{l}: {np.dot(x,np.dot(Q,x))- portfolio_offset}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5e69a34-87dc-47b4-aeb2-3b9a03fd0974",
   "metadata": {},
   "source": [
    "We see that, due to the penalty, the lowest energy solutions correspond to 0111, 1011, 1101, 1110, i.e. the portfolios with only 3 assets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "80c1de6c-d3a4-4ea5-922a-bffb59dd1ea3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial params: [ 0.13778198 -0.75357753 -0.01271329 -0.5461785  -0.1501883   0.36323363]\n",
      "-4.204754\n",
      "-5.681799\n",
      "-5.6837077\n",
      "-5.6837044\n",
      "-5.6837053\n",
      "-5.683704\n",
      "-5.6837063\n",
      "-5.683709\n",
      "-5.6837063\n",
      "-5.683705\n"
     ]
    }
   ],
   "source": [
    "iterations = 1000\n",
    "nlayers = 3\n",
    "final_params = QAOA_solve(portfolio_pauli_terms, portfolio_weights, nlayers, iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8a9064e3-0c61-4d2d-baf9-bdf120c0331d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAHACAYAAACBGTONAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoZElEQVR4nO3df3DU9Z3H8dcmSEJA4g80AZoSQCiCQPghMUixHVOCtZZ4SJHzRKNHPZ1c7aSnbawm+GMMAiI6ZYrHlbO2oqBzMteTRrx4caYQQYlWBamgUJBkkxCE8KPk177vDyerW8KPTTa7n919Pma+Y/Ldz37e33c+u+Hld7+78ZiZCQAAwGEJkT4AAACAsyGwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4r1ekDyAUfD6fampqdP7558vj8UT6cAAAwDkwMx09elSDBg1SQsKZz6HERGCpqalRRkZGpA8DAAB0wf79+/WNb3zjjGNiIrCcf/75kr5suH///hE+GgAAcC6ampqUkZHh/3f8TGIisHS8DNS/f38CCwAAUeZcLufgolsAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA83pF+gAAAJGX+YvXemzuvYuu77G5gxGJHnuqpis/03DiDAsAAHAeZ1gAdEs8/J85gMjjDAsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcF6XAsuKFSuUmZmp5ORkZWdna+vWracdu2rVKn3729/WhRdeqAsvvFC5ubmnjDczlZSUaODAgerTp49yc3O1a9eurhwaAACIQUEHlrVr16qoqEilpaWqrq7W+PHjlZeXp/r6+k7HV1ZWat68efq///s/VVVVKSMjQzNmzNCBAwf8YxYvXqxnnnlGK1eu1JYtW9S3b1/l5eXp5MmTXe8MAADEjKADy7Jly7RgwQIVFBRo9OjRWrlypVJSUrR69epOx7/wwgu65557lJWVpVGjRuk//uM/5PP5VFFRIenLsyvLly/Xgw8+qFmzZmncuHF6/vnnVVNTo/Xr13c6Z3Nzs5qamgI2AAAQu4IKLC0tLdq2bZtyc3O/miAhQbm5uaqqqjqnOU6cOKHW1lZddNFFkqQ9e/bI6/UGzJmamqrs7OzTzllWVqbU1FT/lpGREUwbAAAgygQVWA4ePKj29nalpaUF7E9LS5PX6z2nOX7+859r0KBB/oDScb9g5iwuLtaRI0f82/79+4NpAwAARJle4Sy2aNEivfTSS6qsrFRycnKX50lKSlJSUlIIjwwAALgsqDMsAwYMUGJiourq6gL219XVKT09/Yz3Xbp0qRYtWqSNGzdq3Lhx/v0d9+vKnAAAID4EFVh69+6tSZMm+S+YleS/gDYnJ+e091u8eLEeffRRlZeXa/LkyQG3DR06VOnp6QFzNjU1acuWLWecEwAAxI+gXxIqKirSbbfdpsmTJ2vKlClavny5jh8/roKCAknS/PnzNXjwYJWVlUmSnnjiCZWUlGjNmjXKzMz0X5fSr18/9evXTx6PRz/96U/12GOPacSIERo6dKgeeughDRo0SPn5+aHrFAAARK2gA8vcuXPV0NCgkpISeb1eZWVlqby83H/R7L59+5SQ8NWJm1//+tdqaWnRTTfdFDBPaWmpFi5cKEm6//77dfz4cf34xz/W4cOHNW3aNJWXl3frOhcAABA7unTRbWFhoQoLCzu9rbKyMuD7vXv3nnU+j8ejRx55RI888khXDgcAAMQ4/pYQAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnNelwLJixQplZmYqOTlZ2dnZ2rp162nHbt++XbNnz1ZmZqY8Ho+WL19+ypiFCxfK4/EEbKNGjerKoQEAgBgUdGBZu3atioqKVFpaqurqao0fP155eXmqr6/vdPyJEyc0bNgwLVq0SOnp6aedd8yYMaqtrfVvf/rTn4I9NAAAEKOCDizLli3TggULVFBQoNGjR2vlypVKSUnR6tWrOx1/5ZVXasmSJbr55puVlJR02nl79eql9PR0/zZgwIDTjm1ublZTU1PABgAAYldQgaWlpUXbtm1Tbm7uVxMkJCg3N1dVVVXdOpBdu3Zp0KBBGjZsmG655Rbt27fvtGPLysqUmprq3zIyMrpVGwAAuC2owHLw4EG1t7crLS0tYH9aWpq8Xm+XDyI7O1vPPfecysvL9etf/1p79uzRt7/9bR09erTT8cXFxTpy5Ih/279/f5drAwAA9/WK9AFI0nXXXef/ety4ccrOztaQIUO0bt063XnnnaeMT0pKOuPLSwAAILYEdYZlwIABSkxMVF1dXcD+urq6M15QG6wLLrhAI0eO1O7du0M2JwAAiF5BBZbevXtr0qRJqqio8O/z+XyqqKhQTk5OyA7q2LFj+vTTTzVw4MCQzQkAAKJX0C8JFRUV6bbbbtPkyZM1ZcoULV++XMePH1dBQYEkaf78+Ro8eLDKysokfXmh7o4dO/xfHzhwQO+//7769eunyy67TJL0b//2b7rhhhs0ZMgQ1dTUqLS0VImJiZo3b16o+gQAAFEs6MAyd+5cNTQ0qKSkRF6vV1lZWSovL/dfiLtv3z4lJHx14qampkYTJkzwf7906VItXbpU11xzjSorKyVJn3/+uebNm6fGxkZdcsklmjZtmt5++21dcskl3WwPAADEgi5ddFtYWKjCwsJOb+sIIR0yMzNlZmec76WXXurKYQAAgDjB3xICAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA87oUWFasWKHMzEwlJycrOztbW7duPe3Y7du3a/bs2crMzJTH49Hy5cu7PScAAIgvQQeWtWvXqqioSKWlpaqurtb48eOVl5en+vr6TsefOHFCw4YN06JFi5Senh6SOQEAQHwJOrAsW7ZMCxYsUEFBgUaPHq2VK1cqJSVFq1ev7nT8lVdeqSVLlujmm29WUlJSSOYEAADxJajA0tLSom3btik3N/erCRISlJubq6qqqi4dQFfmbG5uVlNTU8AGAABiV1CB5eDBg2pvb1daWlrA/rS0NHm93i4dQFfmLCsrU2pqqn/LyMjoUm0AABAdovJdQsXFxTpy5Ih/279/f6QPCQAA9KBewQweMGCAEhMTVVdXF7C/rq7utBfU9sScSUlJp70eBgAAxJ6gzrD07t1bkyZNUkVFhX+fz+dTRUWFcnJyunQAPTEnAACILUGdYZGkoqIi3XbbbZo8ebKmTJmi5cuX6/jx4yooKJAkzZ8/X4MHD1ZZWZmkLy+q3bFjh//rAwcO6P3331e/fv102WWXndOcAAAgvgUdWObOnauGhgaVlJTI6/UqKytL5eXl/otm9+3bp4SEr07c1NTUaMKECf7vly5dqqVLl+qaa65RZWXlOc0JAADiW9CBRZIKCwtVWFjY6W0dIaRDZmamzKxbcwIAgPgWle8SAgAA8YXAAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnNcr0gcAAABCI/MXr/XY3HsXXd9jc58LzrAAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAzutSYFmxYoUyMzOVnJys7Oxsbd269YzjX375ZY0aNUrJyckaO3asNmzYEHD77bffLo/HE7DNnDmzK4cGAABiUNCBZe3atSoqKlJpaamqq6s1fvx45eXlqb6+vtPxmzdv1rx583TnnXfqvffeU35+vvLz8/XRRx8FjJs5c6Zqa2v924svvti1jgAAQMwJOrAsW7ZMCxYsUEFBgUaPHq2VK1cqJSVFq1ev7nT8008/rZkzZ+q+++7T5ZdfrkcffVQTJ07Ur371q4BxSUlJSk9P928XXnhh1zoCAAAxJ6jA0tLSom3btik3N/erCRISlJubq6qqqk7vU1VVFTBekvLy8k4ZX1lZqUsvvVTf+ta3dPfdd6uxsfG0x9Hc3KympqaADQAAxK6gAsvBgwfV3t6utLS0gP1paWnyer2d3sfr9Z51/MyZM/X888+roqJCTzzxhN566y1dd911am9v73TOsrIypaam+reMjIxg2gAAAFGmV6QPQJJuvvlm/9djx47VuHHjNHz4cFVWVuraa689ZXxxcbGKior83zc1NRFaAACIYUGdYRkwYIASExNVV1cXsL+urk7p6emd3ic9PT2o8ZI0bNgwDRgwQLt37+709qSkJPXv3z9gAwAAsSuowNK7d29NmjRJFRUV/n0+n08VFRXKycnp9D45OTkB4yXpjTfeOO14Sfr888/V2NiogQMHBnN4AAAgRgX9LqGioiKtWrVKv/3tb/Xxxx/r7rvv1vHjx1VQUCBJmj9/voqLi/3j7733XpWXl+vJJ5/Uzp07tXDhQr377rsqLCyUJB07dkz33Xef3n77be3du1cVFRWaNWuWLrvsMuXl5YWoTQAAEM2CvoZl7ty5amhoUElJibxer7KyslReXu6/sHbfvn1KSPgqB02dOlVr1qzRgw8+qAceeEAjRozQ+vXrdcUVV0iSEhMT9cEHH+i3v/2tDh8+rEGDBmnGjBl69NFHlZSUFKI2AQBANOvSRbeFhYX+MyR/r7Ky8pR9c+bM0Zw5czod36dPH73++utdOQwAABAn+FtCAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJxHYAEAAM4jsAAAAOcRWAAAgPMILAAAwHkEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgPAILAABwHoEFAAA4j8ACAACcR2ABAADOI7AAAADnEVgAAIDzCCwAAMB5BBYAAOA8AgsAAHAegQUAADiPwAIAAJzXK9IHACB0Mn/xWo/NvXfR9T02NwCcDYEFEdNT/7i69A8rASI2hHsd4+FxEw89IrR4SQgAADiPwAIAAJxHYAEAAM4jsAAAAOd1KbCsWLFCmZmZSk5OVnZ2trZu3XrG8S+//LJGjRql5ORkjR07Vhs2bAi43cxUUlKigQMHqk+fPsrNzdWuXbu6cmgAACAGBR1Y1q5dq6KiIpWWlqq6ulrjx49XXl6e6uvrOx2/efNmzZs3T3feeafee+895efnKz8/Xx999JF/zOLFi/XMM89o5cqV2rJli/r27au8vDydPHmy650BAICYEXRgWbZsmRYsWKCCggKNHj1aK1euVEpKilavXt3p+KefflozZ87Ufffdp8svv1yPPvqoJk6cqF/96leSvjy7snz5cj344IOaNWuWxo0bp+eff141NTVav359t5oDAACxIajPYWlpadG2bdtUXFzs35eQkKDc3FxVVVV1ep+qqioVFRUF7MvLy/OHkT179sjr9So3N9d/e2pqqrKzs1VVVaWbb775lDmbm5vV3Nzs//7IkSOSpKampmDaOWdXlL7eI/NK0kcP54W1Zrjrnammr/lEj9Q73eMglnqUOu8z3PUiUZN1DH29SNSMhx57sqZLPYZiTjM7+2ALwoEDB0ySbd68OWD/fffdZ1OmTOn0Puedd56tWbMmYN+KFSvs0ksvNTOzTZs2mSSrqakJGDNnzhz70Y9+1OmcpaWlJomNjY2NjY0tBrb9+/efNYNE5SfdFhcXB5y18fl8OnTokC6++GJ5PJ6IHVdTU5MyMjK0f/9+9e/fn5pRWi8SNekxNmrSY2zUjIceI1Xz75mZjh49qkGDBp11bFCBZcCAAUpMTFRdXV3A/rq6OqWnp3d6n/T09DOO7/hvXV2dBg4cGDAmKyur0zmTkpKUlJQUsO+CCy4IppUe1b9//7AvfjzUpMfYqEmPsVGTHqkZKqmpqec0LqiLbnv37q1JkyapoqLCv8/n86miokI5OTmd3icnJydgvCS98cYb/vFDhw5Venp6wJimpiZt2bLltHMCAID4EvRLQkVFRbrttts0efJkTZkyRcuXL9fx48dVUFAgSZo/f74GDx6ssrIySdK9996ra665Rk8++aSuv/56vfTSS3r33Xf17//+75Ikj8ejn/70p3rsscc0YsQIDR06VA899JAGDRqk/Pz80HUKAACiVtCBZe7cuWpoaFBJSYm8Xq+ysrJUXl6utLQ0SdK+ffuUkPDViZupU6dqzZo1evDBB/XAAw9oxIgRWr9+va644gr/mPvvv1/Hjx/Xj3/8Yx0+fFjTpk1TeXm5kpOTQ9Bi+CQlJam0tPSUl6uoGV31IlGTHmOjJj3GRs146DFSNbvDY3Yu7yUCAACIHP6WEAAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsPcDM5PP5Ar6PtZr0GBs16TE2atIjNaOlXnfwOSxh1N7ersTExJiuSY+xUZMeY6MmPVIzWuqdCwJLCNXX1+vdd9/VX//6V40ZM0YXXHCBamtrlZGRodGjR8dETXqkx2ipSY/0SE136oUCgSWE/vjHP2rNmjX64osv9M4772jQoEHKyMhQdXW1Wlpa9PDDD+v2229Xnz595PP5Av6EQVdt2LBBL774Ythq0iM90iM90mPoa4a7z0j02G2GkKqtrbXDhw+bmdnOnTuturrajh49aitXrrQZM2bY4sWLQ17T6/XakSNHzMxs165d9uc//7lHa9IjPXYVPdJjV8RDj2bh7zMSPXYHgSVEfD6f/+u2tjZrbW0NuL25udleeOEFGzt2rF199dW2efPmHj+mtrY2W7NmTchq0iM99hR6DB49xkaPZu712RM9hgIvCYXQzp07tXjxYp08eVKJiYnKyMjQDTfcoJycHP+YQ4cOqaSkRLt379a9996r6667LqTH0Nmpu0OHDumXv/yl9uzZ0+2a9PhVPXrsHnqkx3MRDz1Kke8zHD12lwMvSkW3jrxXUVGh2bNna/v27WptbVV7e7s2bdqkf/7nf9a8efNUWVkpSbrooov02GOP6Zvf/Ka8Xm/AHMHW/PTTT/Xf//3fqqmpUVtbmyQFPOB8Pp/a2tp00UUX6cknn+xyTXqkR3qkR3oMbY+R6DMSPYZUT5/CiXXt7e1mZnbttdfaXXfd5T+Vd+TIEdu6dastX77cvvvd79qUKVNs+/btAfdta2vrVs3vf//75vF4bNKkSbZw4UJ7++23rbGx0X+7mVl1dbU98sgj3apJj/RIj2evSY9fosfga4arz0j0GEoElhBobW21CRMm2Nq1azu9vba21qZMmWI33nij+Xy+kCx8c3OzDR8+3BYvXmw/+clPbMCAAXbeeefZd7/7XVuxYoV99NFHduzYMZs7d679wz/8g5l17wFHj/TYVfRIj10RDz2ahb/PSPQYKgSWEGhtbbW7777bRo0aZbt37w64gKrDxo0bbcSIEbZnz56Q1Pzggw9s5syZ9sc//tG/r7Ky0mbPnm1JSUmWmppqN954oyUkJFh5ebmZdf+XBz3SY1fQIz12RTz0aBb+PiPRY6gQWEJkx44dNnXqVJs1a5Zt3LjRmpqaAk6vvfrqqzZgwICQ1WtsbLTXXnvNdu3aZWYWUMvM7MUXX7RvfvObNnjw4JDVpEd67Ap6pMeuiocezcLbZ6R6DAUCSwh0LPibb75pOTk55vF4bOzYsXb//ffbs88+a//4j/9oEydOtF/84hdmZqe8ZS0Utc1OfTvcjBkz7K677gpJzY7UX1lZadOmTQtrj1/Xkz2yjqxjd7GOrGOwXFjLnu4xVAgsPaC6utp+8pOf2GWXXWZXXHGFXXvttfbss89aU1OTmVmnp/xC4eun7WpqaiwzM9Pee+89Mzs1RXfX+++/bz/72c9s5MiRPdbj2ebo6R5ZR9axO1hH1rErenotI72O3cHnsHSTz+fTxx9/LK/Xqy+++EKTJ09WZmam//b6+nqdf/756tOnT4/UPHbsmCZOnKiMjIyAMa2trdqzZ49GjhwZsrod8yYmJga8Be7gwYNKSUlRSkpKSGudTXNzs/bu3atvfetb3Z6LdWQdQ1GTdWQdu8qVtQzlOoZcpBNTNOpIqHv37rV/+Zd/sYSEBEtLS7Mrr7zSxo0bZ3PnzrVXXnnFTpw4ETC+J2pOmTLFX3P9+vV28uTJbtf6uhMnTlhdXd0pKdvn81lra2uP/F/N6WqGGuvIOoa6JuvYdfG0jmbhX8twrWNP4gxLF3T82e077rhDn3zyiR5//HENGTJE27Zt0/bt2/XOO++orq5Os2fP1v333x+VNTvqPfzww9q8ebPmzZunq666ShkZGerbt2/A2Lq6OtXX12vs2LFhq1lfXy+v16tx48Z1ux7r+CXW0d2arGNsrOPXa4ZrLcO9jj0q0okpmqWnpwe8NazDzp07rbS01BISEuyJJ56I6poXX3yxDRs2zHr37m0XX3yx3Xrrrfbqq6/a3r17rbm52czMFi5caLfeemvU1mQdWcdoqck6xsY6moX/5xqJx06oEVi6qLGx0XJzc+2Xv/zlacc89NBD9v3vf9//1zCjrea2bdtswoQJ9uGHH5rP57NVq1bZpEmTzOPx2LBhw+xnP/uZvfbaa3bRRRfZqlWrzKz779cPd03WkXWMlpqs45eifR3Nwv9zjcRjpycQWLrhqaeesszMTFuzZo3/Cu6ve/311+2SSy4J6WuR4ay5fft2KykpsU2bNgXsr62ttYcfftiGDh1qHo/HUlJSQtZjJGqyjqxjNNRkHb8U7etoFv6fayTWsScQWLrh0KFDVlBQYMnJyZabm2u/+93v7C9/+Yvt27fPtmzZYnPmzLFZs2aZWejSajhrtre32+eff24tLS1mZtbS0nLKBVszZsywW265xcxC8379SNT84osv7I477gjrOoazZnt7ux04cMD/swrXOoa7Js9Hno/RUjPcz49IPB97AoGli76eQt966y37wQ9+YMnJydavXz/LycmxlJQUy8/P9//BqlBcmR2Jmh3a29v983VcxV5bW2sej8f+8Ic/mFnoTyGGu+amTZssPz/f+vTpE5afaaRqdgjXOoajJs9Hno/RWvPrwv2cjMTvgO7gXULdcPjwYSUnJys5OVmSdPLkSf3v//6vdu/ercmTJysrK0v9+vWL6ppHjx5VUlKSevfu7d/n8/mUkJCghoYG/ed//mfIrp6PRE0z88/t8XgkSX/729/01ltv6ZNPPtHEiRND/jMNd00z89fpTE+sYyRq8nzk+RhNNcP5/IjE87EnEFiC1NLSosrKSi1evFi9e/dWnz59NGTIEP3whz/Ud77znVPGn+2B4mLNzuplZmZq1qxZmj59esDYjrfMdVckav49n88nj8fT7fVyvWZnWltbdd5550VdTZ6PPB+jvebphPs5GYnfAcEisJyjjgS+ZMkSrV69WpdffrnS09PV0tKiTz/9VA0NDcrKytI999yjqVOnRmXNM9X77LPP1NDQoPHjx+tf//VflZ2dHYIOw1+zo15lZaUaGxt11VVX6dJLLz3lidrW1iYzC8kTONw1z7WeFLp/bMJdk+cjz8dYrCmF9vkRzt8BYRHO159iwaBBg+zZZ5/1f3/s2DGrrq62p59+2nJycuyqq66y3bt3R3XNeOhx4MCB5vF4bMyYMVZcXGx/+tOfrKGhIeC123Xr1tmSJUuitmY89BgPj9V46DEeHquRqBmJHnsSgSUI+/fvt1GjRllVVVWnt3u9Xhs1apQVFRVFbc146HH79u02fPhwe/nll+3nP/+5paenW2Jiol199dX21FNP2Z///Gerr6+3yy+/3MrKysys+xfahbtmPPQYD4/VeOgxHh6rkagZiR57GoElCCdOnLDrr7/epk6davv37+/0/eq/+c1vbOzYsf6/dxFtNeOhx40bN9oPf/hD/18jNTN7++237ZZbbrG+fftav3797JprrjGPx2MNDQ1m1v2/WxLumvHQYzw8VuOhx3h4rEaiZiR67GkEliBVVVXZxIkT7YYbbrDXX3/dDh8+HPCe9Yceesiys7PNLHRpNdw1Y73Huro627hxozU2NnY63/r16+3CCy+0mTNnmllo3uIX7prx0KNZ7D9WI1Ev3DXj5bEaD78DehqBJQgdC/7GG2/Y1VdfbR6Px0aOHGmFhYX2+OOP2/Tp023MmDH2yiuvmFnoPmAonDXjocev6/g/ivb29oBfyGlpafbcc8+FvF4kasZ6j2+++aZNnz49rI+bcNeMhx6/LpyP1XiqGYkeQ4l3CXXD9u3b9cILL2jDhg3q37+/hgwZovnz5+t73/tezNSMhx6lwCvlP/zwQ82ePVt/+ctfevTtjeGuGes97tixQy+99JL+53/+R/369QvL4ybcNWOxRzvL27574nETDzUj0WNPI7Cco8OHD+u9995TbW2tWlpaNHXqVI0cOdJ/+9GjR9WvX7+QLny4a8ZTj16vVz6fT1dddZWGDx8eMKa5uVkNDQ36xje+4X97YDTVjIcev66trU2JiYkBj5GjR48qJSWlx96uGe6a8dDj6fztb3/TwYMHlZGREdLHTbzXjESP3Ra5kzvu63ipoqKiwq6//nrzeDw2ZMgQy87OtjFjxtjs2bNt3bp1dvz48YDx0VQz3nu84oor7KabbrJXX33VTp482b3GIlgzHnrseJmpvr7eamtrT3lcdHzMeCjf6RDumvTYM+KhZiR6DDfOsJxBR+rs+GjmhQsXKikpSdu2bdMHH3ygLVu26MCBA7r99tt1zz33RGVNeqTHaKnZUe/GG2/UkSNHNHfuXE2bNk1Dhgw55WPTd+zYoba2No0bNy6qatJjYL2PP/5YLS0tGj9+fJfrxUvNSPQYdpFOTK6rqamxCy64wPbs2XPKbZ988okVFxebx+Ox3//+91Fbkx7pMVpqfv7553beeefZ1KlTLSkpyVJTU+1HP/qRrVu3zj799FM7efKk+Xw+y8/Pt8cee8zMuv9WzXDXpMfY6DESNSPRYzgRWM7is88+s4kTJ9rSpUtPO+auu+6y+fPnW3Nzc1TWpMfQ14tEzVjuseOX6rp162z69OlWX19vZma/+93vbOrUqebxeCwjI8MKCwvtueeeM4/HYx9++GHAfV2vSY+x0WMkakaix0ggsJxBx2uADzzwgI0dO9ZeeeUVO3r06CnjlixZYllZWVFZkx5DXy8SNeOhRzOzrVu32uOPP26fffZZwP7GxkZbvHixjRo1yjwej79eKH4Zh7smPcZGj5GoGYkew4nAcg727dtnc+bMsV69etn06dNt1apV9v7779v27dvtv/7rvywrK8sef/xxM7OA97hHU016pMdoqNna2moHDhwI+P7vPz9i+PDh9tRTT4WkXiRq0mPo68VLzUj0GE4EltNYu3at7dixI2Dfpk2bbPbs2da3b1/r27evTZgwwc4//3y799577dixY2bWvcQa7pr0SI9d5UKP7e3t/l/GPp/P2tra7J133jGPx+M/JR7qHnuyJj3GRo+RqBmJHiOBwNIJn89nN954o2VkZNi0adPsqaee8i9wh4qKCvv9739vO3fuDMknBIa7Jj2Gvl4kasZjj88884z/48Y7tLW1WXt7u+3duzckn9wZ7pr0GPp68VIzEj1GCm9rPo1Dhw5pz5492rhxo8rLy/XXv/5V48aN0z/90z/ppptuOuVDdr7+iZ7RUpMe6TFaanZWb/z48br11ls1e/bsgA83a2lpUe/evbtcK1I16TE2eoxEzUj0GAkElrNoa2tTY2OjPvroI23YsEEVFRU6duyYpk+frjvuuEPTpk3zj+34UXb3U1nDXZMe6TFaav59vTfffFPHjh3Td77zHd1+++26+uqrQ1ovEjXpMTZ6jETNSPQYTgSWIDQ3N8vr9erdd9/VH/7wB1VVVSk5OVk33HCD5s+fH/CR8tFakx7pMVpq0iM9UtOdemERztefYsmxY8dsx44dtnLlSps1a5alpaXZPffcc8pr+dFckx57Bj1Gf71I1KTHnhEPNSPRY0/gDEs3+Xw+NTU1adeuXTpy5IimTZum5OTkmKpJj7FRkx5joyY9UjNa6oUagQUAADgvCv6eNAAAiHcEFgAA4DwCCwAAcB6BBQAAOI/AAgAAnEdgAQAAziOwAAAA5xFYAACA8wgsAADAeQQWAADgvP8HW9a5h3cUX+YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_final = QAOA_from_Ising(\n",
    "    final_params, nlayers, portfolio_pauli_terms, portfolio_weights\n",
    ")\n",
    "print_output(c_final)"
   ]
  }
 ],
 "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.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}