{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# This cell is added by sphinx-gallery\n",
"# It can be customized to whatever you like\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3-qubit Ising model in PyTorch\n",
"\n",
"\n",
"\n",
"The interacting spins with variable coupling strengths of an [Ising model](https://en.wikipedia.org/wiki/Ising_model) can be used to simulate various machine learning concepts like [Hopfield networks](https://en.wikipedia.org/wiki/Hopfield_network) and [Boltzmann machines](https://en.wikipedia.org/wiki/Boltzmann_machine) ([Schuld & Petruccione (2018)](https://www.springer.com/gp/book/9783319964232)). They also closely imitate the underlying mathematics of a subclass of computational problems called [Quadratic Unconstrained Binary Optimization (QUBO)](https://en.wikipedia.org/wiki/Quadratic_unconstrained_binary_optimization) problems.\n",
"\n",
"Ising models are commonly encountered in the subject area of adiabatic quantum computing. Quantum annealing algorithms (for example, as performed on a D-wave system) are often used to find low-energy configurations of Ising problems. The optimization landscape of the Ising model is non-convex, which can make finding global minima challenging. In this tutorial, we get a closer look at this phenomenon by applying gradient descent techniques to a toy Ising model.\n",
"\n",
"## PennyLane implementation\n",
"\n",
"\n",
"This basic tutorial optimizes a 3-qubit Ising model using the PennyLane `default.qubit` device with PyTorch. In the absence of external fields, the Hamiltonian for this system is given by:\n",
"\n",
"$$H=-\\sum_{<i,j>} J_{ij} \\sigma_i \\sigma_{j},$$\n",
"\n",
"where each spin can be in the +1 or -1 spin state and $J_{ij}$ are the nearest-neighbour coupling strengths.\n",
"\n",
"For simplicity, the first spin can be assumed to be in the \\\"up\\\" state (+1 eigenstate of Pauli-Z operator) and the coupling matrix can be set to $J = [1,-1]$. The rotation angles for the other two spins are then\n",
"optimized so that the energy of the system is minimized for the given couplings.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import torch\n",
"from torch.autograd import Variable\n",
"import pennylane as qml\n",
"from pennylane import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A three-qubit quantum circuit is initialized to represent the three\n",
"spins:\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dev = qml.device(\"default.qubit\", wires=3)\n",
"\n",
"@qml.qnode(dev, interface=\"torch\") \n",
"def circuit(p1, p2):\n",
" # We use the general Rot(phi,theta,omega,wires) single-qubit operation\n",
" qml.Rot(p1[0], p1[1], p1[2], wires=1)\n",
" qml.Rot(p2[0], p2[1], p2[2], wires=2)\n",
" return [qml.expval(qml.PauliZ(i)) for i in range(3)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The cost function to be minimized is defined as the energy of the spin\n",
"configuration:\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def cost(var1, var2):\n",
" # the circuit function returns a numpy array of Pauli-Z expectation values\n",
" spins = circuit(var1, var2)\n",
"\n",
" # the expectation value of Pauli-Z is +1 for spin up and -1 for spin down\n",
" energy = -(1 * spins[0] * spins[1]) - (-1 * spins[1] * spins[2])\n",
" return energy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sanity check\n",
"\n",
"\n",
"Let\\'s test the functions above using the\n",
"$[s_1, s_2, s_3] = [1, -1, -1]$ spin configuration and the given\n",
"coupling matrix. The total energy for this Ising model should be:\n",
"\n",
"$$H = -1(J_1 s_1 \\otimes s_2 + J_2 s_2 \\otimes s3) = 2$$\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Energy for [1, -1, -1] spin configuration: tensor(2.0000, dtype=torch.float64)\n"
]
}
],
"source": [
"test1 = torch.tensor([0, np.pi, 0])\n",
"test2 = torch.tensor([0, np.pi, 0])\n",
"\n",
"cost_check = cost(test1, test2)\n",
"print(\"Energy for [1, -1, -1] spin configuration:\", cost_check)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Random initialization\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Randomly initialized angles:\n",
"tensor([1.9632, 2.6022, 2.3277], dtype=torch.float64, requires_grad=True)\n",
"tensor([0.6521, 2.8474, 2.4300], dtype=torch.float64, requires_grad=True)\n",
"Corresponding cost before optimization:\n",
"tensor(1.6792, dtype=torch.float64, grad_fn=<SubBackward0>)\n"
]
}
],
"source": [
"torch.manual_seed(56)\n",
"p1 = Variable((np.pi * torch.rand(3, dtype=torch.float64)), requires_grad=True)\n",
"p2 = Variable((np.pi * torch.rand(3, dtype=torch.float64)), requires_grad=True)\n",
"\n",
"var_init = [p1, p2]\n",
"cost_init = cost(p1, p2)\n",
"\n",
"print(\"Randomly initialized angles:\")\n",
"print(p1)\n",
"print(p2)\n",
"print(\"Corresponding cost before optimization:\")\n",
"print(cost_init)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimization\n",
"\n",
"\n",
"Now we use the PyTorch gradient descent optimizer to minimize the cost:\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Energy after step 5: 0.6846474 | Angles: [array([1.96323939, 1.93604492, 2.32767565]), array([0.65212549, 2.73080219, 2.4299563 ])] \n",
"\n",
"Energy after step 10: -1.0138530 | Angles: [array([1.96323939, 1.0136468 , 2.32767565]), array([0.65212549, 2.73225282, 2.4299563 ])] \n",
"\n",
"Energy after step 15: -1.8171995 | Angles: [array([1.96323939, 0.38483073, 2.32767565]), array([0.65212549, 2.85992571, 2.4299563 ])] \n",
"\n",
"Energy after step 20: -1.9686584 | Angles: [array([1.96323939, 0.13026452, 2.32767565]), array([0.65212549, 2.97097572, 2.4299563 ])] \n",
"\n",
"Energy after step 25: -1.9930403 | Angles: [array([1.96323939, 0.04302756, 2.32767565]), array([0.65212549, 3.04042222, 2.4299563 ])] \n",
"\n",
"Energy after step 30: -1.9980133 | Angles: [array([1.96323939, 0.01413292, 2.32767565]), array([0.65212549, 3.08179844, 2.4299563 ])] \n",
"\n",
"Energy after step 35: -1.9993550 | Angles: [array([1.96323939, 0.00463472, 2.32767565]), array([0.65212549, 3.10627578, 2.4299563 ])] \n",
"\n",
"Energy after step 40: -1.9997802 | Angles: [array([1.96323939e+00, 1.51911413e-03, 2.32767565e+00]), array([0.65212549, 3.12073668, 2.4299563 ])] \n",
"\n",
"Energy after step 45: -1.9999239 | Angles: [array([1.96323939e+00, 4.97829828e-04, 2.32767565e+00]), array([0.65212549, 3.12927707, 2.4299563 ])] \n",
"\n",
"Energy after step 50: -1.9999735 | Angles: [array([1.96323939e+00, 1.63134183e-04, 2.32767565e+00]), array([0.65212549, 3.13432035, 2.4299563 ])] \n",
"\n",
"Energy after step 55: -1.9999908 | Angles: [array([1.96323939e+00, 5.34564150e-05, 2.32767565e+00]), array([0.65212549, 3.13729842, 2.4299563 ])] \n",
"\n",
"Energy after step 60: -1.9999968 | Angles: [array([1.96323939e+00, 1.75166673e-05, 2.32767565e+00]), array([0.65212549, 3.13905695, 2.4299563 ])] \n",
"\n",
"Energy after step 65: -1.9999989 | Angles: [array([1.96323939e+00, 5.73986944e-06, 2.32767565e+00]), array([0.65212549, 3.14009534, 2.4299563 ])] \n",
"\n",
"Energy after step 70: -1.9999996 | Angles: [array([1.96323939e+00, 1.88084132e-06, 2.32767565e+00]), array([0.65212549, 3.14070851, 2.4299563 ])] \n",
"\n",
"Energy after step 75: -1.9999999 | Angles: [array([1.96323939e+00, 6.16314188e-07, 2.32767565e+00]), array([0.65212549, 3.14107057, 2.4299563 ])] \n",
"\n",
"Energy after step 80: -2.0000000 | Angles: [array([1.96323939e+00, 2.01953845e-07, 2.32767565e+00]), array([0.65212549, 3.14128437, 2.4299563 ])] \n",
"\n",
"Energy after step 85: -2.0000000 | Angles: [array([1.96323939e+00, 6.61762372e-08, 2.32767565e+00]), array([0.65212549, 3.14141062, 2.4299563 ])] \n",
"\n",
"Energy after step 90: -2.0000000 | Angles: [array([1.96323939e+00, 2.16846296e-08, 2.32767565e+00]), array([0.65212549, 3.14148516, 2.4299563 ])] \n",
"\n",
"Energy after step 95: -2.0000000 | Angles: [array([1.96323939e+00, 7.10561943e-09, 2.32767565e+00]), array([0.65212549, 3.14152918, 2.4299563 ])] \n",
"\n",
"Energy after step 100: -2.0000000 | Angles: [array([1.96323939e+00, 2.32836938e-09, 2.32767565e+00]), array([0.65212549, 3.14155517, 2.4299563 ])] \n",
"\n"
]
}
],
"source": [
"opt = torch.optim.SGD(var_init, lr=0.1)\n",
"\n",
"def closure():\n",
" opt.zero_grad()\n",
" loss = cost(p1, p2)\n",
" loss.backward()\n",
" return loss\n",
"\n",
"var_pt = [var_init]\n",
"cost_pt = [cost_init]\n",
"x = [0]\n",
"\n",
"for i in range(100):\n",
" opt.step(closure)\n",
" if (i + 1) % 5 == 0:\n",
" x.append(i)\n",
" p1n, p2n = opt.param_groups[0][\"params\"]\n",
" costn = cost(p1n, p2n)\n",
" var_pt.append([p1n, p2n])\n",
" cost_pt.append(costn)\n",
"\n",
" # for clarity, the angles are printed as numpy arrays\n",
" print(\"Energy after step {:5d}: {: .7f} | Angles: {}\".format(\n",
" i+1, costn, [p1n.detach().numpy(), p2n.detach().numpy()]),\"\\n\"\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Note**\n",
"\n",
"When using the *PyTorch* optimizer, keep in mind that:\n",
"\n",
"1. `loss.backward()` computes the gradient of the cost function with\n",
" respect to all parameters with `requires_grad=True`.\n",
"2. `opt.step()` performs the parameter update based on this *current*\n",
" gradient and the learning rate.\n",
"3. `opt.zero_grad()` sets all the gradients back to zero. It\\'s\n",
" important to call this before `loss.backward()` to avoid the\n",
" accumulation of gradients from multiple passes.\n",
"\n",
"Hence, its standard practice to define the `closure()` function that\n",
"clears up the old gradient, evaluates the new gradient and passes it\n",
"onto the optimizer in each step.\n",
"\n",
"The minimum energy is -2 for the spin configuration $[s_1, s_2, s_3] = [1, 1, -1]$ which corresponds to $(\\phi, \\theta, \\omega) = (0, 0, 0)$ for the second spin and $(\\phi, \\theta, \\omega) = (0, \\pi, 0)$ for the third spin. Note that gradient descent optimization might not find this global minimum due to the non-convex cost function, as is shown in the next section.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimized angles:\n",
"tensor([1.9632e+00, 2.3284e-09, 2.3277e+00], dtype=torch.float64,\n",
" requires_grad=True)\n",
"tensor([0.6521, 3.1416, 2.4300], dtype=torch.float64, requires_grad=True)\n",
"Final cost after optimization:\n",
"tensor(-2.0000, dtype=torch.float64, grad_fn=<SubBackward0>)\n"
]
}
],
"source": [
"p1_final, p2_final = opt.param_groups[0][\"params\"]\n",
"print(\"Optimized angles:\")\n",
"print(p1_final)\n",
"print(p2_final)\n",
"print(\"Final cost after optimization:\")\n",
"print(cost(p1_final, p2_final))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": "<Figure size 600x400 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAFzCAYAAADsTAnbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHV0lEQVR4nO3deVhU9f4H8PcZlgEEBhAYQAFRUHBHLC60qFe6qJV661fWtVwyK9NyaZMWS8urt9Qss2vlVbNSy0ozM0vJpdRcUNwFcYOQRVBmWGSb+f7+QCYnAWdghjMD79fznEfmrB9O1rw757tIQggBIiIiIpko5C6AiIiIWjeGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrBhGiIiISFaOchdg6/R6PS5evAgPDw9IkiR3OURERHZDCIHi4mIEBQVBoaj/+QfDyE1cvHgRwcHBcpdBRERkt7KystC+fft6tzOM3ISHhweAmhvp6ekpczVERET2Q6vVIjg42PBdWh+GkZuofTXj6enJMEJERNQIN2vmwAasREREJCuGESIiIpIVwwgRERHJim1GiIjsmBAC1dXV0Ol0cpdCrZCDgwMcHR2bPPQFwwgRkZ2qrKxETk4OysrK5C6FWjE3NzcEBgbC2dm50edgGCEiskN6vR7nzp2Dg4MDgoKC4OzszIEZqVkJIVBZWYlLly7h3LlziIiIaHBgs4YwjBAR2aHKykro9XoEBwfDzc1N7nKolXJ1dYWTkxMuXLiAyspKuLi4NOo8bMBKRGTHGvt/okSWYom/g/xbTERERLJiGGlmVTo9Ui5cxpf7M+UuhYiIyCYwjDQzzdUq3P/fPXjpm6PQXK2SuxwiIpvUoUMHLFy40OT933jjDfTu3bvJ15UkCevXr2/yea63fft2SJKEoqIik48ZM2YMhg8fbtE6bBkbsDYzX3clOrR1w/nCMhzKvIL+XfzlLomIiKwoPj4eOTk5UKlUJh/z3nvvQQhhxapsC5+MyKBPqDcA4GBmkbyFEBGR1Tk7OyMgIMCsrtcqlQpeXl7WK8rGMIzIoE/ItTBy4YrMlRBRSyKEQFlltSyLOf8XX1xcjJEjR6JNmzYIDAzEu+++i/79+2PKlCn1HpOZmYlhw4bB3d0dnp6eePDBB5GXl3fDfh999JGhu/ODDz4IjUZj2LZ//37cdddd8PX1hUqlQr9+/XDw4EGz7nH//v3xzDPPYMqUKfD29oZarcYnn3yC0tJSjB07Fh4eHggPD8ePP/5oOOavr2lWrFgBLy8v/PTTT4iKioK7uzsGDRqEnJwcwzF/fU3TmOvWXud669evNwpFta+3li1bhpCQELi7u+Ppp5+GTqfD22+/jYCAAPj7+2P27Nlm3Sdz8TWNDGKuPRk5lHkFOr2Ag4IDFRFR012t0qHrjJ9kufaJWYlwczbtK2XatGnYtWsXNmzYALVajRkzZuDgwYP1tvnQ6/WGILJjxw5UV1dj4sSJGDFiBLZv327YLyMjA1999RW+//57aLVajBs3Dk8//TS++OILADUhaPTo0Vi0aBGEEJg/fz6GDBmC06dPw8PDw+Tf9dNPP8WLL76Iffv24csvv8SECROwbt06/POf/8TLL7+Md999F48++igyMzPrHQOmrKwM8+bNw2effQaFQoFHHnkEzz//vKFWa123LmfOnMGPP/6IzZs348yZM/i///s/nD17Fp07d8aOHTuwe/duPPbYY0hISEBsbKzJ5zUHn4zIoLPaA+5KR5RW6pCWWyx3OUREzaa4uBiffvop5s2bh4EDB6J79+5Yvnx5g3PrJCcn4+jRo1i1ahViYmIQGxuLlStXYseOHdi/f79hv/LycqxcuRK9e/fGnXfeiUWLFmHNmjXIzc0FAPz973/HI488gsjISERFReHjjz9GWVkZduzYYdbv0KtXL7z66quIiIhAUlISXFxc4Ovri/HjxyMiIgIzZsxAYWEhjhw5Uu85qqqqsGTJEvTt2xd9+vTBpEmTkJycbPXr1kWv12PZsmXo2rUr7r33XgwYMABpaWlYuHAhunTpgrFjx6JLly7Ytm2bWec1B5+MyMBBISE6xAu/ni7Awcwr6BrkKXdJRNQCuDo54MSsRNmubYqzZ8+iqqoKt956q2GdSqVCly5d6j3m5MmTCA4ORnBwsGFd165d4eXlhZMnT+KWW24BAISEhKBdu3aGfeLi4qDX65GWloaAgADk5eXh1Vdfxfbt25Gfnw+dToeysjJkZpo31ELPnj0NPzs4OKBt27bo0aOHYZ1arQYA5Ofn13sONzc3dOrUyfA5MDCwwf0tdd26dOjQwejJkFqthoODg9FgZmq12uzzmoNhRCZ9QrxrwsiFK3jkb6Fyl0NELYAkSSa/KmmNRo8ejcLCQrz33nsIDQ2FUqlEXFwcKisrzTqPk5OT0WdJkozW1bbJ0Ov1Zp3jZu1uzL2uQqG44ZxVVTcOKXGz89aua+j3aSq+ppFJbY+alEw2YiWi1qNjx45wcnIyer2i0WiQnp5e7zFRUVHIyspCVlaWYd2JEydQVFSErl27GtZlZmbi4sWLhs+///47FAqF4anLrl278Oyzz2LIkCHo1q0blEolCgoKLPnr2RQ/Pz8UFxejtLTUsC41NVW+ghrAMCKT3sFekCTgQmEZCkoq5C6HiKhZeHh4YPTo0XjhhRewbds2HD9+HOPGjYNCoai362tCQgJ69OiBkSNH4uDBg9i3bx9GjRqFfv36oW/fvob9XFxcMHr0aBw+fBi//vornn32WTz44IMICAgAAEREROCzzz7DyZMnsXfvXowcORKurq7N8nvLITY2Fm5ubnj55Zdx5swZrFq1CitWrJC7rDrZVRjZuXMn7r33XgQFBZk0Sl5td6q/LrWNmeSkcnVCZ/+ad3Ts4ktErcmCBQsQFxeHe+65BwkJCbjtttsQFRVV74yvkiThu+++g7e3N+68804kJCSgY8eO+PLLL432Cw8Px3333YchQ4bgH//4B3r27IkPP/zQsP1///sfrly5gj59+uDRRx/Fs88+C3//ljvwpI+PDz7//HNs2rQJPXr0wOrVq/HGG2/IXVadJGFHQ7z9+OOP2LVrF2JiYnDfffdh3bp1DQ6Xu337dkOrYE/PPxuJ+vv7mzzLoFarhUqlgkajMTqHJSR9exSr92XiyX4dkTQ4yqLnJqKWrby8HOfOnUNYWFijp223FaWlpWjXrh3mz5+PcePGyV0Omamhv4umfofaVUunwYMHY/DgwWYf5+/vb5Mj2fUJ8cLqfZl8MkJErcqhQ4dw6tQp3HrrrdBoNJg1axYAYNiwYTJXRnKxq9c0jdW7d28EBgbirrvuwq5duxrct6KiAlqt1mixltrBzw7/oUFltfVaKRMR2Zp58+ahV69eSEhIQGlpKX799Vf4+vrKXRbJxK6ejJgrMDDQMKhMRUUFli5div79+2Pv3r3o06dPncfMmTMHM2fObJb6wnzbwNvNCVfKqnD8ogbR14aJJyJqyaKjo5GSkiJ3GWRDWvSTkS5duuDJJ59ETEwM4uPjsWzZMsTHx+Pdd9+t95ikpCRoNBrDcn1XMkuTJMnwdIST5hERUWvVosNIXW699VZkZGTUu12pVMLT09NosaZoTppHRE1gR30QqIWyxN/BVhdGUlNTERgYKHcZBrVPRg5cuMz/qBCRyWpHyCwrK5O5Emrtav8O/nXUVnPYVZuRkpISo6ca586dQ2pqKnx8fBASEoKkpCRkZ2dj5cqVAICFCxciLCwM3bp1Q3l5OZYuXYpffvkFP//8s1y/wg16tfeCg0JCnrYCFzXlaOfVcgfgISLLcXBwgJeXl2G+EDc3t3oHDSOyBiEEysrKkJ+fDy8vLzg4mDY/UV3sKowcOHAAAwYMMHyeNm0agJr5BlasWIGcnByjCY8qKyvx3HPPITs7G25ubujZsye2bt1qdA65uTo7oFuQJ478ocHBC1cYRojIZLUji1pzAjOim/Hy8jL8XWwsuxr0TA7WHPSs1hsbjmPF7vMYE98BbwztZpVrEFHLpdPp6pwAjcjanJycGnwi0iIHPWup+oR6Y8Xu8zjISfOIqBEcHBya9IicSG6trgGrLaptxHriohZXK3UyV0NERNS8GEZsQJDKBQGeLqjWCxz5o0jucoiIiJoVw4gNuH7wsxS+qiEiolaGYcRGRId4AeDgZ0RE1PowjNgIw5ORC1c4+BkREbUqDCM2oluQCs6OClwpq8K5glK5yyEiImo2DCM2wtlRgV7tVQA4aR4REbUuDCM2pM91r2qIiIhaC4YRG9KHM/gSEVErxDBiQ2rDSHp+MbTlHNqZiIhaB4YRG+LnoURoWzcIAaSy3QgREbUSDCM2JiaE7UaIiKh1YRixMdHXGrFy0jwiImotGEZsTO2TkUOZRdDpOfgZERG1fAwjNqZLgAfaODugpKIap/OL5S6HiIjI6hhGbIyDQkI0240QEVErwjBig/pcmzSPYYSIiFoDhhEbVDsSKwc/IyKi1oBhxAbVvqY5X1iGgpIKmashIiKyLoYRG6RydUJntTuAml41RERELRnDiI2K4aR5RETUSjCM2KhoTppHREStBMOIjap9MnL4jyJUVutlroaIiMh6GEZsVEffNvByc0JFtR4nc7Ryl0NERGQ1DCM2SpIkTppHREStAsOIDasdbySFk+YREVELZldhZOfOnbj33nsRFBQESZKwfv36mx6zfft29OnTB0qlEuHh4VixYoXV67SUPmzESkRErYBdhZHS0lL06tULixcvNmn/c+fO4e6778aAAQOQmpqKKVOm4PHHH8dPP/1k5Uoto1ewCg4KCTmaclwsuip3OURERFbhKHcB5hg8eDAGDx5s8v5LlixBWFgY5s+fDwCIiorCb7/9hnfffReJiYnWKtNi3Jwd0TXQE0ezNTiYeQVBXq5yl0RERGRxdvVkxFx79uxBQkKC0brExETs2bOn3mMqKiqg1WqNFjlx8DMiImrpWnQYyc3NhVqtNlqnVquh1Wpx9Wrdrz3mzJkDlUplWIKDg5uj1HpFX5vBl+1GiIiopWrRYaQxkpKSoNFoDEtWVpas9dQ+GTl+UYvyKp2stRAREVmDXbUZMVdAQADy8vKM1uXl5cHT0xOurnW3v1AqlVAqlc1RnknaeblC7alEnrYCR/7Q4NYwH7lLIiIisqgW/WQkLi4OycnJRuu2bNmCuLg4mSoynyRJbDdCREQtml2FkZKSEqSmpiI1NRVATdfd1NRUZGZmAqh5xTJq1CjD/k899RTOnj2LF198EadOncKHH36Ir776ClOnTpWj/Ebrw5FYiYioBbOrMHLgwAFER0cjOjoaADBt2jRER0djxowZAICcnBxDMAGAsLAw/PDDD9iyZQt69eqF+fPnY+nSpXbRrfd6tSOxHsy8AiGEzNUQERFZliT47dYgrVYLlUoFjUYDT09PWWqoqNahxxs/o7Jaj+3P90cH3zay1EFERGQOU79D7erJSGuldHRAz3YqAHxVQ0RELQ/DiJ3gpHlERNRSMYzYCU6aR0RELRXDiJ3oE+oFAEjLK4a2vEreYoiIiCyIYcRO+Hu4IMTHDUIAh7OK5C6HiIjIYhhG7AgHPyMiopaIYcSO9Lk2aR7DCBERtSQMI3aktkdNamYR9HoOD0NERC0Dw4gd6aL2QBtnBxRXVON0fonc5RAREVkEw4gdcXRQoDdf1RARUQvDMGJnOGkeERG1NAwjdub6SfOIiIhaAoYRO9MnuCaMnCsoxeXSSpmrISIiajqGETujcnNChL87AA4NT0RELQPDiB2K4aR5RETUgjCM2CE2YiUiopaEYcQO1TZiPfJHEap0epmrISIiahqGETvU0bcNVK5OKK/SIz2vWO5yiIiImoRhxA4pFBK6qD0AABkciZWIiOwcw4idClfX9Kg5nccwQkRE9o1hxE7Vdu/laxoiIrJ3DCN2qjNf0xARUQvBMGKnap+MnC8sRUW1TuZqiIiIGo9hxE75eSjh6eIIvQDOXiqVuxwiIqJGYxixU5IkGV7VnOarGiIismMMI3Ys4lqPmgw2YiUiIjvGMGLHwv35ZISIiOyf3YWRxYsXo0OHDnBxcUFsbCz27dtX774rVqyAJElGi4uLSzNWa13s3ktERC2BXYWRL7/8EtOmTcPrr7+OgwcPolevXkhMTER+fn69x3h6eiInJ8ewXLhwoRkrtq7aNiPnC8tQWc05aoiIyD7ZVRhZsGABxo8fj7Fjx6Jr165YsmQJ3NzcsGzZsnqPkSQJAQEBhkWtVjdjxdal9lTCQ+kInV7gfCF71BARkX2ymzBSWVmJlJQUJCQkGNYpFAokJCRgz5499R5XUlKC0NBQBAcHY9iwYTh+/HiD16moqIBWqzVabJUkSYZh4fmqhoiI7JXdhJGCggLodLobnmyo1Wrk5ubWeUyXLl2wbNkyfPfdd/j888+h1+sRHx+PP/74o97rzJkzByqVyrAEBwdb9PewtM61jVg5Rw0REdkpuwkjjREXF4dRo0ahd+/e6NevH7799lv4+fnho48+qveYpKQkaDQaw5KVldWMFZvP0L2XPWqIiMhOOcpdgKl8fX3h4OCAvLw8o/V5eXkICAgw6RxOTk6Ijo5GRkZGvfsolUoolcom1dqcwq/1qDmdz9c0RERkn+zmyYizszNiYmKQnJxsWKfX65GcnIy4uDiTzqHT6XD06FEEBgZaq8xmF3GtR825glJU6dijhoiI7I/dhBEAmDZtGj755BN8+umnOHnyJCZMmIDS0lKMHTsWADBq1CgkJSUZ9p81axZ+/vlnnD17FgcPHsQjjzyCCxcu4PHHH5frV7C4IJUL2jg7oEoncIE9aoiIyA6Z/Zpm27ZtGDBggDVquakRI0bg0qVLmDFjBnJzc9G7d29s3rzZ0Kg1MzMTCsWf+erKlSsYP348cnNz4e3tjZiYGOzevRtdu3aVpX5rqOlR44HDWUU4nVdiGJWViIjIXkhCCGHOAUqlEu3bt8fYsWMxevRom+9t0lRarRYqlQoajQaenp5yl1On59cextcpf2BqQmdMToiQuxwiIiIApn+Hmv2aJjs7G5MmTcLXX3+Njh07IjExEV999RUqKyubVDA1XgQbsRIRkR0zO4z4+vpi6tSpSE1Nxd69e9G5c2c8/fTTCAoKwrPPPovDhw9bo05qQO2w8OzeS0RE9qhJDVj79OmDpKQkTJo0CSUlJVi2bBliYmJwxx133HSkU7Kc2u69Zy+Vopo9aoiIyM40KoxUVVXh66+/xpAhQxAaGoqffvoJH3zwAfLy8pCRkYHQ0FA88MADlq6V6tHOyxWuTg6o1Olx4XKZ3OUQERGZxezeNM888wxWr14NIQQeffRRvP322+jevbthe5s2bTBv3jwEBQVZtFCqn0IhIULtjiN/aHA6rwSd/NzlLomIiMhkZoeREydOYNGiRbjvvvvqHanU19cX27Zta3JxZLpw/5owkpFfDMC0EWmJiIhsgdlh5PoRUOs9qaMj+vXr16iCqHEiaifMYyNWIiKyM2aHkQ0bNtS5XpIkuLi4IDw8HGFhYU0ujMxT2703nbP3EhGRnTE7jAwfPhySJOGvY6XVrpMkCbfffjvWr18Pb29vixVKDavt3nvmUgl0egEHhSRzRURERKYxuzfNli1bcMstt2DLli3QaDTQaDTYsmULYmNjsXHjRuzcuROFhYV4/vnnrVEv1aOdtytcnBSorNYjiz1qiIjIjpj9ZGTy5Mn4+OOPER8fb1g3cOBAuLi44IknnsDx48excOFCPPbYYxYtlBrmoJDQyc8dxy9qkZ5XjA6+beQuiYiIyCRmPxk5c+ZMnePLe3p64uzZswCAiIgIFBQUNL06Mkvtqxo2YiUiIntidhiJiYnBCy+8gEuXLhnWXbp0CS+++CJuueUWAMDp06db/AR6tqh2JFYOC09ERPbE7Nc0S5cuxfDhw9G+fXtD4MjKykLHjh3x3XffAQBKSkrw6quvWrZSuilOmEdERPbI7DASGRmJEydO4Oeff0Z6ejoAoEuXLrjrrrugUNQ8aBk+fLhFiyTTRFw3YZ5eL6BgjxoiIrIDZoWRqqoquLq6IjU1FYMGDcKgQYOsVRc1QoiPG5wdFSiv0uOPK1cR0tZN7pKIiIhuyqw2I05OTggJCYFOp7NWPdQEtT1qAL6qISIi+2F2A9ZXXnkFL7/8Mi5fvmyNeqiJOBIrERHZG7PbjHzwwQfIyMhAUFAQQkND0aaN8XgWBw8etFhxZL7Oaj4ZISIi+9Ko4eDJdoX7/9mIlYiIyB6YHUZef/11a9RBFhKh/nOsEfaoISIie2B2mxEAKCoqwtKlS5GUlGRoO3Lw4EFkZ2dbtDgyX6iPG5wdFCir1CG76Krc5RAREd2U2U9Gjhw5goSEBKhUKpw/fx7jx4+Hj48Pvv32W2RmZmLlypXWqJNM5OigQEe/NjiVW4yM/BIE+7B7LxER2Tazn4xMmzYNY8aMwenTp+Hi4mJYP2TIEOzcudOixVHjhHMkViIisiNmh5H9+/fjySefvGF9u3btkJuba5GiqGkirjViZfdeIiKyB2aHEaVSCa1We8P69PR0+Pn5WaQoapoIQ/dehhEiIrJ9ZoeRoUOHYtasWaiqqgIASJKEzMxMvPTSS7j//vstXiCZr3askYy8YgghZK6GiIioYWaHkfnz56OkpAT+/v64evUq+vXrh/DwcHh4eGD27NnWqNHI4sWL0aFDB7i4uCA2Nhb79u1rcP+1a9ciMjISLi4u6NGjBzZt2mT1GuUW2rYNHBUSSit1yNGUy10OERFRg8wOIyqVClu2bMH333+P999/H5MmTcKmTZuwY8eOG0ZjtbQvv/wS06ZNw+uvv46DBw+iV69eSExMRH5+fp377969Gw8//DDGjRuHQ4cOYfjw4Rg+fDiOHTtm1Trl5uSgQJhvzT+L9Dw2YiUiItsmCTt6jh8bG4tbbrkFH3zwAQBAr9cjODgYzzzzDKZPn37D/iNGjEBpaSk2btxoWPe3v/0NvXv3xpIlS0y6plarhUqlgkajgaenp2V+kWYw8YuD+OFoDl69OwqP39FR7nKIiKgVMvU71OxxRgAgOTkZycnJyM/Ph16vN9q2bNmyxpzypiorK5GSkoKkpCTDOoVCgYSEBOzZs6fOY/bs2YNp06YZrUtMTMT69eutUqMtMXTvZY8aIiKycWaHkZkzZ2LWrFno27cvAgMDIUnNM9x4QUEBdDod1Gq10Xq1Wo1Tp07VeUxubm6d+zfUBbmiogIVFRWGz3X1HLIHtT1q0jnWCBER2Tizw8iSJUuwYsUKPProo9aoR3Zz5szBzJkz5S6jyWrHGsnIK4EQotlCIxERkbnMbsBaWVmJ+Ph4a9TSIF9fXzg4OCAvL89ofV5eHgICAuo8JiAgwKz9ASApKQkajcawZGVlNb14GYT5toGDQkJxRTXytBU3P4CIiEgmZoeRxx9/HKtWrbJGLQ1ydnZGTEwMkpOTDev0ej2Sk5MRFxdX5zFxcXFG+wPAli1b6t0fqBnUzdPT02ixR86OCnRoWzMvDYeFJyIiW2b2a5ry8nJ8/PHH2Lp1K3r27AknJyej7QsWLLBYcX81bdo0jB49Gn379sWtt96KhQsXorS0FGPHjgUAjBo1Cu3atcOcOXMAAJMnT0a/fv0wf/583H333VizZg0OHDiAjz/+2Go12pIIfw+cuVSK9LwS3BHB0XGJiMg2NWrW3t69ewPADeN1WLtdwogRI3Dp0iXMmDEDubm56N27NzZv3mxopJqZmQmF4s+HPfHx8Vi1ahVeffVVvPzyy4iIiMD69evRvXt3q9ZpKzqr3bH5OJDBJyNERGTD7GqcETnY6zgjALDh8EU8u/oQ+oZ64+sJzd/Oh4iIWjdTv0PNbjPSkPpGQiV5RFwbaySdc9QQEZENMzmMuLm54dKlS4bPd999N3Jycgyf8/LyEBgYaNnqqEnCfNtAIQHa8mpcKmaPGiIisk0mh5Hy8nKj/7veuXMnrl69arQP/+/btrg4OaBD25o5ak7ncyRWIiKyTRZ9TcOBtWzPn8PCsxErERHZJouGEbI9fw4LzycjRERkm0wOI5IkGT35+Otnsk2d1X8OC09ERGSLTB5nRAiBzp07GwJISUkJoqOjDeN6sL2Ibap9TZOeX8w5aoiIyCaZHEaWL19uzTrISjr5uUOSgKKyKhSWVsLXXSl3SUREREZMDiOjR4+2Zh1kJS5ODgjxccOFwjKk5xUzjBARkc1hA9ZWIML/WrsRNmIlIiIbxDDSCtT2qDnNRqxERGSDGEZageuHhSciIrI1DCOtgKF7L1/TEBGRDTI5jNxxxx2YN28e0tPTrVkPWUFtj5rC0koUlnCOGiIisi0mh5Hx48djz549iImJQVRUFF566SXs2rWL44vYAVdnB7T3dgXApyNERGR7TA4jo0aNwjfffIOCggLMnz8fRUVFeOCBBxAQEIDHHnsM69evv2HiPLIdtT1qOCw8ERHZGrPbjCiVSgwZMgQfffQRLl68iA0bNiAwMBCvvfYa2rZti3vuuQe7du2yRq3UBLU9ajLYiJWIiGxMkxuwxsbGYvbs2Th69CiOHj2KgQMHIicnxxK1kQXVPhk5zScjRERkY0wegdUUnTp1wtSpUy15SrKQP7v3MowQEZFtYdfeVqJ2wryCkgpcKa2UuRoiIqI/MYy0Em2Ujmjnda1HzSU+HSEiItvBMNKKcFh4IiKyRWaHkVmzZqGsrOyG9VevXsWsWbMsUhRZB4eFJyIiW2R2GJk5cyZKSm78P+uysjLMnDnTIkWRdURwWHgiIrJBZocRIQQkSbph/eHDh+Hj42ORosg6ap+MnM7nkxEiIrIdJnft9fb2hiRJkCQJnTt3NgokOp0OJSUleOqpp6xSJFlGbY+aPG0FNFeroHJ1krkiIiIiM8LIwoULIYTAY489hpkzZ0KlUhm2OTs7o0OHDoiLi7NKkWQZHi5OCFS5IEdTjoz8YsSE8kkWERHJz+QwMnr0aABAWFgYbrvtNjg6WnS8tJu6fPkynnnmGXz//fdQKBS4//778d5778Hd3b3eY/r3748dO3YYrXvyySexZMkSa5drsyLUHsjRlON0XgnDCBER2QSz24x4eHjg5MmThs/fffcdhg8fjpdffhmVldYbTGvkyJE4fvw4tmzZgo0bN2Lnzp144oknbnrc+PHjkZOTY1jefvttq9VoD/5sN8JGrEREZBvMDiNPPvkk0tPTAQBnz57FiBEj4ObmhrVr1+LFF1+0eIEAcPLkSWzevBlLly5FbGwsbr/9dixatAhr1qzBxYsXGzzWzc0NAQEBhsXT09MqNdoLdu8lIiJbY3YYSU9PR+/evQEAa9euRb9+/bBq1SqsWLEC33zzjaXrAwDs2bMHXl5e6Nu3r2FdQkICFAoF9u7d2+CxX3zxBXx9fdG9e3ckJSXVOUZKa8LuvUREZGvMbvghhIBerwcAbN26Fffccw8AIDg4GAUFBZat7prc3Fz4+/sbrXN0dISPjw9yc3PrPe5f//oXQkNDERQUhCNHjuCll15CWloavv3223qPqaioQEVFheGzVqtt+i9gQ2p71ORoylFcXgUPF/aoISIieZn9ZKRv375466238Nlnn2HHjh24++67AQDnzp2DWq0261zTp083dBeubzl16pS5JRo88cQTSExMRI8ePTBy5EisXLkS69atw5kzZ+o9Zs6cOVCpVIYlODi40de3RSpXJ6g9lQDYboSIiGyD2U9GFi5ciJEjR2L9+vV45ZVXEB4eDgD4+uuvER8fb9a5nnvuOYwZM6bBfTp27IiAgADk5+cbra+ursbly5cREBBg8vViY2MBABkZGejUqVOd+yQlJWHatGmGz1qttsUFkgh/D+RpK5CRV4I+Id5yl0NERK2c2WGkZ8+eOHr06A3r33nnHTg4OJh1Lj8/P/j5+d10v7i4OBQVFSElJQUxMTEAgF9++QV6vd4QMEyRmpoKAAgMDKx3H6VSCaVSafI57VGE2h2/ZRRwJFYiIrIJjR4sJCUlxdDFt2vXrujTp4/FivqrqKgoDBo0COPHj8eSJUtQVVWFSZMm4aGHHkJQUBAAIDs7GwMHDsTKlStx66234syZM1i1ahWGDBmCtm3b4siRI5g6dSruvPNO9OzZ02q12oMI/5pGrHxNQ0REtsDsMJKfn48RI0Zgx44d8PLyAgAUFRVhwIABWLNmjUlPOhrjiy++wKRJkzBw4EDDoGfvv/++YXtVVRXS0tIMvWWcnZ2xdetWLFy4EKWlpQgODsb999+PV1991Sr12ZMI9bWxRvIYRoiISH6SEEKYc8CIESNw9uxZrFy5ElFRUQCAEydOYPTo0QgPD8fq1autUqhctFotVCoVNBpNixmjpKisEr1nbQEAHJuZCHdl846mS0RErYOp36Fm96bZvHkzPvzwQ0MQAWpe0yxevBg//vhj46qlZuXl5gw/j5p2MWf4qoaIiGRmdhjR6/VwcrpxbAonJyfD+CNk+zgSKxER2Qqzw8jf//53TJ482WgY9uzsbEydOhUDBw60aHFkPbVhhCOxEhGR3MwOIx988AG0Wi06dOiATp06oVOnTggLC4NWq8WiRYusUSNZQe2w8OxRQ0REcjO75WJwcDAOHjyIrVu3GkZHjYqKQkJCgsWLI+v5c/ZevqYhIiJ5NaobhSRJuOuuu3DXXXdZuh5qJrVPRrIuX0VZZTXcnNmjhoiI5GHya5pffvkFXbt2rXPiOI1Gg27duuHXX3+1aHFkPT5tnOHr7gwAOJNfKnM1RETUmpkcRhYuXIjx48fX2U9YpVLhySefxIIFCyxaHFlXOF/VEBGRDTA5jBw+fBiDBg2qd/s//vEPpKSkWKQoah61w8KncyRWIiKSkclhJC8vr87xRWo5Ojri0qVLFimKmkftsPAZfDJCREQyMjmMtGvXDseOHat3+5EjRxqcDZdsDyfMIyIiW2ByGBkyZAhee+01lJeX37Dt6tWreP3113HPPfdYtDiyrtonI5mXy1BepZO5GiIiaq1M7s/56quv4ttvv0Xnzp0xadIkdOnSBQBw6tQpLF68GDqdDq+88orVCiXLa9vGGd5uTrhSVoWM/BJ0b6eSuyQiImqFTA4jarUau3fvxoQJE5CUlITayX4lSUJiYiIWL14MtVpttULJ8iRJQoS/B/adv8wwQkREsjFrpKvQ0FBs2rQJV65cQUZGBoQQiIiIgLe3t7XqIyuLULtj3/nL7N5LRESyadSwm97e3rjlllssXQvJwDAsPLv3EhGRTMyeKI9aFk6YR0REcmMYaeVqe9RcKCxljxoiIpIFw0gr5+euhMrVCXoBnCvgHDVERNT8GEZauZoeNTVPR9Lz2IiViIiaH8MIXTcsPNuNEBFR82MYoT+HhWePGiIikgHDCBmejHCsESIikgPDCBmejJwvLENFNXvUEBFR82IYIag9lfBwcYROL3C+oEzucoiIqJVhGCGjHjV8VUNERM2NYYQA/PmqJp2NWImIqJnZTRiZPXs24uPj4ebmBi8vL5OOEUJgxowZCAwMhKurKxISEnD69GnrFmqn/uzeyycjRETUvOwmjFRWVuKBBx7AhAkTTD7m7bffxvvvv48lS5Zg7969aNOmDRITE1FeXm7FSu2TYY4aPhkhIqJm1qhZe+Uwc+ZMAMCKFStM2l8IgYULF+LVV1/FsGHDAAArV66EWq3G+vXr8dBDD1mrVLtU22bkXEEpqnR6ODnYTU4lIiI712K/cc6dO4fc3FwkJCQY1qlUKsTGxmLPnj0yVmabAlUucFc6olovcJ5z1BARUTNqsWEkNzcXAKBWq43Wq9Vqw7a6VFRUQKvVGi2tgSRJCDf0qOGrGiIiaj6yhpHp06dDkqQGl1OnTjVrTXPmzIFKpTIswcHBzXp9OdW+qjmarZG5EiIiak1kDSPPPfccTp482eDSsWPHRp07ICAAAJCXl2e0Pi8vz7CtLklJSdBoNIYlKyurUde3R7dH+AIANh3NgRBC5mqIiKi1kLUBq5+fH/z8/Kxy7rCwMAQEBCA5ORm9e/cGAGi1Wuzdu7fBHjlKpRJKpdIqNdm6u7qq4erkgAuFZTjyhwa9gr3kLomIiFoBu2kzkpmZidTUVGRmZkKn0yE1NRWpqakoKfmzfUNkZCTWrVsHoKYNxJQpU/DWW29hw4YNOHr0KEaNGoWgoCAMHz5cpt/Ctrk5OyKha00bm+9SL8pcDRERtRZ207V3xowZ+PTTTw2fo6OjAQDbtm1D//79AQBpaWnQaP5s7/Diiy+itLQUTzzxBIqKinD77bdj8+bNcHFxadba7cmwXkH4/vBFbDxyEa/cHQUHhSR3SURE1MJJgo0DGqTVaqFSqaDRaODp6Sl3OVZXWa3HLbO3QnO1CqvGxyK+k6/cJRERkZ0y9TvUbl7TUPNwdlRgSI+aBr4b+KqGiIiaAcMI3eDeXkEAgB+P5aKiWidzNURE1NIxjNANYsPawt9DCc3VKuxML5C7HCIiauEYRugGDgrJ8HRkw2G+qiEiIutiGKE6Db0WRraeyENpRbXM1RARUUvGMEJ16tlehQ5t3XC1SoetJ/NufgAREVEjMYxQnSRJMjwdYa8aIiKyJoYRqtfQ3jVhZEf6JVwprZS5GiIiaqkYRqhe4f4e6BroiWq9wI/HcuUuh4iIWiiGEWpQ7dORDYezZa6EiIhaKoYRalBtF9+95y4jV1MuczVERNQSMYxQg9p5ueKWDt4QAth4hA1ZiYjI8hhG6KZqe9V8x141RERkBQwjdFNDegTCQSHhaLYG5wpK5S6HiIhaGIYRuqm27krcHu4LgGOOEBGR5TGMkEmGXetV893hbAghZK6GiIhaEoYRMsk/ugVA6ajA2UulOH5RK3c5RETUgjCMkEnclY4YGOUPgDP5EhGRZTGMkMmG9moHAPj+8EXo9XxVQ0RElsEwQibr38UPHkpH5GjKceDCFbnLISKiFoJhhEzm4uSAQd0DAADfpXJ4eCIisgyGETJL7Vw1m47moEqnl7kaIiJqCRhGyCxxHdvC190ZV8qq8NvpArnLISKiFoBhhMzi6KDAPT1rZ/JlrxoiImo6hhEyW+1Mvj8fz8XVSp3M1RARkb1jGCGz9QnxQntvV5RW6pB8Kk/ucoiIyM4xjJDZJEkyzOTLuWqIiKip7CaMzJ49G/Hx8XBzc4OXl5dJx4wZMwaSJBktgwYNsm6hrcSw3jUDoG1PuwTN1SqZqyEiIntmN2GksrISDzzwACZMmGDWcYMGDUJOTo5hWb16tZUqbF26BHigi9oDlTo9fjqWK3c5RERkx+wmjMycORNTp05Fjx49zDpOqVQiICDAsHh7e1upwtandswR9qohIqKmsJsw0ljbt2+Hv78/unTpggkTJqCwsFDuklqM2nYju88UIL+4XOZqiIjIXrXoMDJo0CCsXLkSycnJ+M9//oMdO3Zg8ODB0Onq745aUVEBrVZrtFDdgn3cEB3iBb0AfjiSI3c5RERkp2QNI9OnT7+hgelfl1OnTjX6/A899BCGDh2KHj16YPjw4di4cSP279+P7du313vMnDlzoFKpDEtwcHCjr98aDLv2dOQ79qohIqJGkjWMPPfcczh58mSDS8eOHS12vY4dO8LX1xcZGRn17pOUlASNRmNYsrKyLHb9lujunkFQSEBqVhEyC8vkLoeIiOyQo5wX9/Pzg5+fX7Nd748//kBhYSECAwPr3UepVEKpVDZbTfbOz0OJ+E6++C2jAN8fuYiJA8LlLomIiOyM3bQZyczMRGpqKjIzM6HT6ZCamorU1FSUlJQY9omMjMS6desAACUlJXjhhRfw+++/4/z580hOTsawYcMQHh6OxMREuX6NFqm2V813qdkyV0JERPbIbsLIjBkzEB0djddffx0lJSWIjo5GdHQ0Dhw4YNgnLS0NGo0GAODg4IAjR45g6NCh6Ny5M8aNG4eYmBj8+uuvfPJhYYndAuDsoEB6XglO5bLBLxERmUcSQgi5i7BlWq0WKpUKGo0Gnp6ecpdjs5787AB+Op6HCf074aVBkXKXQ0RENsDU71C7eTJCtm1or5rh4b8/fBHMt0REZA6GEbKIgVH+aOPsgD+uXMXBzCK5yyEiIjvCMEIW4eLkgMRuAQCADWzISkREZmAYIYu591qvmh+O5qBap5e5GiIishcMI2Qxt4f7wqeNMwpKKrH7DOcAIiIi0zCMkMU4OSgwpMe1VzWcyZeIiEzEMEIWVdur5qdjuSivqn9CQiIioloMI2RRfUO9EaRyQXFFNban5ctdDhER2QGGEbIohULCvddm8uWrGiIiMgXDCFlc7Vw1W0/mo7i8SuZqiIjI1jGMkMV1DfREJ782qKzW4+fjeXKXQ0RENo5hhCxOkiRDQ1a+qiEiopthGCGrqH1V81tGAQpLKmSuhoiIbBnDCFlFmG8b9Gyvgk4vsOlojtzlEBGRDWMYIasZeq1XzXepfFVDRET1Yxghq7m3VxAkCThw4Qqyi67KXQ4REdkohhGyGrWnC2LDfAAA37MhKxER1YNhhKxqWO+aXjWLf8nAh9szOEQ8ERHdgGGErGporyD0DvZCcUU13t6chv7vbMeX+zNRrdPLXRoREdkISQgh5C7Clmm1WqhUKmg0Gnh6espdjl3S6wXWp2Zj/s/phrYjEf7ueHFQJBKi/CFJkswVEhGRNZj6HcowchMMI5ZTXqXD579fwAfbMlBUVjNMfN9Qb0wfHIm+HXxkro6IiCyNYcRCGEYsT1tehSXbz2DZrnMor6p5XXNXVzVeTOyCCLWHzNUREZGlMIxYCMOI9eRqyvFecjq+3J8FvQAUEvBATDCm3tUZASoXucsjIqImYhixEIYR68vIL8Y7P6Xhp2uT6ikdFXjs9jA81a8TVK5OMldHRESNxTBiIQwjzSflwhXM/fEk9p+/AgBQuTph0oBwPBoXChcnB5mrIyIiczGMWAjDSPMSQiD5ZD7+s/kUTueXAACCVC6Y9o8u+Gd0Ozgo2POGiMheMIxYCMOIPHR6gW8O/oF3t6QjR1MOAIgM8MBLgyLRv4sfuwMTEdkBhhELYRiRV3mVDit2n8eH2zKgLa8GAMSG+WD64EhEh3jLXB0RETWkRYWR8+fP480338Qvv/yC3NxcBAUF4ZFHHsErr7wCZ2fneo8rLy/Hc889hzVr1qCiogKJiYn48MMPoVarTb42w4htKCqrxH+3n8Hy3edRWV3THfi28LYI8WkDP3dntHVXoq27M9q2UcLPo+ZPlasTFHytQ0QkG1O/Qx2bsaZGO3XqFPR6PT766COEh4fj2LFjGD9+PEpLSzFv3rx6j5s6dSp++OEHrF27FiqVCpMmTcJ9992HXbt2NWP1ZAlebs5IGhKFUfEd8O6WdHxz8A/syijELhTWe4yDQoJPG2f4uivh6+6Mttd+rg0uvu7XfW7jzEayREQysYsnI3V555138N///hdnz56tc7tGo4Gfnx9WrVqF//u//wNQE2qioqKwZ88e/O1vfzPpOnwyYpsy8oux5+xlFJZUoKCkAoUllSgsqURBaQUKiisMr3TM4a50hE8bZzg7KuCokODkoICDQoKTg3Ttz5rPjoqa7Y4ON+7z57Y/93FUSFAoJEiQIEk146nU/gwAkiRdW1fzs3TDz8bH4frtqD3Hn7+H0c/X9qiviU1t2xvJaF3990hC/RsbPs462HSIyLI8XJxwW7ivxc7Xop6M1EWj0cDHp/4hxFNSUlBVVYWEhATDusjISISEhDQYRioqKlBRUWH4rNVqLVc0WUy4vwfC/esfrbWyWo/LpZUouD6slNb8eekvnwtKKlClEyipqEZJhfkhhoiopega6IlNk+9o9uvaZRjJyMjAokWLGnxFk5ubC2dnZ3h5eRmtV6vVyM3Nrfe4OXPmYObMmZYqlWTi7KhAgMrFpJFchRDQllejsKQCV8oqUVktUK3Xo1onUK0XqNbpa/68ft219Tq9QJVOQKfXX/tToOrafjXban7WCwEBQAhAQNT8abQO0F/7oXa7XgjDtprnlwL6vxxXs9b4d7nx97vu5+v2Nhxfz/aGznPDtvo33WRj4zVUa2tgn8+0ydZ18G0jy3VlDSPTp0/Hf/7znwb3OXnyJCIjIw2fs7OzMWjQIDzwwAMYP368xWtKSkrCtGnTDJ+1Wi2Cg4Mtfh2yHZIkQeXqxNFeiYhkImsYee655zBmzJgG9+nYsaPh54sXL2LAgAGIj4/Hxx9/3OBxAQEBqKysRFFRkdHTkby8PAQEBNR7nFKphFKpNKl+IiIiajpZw4ifnx/8/PxM2jc7OxsDBgxATEwMli9fDoVC0eD+MTExcHJyQnJyMu6//34AQFpaGjIzMxEXF9fk2omIiMgyGv5GtxHZ2dno378/QkJCMG/ePFy6dAm5ublGbT+ys7MRGRmJffv2AQBUKhXGjRuHadOmYdu2bUhJScHYsWMRFxdnck8aIiIisj67aMC6ZcsWZGRkICMjA+3btzfaVttgr6qqCmlpaSgrKzNse/fdd6FQKHD//fcbDXpGREREtsNuxxlpLhxnhIiIqHFM/Q61i9c0RERE1HIxjBAREZGsGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJyi7GGZFTbc9nzt5LRERkntrvzpuNIsIwchPFxcUAwMnyiIiIGqm4uBgqlare7Rz07Cb0ej0uXrwIDw8PSJJkkXPWzgSclZXFgdQshPfUOnhfLY/31Dp4Xy3PEvdUCIHi4mIEBQU1OKccn4zchEKhuGEIekvx9PTkvzQWxntqHbyvlsd7ah28r5bX1Hva0BORWmzASkRERLJiGCEiIiJZMYzIQKlU4vXXX4dSqZS7lBaD99Q6eF8tj/fUOnhfLa857ykbsBIREZGs+GSEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrBhGZLB48WJ06NABLi4uiI2Nxb59++QuyW7MmTMHt9xyCzw8PODv74/hw4cjLS3NaJ/y8nJMnDgRbdu2hbu7O+6//37k5eXJVLH9mTt3LiRJwpQpUwzreE/Nl52djUceeQRt27aFq6srevTogQMHDhi2CyEwY8YMBAYGwtXVFQkJCTh9+rSMFds+nU6H1157DWFhYXB1dUWnTp3w5ptvGs17wvvasJ07d+Lee+9FUFAQJEnC+vXrjbabcv8uX76MkSNHwtPTE15eXhg3bhxKSkqaVpigZrVmzRrh7Owsli1bJo4fPy7Gjx8vvLy8RF5entyl2YXExESxfPlycezYMZGamiqGDBkiQkJCRElJiWGfp556SgQHB4vk5GRx4MAB8be//U3Ex8fLWLX92Ldvn+jQoYPo2bOnmDx5smE976l5Ll++LEJDQ8WYMWPE3r17xdmzZ8VPP/0kMjIyDPvMnTtXqFQqsX79enH48GExdOhQERYWJq5evSpj5bZt9uzZom3btmLjxo3i3LlzYu3atcLd3V289957hn14Xxu2adMm8corr4hvv/1WABDr1q0z2m7K/Rs0aJDo1auX+P3338Wvv/4qwsPDxcMPP9ykuhhGmtmtt94qJk6caPis0+lEUFCQmDNnjoxV2a/8/HwBQOzYsUMIIURRUZFwcnISa9euNexz8uRJAUDs2bNHrjLtQnFxsYiIiBBbtmwR/fr1M4QR3lPzvfTSS+L222+vd7terxcBAQHinXfeMawrKioSSqVSrF69ujlKtEt33323eOyxx4zW3XfffWLkyJFCCN5Xc/01jJhy/06cOCEAiP379xv2+fHHH4UkSSI7O7vRtfA1TTOqrKxESkoKEhISDOsUCgUSEhKwZ88eGSuzXxqNBgDg4+MDAEhJSUFVVZXRPY6MjERISAjv8U1MnDgRd999t9G9A3hPG2PDhg3o27cvHnjgAfj7+yM6OhqffPKJYfu5c+eQm5trdE9VKhViY2N5TxsQHx+P5ORkpKenAwAOHz6M3377DYMHDwbA+9pUpty/PXv2wMvLC3379jXsk5CQAIVCgb179zb62pworxkVFBRAp9NBrVYbrVer1Th16pRMVdkvvV6PKVOm4LbbbkP37t0BALm5uXB2doaXl5fRvmq1Grm5uTJUaR/WrFmDgwcPYv/+/Tds4z0139mzZ/Hf//4X06ZNw8svv4z9+/fj2WefhbOzM0aPHm24b3X9t4D3tH7Tp0+HVqtFZGQkHBwcoNPpMHv2bIwcORIAeF+byJT7l5ubC39/f6Ptjo6O8PHxadI9ZhghuzVx4kQcO3YMv/32m9yl2LWsrCxMnjwZW7ZsgYuLi9zltAh6vR59+/bFv//9bwBAdHQ0jh07hiVLlmD06NEyV2e/vvrqK3zxxRdYtWoVunXrhtTUVEyZMgVBQUG8r3aOr2maka+vLxwcHG7ohZCXl4eAgACZqrJPkyZNwsaNG7Ft2za0b9/esD4gIACVlZUoKioy2p/3uH4pKSnIz89Hnz594OjoCEdHR+zYsQPvv/8+HB0doVareU/NFBgYiK5duxqti4qKQmZmJgAY7hv/W2CeF154AdOnT8dDDz2EHj164NFHH8XUqVMxZ84cALyvTWXK/QsICEB+fr7R9urqaly+fLlJ95hhpBk5OzsjJiYGycnJhnV6vR7JycmIi4uTsTL7IYTApEmTsG7dOvzyyy8ICwsz2h4TEwMnJyeje5yWlobMzEze43oMHDgQR48eRWpqqmHp27cvRo4cafiZ99Q8t9122w1dztPT0xEaGgoACAsLQ0BAgNE91Wq12Lt3L+9pA8rKyqBQGH9tOTg4QK/XA+B9bSpT7l9cXByKioqQkpJi2OeXX36BXq9HbGxs4y/e6Kav1Chr1qwRSqVSrFixQpw4cUI88cQTwsvLS+Tm5spdml2YMGGCUKlUYvv27SInJ8ewlJWVGfZ56qmnREhIiPjll1/EgQMHRFxcnIiLi5OxavtzfW8aIXhPzbVv3z7h6OgoZs+eLU6fPi2++OIL4ebmJj7//HPDPnPnzhVeXl7iu+++E0eOHBHDhg1jF9SbGD16tGjXrp2ha++3334rfH19xYsvvmjYh/e1YcXFxeLQoUPi0KFDAoBYsGCBOHTokLhw4YIQwrT7N2jQIBEdHS327t0rfvvtNxEREcGuvfZo0aJFIiQkRDg7O4tbb71V/P7773KXZDcA1LksX77csM/Vq1fF008/Lby9vYWbm5v45z//KXJycuQr2g79NYzwnprv+++/F927dxdKpVJERkaKjz/+2Gi7Xq8Xr732mlCr1UKpVIqBAweKtLQ0maq1D1qtVkyePFmEhIQIFxcX0bFjR/HKK6+IiooKwz68rw3btm1bnf8NHT16tBDCtPtXWFgoHn74YeHu7i48PT3F2LFjRXFxcZPqkoS4bug6IiIiombGNiNEREQkK4YRIiIikhXDCBEREcmKYYSIiIhkxTBCREREsmIYISIiIlkxjBAREZGsGEaI6KbeeOMN9O7du0nnOH/+PCRJQmpqqkVqqk///v0xZcoUq16DiCyLYYSoBcjKysJjjz2GoKAgODs7IzQ0FJMnT0ZhYaHZ55IkCevXrzda9/zzzxvNV9EYwcHByMnJQffu3Zt0nlrbt2+HJEk3TOD37bff4s0337TINRqjuUIXUUvCMEJk586ePYu+ffvi9OnTWL16NTIyMrBkyRLDBIyXL19u8jXc3d3Rtm3bJp3DwcEBAQEBcHR0bHI9DfHx8YGHh4dVr0FElsUwQmTnJk6cCGdnZ/z888/o168fQkJCMHjwYGzduhXZ2dl45ZVXDPt26NABb775Jh5++GG0adMG7dq1w+LFi422A8A///lPSJJk+PzX1zRjxozB8OHD8e9//xtqtRpeXl6YNWsWqqur8cILL8DHxwft27fH8uXLDcf89YnBmDFjIEnSDcv27dsBAJ999hn69u0LDw8PBAQE4F//+pdh6vLz589jwIABAABvb29IkoQxY8YAuPE1zZUrVzBq1Ch4e3vDzc0NgwcPxunTpw3bV6xYAS8vL/z000+IioqCu7s7Bg0ahJycnHrv+ZUrVzBy5Ej4+fnB1dUVERERht+1dibp6OhoSJKE/v37G45bunQpoqKi4OLigsjISHz44Yc33J81a9YgPj4eLi4u6N69O3bs2GHSdYnsWpNmtiEiWRUWFgpJksS///3vOrePHz9eeHt7C71eL4QQIjQ0VHh4eIg5c+aItLQ08f777wsHBwfx888/CyGEyM/PN0w8mJOTI/Lz84UQQrz++uuiV69ehvOOHj1aeHh4iIkTJ4pTp06J//3vfwKASExMFLNnzxbp6enizTffFE5OTiIrK0sIIcS5c+cEAHHo0CEhhBBFRUVGMy9PnjxZ+Pv7Gybg+9///ic2bdokzpw5I/bs2SPi4uLE4MGDhRBCVFdXi2+++UYAEGlpaSInJ0cUFRUJIW6c5G/o0KEiKipK7Ny5U6SmporExEQRHh4uKisrhRBCLF++XDg5OYmEhASxf/9+kZKSIqKiosS//vWveu/7xIkTRe/evcX+/fvFuXPnxJYtW8SGDRuEEDUz9gIQW7duFTk5OaKwsFAIIcTnn38uAgMDxTfffCPOnj0rvvnmG+Hj4yNWrFhhdH/at28vvv76a3HixAnx+OOPCw8PD1FQUHDT6xLZM4YRIjv2+++/CwBi3bp1dW5fsGCBACDy8vKEEDVhZNCgQUb7jBgxwvAlL4So83x1hZHQ0FCh0+kM67p06SLuuOMOw+fq6mrRpk0bsXr1aiHEjWHket98841wcXERv/32W72/6/79+wUAw+ygtbOPXrlyxWi/68NIenq6ACB27dpl2F5QUCBcXV3FV199JYSoCSMAREZGhmGfxYsXC7VaXW8t9957rxg7dmyd2+r7PTt16iRWrVpltO7NN98UcXFxRsfNnTvXsL2qqkq0b99e/Oc//7npdYnsGV/TELUAwozJt+Pi4m74fPLkSbOv2a1bNygUf/4nRK1Wo0ePHobPDg4OaNu2reHVSn0OHTqERx99FB988AFuu+02w/qUlBTce++9CAkJgYeHB/r16wcAyMzMNLnGkydPwtHREbGxsYZ1bdu2RZcuXYx+Zzc3N3Tq1MnwOTAwsMG6J0yYgDVr1qB379548cUXsXv37gbrKC0txZkzZzBu3Di4u7sblrfeegtnzpwx2vf6fz6Ojo7o27evoVZzr0tkLxhGiOxYeHg4JEmqN0ycPHkS3t7e8PPzs/i1nZycjD5LklTnOr1eX+85cnNzMXToUDz++OMYN26cYX1paSkSExPh6emJL774Avv378e6desAAJWVlRb8LWrUVXdDAW/w4MG4cOECpk6diosXL2LgwIF4/vnn692/pKQEAPDJJ58gNTXVsBw7dgy///67yXWae10ie8EwQmTH2rZti7vuugsffvghrl69arQtNzcXX3zxBUaMGAFJkgzr//rl9/vvvyMqKsrw2cnJCTqdzrqFAygvL8ewYcMQGRmJBQsWGG07deoUCgsLMXfuXNxxxx2IjIy84UmFs7MzADRYa1RUFKqrq7F3717DusLCQqSlpaFr165Nqt/Pzw+jR4/G559/joULF+Ljjz+uty61Wo2goCCcPXsW4eHhRkttg9da1//zqa6uRkpKitE/n/quS2TPrNvHjois7oMPPkB8fDwSExPx1ltvISwsDMePH8cLL7yAdu3aYfbs2Ub779q1C2+//TaGDx+OLVu2YO3atfjhhx8M2zt06IDk5GTcdtttUCqV8Pb2tkrdTz75JLKyspCcnIxLly4Z1vv4+CAkJATOzs5YtGgRnnrqKRw7duyGsUNCQ0MhSRI2btyIIUOGwNXVFe7u7kb7REREYNiwYRg/fjw++ugjeHh4YPr06WjXrh2GDRvW6NpnzJiBmJgYdOvWDRUVFdi4caMhMPj7+8PV1RWbN29G+/bt4eLiApVKhZkzZ+LZZ5+FSqXCoEGDUFFRgQMHDuDKlSuYNm2a4dyLFy9GREQEoqKi8O677+LKlSt47LHHbnpdInvGJyNEdi4iIgIHDhxAx44d8eCDD6JTp0544oknMGDAAOzZswc+Pj5G+z/33HM4cOAAoqOj8dZbb2HBggVITEw0bJ8/fz62bNmC4OBgREdHW63uHTt2ICcnB127dkVgYKBh2b17N/z8/LBixQqsXbsWXbt2xdy5czFv3jyj49u1a4eZM2di+vTpUKvVmDRpUp3XWb58OWJiYnDPPfcgLi4OQghs2rTphlcz5nB2dkZSUhJ69uyJO++8Ew4ODlizZg2AmnYe77//Pj766CMEBQUZQs/jjz+OpUuXYvny5ejRowf69euHFStW3PBkZO7cuZg7dy569eqF3377DRs2bICvr+9Nr0tkzyRhTss3IrJrHTp0wJQpUzhcug06f/48wsLCcOjQoSYPvU9kb/hkhIiIiGTFMEJERESy4msaIiIikhWfjBAREZGsGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrP4fLp2NEM3aZUkAAAAASUVORK5CYII=\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"\n",
"fig = plt.figure(figsize=(6, 4))\n",
"\n",
"# Enable processing the Torch trainable tensors\n",
"with torch.no_grad():\n",
" plt.plot(x, cost_pt, label = 'global minimum')\n",
" plt.xlabel(\"Optimization steps\")\n",
" plt.ylabel(\"Cost / Energy\")\n",
" plt.legend()\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local minimum\n",
"\n",
"If the spins are initialized close to the local minimum of zero energy, the optimizer is likely to get stuck here and never find the global minimum at -2.\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Corresponding cost before optimization:\n",
"tensor(0.0082, dtype=torch.float64, grad_fn=<SubBackward0>)\n"
]
}
],
"source": [
"torch.manual_seed(9)\n",
"p3 = Variable((np.pi*torch.rand(3, dtype = torch.float64)), requires_grad = True)\n",
"p4 = Variable((np.pi*torch.rand(3, dtype = torch.float64)), requires_grad = True)\n",
"\n",
"var_init_loc = [p3, p4]\n",
"cost_init_loc = cost(p3, p4)\n",
"\n",
"print(\"Corresponding cost before optimization:\")\n",
"print(cost_init_loc)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Energy after step 5: 0.0032761 | Angles: [array([0.77369911, 2.63471297, 1.07981163]), array([0.26038622, 0.08659858, 1.91060734])] \n",
"\n",
"Energy after step 10: 0.0013137 | Angles: [array([0.77369911, 2.63406019, 1.07981163]), array([0.26038622, 0.05483683, 1.91060734])] \n",
"\n",
"Energy after step 15: 0.0005266 | Angles: [array([0.77369911, 2.63379816, 1.07981163]), array([0.26038622, 0.03471974, 1.91060734])] \n",
"\n",
"Energy after step 20: 0.0002111 | Angles: [array([0.77369911, 2.63369307, 1.07981163]), array([0.26038622, 0.02198151, 1.91060734])] \n",
"\n",
"Energy after step 25: 0.0000846 | Angles: [array([0.77369911, 2.63365094, 1.07981163]), array([0.26038622, 0.01391648, 1.91060734])] \n",
"\n",
"Energy after step 30: 0.0000339 | Angles: [array([0.77369911, 2.63363405, 1.07981163]), array([0.26038622, 0.00881044, 1.91060734])] \n",
"\n",
"Energy after step 35: 0.0000136 | Angles: [array([0.77369911, 2.63362729, 1.07981163]), array([0.26038622, 0.00557782, 1.91060734])] \n",
"\n",
"Energy after step 40: 0.0000054 | Angles: [array([0.77369911, 2.63362457, 1.07981163]), array([0.26038622, 0.00353126, 1.91060734])] \n",
"\n",
"Energy after step 45: 0.0000022 | Angles: [array([0.77369911, 2.63362348, 1.07981163]), array([0.26038622, 0.00223561, 1.91060734])] \n",
"\n",
"Energy after step 50: 0.0000009 | Angles: [array([0.77369911, 2.63362305, 1.07981163]), array([2.60386222e-01, 1.41534339e-03, 1.91060734e+00])] \n",
"\n",
"Energy after step 55: 0.0000004 | Angles: [array([0.77369911, 2.63362287, 1.07981163]), array([2.60386222e-01, 8.96040252e-04, 1.91060734e+00])] \n",
"\n",
"Energy after step 60: 0.0000001 | Angles: [array([0.77369911, 2.6336228 , 1.07981163]), array([2.60386222e-01, 5.67274421e-04, 1.91060734e+00])] \n",
"\n",
"Energy after step 65: 0.0000001 | Angles: [array([0.77369911, 2.63362278, 1.07981163]), array([2.60386222e-01, 3.59135947e-04, 1.91060734e+00])] \n",
"\n",
"Energy after step 70: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 2.27365491e-04, 1.91060734e+00])] \n",
"\n",
"Energy after step 75: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 1.43942891e-04, 1.91060734e+00])] \n",
"\n",
"Energy after step 80: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 9.11288509e-05, 1.91060734e+00])] \n",
"\n",
"Energy after step 85: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 5.76927932e-05, 1.91060734e+00])] \n",
"\n",
"Energy after step 90: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 3.65247488e-05, 1.91060734e+00])] \n",
"\n",
"Energy after step 95: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 2.31234648e-05, 1.91060734e+00])] \n",
"\n",
"Energy after step 100: 0.0000000 | Angles: [array([0.77369911, 2.63362276, 1.07981163]), array([2.60386222e-01, 1.46392417e-05, 1.91060734e+00])] \n",
"\n"
]
}
],
"source": [
"opt = torch.optim.SGD(var_init_loc, lr = 0.1)\n",
"\n",
"def closure():\n",
" opt.zero_grad()\n",
" loss = cost(p3, p4)\n",
" loss.backward()\n",
" return loss\n",
"\n",
"var_pt_loc = [var_init_loc]\n",
"cost_pt_loc = [cost_init_loc]\n",
"\n",
"for j in range(100):\n",
" opt.step(closure)\n",
" if (j + 1) % 5 == 0:\n",
" p3n, p4n = opt.param_groups[0]['params']\n",
" costn = cost(p3n, p4n)\n",
" var_pt_loc.append([p3n, p4n])\n",
" cost_pt_loc.append(costn)\n",
"\n",
" # for clarity, the angles are printed as numpy arrays\n",
" print('Energy after step {:5d}: {: .7f} | Angles: {}'.format(\n",
" j+1, costn, [p3n.detach().numpy(), p4n.detach().numpy()]),\"\\n\"\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": "<Figure size 600x400 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAFzCAYAAAD7bpkSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABR9ElEQVR4nO3de1xUZf4H8M9wHVAYUIQBREEl0bxgEBNka/6kHZNSqp+ZXUAztVZNY9MNU7TScDXLtehHVqttpRhp5JpLEXZZE1EuWpZ3MVxluKgMisptnt8fs5waBZxB4AzM5/16ndfMnPOcc75z9OV8POc5z1EIIQSIiIiIrIyd3AUQERERNYUhhYiIiKwSQwoRERFZJYYUIiIiskoMKURERGSVGFKIiIjIKjGkEBERkVViSCEiIiKr5CB3AZ2VwWDA2bNn4ebmBoVCIXc5REREnYYQAhcvXoSfnx/s7Jo/X8KQ0kpnz55FQECA3GUQERF1WqdPn0bv3r2bXc6Q0kpubm4AjAfY3d1d5mqIiIg6j6qqKgQEBEi/pc1hSGmlxks87u7uDClEREStcKPuEuw4S0RERFaJIYWIiIisEkMKERERWSX2SSEismFCCNTX16OhoUHuUqgLsbe3h4ODw00P0cGQQkRko2pra1FSUoLLly/LXQp1Qa6urvD19YWTk1Ort8GQQkRkgwwGA4qKimBvbw8/Pz84OTlxYEpqE0II1NbWory8HEVFRQgODm5xwLaWMKQQEdmg2tpaGAwGBAQEwNXVVe5yqItxcXGBo6Mjfv31V9TW1kKpVLZqO+w4S0Rkw1r7P1yiG2mLv1v820lERERWiZd7rMXFi8APPwCXLgH/+79yV0NERCQ7nkmxFmfOAPfeC0yfLnclRERW7e6778a8efM6bH+nTp2CQqHA/v3723S7GzZsgIeHh0XrdPR3lxvPpFgLPz/ja2UlcPkywI5sRERd2qRJkzBu3DiL1tm6dSscHR3bqSLrI/uZlJSUFAQGBkKpVEKj0WDv3r0ttk9PT0dISAiUSiWGDh2KHTt2mCwXQiApKQm+vr5wcXFBdHQ0jh07ZtLm6NGjmDBhAry8vODu7o6RI0fim2++afPvZhE3N6BbN+P7s2flrYWIiNqdi4sLvL29LVqnR48eN3xycFcia0jZvHkzEhISsGTJEhQUFGD48OHQarUoKytrsv3u3bsxefJkTJs2DYWFhYiNjUVsbCwOHjwotVm5ciXWrl2L1NRU5Obmolu3btBqtbh69arU5r777kN9fT127tyJ/Px8DB8+HPfddx90Ol27f+dmKRS/nU1hSCEiOQgBVFfLMwnR6rIvXLiAuLg4eHp6wtXVFffee+91/zn94YcfcPfdd8PV1RWenp7QarW4cOECACAzMxMjR46Eh4cHevbsifvuuw8nTpywqIbAwEAsW7YMcXFx6N69O/r27Ytt27ahvLwcEyZMQPfu3TFs2DDk5eVJ61x7uWfp0qUIDQ3Fhx9+iMDAQKhUKjzyyCO4ePGi1Obayz2t2W/jfn5vzZo1CAwMlD5PmTIFsbGxePXVV+Hj4wMPDw+8/PLLqK+vx/z589GjRw/07t0b69evt+g4WUrWkPL6669j+vTpmDp1KgYPHozU1FS4urri73//e5Pt//a3v2Hs2LGYP38+Bg0ahFdeeQW33XYb3nrrLQDGsyhr1qzBokWLMGHCBAwbNgz/+Mc/cPbsWWRkZAAAKioqcOzYMbzwwgsYNmwYgoODsWLFCly+fNkk7MiCIYWI5HT5MtC9uzzTTYx6O2XKFOTl5WHbtm3IycmBEALjxo1DXV0dAGD//v0YM2YMBg8ejJycHOzatQv333+/9CiA6upqJCQkIC8vD9nZ2bCzs8MDDzwAg8FgUR1vvPEG7rzzThQWFiImJgZPPPEE4uLi8Pjjj6OgoAD9+/dHXFwcRAuB7MSJE8jIyMD27duxfft2fPfdd1ixYkW777cpO3fuxNmzZ/H999/j9ddfx5IlS3DffffB09MTubm5ePrppzFz5kz85z//sWi7FhEyqampEfb29uKzzz4zmR8XFyfGjx/f5DoBAQHijTfeMJmXlJQkhg0bJoQQ4sSJEwKAKCwsNGnzhz/8QTz77LNCCCEMBoMYOHCgeOqpp8SlS5dEXV2dWLVqlfD29hbnz59vtt6rV68KvV4vTadPnxYAhF6vt+yLt2TyZCEAIVavbrttEhE14cqVK+KXX34RV65c+W3mpUvGf4PkmC5dMrv2UaNGiblz5wohhDh69KgAIH744QdpeUVFhXBxcRGffPKJEEKIyZMnizvvvNPs7ZeXlwsA4qeffhJCCFFUVNTkb8vv9e3bVzz++OPS55KSEgFALF68WJqXk5MjAIiSkhIhhBDr168XKpVKWr5kyRLh6uoqqqqqpHnz588XGo2mye/e2v0uWbJEDB8+3KT+N954Q/Tt21f6HB8fL/r27SsaGhqkeQMHDhR33XWX9Lm+vl5069ZNbNq0qclj0uTfsf/S6/Vm/YbK1nG2oqICDQ0N8PHxMZnv4+ODw4cPN7mOTqdrsn3jZZrG15baKBQKfP3114iNjYWbmxvs7Ozg7e2NzMxMeHp6NltvcnIyXnrpJcu+pKV4JoWI5OTqahwGQa59t8KhQ4fg4OAAjUYjzevZsycGDhyIQ4cOATCeSZk4cWKz2zh27BiSkpKQm5uLiooK6QxKcXExhgwZYnYtw4YNk943/g4NHTr0unllZWVQq9VNbiMwMNCkz4mvr2+zXSDacr9NufXWW00GZPPx8TE5Hvb29ujZs+cN67sZNnd3jxACs2bNgre3N/7973/DxcUF7733Hu6//37s27cPvr6+Ta6XmJiIhIQE6XNVVRUCAgLatjiGFCKSk0LxWwf+LsTFxaXF5ffffz/69u2Ld999F35+fjAYDBgyZAhqa2st2s/v77ppfA5SU/Nauox07Z07CoXihpedLN2vnZ3ddZd+Gi+N3aiW1tR3M2Trk+Ll5QV7e3uUlpaazC8tLW026anV6hbbN7621Gbnzp3Yvn070tLScOedd+K2227D22+/DRcXF3zwwQfN1uvs7Ax3d3eTqc35+xtfGVKIiMwyaNAg1NfXIzc3V5p37tw5HDlyBIMHDwZgPNOQnZ3d5PqNbRctWoQxY8Zg0KBBUofarqpXr17Q6XQmQaWtx4BpK7KFFCcnJ4SFhZn8xTEYDMjOzkZkZGST60RGRl73Fy0rK0tqHxQUBLVabdKmqqoKubm5UpvGR5Jf+0wBOzu7dk2DZuGZFCIiiwQHB2PChAmYPn06du3ahQMHDuDxxx+Hv78/JkyYAMB4Jnzfvn3405/+hB9//BGHDx/G//3f/6GiogKenp7o2bMn1q1bh+PHj2Pnzp0mZ827orvvvhvl5eVYuXIlTpw4gZSUFPzrX/+Su6wmyXp3T0JCAt5991188MEHOHToEJ555hlUV1dj6tSpAIC4uDgkJiZK7efOnYvMzEysXr0ahw8fxtKlS5GXl4fZs2cDMJ52mjdvHpYtW4Zt27bhp59+QlxcHPz8/BAbGwvAGHQ8PT0RHx+PAwcO4OjRo5g/fz6KiooQExPT4cfAxO9Dyk3cjkdEZEvWr1+PsLAw3HfffYiMjIQQAjt27JAuTdxyyy346quvcODAAURERCAyMhKff/45HBwcYGdnh7S0NOTn52PIkCF47rnnsGrVKpm/UfsaNGgQ3n77baSkpGD48OHYu3cvnn/+ebnLalqL3Wo7wJtvvin69OkjnJycREREhNizZ4+0bNSoUSI+Pt6k/SeffCJuueUW4eTkJG699VbxxRdfmCw3GAxi8eLFwsfHRzg7O4sxY8aII0eOmLTZt2+f+OMf/yh69Ogh3NzcxB133CF27NhhUd3m9ky2SHX1bz3d23K7RETXaOnOC6K20BZ39yiE4H/ZW6OqqgoqlQp6vb5t+6d4ehqHxj90CAgJabvtEhH9ztWrV1FUVISgoCAolUq5y6EuqKW/Y+b+hso+LD5do/GSz5kz8tZBREQkM4YUa8POs0RERAAYUqwPQwoREREAhhTrw5BCRB2I3RKpvbTF3y2GFGvDkEJEHaDx9tzLN/FgP6KWNP7dunaUWkvY3LD4Vo8hhYg6gL29PTw8PKTnrri6ukrDpxPdDCEELl++jLKyMnh4eMDe3r7V22JIsTYMKUTUQRofF9KeD4gj2+Xh4WHRAw2bwpBiba4ddZb/syGidqJQKODr6wtvb+8mHzBH1FqOjo43dQalEUOKtWl8CnNtLXD+PNCzp7z1EFGXZ29v3yY/KERtjR1nrY2TE9Crl/E9L/kQEZENY0ixRuyXQkRExJBilRhSiIiIGFKsEp/fQ0RExJBilXgmhYiIiCHFKjGkEBERMaRYJYYUIiIihhSrxJBCRETEkGKVGkOKTgc0NMhbCxERkUwYUqyRtzdgZ2cMKOXlcldDREQkC4YUa+TgAPj4GN/zkg8REdkohhRrxX4pRERk4xhSrJW/v/GVIYWIiGwUQ4q14pkUIiKycVYRUlJSUhAYGAilUgmNRoO9e/e22D49PR0hISFQKpUYOnQoduzYYbJcCIGkpCT4+vrCxcUF0dHROHbsmLT822+/hUKhaHLat29fu3xHizGkEBGRjZM9pGzevBkJCQlYsmQJCgoKMHz4cGi1WpSVlTXZfvfu3Zg8eTKmTZuGwsJCxMbGIjY2FgcPHpTarFy5EmvXrkVqaipyc3PRrVs3aLVaXL16FQAQFRWFkpISk+mpp55CUFAQwsPDO+R73xCf30NERDZOIYQQchag0Whw++2346233gIAGAwGBAQEYM6cOXjhhReuaz9p0iRUV1dj+/bt0rw77rgDoaGhSE1NhRACfn5++POf/4znn38eAKDX6+Hj44MNGzbgkUceuW6bdXV18Pf3x5w5c7B48WKz6q6qqoJKpYJer4e7u3trvnrL/vUvYNw4IDQUKCxs++0TERHJxNzfUFnPpNTW1iI/Px/R0dHSPDs7O0RHRyMnJ6fJdXJyckzaA4BWq5XaFxUVQafTmbRRqVTQaDTNbnPbtm04d+4cpk6derNfqe3wcg8REdk4Bzl3XlFRgYaGBvg0jgnyXz4+Pjh8+HCT6+h0uibb63Q6aXnjvObaXOv999+HVqtF7969m621pqYGNTU10ueqqqpm27aJxpBSVgbU1QGOju27PyIiIisje58Uuf3nP//Bl19+iWnTprXYLjk5GSqVSpoCAgLat7CePX8LJs2EKyIioq5M1pDi5eUFe3t7lJaWmswvLS2FWq1uch21Wt1i+8ZXc7e5fv169OzZE+PHj2+x1sTEROj1emk6ffp0y1/uZtnZAb6+xve85ENERDZI1pDi5OSEsLAwZGdnS/MMBgOys7MRGRnZ5DqRkZEm7QEgKytLah8UFAS1Wm3SpqqqCrm5uddtUwiB9evXIy4uDo43uJzi7OwMd3d3k6ndsV8KERHZMFn7pABAQkIC4uPjER4ejoiICKxZswbV1dVSJ9a4uDj4+/sjOTkZADB37lyMGjUKq1evRkxMDNLS0pCXl4d169YBABQKBebNm4dly5YhODgYQUFBWLx4Mfz8/BAbG2uy7507d6KoqAhPPfVUh35nszGkEBGRDZM9pEyaNAnl5eVISkqCTqdDaGgoMjMzpY6vxcXFsLP77YRPVFQUNm7ciEWLFmHhwoUIDg5GRkYGhgwZIrVZsGABqqurMWPGDFRWVmLkyJHIzMyEUqk02ff777+PqKgohISEdMyXtRRDChER2TDZx0nprNp9nBQAWLECSEwEpkwB1q9vn30QERF1sE4xTgrdAM+kEBGRDWNIsWYMKUREZMMYUqwZn99DREQ2jCHFmjWGlAsXgCtX5K2FiIiogzGkWDOVCnBxMb4vKZG3FiIiog7GkGLNFAr2SyEiIpvFkGLtGFKIiMhGMaRYO4YUIiKyUQwp1o4hhYiIbBRDirVjSCEiIhvFkGLtGFKIiMhGMaRYO4YUIiKyUQwp1s7f3/jKkEJERDaGIcXa+foaXy9eNE5EREQ2giHF2nXvDjQ+xppnU4iIyIYwpHQG7JdCREQ2iCGlM2BIISIiG8SQ0hkwpBARkQ1iSOkMGFKIiMgGMaR0BgwpRERkgxhSOgOGFCIiskEMKZ0BQwoREdkghpTO4PchRQh5ayEiIuogDCmdQeOos1evApWVspZCRETUUWQPKSkpKQgMDIRSqYRGo8HevXtbbJ+eno6QkBAolUoMHToUO3bsMFkuhEBSUhJ8fX3h4uKC6OhoHDt27LrtfPHFF9BoNHBxcYGnpydiY2Pb8mu1LaUS6NHD+J6XfIiIyEbIGlI2b96MhIQELFmyBAUFBRg+fDi0Wi3KysqabL97925MnjwZ06ZNQ2FhIWJjYxEbG4uDBw9KbVauXIm1a9ciNTUVubm56NatG7RaLa5evSq12bJlC5544glMnToVBw4cwA8//IBHH3203b/vTeGDBomIyNYIGUVERIhZs2ZJnxsaGoSfn59ITk5usv3DDz8sYmJiTOZpNBoxc+ZMIYQQBoNBqNVqsWrVKml5ZWWlcHZ2Fps2bRJCCFFXVyf8/f3Fe++9d1O16/V6AUDo9fqb2o7ZtFohACHWr++Y/REREbUTc39DZTuTUltbi/z8fERHR0vz7OzsEB0djZycnCbXycnJMWkPAFqtVmpfVFQEnU5n0kalUkGj0UhtCgoKcObMGdjZ2WHEiBHw9fXFvffea3I2pik1NTWoqqoymToU7/AhIiIbI1tIqaioQENDA3x8fEzm+/j4QKfTNbmOTqdrsX3ja0ttTp48CQBYunQpFi1ahO3bt8PT0xN33303zp8/32y9ycnJUKlU0hQQEGDBt20DDClERGRjZO8429EMBgMA4MUXX8RDDz2EsLAwrF+/HgqFAunp6c2ul5iYCL1eL02nT5/uqJKNGFKIiMjGyBZSvLy8YG9vj9LSUpP5paWlUKvVTa6jVqtbbN/42lIb3//ezjt48GBpubOzM/r164fi4uJm63V2doa7u7vJ1KEYUoiIyMbIFlKcnJwQFhaG7OxsaZ7BYEB2djYiIyObXCcyMtKkPQBkZWVJ7YOCgqBWq03aVFVVITc3V2oTFhYGZ2dnHDlyRGpTV1eHU6dOoW/fvm32/docQwoREdkYBzl3npCQgPj4eISHhyMiIgJr1qxBdXU1pk6dCgCIi4uDv78/kpOTAQBz587FqFGjsHr1asTExCAtLQ15eXlYt24dAEChUGDevHlYtmwZgoODERQUhMWLF8PPz08aB8Xd3R1PP/00lixZgoCAAPTt2xerVq0CAEycOLHjD4K5GkNKSQlgMAB2NneljoiIbIysIWXSpEkoLy9HUlISdDodQkNDkZmZKXV8LS4uht3vfoyjoqKwceNGLFq0CAsXLkRwcDAyMjIwZMgQqc2CBQtQXV2NGTNmoLKyEiNHjkRmZiaUSqXUZtWqVXBwcMATTzyBK1euQKPRYOfOnfD09Oy4L28pHx9AoQDq64GKCsDbW+6KiIiI2pVCCD4MpjWqqqqgUqmg1+s7rn+KWg2UlgKFhUBoaMfsk4iIqI2Z+xvKawadCfulEBGRDWFI6UwYUoiIyIYwpHQmDClERGRDGFI6k8aHDJ45I28dREREHYAhpTPhmRQiIrIhDCmdCUMKERHZEIaUzoQhhYiIbAhDSmfSGFJKS42DuhEREXVhDCmdSa9egL09IIQxqBAREXVhDCmdiZ0d8N+nOPOSDxERdXUMKZ0N+6UQEZGNYEjpbBhSiIjIRjCkdDYMKUREZCMsDinffPNNe9RB5mJIISIiG2FxSBk7diz69++PZcuW4fTp0+1RE7WEIYWIiGyExSHlzJkzmD17Nj799FP069cPWq0Wn3zyCWpra9ujProWn99DREQ2wuKQ4uXlheeeew779+9Hbm4ubrnlFvzpT3+Cn58fnn32WRw4cKA96qRGPJNCREQ24qY6zt52221ITEzE7NmzcenSJfz9739HWFgY7rrrLvz8889tVSP9XmNIOXcOqKmRtxYiIqJ21KqQUldXh08//RTjxo1D37598eWXX+Ktt95CaWkpjh8/jr59+2LixIltXSsBgKcn4OxsfF9SIm8tRERE7cjB0hXmzJmDTZs2QQiBJ554AitXrsSQIUOk5d26dcNrr70Gv8b/8VPbUiiMZ1OKioyXfAID5a6IiIioXVgcUn755Re8+eabePDBB+Hc+D/6a3h5efFW5fb0+5BCRETURVkcUrKzs2+8UQcHjBo1qlUFkRnYeZaIiGyAxSFl27ZtTc5XKBRQKpUYMGAAgoKCbrowagFDChER2QCLO87GxsbigQceQGxs7HWTVqvFgAEDMGrUKFy4cMHsbaakpCAwMBBKpRIajQZ79+5tsX16ejpCQkKgVCoxdOhQ7Nixw2S5EAJJSUnw9fWFi4sLoqOjcezYMZM2gYGBUCgUJtOKFSvMPxByYkghIiIbYHFIycrKwu23346srCzo9Xro9XpkZWVBo9Fg+/bt+P7773Hu3Dk8//zzZm1v8+bNSEhIwJIlS1BQUIDhw4dDq9WirKysyfa7d+/G5MmTMW3aNBQWFkoB6eDBg1KblStXYu3atUhNTUVubi66desGrVaLq1evmmzr5ZdfRklJiTTNmTPH0sMhD4YUIiKyBcJCt956q/jhhx+um79r1y4xePBgIYQQWVlZIiAgwKztRUREiFmzZkmfGxoahJ+fn0hOTm6y/cMPPyxiYmJM5mk0GjFz5kwhhBAGg0Go1WqxatUqaXllZaVwdnYWmzZtkub17dtXvPHGG2bV2BS9Xi8ACL1e3+pttFp2thCAEIMGdfy+iYiIbpK5v6EWn0k5ceIE3N3dr5vv7u6OkydPAgCCg4NRUVFxw23V1tYiPz8f0dHR0jw7OztER0cjJyenyXVycnJM2gOAVquV2hcVFUGn05m0UalU0Gg0121zxYoV6NmzJ0aMGIFVq1ahvr6+2VprampQVVVlMsmGZ1KIiMgGWBxSwsLCMH/+fJSXl0vzysvLsWDBAtx+++0AgGPHjiEgIOCG26qoqEBDQwN8fHxM5vv4+ECn0zW5jk6na7F94+uNtvnss88iLS0N33zzDWbOnIlXX30VCxYsaLbW5ORkqFQqaTLn+7WbxpCi1wPV1fLVQURE1I4svrvnvffeQ2xsLHr37i39UJ8+fRr9+vXD559/DgC4dOkSFi1a1LaVtrGEhATp/bBhw+Dk5ISZM2ciOTm5yfFfEhMTTdapqqqSL6i4uwPduwOXLhnPpgQHy1MHERFRO7I4pISEhOCXX37BV199haNHjwIABg4ciHvuuQd2dsYTM7GxsWZty8vLC/b29igtLTWZX1paCrVa3eQ6arW6xfaNr6WlpfD19TVpExoa2mwtGo0G9fX1OHXqFAYOHHjdcmdn52YHr5OFnx9w9ChDChERdVkWXe6pq6uDg4MDfvnlF4wdOxbPPvssnn32WWi1WimgWMLJyQlhYWEmA8QZDAZkZ2cjMjKyyXUiIyOvG1AuKytLah8UFAS1Wm3SpqqqCrm5uc1uEwD2798POzs7eHt7W/w9ZMF+KURE1MVZdCbF0dERffr0QUNDQ5sVkJCQgPj4eISHhyMiIgJr1qxBdXU1pk6dCgCIi4uDv78/kpOTAQBz587FqFGjsHr1asTExCAtLQ15eXlYt24dAOOgcvPmzcOyZcsQHByMoKAgLF68GH5+ftIZnpycHOTm5mL06NFwc3NDTk4OnnvuOTz++OPw9PRss+/WrhhSiIioi7P4cs+LL76IhQsX4sMPP0SPHj1uuoBJkyahvLwcSUlJ0Ol0CA0NRWZmptTxtbi42OQsTVRUFDZu3IhFixZh4cKFCA4ORkZGhslDDhcsWIDq6mrMmDEDlZWVGDlyJDIzM6FUKgEYL92kpaVh6dKlqKmpQVBQEJ577jmTPidWjyGFiIi6OIUQQliywogRI3D8+HHU1dWhb9++6Natm8nygoKCNi3QWlVVVUGlUkGv1zd5S3a7e+MNICEBeOQRYNOmjt8/ERFRK5n7G2rxmRRzO8VSO+OZFCIi6uIsDilLlixpjzrIUgwpRETUxVl+Sw6AyspKvPfee0hMTMT58+cBGC/znDlzpk2Loxb8PqRYdsWOiIioU7D4TMqPP/6I6OhoqFQqnDp1CtOnT0ePHj2wdetWFBcX4x//+Ed71EnXahwD5vJloKoKUKnkrYeIiKiNWXwmJSEhAVOmTMGxY8eku2UAYNy4cfj+++/btDhqgasr4OFhfM8zWERE1AVZHFL27duHmTNnXjff39+/2eftUDthvxQiIurCLA4pzs7OTT4B+OjRo+jVq1ebFEVmYkghIqIuzOKQMn78eLz88suoq6sDYBzhtbi4GH/5y1/w0EMPtXmB1AJ/f+MrQwoREXVBFoeU1atX49KlS/D29saVK1cwatQoDBgwAG5ubli+fHl71EjN4ZkUIiLqwiy+u0elUiErKwu7du3Cjz/+iEuXLuG2225DdHR0e9RHLWFIISKiLszikNJo5MiRGDlyZFvWQpZiSCEioi6sVSElOzsb2dnZKCsrg8FgMFn297//vU0KIzMwpBARURdmcUh56aWX8PLLLyM8PBy+vr5QKBTtUReZ49pRZ/lnQUREXYjFISU1NRUbNmzAE0880R71kCXUauNrXR1w7hzg5SVvPURERG3I4rt7amtrERUV1R61kKWcnIDGsWl4yYeIiLoYi0PKU089hY0bN7ZHLdQa7JdCRERdlMWXe65evYp169bh66+/xrBhw+Do6Giy/PXXX2+z4sgMfn7AgQN8fg8REXU5rXoKcmhoKADg4MGDJsvYiVYGPJNCRERdlMUh5ZtvvmmPOqi1GFKIiKiLsrhPSkvKysracnNkDoYUIiLqoswOKa6urigvL5c+x8TEoKSkRPpcWloKX1/ftq2ObowPGSQioi7K7JBy9epVCCGkz99//z2uXLli0ub3y6mD8EwKERF1UW16uYcdZ2XQGFJ0OqChQd5aiIiI2lCbhpTWSklJQWBgIJRKJTQaDfbu3dti+/T0dISEhECpVGLo0KHYsWOHyXIhBJKSkuDr6wsXFxdER0fj2LFjTW6rpqYGoaGhUCgU2L9/f1t9pY7j7Q3Y2QEGA8A+QURE1IWYHVIUCoXJmZJrP7fW5s2bkZCQgCVLlqCgoADDhw+HVqttthPu7t27MXnyZEybNg2FhYWIjY1FbGysye3QK1euxNq1a5Gamorc3Fx069YNWq0WV69evW57CxYsgF/j2YjOyN7+t+HxecmHiIi6EmEmhUIhPDw8hKenp/D09BQKhUKoVCrps4eHh7CzszN3c5KIiAgxa9Ys6XNDQ4Pw8/MTycnJTbZ/+OGHRUxMjMk8jUYjZs6cKYQQwmAwCLVaLVatWiUtr6ysFM7OzmLTpk0m6+3YsUOEhISIn3/+WQAQhYWFZtet1+sFAKHX681ep92EhwsBCLFtm9yVEBER3ZC5v6Fmj5Oyfv36Ng9ItbW1yM/PR2JiojTPzs4O0dHRyMnJaXKdnJwcJCQkmMzTarXIyMgAABQVFUGn0yE6OlparlKpoNFokJOTg0ceeQSA8W6k6dOnIyMjA66urm38zToYO88SEVEXZHZIiY+Pb/OdV1RUoKGhAT4+PibzfXx8cPjw4SbX0el0TbbX6XTS8sZ5zbURQmDKlCl4+umnER4ejlOnTt2w1pqaGtTU1Eifq6qqbrhOh2FIISKiLsgqOs52tDfffBMXL140OYNzI8nJyVCpVNIUEBDQjhVaqDGk8Pk9RETUhcgaUry8vGBvb4/S0lKT+aWlpVA3dga9hlqtbrF942tLbXbu3ImcnBw4OzvDwcEBAwYMAACEh4c3e8YoMTERer1emk6fPm3ht21HPJNCRERdkKwhxcnJCWFhYcjOzpbmGQwGZGdnIzIyssl1IiMjTdoDQFZWltQ+KCgIarXapE1VVRVyc3OlNmvXrsWBAwewf/9+7N+/X7qFefPmzVi+fHmT+3V2doa7u7vJZDUYUoiIqAuy+AGDbS0hIQHx8fEIDw9HREQE1qxZg+rqakydOhUAEBcXB39/fyQnJwMA5s6di1GjRmH16tWIiYlBWloa8vLysG7dOgDGW6PnzZuHZcuWITg4GEFBQVi8eDH8/PwQGxsLAOjTp49JDd27dwcA9O/fH7179+6gb96GGFKIiKgLMjuk3HXXXZgwYQLGjx+PW265pc0KmDRpEsrLy5GUlASdTofQ0FBkZmZKHV+Li4thZ/fbCZ+oqChs3LgRixYtwsKFCxEcHIyMjAwMGTJEarNgwQJUV1djxowZqKysxMiRI5GZmQmlUtlmdVuVxuf3lJcDtbWAk5O89RAREbUBhRDmPXDnH//4Bz7//HN89dVX6N27N8aPH4/x48cjKirKJofDr6qqgkqlgl6vl//SjxCAszNQVwf8+itwzZkiIiIia2Lub6jZfVLi4uKwZcsWVFRUYPXq1aisrMTEiROhVqvx5JNPIiMj47oHDlIHUSh4yYeIiLocizvOOjs7Y9y4cXjnnXdw9uxZbNu2Db6+vli8eDF69uyJ++67Dz/88EN71EotYUghIqIu5qbv7tFoNFi+fDl++ukn/PTTTxgzZgxKSkraojayBEMKERF1MW16d0///v3x3HPPteUmyVwMKURE1MXY5IizXRJDChERdTEMKV0FQwoREXUxDCldBZ/fQ0REXYzFIeXll1/G5cuXr5t/5coVvPzyy21SFLUCz6QQEVEXY/Zgbo3s7e1RUlICb29vk/nnzp2Dt7c3Ghoa2rRAa2VVg7kBQGUl4OlpfF9dDbi6yloOERFRc9p8MLdGQogmR5g9cOAAevToYenmqK2oVICLi/E9bwEnIqIuwOxbkD09PaFQKKBQKHDLLbeYBJWGhgZcunQJTz/9dLsUSWZoHHX2xAnjJZ/+/eWuiIiI6KaYHVLWrFkDIQSefPJJvPTSS1CpVNIyJycnBAYGIjIysl2KJDP5+/8WUoiIiDo5s0NKfHw8ACAoKAh33nknHBzadBw4agvsPEtERF2IxX1S3NzccOjQIenz559/jtjYWCxcuBC1tbVtWhxZiCGFiIi6EItDysyZM3H06FEAwMmTJzFp0iS4uroiPT0dCxYsaPMCyQIMKURE1IVYHFKOHj2K0NBQAEB6ejpGjRqFjRs3YsOGDdiyZUtb10eWYEghIqIupFW3IBsMBgDA119/jXHjxgEAAgICUFFR0bbVkWUYUoiIqAuxOKSEh4dj2bJl+PDDD/Hdd98hJiYGAFBUVAQfH582L5AswJBCRERdiMUhZc2aNSgoKMDs2bPx4osvYsCAAQCATz/9FFFRUW1eIFnA3x+wswMuXQJOnpS7GiIiopti8bD4zbl69Srs7e3h6OjYFpuzelY3LH6j0aOBb78FVq8GEhLkroaIiOg67TYsfqP8/Hx89NFH+Oijj1BQUAClUmkzAcWqPfig8fWzz+Stg4iI6CZZfCalrKwMkyZNwnfffQcPDw8AQGVlJUaPHo20tDT06tWrPeq0OlZ7JuX0aaBPH+Mw+SUlAPsJERGRlWm3Mylz5szBpUuX8PPPP+P8+fM4f/48Dh48iKqqKjz77LM3VTS1gYAAIDwcEAL4/HO5qyEiImo1i0NKZmYm3n77bQwaNEiaN3jwYKSkpOBf//pXmxZHrfTAA8ZXXvIhIqJOzOKQYjAYmux74ujoKI2fYqmUlBQEBgZCqVRCo9Fg7969LbZPT09HSEgIlEolhg4dih07dpgsF0IgKSkJvr6+cHFxQXR0NI4dO2bSZvz48ejTpw+USiV8fX3xxBNP4GxXuXW3sV9Kdjag18tbCxERUStZHFL+53/+B3PnzjX5QT9z5gyee+45jBkzxuICNm/ejISEBCxZsgQFBQUYPnw4tFotysrKmmy/e/duTJ48GdOmTUNhYSFiY2MRGxuLgwcPSm1WrlyJtWvXIjU1Fbm5uejWrRu0Wi2uXr0qtRk9ejQ++eQTHDlyBFu2bMGJEyfwv//7vxbXb5VCQoxTXR1wTYAjIiLqNISFiouLRWhoqHB0dBT9+vUT/fr1E46OjmLEiBHi9OnTlm5OREREiFmzZkmfGxoahJ+fn0hOTm6y/cMPPyxiYmJM5mk0GjFz5kwhhBAGg0Go1WqxatUqaXllZaVwdnYWmzZtaraOzz//XCgUClFbW2tW3Xq9XgAQer3erPYdLjFRCECIiRPlroSIiMiEub+hFp9JCQgIQEFBAb744gvMmzcP8+bNw44dO1BQUIDevXtbtK3a2lrk5+cjOjpammdnZ4fo6Gjk5OQ0uU5OTo5JewDQarVS+6KiIuh0OpM2KpUKGo2m2W2eP38eH3/8MaKiopq9jbqmpgZVVVUmk1Vr7JeyYwdw5Yq8tRAREbVCq8ZJUSgUuOeeezBnzhzMmTPnutBgroqKCjQ0NFw3nL6Pjw90Ol2T6+h0uhbbN76as82//OUv6NatG3r27Ini4mJ83sLdMMnJyVCpVNIUEBBg3peUS3g40Ls3UF0NfP213NUQERFZzOyQsnPnTgwePLjJMwh6vR633nor/v3vf7dpce1t/vz5KCwsxFdffQV7e3vExcVBNDNsTGJiIvR6vTSdPn26g6u1kELBu3yIiKhTMzukrFmzBtOnT29y0BWVSoWZM2fi9ddft2jnXl5esLe3R2lpqcn80tJSqNXqJtdRq9Uttm98NWebXl5euOWWW3DPPfcgLS0NO3bswJ49e5rcr7OzM9zd3U0mq9cYUrZtA+rr5a2FiIjIQmaHlAMHDmDs2LHNLv/jH/+I/Px8i3bu5OSEsLAwZGdnS/MMBgOys7MRGRnZ5DqRkZEm7QEgKytLah8UFAS1Wm3SpqqqCrm5uc1us3G/gLHvSZdx111Az57AuXPArl1yV0NERGQRs0NKaWlpi8/mcXBwQHl5ucUFJCQk4N1338UHH3yAQ4cO4ZlnnkF1dTWmTp0KAIiLi0NiYqLUfu7cucjMzMTq1atx+PBhLF26FHl5eZg9ezYAY3+ZefPmYdmyZdi2bRt++uknxMXFwc/PD7GxsQCA3NxcvPXWW9i/fz9+/fVX7Ny5E5MnT0b//v1bDDKdjoMDcP/9xvdbt8pbCxERkYXMDin+/v4mY5Fc68cff4Svr6/FBUyaNAmvvfYakpKSEBoaiv379yMzM1Pq+FpcXIySkhKpfVRUFDZu3Ih169Zh+PDh+PTTT5GRkYEhQ4ZIbRYsWIA5c+ZgxowZuP3223Hp0iVkZmZCqVQCAFxdXbF161aMGTMGAwcOxLRp0zBs2DB89913cHZ2tvg7WLXGSz4ZGcah8omIiDoJsx8wOGfOHHz77bfYt2+f9GPf6MqVK4iIiMDo0aOxdu3adinU2ljtAwavdeUK0KuX8S6fffuMd/0QERHJyNzfULNDSmlpKW677TbY29tj9uzZGDhwIADg8OHDSElJQUNDAwoKCq679ber6jQhBQAmTgQ+/RRYuBBYvlzuaoiIyMa1eUgBgF9//RXPPPMMvvzyS+lWXYVCAa1Wi5SUFAQFBd185Z1EpwopGzcCjz1mHCr/0CG5qyEiIhvXLiGl0YULF3D8+HEIIRAcHAxPT8+bKrYz6lQhRa83XvKpqzOGlJAQuSsiIiIbZu5vaKtGnPX09MTtt9+OiIgImwwonY5KBTQ+/JEDuxERUSfRqpBCnRBHnyUiok6GIcVWTJhgHCp/3z7A2of0JyIiAkOK7fDxAaKijO8zMmQthYiIyBwMKbaEl3yIiKgTYUixJY0h5fvvjc/zISIismIMKbakXz9g+HCgoQH45z/lroaIiKhFDCm2hpd8iIiok2BIsTWNIeXLL4FLl+SthYiIqAUMKbZm6FCgf3+gpgbIzJS7GiIiomYxpNgahYKXfIiIqFNgSLFFjSHliy+A2lp5ayEiImoGQ4otuuMOQK02Pnjwm2/kroaIiKhJDCm2yM7OOEw+AGzdKm8tREREzWBIsVUPPmh8/fxz47gpREREVoYhxVbdfTegUgGlpcCePXJXQ0REdB2GFFvl5ATcd5/xPe/yISIiK8SQYst+fyuyEPLWQkREdA2GFFs2diygVAInTwI//ih3NURERCYYUmxZt26AVmt8z0s+RERkZawipKSkpCAwMBBKpRIajQZ79+5tsX16ejpCQkKgVCoxdOhQ7Nixw2S5EAJJSUnw9fWFi4sLoqOjcezYMWn5qVOnMG3aNAQFBcHFxQX9+/fHkiVLUGuLA5tx9FkiIrJSsoeUzZs3IyEhAUuWLEFBQQGGDx8OrVaLsrKyJtvv3r0bkydPxrRp01BYWIjY2FjExsbi4MGDUpuVK1di7dq1SE1NRW5uLrp16watVourV68CAA4fPgyDwYB33nkHP//8M9544w2kpqZi4cKFHfKdrcr99wP29sbLPSdPyl0NERHRb4TMIiIixKxZs6TPDQ0Nws/PTyQnJzfZ/uGHHxYxMTEm8zQajZg5c6YQQgiDwSDUarVYtWqVtLyyslI4OzuLTZs2NVvHypUrRVBQkNl16/V6AUDo9Xqz17Fa//M/QgBCvPaa3JUQEZENMPc3VNYzKbW1tcjPz0d0dLQ0z87ODtHR0cjJyWlynZycHJP2AKDVaqX2RUVF0Ol0Jm1UKhU0Gk2z2wQAvV6PHj16NLu8pqYGVVVVJlOX0TiwG0efJSIiKyJrSKmoqEBDQwN8fHxM5vv4+ECn0zW5jk6na7F946sl2zx+/DjefPNNzJw5s9lak5OToVKppCkgIKDlL9eZxMYaX3NygGaOERERUUeTvU+K3M6cOYOxY8di4sSJmD59erPtEhMTodfrpen06dMdWGU78/cHIiKMY6V8/rnc1RAREQGQOaR4eXnB3t4epaWlJvNLS0uhVqubXEetVrfYvvHVnG2ePXsWo0ePRlRUFNatW9dirc7OznB3dzeZuhTe5UNERFZG1pDi5OSEsLAwZGdnS/MMBgOys7MRGRnZ5DqRkZEm7QEgKytLah8UFAS1Wm3SpqqqCrm5uSbbPHPmDO6++26EhYVh/fr1sLOz8ZNKjSFl506gslLWUoiIiAAruNyTkJCAd999Fx988AEOHTqEZ555BtXV1Zg6dSoAIC4uDomJiVL7uXPnIjMzE6tXr8bhw4exdOlS5OXlYfbs2QAAhUKBefPmYdmyZdi2bRt++uknxMXFwc/PD7H/7XvRGFD69OmD1157DeXl5dDpdM32WbEJAwcCgwcDdXXAF1/IXQ0REREc5C5g0qRJKC8vR1JSEnQ6HUJDQ5GZmSl1fC0uLjY5yxEVFYWNGzdi0aJFWLhwIYKDg5GRkYEhQ4ZIbRYsWIDq6mrMmDEDlZWVGDlyJDIzM6FUKgEYz7wcP34cx48fR+/evU3qEbb8DJsHHgB++cV4yeexx+SuhoiIbJxC2PSvcutVVVVBpVJBr9d3nf4p+flAeDjg6gpUVAAuLnJXREREXZC5v6GyX+4hK3LbbUCfPsDly0BWltzVEBGRjWNIod8oFLzLh4iIrAZDCplqDCnbtgH19fLWQkRENo0hhUyNHAl4eQHnzwPffy93NUREZMMYUsiUvT0wfrzxPS/5EBGRjBhS6HqNl3wyMoxD5RMREcmAIYWuFx0NdO8O/Oc/QF6e3NUQEZGNYkih6ymVwLhxxvdbt8pbCxER2SyGFGoab0UmIiKZMaRQ08aNA5ycgCNHgEOH5K6GiIhsEEMKNc3dHRgzxvieZ1OIiEgGDCnUvAcfNL6yXwoREcmAIYWaN348YGdnfPBgcbHc1RARkY1hSKHmeXsDd95pfJ+RIWspRERkexhSqGW8y4eIiGTCkEItawwp338PVFTIWwsREdkUhhRqWWAgMGIEYDAAf/2r3NUQEZENYUihG1u0yPj62mtAerq8tRARkc1gSKEbe/BBYP584/upU4Gff5a3HiIisgkMKWSeV181Du5WXW3sp1JZKXdFRETUxTGkkHkcHIBNm4A+fYBjx4C4OGM/FSIionbCkELm69XLOPqsszPwz38Cy5bJXREREXVhDClkmbAwIDXV+H7pUuCLL2Qth4iIui7ZQ0pKSgoCAwOhVCqh0Wiwd+/eFtunp6cjJCQESqUSQ4cOxY4dO0yWCyGQlJQEX19fuLi4IDo6GseOHTNps3z5ckRFRcHV1RUeHh5t/ZW6vilTgGeeAYQAHnsMOH5c7oqIiKgLkjWkbN68GQkJCViyZAkKCgowfPhwaLValJWVNdl+9+7dmDx5MqZNm4bCwkLExsYiNjYWBw8elNqsXLkSa9euRWpqKnJzc9GtWzdotVpcvXpValNbW4uJEyfimWeeaffv2GWtWQNERgJ6vbEjbXW13BUREVFXI2QUEREhZs2aJX1uaGgQfn5+Ijk5ucn2Dz/8sIiJiTGZp9FoxMyZM4UQQhgMBqFWq8WqVauk5ZWVlcLZ2Vls2rTpuu2tX79eqFSqVtWu1+sFAKHX61u1fpdw5owQarUQgBCTJglhMMhdERERdQLm/obKdialtrYW+fn5iI6OlubZ2dkhOjoaOTk5Ta6Tk5Nj0h4AtFqt1L6oqAg6nc6kjUqlgkajaXab5qqpqUFVVZXJZPP8/IyDuzk4AJs3A2+8IXdFRETUhcgWUioqKtDQ0AAfHx+T+T4+PtDpdE2uo9PpWmzf+GrJNs2VnJwMlUolTQEBATe1vS5j5MjfwsmCBcA338hbDxERdRmyd5ztLBITE6HX66Xp9OnTcpdkPWbNMo6b0tAATJoE8NgQEVEbkC2keHl5wd7eHqWlpSbzS0tLoVarm1xHrVa32L7x1ZJtmsvZ2Rnu7u4mE/2XQmG8LTk0FCgvBx56CPhdR2UiIqLWkC2kODk5ISwsDNnZ2dI8g8GA7OxsREZGNrlOZGSkSXsAyMrKktoHBQVBrVabtKmqqkJubm6z26Q24uICfPYZ0KMHsG8fMHu28RZlIiKiVnKQc+cJCQmIj49HeHg4IiIisGbNGlRXV2Pq1KkAgLi4OPj7+yM5ORkAMHfuXIwaNQqrV69GTEwM0tLSkJeXh3Xr1gEAFAoF5s2bh2XLliE4OBhBQUFYvHgx/Pz8EBsbK+23uLgY58+fR3FxMRoaGrB//34AwIABA9C9e/cOPQZdSmAgkJYGjB0LvP8+EBEBzJghd1VERNRZddDdRs168803RZ8+fYSTk5OIiIgQe/bskZaNGjVKxMfHm7T/5JNPxC233CKcnJzErbfeKr744guT5QaDQSxevFj4+PgIZ2dnMWbMGHHkyBGTNvHx8QLAddM333xjdt28BbkFycnG25IdHYXIyZG7GiIisjLm/oYqhOA5+daoqqqCSqWCXq9n/5RrCQFMnAhs2WK8TbmgALjmjisiIrJd5v6G8u4eansKBbB+PTBoEHD2rDGw1NXJXRUREXUyDCnUPtzcjB1p3dyAf/8bmD9f7oqIiKiTYUih9jNwIPDhh8b3f/sb8NFH8tZDRESdCkMKta8JE4BFi4zvZ8wA/nsnFRER0Y0wpFD7W7rUeFvylSvAgw8C58/LXREREXUCDCnU/uztgY8/Bvr1A4qKgEcfNQ6hT0RE1AKGFOoYPXoAW7caR6b98ksgKUnuioiIyMoxpFDHGT4ceO894/tXXzXe/UNERNQMhhTqWI8+CsybZ3wfHw8cPChrOUREZL0YUqjjrVwJjBoFXLwI3HYb8MwzwOnTcldFRERWhiGFOp6jI/DJJ8Af/2gciTY1FRgwwPjk5DNn5K6OiIisBEMKycPb29iB9vvvgdGjgdpaICUF6N/feDmopETuComISGYMKSSvu+4Cdu40TnfdBdTUGEen7dcP+POfgdJSuSskIiKZMKSQdRg9GvjuOyArC4iMBK5eBV5/3RhW/vIXoKJC7gqJiKiDMaSQ9VAogOho4IcfgMxMICICuHzZ2NE2MBBYuBA4d07uKomIqIMwpJD1USgArRbYswfYvh0ICwOqq4HkZCAoyDgQ3IULcldJRETtjCGFrJdCAcTEAPv2AZ9/bhwM7uJF4JVXjGHlpZcAvV7uKomIqJ0wpJD1UyiA8eOBggJgyxZgyBBjOFm61HgZaPlyY3ghIqIuhSGFOg87O+NTlA8cADZvBgYNAiorgUWLjGFlxQrg0iW5qyQiojbCkEKdj50d8PDDwE8/ARs3ArfcApw/DyQmGi8DLV4M/OtfQHm53JUSEdFNUAghhNxFdEZVVVVQqVTQ6/Vwd3eXuxzbVl8PbNpk7KNy4oTpsj59gPDw36awMOMTmYmISDbm/oYypLQSQ4oVqq8H0tKMI9nm5QFHjgBN/fXu1+/64MI/QyKiDsOQ0s4YUjqBqiqgsNAYWBqn48ebbjtwoGlwCQ0Funfv0HKJiGwFQ0o7Y0jppC5cMN4l9PvgcurU9e3s7IwdcxvPtIwYAfj6Ar16AW5uxjuOiIioVcz9DbWKjrMpKSkIDAyEUqmERqPB3r17W2yfnp6OkJAQKJVKDB06FDt27DBZLoRAUlISfH194eLigujoaBw7dsykzfnz5/HYY4/B3d0dHh4emDZtGi7xzpCuz9MTGDPGONR+ejpQVGTsYJuZCSxbBsTGAr17AwYD8PPPwAcfAM8+a3yu0IABgEoFKJXGNqGhwD33AI8+Csyda1z/nXeArVuND048dMg4nL/BIPe3JiLqlGQ/k7J582bExcUhNTUVGo0Ga9asQXp6Oo4cOQJvb+/r2u/evRt/+MMfkJycjPvuuw8bN27EX//6VxQUFGDIkCEAgL/+9a9ITk7GBx98gKCgICxevBg//fQTfvnlFyiVSgDAvffei5KSErzzzjuoq6vD1KlTcfvtt2Pjxo1m1c0zKV2cTgfk5/92tuXgQaCszDhMv6Xs7ICePY1nYby8jK+NU/fugIuLMfgolea/d3Y2bpeIqBPqNJd7NBoNbr/9drz11lsAAIPBgICAAMyZMwcvvPDCde0nTZqE6upqbN++XZp3xx13IDQ0FKmpqRBCwM/PD3/+85/x/PPPAwD0ej18fHywYcMGPPLIIzh06BAGDx6Mffv2ITw8HACQmZmJcePG4T//+Q/8/PxuWDdDio26fNl4dqS83HRqal55uXEcl/bi5NR0kHFyAhwcAHt749RW7+3sjJe5FIobv7ekLdD0a0vLzG3ze+bOs7StuevfTDtL8XIktTVfX+PDX9uIub+hDm22x1aora1Ffn4+EhMTpXl2dnaIjo5GTk5Ok+vk5OQgISHBZJ5Wq0VGRgYAoKioCDqdDtHR0dJylUoFjUaDnJwcPPLII8jJyYGHh4cUUAAgOjoadnZ2yM3NxQMPPHDdfmtqalBTUyN9rqqqatV3pk7O1dV4W3OfPua1r6szPhSxuUBTXW184vPVq8CVKy2/v3IFaGj4bdu1tcaJjwYgovYWGwt89lmH71bWkFJRUYGGhgb4+PiYzPfx8cHhw4ebXEen0zXZXqfTScsb57XU5tpLSQ4ODujRo4fU5lrJycl46aWXzPxmRP/l6Aio1capLdTXtxxirl4FamqMYaZxqq+//n1T8260XAhj/xohfpt+/7k174HfbhP//WtT8yxt83vmzrO0rbnr30w7S9n6vRC2/v3by6BBsuxW1pDSmSQmJpqcwamqqkJAQICMFZFNcnAw9mPh7dFEZANk7Xnn5eUFe3t7lJaWmswvLS2Fupn/earV6hbbN77eqE1ZWZnJ8vr6epw/f77Z/To7O8Pd3d1kIiIiovYja0hxcnJCWFgYsrOzpXkGgwHZ2dmIbKaDTmRkpEl7AMjKypLaBwUFQa1Wm7SpqqpCbm6u1CYyMhKVlZXIz8+X2uzcuRMGgwEajabNvh8RERG1nuyXexISEhAfH4/w8HBERERgzZo1qK6uxtSpUwEAcXFx8Pf3R3JyMgBg7ty5GDVqFFavXo2YmBikpaUhLy8P69atAwAoFArMmzcPy5YtQ3BwsHQLsp+fH2JjYwEAgwYNwtixYzF9+nSkpqairq4Os2fPxiOPPGLWnT1ERETU/mQPKZMmTUJ5eTmSkpKg0+kQGhqKzMxMqeNrcXEx7H43HkRUVBQ2btyIRYsWYeHChQgODkZGRoY0RgoALFiwANXV1ZgxYwYqKysxcuRIZGZmSmOkAMDHH3+M2bNnY8yYMbCzs8NDDz2EtWvXdtwXJyIiohbJPk5KZ8VxUoiIiFqnUw2LT0RERHQthhQiIiKySgwpREREZJUYUoiIiMgqMaQQERGRVWJIISIiIqsk+zgpnVXjndt8GjIREZFlGn87bzQKCkNKK128eBEA+JBBIiKiVrp48SJUKlWzyzmYWysZDAacPXsWbm5uUCgUbbLNxicrnz59mgPEtREe0/bB49r2eEzbB49r22uLYyqEwMWLF+Hn52cyqvy1eCallezs7NC7d+922Tafstz2eEzbB49r2+MxbR88rm3vZo9pS2dQGrHjLBEREVklhhQiIiKySgwpVsTZ2RlLliyBs7Oz3KV0GTym7YPHte3xmLYPHte215HHlB1niYiIyCrxTAoRERFZJYYUIiIiskoMKURERGSVGFKIiIjIKjGkWImUlBQEBgZCqVRCo9Fg7969cpfUqSQnJ+P222+Hm5sbvL29ERsbiyNHjpi0uXr1KmbNmoWePXuie/fueOihh1BaWipTxZ3PihUroFAoMG/ePGkej6nlzpw5g8cffxw9e/aEi4sLhg4diry8PGm5EAJJSUnw9fWFi4sLoqOjcezYMRkrtn4NDQ1YvHgxgoKC4OLigv79++OVV14xeS4Mj2vLvv/+e9x///3w8/ODQqFARkaGyXJzjt/58+fx2GOPwd3dHR4eHpg2bRouXbp0c4UJkl1aWppwcnISf//738XPP/8spk+fLjw8PERpaancpXUaWq1WrF+/Xhw8eFDs379fjBs3TvTp00dcunRJavP000+LgIAAkZ2dLfLy8sQdd9whoqKiZKy689i7d68IDAwUw4YNE3PnzpXm85ha5vz586Jv375iypQpIjc3V5w8eVJ8+eWX4vjx41KbFStWCJVKJTIyMsSBAwfE+PHjRVBQkLhy5YqMlVu35cuXi549e4rt27eLoqIikZ6eLrp37y7+9re/SW14XFu2Y8cO8eKLL4qtW7cKAOKzzz4zWW7O8Rs7dqwYPny42LNnj/j3v/8tBgwYICZPnnxTdTGkWIGIiAgxa9Ys6XNDQ4Pw8/MTycnJMlbVuZWVlQkA4rvvvhNCCFFZWSkcHR1Fenq61ObQoUMCgMjJyZGrzE7h4sWLIjg4WGRlZYlRo0ZJIYXH1HJ/+ctfxMiRI5tdbjAYhFqtFqtWrZLmVVZWCmdnZ7Fp06aOKLFTiomJEU8++aTJvAcffFA89thjQggeV0tdG1LMOX6//PKLACD27dsntfnXv/4lFAqFOHPmTKtr4eUemdXW1iI/Px/R0dHSPDs7O0RHRyMnJ0fGyjo3vV4PAOjRowcAID8/H3V1dSbHOSQkBH369OFxvoFZs2YhJibG5NgBPKatsW3bNoSHh2PixInw9vbGiBEj8O6770rLi4qKoNPpTI6pSqWCRqPhMW1BVFQUsrOzcfToUQDAgQMHsGvXLtx7770AeFxvljnHLycnBx4eHggPD5faREdHw87ODrm5ua3eNx8wKLOKigo0NDTAx8fHZL6Pjw8OHz4sU1Wdm8FgwLx583DnnXdiyJAhAACdTgcnJyd4eHiYtPXx8YFOp5Ohys4hLS0NBQUF2Ldv33XLeEwtd/LkSfzf//0fEhISsHDhQuzbtw/PPvssnJycEB8fLx23pv494DFt3gsvvICqqiqEhITA3t4eDQ0NWL58OR577DEA4HG9SeYcP51OB29vb5PlDg4O6NGjx00dY4YU6nJmzZqFgwcPYteuXXKX0qmdPn0ac+fORVZWFpRKpdzldAkGgwHh4eF49dVXAQAjRozAwYMHkZqaivj4eJmr67w++eQTfPzxx9i4cSNuvfVW7N+/H/PmzYOfnx+PayfHyz0y8/Lygr29/XV3RJSWlkKtVstUVec1e/ZsbN++Hd988w169+4tzVer1aitrUVlZaVJex7n5uXn56OsrAy33XYbHBwc4ODggO+++w5r166Fg4MDfHx8eEwt5Ovri8GDB5vMGzRoEIqLiwFAOm7898Ay8+fPxwsvvIBHHnkEQ4cOxRNPPIHnnnsOycnJAHhcb5Y5x0+tVqOsrMxkeX19Pc6fP39Tx5ghRWZOTk4ICwtDdna2NM9gMCA7OxuRkZEyVta5CCEwe/ZsfPbZZ9i5cyeCgoJMloeFhcHR0dHkOB85cgTFxcU8zs0YM2YMfvrpJ+zfv1+awsPD8dhjj0nveUwtc+edd153a/zRo0fRt29fAEBQUBDUarXJMa2qqkJubi6PaQsuX74MOzvTnzN7e3sYDAYAPK43y5zjFxkZicrKSuTn50ttdu7cCYPBAI1G0/qdt7rLLbWZtLQ04ezsLDZs2CB++eUXMWPGDOHh4SF0Op3cpXUazzzzjFCpVOLbb78VJSUl0nT58mWpzdNPPy369Okjdu7cKfLy8kRkZKSIjIyUserO5/d39wjBY2qpvXv3CgcHB7F8+XJx7Ngx8fHHHwtXV1fx0UcfSW1WrFghPDw8xOeffy5+/PFHMWHCBN4qewPx8fHC399fugV569atwsvLSyxYsEBqw+PasosXL4rCwkJRWFgoAIjXX39dFBYWil9//VUIYd7xGzt2rBgxYoTIzc0Vu3btEsHBwbwFuat48803RZ8+fYSTk5OIiIgQe/bskbukTgVAk9P69eulNleuXBF/+tOfhKenp3B1dRUPPPCAKCkpka/oTujakMJjarl//vOfYsiQIcLZ2VmEhISIdevWmSw3GAxi8eLFwsfHRzg7O4sxY8aII0eOyFRt51BVVSXmzp0r+vTpI5RKpejXr5948cUXRU1NjdSGx7Vl33zzTZP/hsbHxwshzDt+586dE5MnTxbdu3cX7u7uYurUqeLixYs3VZdCiN8NyUdERERkJdgnhYiIiKwSQwoRERFZJYYUIiIiskoMKURERGSVGFKIiIjIKjGkEBERkVViSCEiIiKrxJBCRK22dOlShIaG3tQ2Tp06BYVCgf3797dJTc25++67MW/evHbdBxG1LYYUoi7s9OnTePLJJ+Hn5wcnJyf07dsXc+fOxblz5yzelkKhQEZGhsm8559/3uR5Hq0REBCAkpISDBky5Ka20+jbb7+FQqG47sGHW7duxSuvvNIm+2iNjgpjRF0JQwpRF3Xy5EmEh4fj2LFj2LRpE44fP47U1FTp4ZXnz5+/6X10794dPXv2vKlt2NvbQ61Ww8HB4abraUmPHj3g5ubWrvsgorbFkELURc2aNQtOTk746quvMGrUKPTp0wf33nsvvv76a5w5cwYvvvii1DYwMBCvvPIKJk+ejG7dusHf3x8pKSkmywHggQcegEKhkD5fe7lnypQpiI2NxauvvgofHx94eHjg5ZdfRn19PebPn48ePXqgd+/eWL9+vbTOtWcYpkyZAoVCcd307bffAgA+/PBDhIeHw83NDWq1Go8++qj0iPhTp05h9OjRAABPT08oFApMmTIFwPWXey5cuIC4uDh4enrC1dUV9957L44dOyYt37BhAzw8PPDll19i0KBB6N69O8aOHYuSkpJmj/mFCxfw2GOPoVevXnBxcUFwcLD0XRufzD1ixAgoFArcfffd0nrvvfceBg0aBKVSiZCQELz99tvXHZ+0tDRERUVBqVRiyJAh+O6778zaL1GndlNP/iEiq3Tu3DmhUCjEq6++2uTy6dOnC09PT2EwGIQQQvTt21e4ubmJ5ORkceTIEbF27Vphb28vvvrqKyGEEGVlZdIDG0tKSkRZWZkQQoglS5aI4cOHS9uNj48Xbm5uYtasWeLw4cPi/fffFwCEVqsVy5cvF0ePHhWvvPKKcHR0FKdPnxZCCFFUVCQAiMLCQiGEEJWVlSZPsp47d67w9vaWHlz4/vvvix07dogTJ06InJwcERkZKe69914hhBD19fViy5YtAoA4cuSIKCkpEZWVlUKI6x+OOH78eDFo0CDx/fffi/379wutVisGDBggamtrhRBCrF+/Xjg6Ooro6Gixb98+kZ+fLwYNGiQeffTRZo/7rFmzRGhoqNi3b58oKioSWVlZYtu2bUII4xOQAYivv/5alJSUiHPnzgkhhPjoo4+Er6+v2LJlizh58qTYsmWL6NGjh9iwYYPJ8endu7f49NNPxS+//CKeeuop4ebmJioqKm64X6LOjCGFqAvas2ePACA+++yzJpe//vrrAoAoLS0VQhhDytixY03aTJo0SfrxF0I0ub2mQkrfvn1FQ0ODNG/gwIHirrvukj7X19eLbt26iU2bNgkhrg8pv7dlyxahVCrFrl27mv2u+/btEwCkp602Ps31woULJu1+H1KOHj0qAIgffvhBWl5RUSFcXFzEJ598IoQwhhQA4vjx41KblJQU4ePj02wt999/v5g6dWqTy5r7nv379xcbN240mffKK6+IyMhIk/VWrFghLa+rqxO9e/cWf/3rX2+4X6LOjJd7iLowYcFDziMjI6/7fOjQIYv3eeutt8LO7rd/Wnx8fDB06FDps729PXr27CldomlOYWEhnnjiCbz11lu48847pfn5+fm4//770adPH7i5uWHUqFEAgOLiYrNrPHToEBwcHKDRaKR5PXv2xMCBA02+s6urK/r37y999vX1bbHuZ555BmlpaQgNDcWCBQuwe/fuFuuorq7GiRMnMG3aNHTv3l2ali1bhhMnTpi0/f2fj4ODA8LDw6VaLd0vUWfBkELUBQ0YMAAKhaLZkHHo0CF4enqiV69ebb5vR0dHk88KhaLJeQaDodlt6HQ6jB8/Hk899RSmTZsmza+uroZWq4W7uzs+/vhj7Nu3D5999hkAoLa2tg2/hVFTdbcU/O699178+uuveO6553D27FmMGTMGzz//fLPtL126BAB49913sX//fmk6ePAg9uzZY3adlu6XqLNgSCHqgnr27Il77rkHb7/9Nq5cuWKyTKfT4eOPP8akSZOgUCik+df+KO7ZsweDBg2SPjs6OqKhoaF9Cwdw9epVTJgwASEhIXj99ddNlh0+fBjnzp3DihUrcNdddyEkJOS6MxtOTk4A0GKtgwYNQn19PXJzc6V5586dw5EjRzB48OCbqr9Xr16Ij4/HRx99hDVr1mDdunXN1uXj4wM/Pz+cPHkSAwYMMJkaO9o2+v2fT319PfLz803+fJrbL1Fn1r73/BGRbN566y1ERUVBq9Vi2bJlCAoKws8//4z58+fD398fy5cvN2n/ww8/YOXKlYiNjUVWVhbS09PxxRdfSMsDAwORnZ2NO++8E87OzvD09GyXumfOnInTp08jOzsb5eXl0vwePXqgT58+cHJywptvvomnn34aBw8evG7sk759+0KhUGD79u0YN24cXFxc0L17d5M2wcHBmDBhAqZPn4533nkHbm5ueOGFF+Dv748JEya0uvakpCSEhYXh1ltvRU1NDbZv3y4FCW9vb7i4uCAzMxO9e/eGUqmESqXCSy+9hGeffRYqlQpjx45FTU0N8vLycOHCBSQkJEjbTklJQXBwMAYNGoQ33ngDFy5cwJNPPnnD/RJ1ZjyTQtRFBQcHIy8vD/369cPDDz+M/v37Y8aMGRg9ejRycnLQo0cPk/Z//vOfkZeXhxEjRmDZsmV4/fXXodVqpeWrV69GVlYWAgICMGLEiHar+7vvvkNJSQkGDx4MX19fadq9ezd69eqFDRs2ID09HYMHD8aKFSvw2muvmazv7++Pl156CS+88AJ8fHwwe/bsJvezfv16hIWF4b777kNkZCSEENixY8d1l3gs4eTkhMTERAwbNgx/+MMfYG9vj7S0NADGfiRr167FO++8Az8/PykMPfXUU3jvvfewfv16DB06FKNGjcKGDRuuO5OyYsUKrFixAsOHD8euXbuwbds2eHl53XC/RJ2ZQljSs46IuqTAwEDMmzePw8ZboVOnTiEoKAiFhYU3/QgCos6GZ1KIiIjIKjGkEBERkVXi5R4iIiKySjyTQkRERFaJIYWIiIisEkMKERERWSWGFCIiIrJKDClERERklRhSiIiIyCoxpBAREZFVYkghIiIiq8SQQkRERFbp/wFo8pjFfsMLowAAAABJRU5ErkJggg==\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(6, 4))\n",
"\n",
"# Enable processing the Torch trainable tensors\n",
"with torch.no_grad():\n",
" plt.plot(x, cost_pt_loc, 'r', label = 'local minimum')\n",
" plt.xlabel(\"Optimization steps\")\n",
" plt.ylabel(\"Cost / Energy\")\n",
" plt.legend()\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Send it after class |\n",
"\n",
"Try with different initialization parameters and see how the results change.\n",
"\n"
]
},
{
"cell_type": "markdown",
"source": [
"# Two Types of Quantum Computers\n",
"\n",
"It is not well-known what problems quantum computers can solve. We can classify quantum computers into “gate model” and “quantum annealing”. They have different use cases and architecture.\n",
"\n",
"## Gate Model\n",
"\n",
"The quantum gate model computers are designed for universal purposes by applying gates to build complex algorithms. They can be considered the upward compatible machines of classical computers but exponentially fast if appropriate quantum algorithms like Shor and Grover are used. The applications include prime factorisation, which becomes a threat to the existing secure data transmission methods over the internet using RSA (Rivest–Shamir–Adleman) keys, quantum simulations which have the potential for finding new medicines, and machine learning. However, the gate model requires high-quality qubit processors; the number of qubits is limited to two digits (except for IBM’s exploratory 127-qubit Eagle system.)\n",
"\n",
"## Quantum Annealing\n",
"\n",
"On the other hand, quantum annealers shine when they are used to search for the best from a vast number of possible combinations. Such challenges are called combinatorial optimisation problems. Many applications of quantum annealing have been reported recently. There are also researches to develop novel machine learning algorithms using quantum annealers. Amin et al. showed that there were possibilities to use quantum annealing hardware as a sampler for Boltzmann Machine by exploiting its quantum nature.\n",
"\n",
"In physics, the minimum energy principle states that the internal energy decreases and approaches the minimum values. If we can formulate our problems in the form of an energy minimisation equation, the quantum annealer will be able to search for the best possible answer.\n",
"\n",
"Quantum annealer is relatively noise-tolerant compared to the gate model. At the time of this writing, the most considerable number of qubits on a quantum processing unit (QPU) is at least five thousand provided by D-Wave’s Advantage solver. It is expected to increase more in the coming years. Five thousand qubits is a large number compared to what is achieved by the gate model. However, what we obtain from the annealer is as limited as a set of 5k-bit values per sampling. Furthermore, not all qubits are connected. Therefore, a hybrid system of classical computing and the quantum annealer is being developed to complement the capabilities of the two systems.\n",
"\n",
"Annealing is a technique used to make quality iron products by heating the material to a certain degree of temperature and cooling it down slowly. The iron’s mechanical characteristics are altered through the heat treatment, which reduces the hardness and increases the elasticity.\n",
"\n",
"Simulated annealing (SA) is a probabilistic technique whose name and idea are derived from annealing in material science. SA is a meta-heuristic in which algorithms approximate the global minimum or maximum of a given function. It uses a parameter called temperature that is initialised with a high value and decreases as the algorithm proceeds, like the physical annealing. High-temperature value promotes the algorithm to look for the various parts of the problem space, which helps to avoid being stuck in local minima or maxima. When the temperature goes down to zero, the state reaches the optimal point if successful. Going through the SA algorithm is outside the scope of this article, and the details on how it works can be found, for example, on the MathWorks website.\n",
"\n",
"Quantum annealing (QA) is an algorithm to solve combinatorial optimisation problems mainly. Instead of using temperature to explore the problem space, QA uses the laws of quantum mechanics to measure the energy state. Starting from the qubits’ dual state, where all the possible combinations of the values are equally likely to happen, we can gradually reduce such quantum mechanical effects by the process called quantum fluctuation to reach the most optimal point where the lowest energy state is realised.\n",
"\n",
"Kadowaki and Nishimori proved that QA converges to the optimal state with a much larger probability than SA in almost all cases on classical computers. The actual quantum annealing machines were developed by D-Wave and built on the ground of their theoretical framework.\n",
"\n",
"## Formulating Problem for QA\n",
"\n",
"QA machines are specialised hardware to solve combinatorial optimisation problems. These problems can be found in many places in our life. For example, an Amazon delivery person needs to find the most optimal route to deliver all parcels given for the day. The objective is to find a route with the shortest distance to visit all the delivery addresses. We can treat the distance as the energy in the QA context and find the path that reaches the ground state. We can also think of a futuristic example where multiple autonomous cars are navigated to optimise the traffic in a specific city area with heavy traffic jams. The objective is to reduce the traffic by minimising the number of overlaps of cars taking the same route simultaneously while each car aims to optimise the travel time to a destination.\n",
"\n",
"Quantum annealers require an objective function to be defined in the Ising model or Quadratic Unconstrained Binary Optimisation (QUBO). Ising model is a mathematical model in statistical mechanics. QUBO is mathematically equivalent to the Ising model and is much simpler to formulate a problem.\n",
"\n",
"QUBO is a minimisation objective function expressed as follows:\n",
"\n",
"$$ H_{QUBO}(x)=\\sum_{i=1}^N \\sum_{j=1}^N Q_{ij}x_ix_j $$\n",
"\n",
"\n",
"where $x_i$ and $x_j\\_ {i,j =0,1,…,N}$ takes 0 or 1. $Q_{ij}$ is a coefficient matrix of size $n \\times n$, which we need to prepare for a QA machine. In case of a problem to search for a global maximum, we can flip the sign to minus to convert the problem into a minimisation problem.\n",
"\n",
"In reality, a problem must be satisfied with one or more constraints. For example, a knapsack problem has a constraint of how many weights you can pack things in the bag. Such constraints are multiplied by weight parameters and are added to the objective function as penalty terms.\n",
"\n",
"## Implementation Using Quantum Annealer\n",
"\n",
"In this section, we would like to code and solve a simple combinatorial optimisation problem in Python and the real Quantum Annealer by D-Wave.\n",
"\n",
"## Preparation\n",
"\n",
"Before starting, we need to sign up for a cloud computing service. D-wave offers trial access to QA solvers with a sixty-second computational time allowance for one month. Please go to D-Wave Leap to find out more. After successful registration, the API token is available on the dashboard."
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "markdown",
"source": [
"## Problem Setting\n",
"\n",
"We would like to challenge a simplified portfolio optimisation problem. Suppose there are twenty candidate stocks from FTSE100 to invest in. Our objective is to maximise an expected return rate at the minimum risk while keeping more or less ten stocks in our portfolio.\n",
"\n",
"We will take the following six steps to solve the problem:\n",
"\n",
"1. Import historical stock price data from the internet\n",
"2. Prepare data by cleansing and feature engineering\n",
"3. Formulate the cost functions and constraints to obtain a QUBO matrix\n",
"4. Sample optimisation results from a D-Wave solver\n",
"5. Sample optimisation results from a simulated quantum annealing solver\n",
"6. Evaluate the outcomes\n",
"\n",
"Each of the above steps is demonstrated in a Colab notebook found in Python Coding in Colab Notebook section.\n",
"\n",
"## Dependencies\n",
"\n",
"Some notes on the Python requirements:\n",
"\n",
"* yfinance is used to retrieve stock price data from Yahoo! Finance.\n",
"* Prophet is a python package used to make time-series forecasting. It is possible to do multivariate time series forecasting and add holiday effects easily.\n",
"* D-Wave Ocean SDK needs to be installed. The SDK includes packages that are necessary to communicate with D-Wave resources.\n",
"* PyQUBO is a library that can construct QUBO from cost and constraints functions.\n",
"* OpenJIJ is an open-source meta-heuristics library for the Ising model and QUBO. It simulates the quantum annealers and can be used to experiment without the D-Wave computers. It accepts QUBOs stored in a Python dictionary type.\n",
"\n",
"## Cost and Constraint Function\n",
"\n",
"Our objective in this example is to minimise the sum of two cost functions and one constraint term, namely:\n",
"\n",
"1. the maximisation of return rates\n",
"2. the minimisation of risks by lowering portfolio variance\n",
"3. the minimisation of the penalty term\n",
"\n",
"1. The maximisation of the return rates can be formulated as a minimisation cost function:\n",
"\n",
"$$ H_{cost1}=-\\sum_{i=1}^N R_ix_i $$ (cost function 1)\n",
"\n",
"\n",
"where $R_i$ is the return rate of the $i$th stock $x_i$. Because $x_i$ only takes a binary value $0$ or $1$, $x_i$ is equivalent with the quadratic term $x_i \\times x_i$. Therefore, $H_{cost1}$ function can be equivalent to the following QUBO:\n",
"\n",
"$$ H_{cost1_{QUBO}} = -\\sum_{i=j\\, i,j=1}^NR_{ij}x_ix_j $$ (cost function as QUBO)\n",
"\n",
"\n",
"2. The minimisation of portfolio variance is represented as :\n",
"\n",
"$$ H_{cost2}=\\sum_{i=1}^N\\sum_{j=1}^NW_{ij}x_ix_j $$ (cost function 2)\n",
"\n",
"where W denotes a covariance matrix.\n",
"\n",
"3. The minimisation of the penalty can be formed in mean squared error (MSE) as a constraint term:\n",
"\n",
"$$ H_{const}=(\\sum_{i=1}^N x_i- K)^2 $$ (constraint term)\n",
"\n",
"\n",
"where K is the number of stocks to be selected for the portfolio.\n",
"\n",
"$H_{const}$ can also be transformed to the equivalent QUBOs by expanding the formula:\n",
"\n",
"$$ H_{const}=(\\sum_{i=1}^N x_i- K)^2\\\\ =(\\sum_{i=1}^N x_i- K)(\\sum_{j=1}^N x_j- K)\\\\=\\sum_{i,j=1}^N x_ix_j-2K\\sum_{i=1}^Nx_i + K^2 $$\n",
"(Constraint term expanded for QUBO formulation)\n",
"\n",
"$-2K\\sum_{i=1}^Nx_i$ part can also be equivalent to the quadratic form, as seen in the above $H_{cost1}$ function. $K^2$ is a constant value called offset which is not included in the QUBO matrix.\n",
"\n",
"We will use the `PyQUBO` library in Python code to generate the QUBO matrix from functions for simplicity.\n",
"\n",
"Thus, we can solve the problem by finding such x that minimises the function H(x) = (H_cost1 + H_cost2 + lambda*H_const), where lambda is a weight parameter for the constraint"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "aQbHUjMPD3GZ"
},
"source": [
"# 1. Import historical stock price data from the internet\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "slZHyufhIy_m"
},
"source": [
"First of all, Install Yahoo finance library to retrieve stock data.\n",
"We will import randomly selected 20 stocks from FTSE100 and extract the historical closing prices for the past 2 years from July 2020."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"id": "mnsygle_Dpw6"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collecting yfinance\r\n",
" Downloading yfinance-0.2.3-py2.py3-none-any.whl (50 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m50.4/50.4 kB\u001B[0m \u001B[31m674.3 kB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: requests>=2.26 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (2.28.1)\r\n",
"Collecting pytz>=2022.5\r\n",
" Downloading pytz-2022.7-py2.py3-none-any.whl (499 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m499.4/499.4 kB\u001B[0m \u001B[31m3.5 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hCollecting html5lib>=1.1\r\n",
" Downloading html5lib-1.1-py2.py3-none-any.whl (112 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m112.2/112.2 kB\u001B[0m \u001B[31m26.3 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: pandas>=1.3.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (1.5.1)\r\n",
"Requirement already satisfied: cryptography>=3.3.2 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (37.0.1)\r\n",
"Collecting lxml>=4.9.1\r\n",
" Downloading lxml-4.9.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl (7.1 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m7.1/7.1 MB\u001B[0m \u001B[31m29.2 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0ma \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: beautifulsoup4>=4.11.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (4.11.1)\r\n",
"Requirement already satisfied: frozendict>=2.3.4 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (2.3.4)\r\n",
"Collecting multitasking>=0.0.7\r\n",
" Downloading multitasking-0.0.11-py3-none-any.whl (8.5 kB)\r\n",
"Requirement already satisfied: appdirs>=1.4.4 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (1.4.4)\r\n",
"Requirement already satisfied: numpy>=1.16.5 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from yfinance) (1.23.4)\r\n",
"Requirement already satisfied: soupsieve>1.2 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from beautifulsoup4>=4.11.1->yfinance) (2.3.1)\r\n",
"Requirement already satisfied: cffi>=1.12 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from cryptography>=3.3.2->yfinance) (1.15.1)\r\n",
"Requirement already satisfied: webencodings in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from html5lib>=1.1->yfinance) (0.5.1)\r\n",
"Requirement already satisfied: six>=1.9 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from html5lib>=1.1->yfinance) (1.16.0)\r\n",
"Requirement already satisfied: python-dateutil>=2.8.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pandas>=1.3.0->yfinance) (2.8.2)\r\n",
"Requirement already satisfied: urllib3<1.27,>=1.21.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests>=2.26->yfinance) (1.26.12)\r\n",
"Requirement already satisfied: certifi>=2017.4.17 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests>=2.26->yfinance) (2022.9.24)\r\n",
"Requirement already satisfied: idna<4,>=2.5 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests>=2.26->yfinance) (3.4)\r\n",
"Requirement already satisfied: charset-normalizer<3,>=2 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests>=2.26->yfinance) (2.1.1)\r\n",
"Requirement already satisfied: pycparser in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from cffi>=1.12->cryptography>=3.3.2->yfinance) (2.21)\r\n",
"Installing collected packages: pytz, multitasking, lxml, html5lib, yfinance\r\n",
" Attempting uninstall: pytz\r\n",
" Found existing installation: pytz 2022.1\r\n",
" Uninstalling pytz-2022.1:\r\n",
" Successfully uninstalled pytz-2022.1\r\n",
"Successfully installed html5lib-1.1 lxml-4.9.2 multitasking-0.0.11 pytz-2022.7 yfinance-0.2.3\r\n",
"\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m A new release of pip available: \u001B[0m\u001B[31;49m22.3\u001B[0m\u001B[39;49m -> \u001B[0m\u001B[32;49m22.3.1\u001B[0m\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m To update, run: \u001B[0m\u001B[32;49mpip install --upgrade pip\u001B[0m\r\n"
]
}
],
"source": [
"!pip install yfinance"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"id": "VXCVsQj_uTcO"
},
"outputs": [],
"source": [
"import yfinance as yf \n",
"import numpy as np\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "B8fQ4MloD_zP",
"outputId": "20fffb52-6657-48f5-d747-d45588eac329"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[*********************100%***********************] 20 of 20 completed\n",
"[*********************100%***********************] 20 of 20 completed\n"
]
}
],
"source": [
"# list of FTSE100 stocks available at \n",
"# https://www.londonstockexchange.com/indices/ftse-100/constituents/table\n",
"# below are randomly selected stocks\n",
"\n",
"stocks_dict = {\n",
" 'FLTR.L':'FLUTTER ENTERTAINMENT PLC',\n",
" 'SVT.L':'SEVERN TRENT PLC', \n",
" 'LSEG.L':'LONDON STOCK EXCHANGE GROUP PLC', \n",
" 'BNZL.L':'BUNZL PLC', \n",
" 'CPG.L':'COMPASS GROUP PLC', \n",
" 'NXT.L':'NEXT PLC', \n",
" 'ULVR.L':'UNILEVER PLC', \n",
" 'SHEL.L':'SHELL PLC', \n",
" 'CCH.L':'COCA-COLA HBC A', \n",
" 'BRBY.L':'BURBERRY GROUP PLC', \n",
" 'SSE.L':'SSE PLC', \n",
" 'REL.L':'RELX PLC', \n",
" 'EXPN.L':'EXPERIAN PLC', \n",
" 'DCC.L':'DCC PLC ORD', \n",
" 'AZN.L':'ASTRAZENECA PLC', \n",
" 'GSK.L':'GSK PLC ORD', \n",
" 'ADM.L':'ADMIRAL GROUP PLC', \n",
" 'SDR.L':'SCHRODERS PLC', \n",
" 'BATS.L':'BRITISH AMERICAN TOBACCO PLC', \n",
" 'IHG.L':'INTERCONTINENTAL HOTELS GROUP PLC', \n",
" }\n",
" \n",
"stocks = list(stocks_dict.keys()) # list of candidate stocks\n",
"\n",
"data = yf.download(stocks,'2020-07-01','2022-04-30')['Close'] # main data\n",
"data_test = yf.download(stocks,'2022-05-01','2022-07-24')['Close'] # testing data"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 364
},
"id": "kwMT3Yb_X1V5",
"outputId": "3bd6d8ef-b5ec-4975-cb3b-3c96fe2e91d0"
},
"outputs": [
{
"data": {
"text/plain": " ADM.L AZN.L BATS.L BNZL.L BRBY.L \\\ncount 464.000000 464.000000 463.000000 464.000000 464.000000 \nmean 2949.351293 8366.258621 2787.561555 2537.961207 1804.203375 \nstd 298.590229 753.326915 223.859266 228.969271 228.739954 \nmin 2240.000000 6794.000000 2448.000000 2131.000000 1252.500000 \n25% 2756.750000 7931.500000 2643.750000 2363.000000 1630.000000 \n50% 2962.000000 8413.000000 2734.000000 2479.500000 1820.250000 \n75% 3116.250000 8678.250000 2825.750000 2712.750000 1965.250000 \nmax 3688.000000 10930.000000 3439.000000 3163.000000 2264.000000 \n\n CCH.L CPG.L DCC.L EXPN.L FLTR.L \\\ncount 464.000000 463.000000 463.000000 463.000000 464.000000 \nmean 2294.758621 1466.997840 6042.596112 2933.794816 12754.031149 \nstd 314.297710 176.851449 420.743959 288.931273 2003.651264 \nmin 1460.500000 1050.500000 5006.000000 2273.000000 7922.000000 \n25% 2073.000000 1393.250000 5785.000000 2738.000000 11400.000000 \n50% 2372.500000 1497.500000 6030.000000 2893.000000 12890.000000 \n75% 2532.250000 1586.250000 6268.000000 3109.000000 14282.500000 \nmax 2784.000000 1820.500000 7204.000000 3667.000000 16915.000000 \n\n GSK.L IHG.L LSEG.L NXT.L REL.L \\\ncount 464.000000 463.000000 464.000000 464.000000 464.000000 \nmean 1475.618385 4722.868251 7975.392241 7204.778017 1997.713362 \nstd 128.447452 399.235132 794.814741 952.486937 248.250476 \nmin 1199.561523 3516.000000 6370.000000 4673.000000 1527.500000 \n25% 1385.821991 4549.000000 7371.500000 6282.000000 1788.625000 \n50% 1446.565674 4812.000000 7939.000000 7660.000000 1918.250000 \n75% 1579.234680 5026.000000 8596.000000 7980.500000 2227.250000 \nmax 1823.720459 5338.000000 9910.000000 8426.000000 2449.000000 \n\n SDR.L SHEL.L SSE.L SVT.L ULVR.L \ncount 464.000000 463.000000 464.000000 464.000000 464.000000 \nmean 568.022027 1489.860042 1503.168103 2601.586207 4142.657884 \nstd 52.854523 299.833511 144.853407 241.375319 357.982299 \nmin 442.510010 900.000000 1171.000000 2168.000000 3328.000000 \n25% 519.009995 1333.299988 1384.875000 2413.250000 3911.125000 \n50% 587.179993 1440.400024 1524.250000 2511.000000 4124.500000 \n75% 606.942520 1660.300049 1612.000000 2823.000000 4358.500000 \nmax 658.070007 2228.000000 1868.500000 3211.000000 4892.000000 ",
"text/html": "<div>\n<style scoped>\n .dataframe tbody tr th:only-of-type {\n vertical-align: middle;\n }\n\n .dataframe tbody tr th {\n vertical-align: top;\n }\n\n .dataframe thead th {\n text-align: right;\n }\n</style>\n<table border=\"1\" class=\"dataframe\">\n <thead>\n <tr style=\"text-align: right;\">\n <th></th>\n <th>ADM.L</th>\n <th>AZN.L</th>\n <th>BATS.L</th>\n <th>BNZL.L</th>\n <th>BRBY.L</th>\n <th>CCH.L</th>\n <th>CPG.L</th>\n <th>DCC.L</th>\n <th>EXPN.L</th>\n <th>FLTR.L</th>\n <th>GSK.L</th>\n <th>IHG.L</th>\n <th>LSEG.L</th>\n <th>NXT.L</th>\n <th>REL.L</th>\n <th>SDR.L</th>\n <th>SHEL.L</th>\n <th>SSE.L</th>\n <th>SVT.L</th>\n <th>ULVR.L</th>\n </tr>\n </thead>\n <tbody>\n <tr>\n <th>count</th>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>463.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>463.000000</td>\n <td>463.000000</td>\n <td>463.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>463.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>463.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n <td>464.000000</td>\n </tr>\n <tr>\n <th>mean</th>\n <td>2949.351293</td>\n <td>8366.258621</td>\n <td>2787.561555</td>\n <td>2537.961207</td>\n <td>1804.203375</td>\n <td>2294.758621</td>\n <td>1466.997840</td>\n <td>6042.596112</td>\n <td>2933.794816</td>\n <td>12754.031149</td>\n <td>1475.618385</td>\n <td>4722.868251</td>\n <td>7975.392241</td>\n <td>7204.778017</td>\n <td>1997.713362</td>\n <td>568.022027</td>\n <td>1489.860042</td>\n <td>1503.168103</td>\n <td>2601.586207</td>\n <td>4142.657884</td>\n </tr>\n <tr>\n <th>std</th>\n <td>298.590229</td>\n <td>753.326915</td>\n <td>223.859266</td>\n <td>228.969271</td>\n <td>228.739954</td>\n <td>314.297710</td>\n <td>176.851449</td>\n <td>420.743959</td>\n <td>288.931273</td>\n <td>2003.651264</td>\n <td>128.447452</td>\n <td>399.235132</td>\n <td>794.814741</td>\n <td>952.486937</td>\n <td>248.250476</td>\n <td>52.854523</td>\n <td>299.833511</td>\n <td>144.853407</td>\n <td>241.375319</td>\n <td>357.982299</td>\n </tr>\n <tr>\n <th>min</th>\n <td>2240.000000</td>\n <td>6794.000000</td>\n <td>2448.000000</td>\n <td>2131.000000</td>\n <td>1252.500000</td>\n <td>1460.500000</td>\n <td>1050.500000</td>\n <td>5006.000000</td>\n <td>2273.000000</td>\n <td>7922.000000</td>\n <td>1199.561523</td>\n <td>3516.000000</td>\n <td>6370.000000</td>\n <td>4673.000000</td>\n <td>1527.500000</td>\n <td>442.510010</td>\n <td>900.000000</td>\n <td>1171.000000</td>\n <td>2168.000000</td>\n <td>3328.000000</td>\n </tr>\n <tr>\n <th>25%</th>\n <td>2756.750000</td>\n <td>7931.500000</td>\n <td>2643.750000</td>\n <td>2363.000000</td>\n <td>1630.000000</td>\n <td>2073.000000</td>\n <td>1393.250000</td>\n <td>5785.000000</td>\n <td>2738.000000</td>\n <td>11400.000000</td>\n <td>1385.821991</td>\n <td>4549.000000</td>\n <td>7371.500000</td>\n <td>6282.000000</td>\n <td>1788.625000</td>\n <td>519.009995</td>\n <td>1333.299988</td>\n <td>1384.875000</td>\n <td>2413.250000</td>\n <td>3911.125000</td>\n </tr>\n <tr>\n <th>50%</th>\n <td>2962.000000</td>\n <td>8413.000000</td>\n <td>2734.000000</td>\n <td>2479.500000</td>\n <td>1820.250000</td>\n <td>2372.500000</td>\n <td>1497.500000</td>\n <td>6030.000000</td>\n <td>2893.000000</td>\n <td>12890.000000</td>\n <td>1446.565674</td>\n <td>4812.000000</td>\n <td>7939.000000</td>\n <td>7660.000000</td>\n <td>1918.250000</td>\n <td>587.179993</td>\n <td>1440.400024</td>\n <td>1524.250000</td>\n <td>2511.000000</td>\n <td>4124.500000</td>\n </tr>\n <tr>\n <th>75%</th>\n <td>3116.250000</td>\n <td>8678.250000</td>\n <td>2825.750000</td>\n <td>2712.750000</td>\n <td>1965.250000</td>\n <td>2532.250000</td>\n <td>1586.250000</td>\n <td>6268.000000</td>\n <td>3109.000000</td>\n <td>14282.500000</td>\n <td>1579.234680</td>\n <td>5026.000000</td>\n <td>8596.000000</td>\n <td>7980.500000</td>\n <td>2227.250000</td>\n <td>606.942520</td>\n <td>1660.300049</td>\n <td>1612.000000</td>\n <td>2823.000000</td>\n <td>4358.500000</td>\n </tr>\n <tr>\n <th>max</th>\n <td>3688.000000</td>\n <td>10930.000000</td>\n <td>3439.000000</td>\n <td>3163.000000</td>\n <td>2264.000000</td>\n <td>2784.000000</td>\n <td>1820.500000</td>\n <td>7204.000000</td>\n <td>3667.000000</td>\n <td>16915.000000</td>\n <td>1823.720459</td>\n <td>5338.000000</td>\n <td>9910.000000</td>\n <td>8426.000000</td>\n <td>2449.000000</td>\n <td>658.070007</td>\n <td>2228.000000</td>\n <td>1868.500000</td>\n <td>3211.000000</td>\n <td>4892.000000</td>\n </tr>\n </tbody>\n</table>\n</div>"
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data.describe()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 364
},
"id": "qP_ey5ncNfmv",
"outputId": "15310f68-22a7-48ed-e0f7-c56d0ef7daf2"
},
"outputs": [
{
"data": {
"text/plain": " ADM.L AZN.L BATS.L BNZL.L BRBY.L \\\ncount 57.000000 57.000000 57.000000 57.000000 57.000000 \nmean 2185.973684 10508.228070 3466.447368 2814.035088 1629.964912 \nstd 173.503948 403.405858 81.319788 146.057225 60.288581 \nmin 1729.000000 9750.000000 3300.000000 2575.000000 1482.000000 \n25% 2157.000000 10244.000000 3426.000000 2671.000000 1586.000000 \n50% 2228.000000 10506.000000 3478.000000 2837.000000 1636.500000 \n75% 2267.000000 10800.000000 3528.500000 2919.000000 1678.500000 \nmax 2561.000000 11232.000000 3628.000000 3108.000000 1743.000000 \n\n CCH.L CPG.L DCC.L EXPN.L FLTR.L \\\ncount 57.000000 57.000000 57.000000 57.000000 57.000000 \nmean 1777.026316 1710.345614 5494.614035 2555.052632 8403.048070 \nstd 101.134432 237.407090 409.703802 152.041397 1243.157730 \nmin 1548.000000 16.700001 4889.000000 2285.000000 81.739998 \n25% 1712.000000 1691.500000 5162.000000 2424.000000 8200.000000 \n50% 1775.500000 1730.000000 5308.000000 2578.000000 8468.000000 \n75% 1837.500000 1796.500000 5710.000000 2661.000000 8944.000000 \nmax 1945.500000 1854.000000 6268.000000 2835.000000 9760.000000 \n\n GSK.L IHG.L LSEG.L NXT.L REL.L \\\ncount 57.000000 57.000000 57.000000 57.000000 57.000000 \nmean 1758.225729 4650.421053 7403.684211 6209.087719 2244.596491 \nstd 34.193845 256.656103 264.537545 239.741846 75.721260 \nmin 1684.302124 4193.000000 6738.000000 5764.000000 2081.000000 \n25% 1732.252319 4404.000000 7220.000000 5986.000000 2210.000000 \n50% 1758.242188 4695.000000 7390.000000 6148.000000 2265.000000 \n75% 1783.400024 4841.000000 7632.000000 6428.000000 2299.000000 \nmax 1815.863037 5150.000000 7894.000000 6686.000000 2375.000000 \n\n SDR.L SHEL.L SSE.L SVT.L ULVR.L \ncount 57.000000 57.000000 57.000000 57.000000 57.000000 \nmean 471.872281 2205.994736 1753.149123 2853.773158 3699.675439 \nstd 16.497939 148.555222 89.980502 401.790333 129.866329 \nmin 438.260010 1936.400024 1587.500000 28.070000 3451.500000 \n25% 458.320007 2044.000000 1677.500000 2794.000000 3609.500000 \n50% 469.200012 2225.000000 1761.000000 2872.000000 3690.500000 \n75% 486.200012 2338.000000 1821.000000 2996.000000 3802.500000 \nmax 505.239990 2440.000000 1920.000000 3148.000000 3943.000000 ",
"text/html": "<div>\n<style scoped>\n .dataframe tbody tr th:only-of-type {\n vertical-align: middle;\n }\n\n .dataframe tbody tr th {\n vertical-align: top;\n }\n\n .dataframe thead th {\n text-align: right;\n }\n</style>\n<table border=\"1\" class=\"dataframe\">\n <thead>\n <tr style=\"text-align: right;\">\n <th></th>\n <th>ADM.L</th>\n <th>AZN.L</th>\n <th>BATS.L</th>\n <th>BNZL.L</th>\n <th>BRBY.L</th>\n <th>CCH.L</th>\n <th>CPG.L</th>\n <th>DCC.L</th>\n <th>EXPN.L</th>\n <th>FLTR.L</th>\n <th>GSK.L</th>\n <th>IHG.L</th>\n <th>LSEG.L</th>\n <th>NXT.L</th>\n <th>REL.L</th>\n <th>SDR.L</th>\n <th>SHEL.L</th>\n <th>SSE.L</th>\n <th>SVT.L</th>\n <th>ULVR.L</th>\n </tr>\n </thead>\n <tbody>\n <tr>\n <th>count</th>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n <td>57.000000</td>\n </tr>\n <tr>\n <th>mean</th>\n <td>2185.973684</td>\n <td>10508.228070</td>\n <td>3466.447368</td>\n <td>2814.035088</td>\n <td>1629.964912</td>\n <td>1777.026316</td>\n <td>1710.345614</td>\n <td>5494.614035</td>\n <td>2555.052632</td>\n <td>8403.048070</td>\n <td>1758.225729</td>\n <td>4650.421053</td>\n <td>7403.684211</td>\n <td>6209.087719</td>\n <td>2244.596491</td>\n <td>471.872281</td>\n <td>2205.994736</td>\n <td>1753.149123</td>\n <td>2853.773158</td>\n <td>3699.675439</td>\n </tr>\n <tr>\n <th>std</th>\n <td>173.503948</td>\n <td>403.405858</td>\n <td>81.319788</td>\n <td>146.057225</td>\n <td>60.288581</td>\n <td>101.134432</td>\n <td>237.407090</td>\n <td>409.703802</td>\n <td>152.041397</td>\n <td>1243.157730</td>\n <td>34.193845</td>\n <td>256.656103</td>\n <td>264.537545</td>\n <td>239.741846</td>\n <td>75.721260</td>\n <td>16.497939</td>\n <td>148.555222</td>\n <td>89.980502</td>\n <td>401.790333</td>\n <td>129.866329</td>\n </tr>\n <tr>\n <th>min</th>\n <td>1729.000000</td>\n <td>9750.000000</td>\n <td>3300.000000</td>\n <td>2575.000000</td>\n <td>1482.000000</td>\n <td>1548.000000</td>\n <td>16.700001</td>\n <td>4889.000000</td>\n <td>2285.000000</td>\n <td>81.739998</td>\n <td>1684.302124</td>\n <td>4193.000000</td>\n <td>6738.000000</td>\n <td>5764.000000</td>\n <td>2081.000000</td>\n <td>438.260010</td>\n <td>1936.400024</td>\n <td>1587.500000</td>\n <td>28.070000</td>\n <td>3451.500000</td>\n </tr>\n <tr>\n <th>25%</th>\n <td>2157.000000</td>\n <td>10244.000000</td>\n <td>3426.000000</td>\n <td>2671.000000</td>\n <td>1586.000000</td>\n <td>1712.000000</td>\n <td>1691.500000</td>\n <td>5162.000000</td>\n <td>2424.000000</td>\n <td>8200.000000</td>\n <td>1732.252319</td>\n <td>4404.000000</td>\n <td>7220.000000</td>\n <td>5986.000000</td>\n <td>2210.000000</td>\n <td>458.320007</td>\n <td>2044.000000</td>\n <td>1677.500000</td>\n <td>2794.000000</td>\n <td>3609.500000</td>\n </tr>\n <tr>\n <th>50%</th>\n <td>2228.000000</td>\n <td>10506.000000</td>\n <td>3478.000000</td>\n <td>2837.000000</td>\n <td>1636.500000</td>\n <td>1775.500000</td>\n <td>1730.000000</td>\n <td>5308.000000</td>\n <td>2578.000000</td>\n <td>8468.000000</td>\n <td>1758.242188</td>\n <td>4695.000000</td>\n <td>7390.000000</td>\n <td>6148.000000</td>\n <td>2265.000000</td>\n <td>469.200012</td>\n <td>2225.000000</td>\n <td>1761.000000</td>\n <td>2872.000000</td>\n <td>3690.500000</td>\n </tr>\n <tr>\n <th>75%</th>\n <td>2267.000000</td>\n <td>10800.000000</td>\n <td>3528.500000</td>\n <td>2919.000000</td>\n <td>1678.500000</td>\n <td>1837.500000</td>\n <td>1796.500000</td>\n <td>5710.000000</td>\n <td>2661.000000</td>\n <td>8944.000000</td>\n <td>1783.400024</td>\n <td>4841.000000</td>\n <td>7632.000000</td>\n <td>6428.000000</td>\n <td>2299.000000</td>\n <td>486.200012</td>\n <td>2338.000000</td>\n <td>1821.000000</td>\n <td>2996.000000</td>\n <td>3802.500000</td>\n </tr>\n <tr>\n <th>max</th>\n <td>2561.000000</td>\n <td>11232.000000</td>\n <td>3628.000000</td>\n <td>3108.000000</td>\n <td>1743.000000</td>\n <td>1945.500000</td>\n <td>1854.000000</td>\n <td>6268.000000</td>\n <td>2835.000000</td>\n <td>9760.000000</td>\n <td>1815.863037</td>\n <td>5150.000000</td>\n <td>7894.000000</td>\n <td>6686.000000</td>\n <td>2375.000000</td>\n <td>505.239990</td>\n <td>2440.000000</td>\n <td>1920.000000</td>\n <td>3148.000000</td>\n <td>3943.000000</td>\n </tr>\n </tbody>\n</table>\n</div>"
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_test.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "W2Wy0fxgsPrz"
},
"source": [
"# 2. Prepare data by cleansing and feature engineering\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "VIbKQ2AzvIAl"
},
"source": [
"Check null values. We will drop a row with null in this example for simplicity but please cosnider replacing the values if possible."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Q7UlzUv7ZpHS",
"outputId": "7d52ff9e-99e9-4458-bebb-1227b8801db6"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ADM.L 0\n",
"AZN.L 0\n",
"BATS.L 1\n",
"BNZL.L 0\n",
"BRBY.L 0\n",
"CCH.L 0\n",
"CPG.L 1\n",
"DCC.L 1\n",
"EXPN.L 1\n",
"FLTR.L 0\n",
"GSK.L 0\n",
"IHG.L 1\n",
"LSEG.L 0\n",
"NXT.L 0\n",
"REL.L 0\n",
"SDR.L 0\n",
"SHEL.L 1\n",
"SSE.L 0\n",
"SVT.L 0\n",
"ULVR.L 0\n",
"dtype: int64\n"
]
}
],
"source": [
"print(data.isna().sum()) # null row for each stock"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "kF5eB8_LN2Hz",
"outputId": "aeb3f233-7292-4132-c4a5-3abc62dd112a"
},
"outputs": [
{
"data": {
"text/plain": "ADM.L 0\nAZN.L 0\nBATS.L 0\nBNZL.L 0\nBRBY.L 0\nCCH.L 0\nCPG.L 0\nDCC.L 0\nEXPN.L 0\nFLTR.L 0\nGSK.L 0\nIHG.L 0\nLSEG.L 0\nNXT.L 0\nREL.L 0\nSDR.L 0\nSHEL.L 0\nSSE.L 0\nSVT.L 0\nULVR.L 0\ndtype: int64"
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_test.isna().sum() # check for test data"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "TrqWiGbzcmuC",
"outputId": "0d5d9a6a-7c93-4857-87a5-89da4e033616"
},
"outputs": [
{
"data": {
"text/plain": "(459, 20)"
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data.dropna(inplace=True) # drop the rows\n",
"data.shape"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "PFC3Fv6-c1ii",
"outputId": "88282590-12a5-4cf7-deba-f4367ee770fd"
},
"outputs": [
{
"data": {
"text/plain": "(57, 20)"
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_test.shape"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "XYuopxu2M16X"
},
"source": [
"## Feature Engineering"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "7GSqM12zuK8s"
},
"source": [
"To make prediction, we use Prophet. See more at https://facebook.github.io/prophet/\n",
"\n",
"We make time series prediction based on the historical data between July 2020 and April 2022. The period to predict is May 2022 till late July 2022, which is the same period of the data in test set.\n",
"\n",
"The predicted values are transformed to be the average return rate for each stock. We can use it for one of the objectives to maximise the sum of return rates."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"id": "p0bAaw8huJxJ"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collecting prophet\r\n",
" Downloading prophet-1.1.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.9 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m8.9/8.9 MB\u001B[0m \u001B[31m18.5 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hCollecting convertdate>=2.1.2\r\n",
" Downloading convertdate-2.4.0-py3-none-any.whl (47 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m47.9/47.9 kB\u001B[0m \u001B[31m13.5 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: pandas>=1.0.4 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (1.5.1)\r\n",
"Collecting holidays>=0.14.2\r\n",
" Downloading holidays-0.17.2-py3-none-any.whl (187 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m187.5/187.5 kB\u001B[0m \u001B[31m42.2 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting tqdm>=4.36.1\r\n",
" Downloading tqdm-4.64.1-py2.py3-none-any.whl (78 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m78.5/78.5 kB\u001B[0m \u001B[31m6.9 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: numpy>=1.15.4 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (1.23.4)\r\n",
"Requirement already satisfied: matplotlib>=2.0.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (3.6.0)\r\n",
"Collecting LunarCalendar>=0.0.9\r\n",
" Downloading LunarCalendar-0.0.9-py2.py3-none-any.whl (18 kB)\r\n",
"Requirement already satisfied: setuptools>=42 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (63.4.1)\r\n",
"Collecting setuptools-git>=1.2\r\n",
" Downloading setuptools_git-1.2-py2.py3-none-any.whl (10 kB)\r\n",
"Requirement already satisfied: wheel>=0.37.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (0.37.1)\r\n",
"Requirement already satisfied: python-dateutil>=2.8.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from prophet) (2.8.2)\r\n",
"Collecting cmdstanpy>=1.0.4\r\n",
" Downloading cmdstanpy-1.0.8-py3-none-any.whl (81 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m81.2/81.2 kB\u001B[0m \u001B[31m23.5 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting pymeeus<=1,>=0.3.13\r\n",
" Downloading PyMeeus-0.5.12.tar.gz (5.8 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m5.8/5.8 MB\u001B[0m \u001B[31m25.7 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25h Preparing metadata (setup.py) ... \u001B[?25ldone\r\n",
"\u001B[?25hCollecting korean-lunar-calendar\r\n",
" Downloading korean_lunar_calendar-0.3.1-py3-none-any.whl (9.0 kB)\r\n",
"Collecting hijri-converter\r\n",
" Downloading hijri_converter-2.2.4-py3-none-any.whl (14 kB)\r\n",
"Requirement already satisfied: pytz in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from LunarCalendar>=0.0.9->prophet) (2022.7)\r\n",
"Collecting ephem>=3.7.5.3\r\n",
" Downloading ephem-4.1.4-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (1.8 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m1.8/1.8 MB\u001B[0m \u001B[31m36.6 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: contourpy>=1.0.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (1.0.5)\r\n",
"Requirement already satisfied: pyparsing>=2.2.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (3.0.9)\r\n",
"Requirement already satisfied: pillow>=6.2.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (9.2.0)\r\n",
"Requirement already satisfied: cycler>=0.10 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (0.11.0)\r\n",
"Requirement already satisfied: kiwisolver>=1.0.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (1.4.4)\r\n",
"Requirement already satisfied: packaging>=20.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (21.3)\r\n",
"Requirement already satisfied: fonttools>=4.22.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from matplotlib>=2.0.0->prophet) (4.37.4)\r\n",
"Requirement already satisfied: six>=1.5 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from python-dateutil>=2.8.0->prophet) (1.16.0)\r\n",
"Building wheels for collected packages: pymeeus\r\n",
" Building wheel for pymeeus (setup.py) ... \u001B[?25ldone\r\n",
"\u001B[?25h Created wheel for pymeeus: filename=PyMeeus-0.5.12-py3-none-any.whl size=732001 sha256=d1adb90fefb904067f35650132f35c1462d0335672a32dac671aa903f3d5a864\r\n",
" Stored in directory: /home/obm/.cache/pip/wheels/c2/3a/3d/11734e652782d3f823a08aae1c452e887eb16349750cca3f8a\r\n",
"Successfully built pymeeus\r\n",
"Installing collected packages: setuptools-git, pymeeus, korean-lunar-calendar, ephem, tqdm, hijri-converter, convertdate, LunarCalendar, holidays, cmdstanpy, prophet\r\n",
"Successfully installed LunarCalendar-0.0.9 cmdstanpy-1.0.8 convertdate-2.4.0 ephem-4.1.4 hijri-converter-2.2.4 holidays-0.17.2 korean-lunar-calendar-0.3.1 prophet-1.1.1 pymeeus-0.5.12 setuptools-git-1.2 tqdm-4.64.1\r\n",
"\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m A new release of pip available: \u001B[0m\u001B[31;49m22.3\u001B[0m\u001B[39;49m -> \u001B[0m\u001B[32;49m22.3.1\u001B[0m\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m To update, run: \u001B[0m\u001B[32;49mpip install --upgrade pip\u001B[0m\r\n"
]
}
],
"source": [
"!pip install prophet"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"id": "9kmSFx0mzU4x"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR:prophet.plot:Importing plotly failed. Interactive plots will not work.\n"
]
}
],
"source": [
"from prophet import Prophet"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"id": "rgmtGCvLzb7B"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/04vgcbiz.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/gc74709l.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=76984', 'data', 'file=/tmp/tmpygu7z4jl/04vgcbiz.json', 'init=/tmp/tmpygu7z4jl/gc74709l.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model6kmbl5_0/prophet_model-20221222121305.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:05 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:05 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/egnks3q5.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/zc53ldjl.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=4851', 'data', 'file=/tmp/tmpygu7z4jl/egnks3q5.json', 'init=/tmp/tmpygu7z4jl/zc53ldjl.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelt3ce1t5_/prophet_model-20221222121305.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:05 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:05 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/4jh1cdyq.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/2ophhnw5.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=74384', 'data', 'file=/tmp/tmpygu7z4jl/4jh1cdyq.json', 'init=/tmp/tmpygu7z4jl/2ophhnw5.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelrndl262l/prophet_model-20221222121305.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:05 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:06 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/ujh1q3dy.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/q6mmlqra.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=44191', 'data', 'file=/tmp/tmpygu7z4jl/ujh1q3dy.json', 'init=/tmp/tmpygu7z4jl/q6mmlqra.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelp3ozgbn3/prophet_model-20221222121306.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:06 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:06 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/v5dtfk74.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/gimcyzzu.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=37452', 'data', 'file=/tmp/tmpygu7z4jl/v5dtfk74.json', 'init=/tmp/tmpygu7z4jl/gimcyzzu.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model061sru16/prophet_model-20221222121306.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:06 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:06 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/46qo06w9.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/ydyysnd7.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=34946', 'data', 'file=/tmp/tmpygu7z4jl/46qo06w9.json', 'init=/tmp/tmpygu7z4jl/ydyysnd7.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelrb2ph4h_/prophet_model-20221222121307.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/rc02zqjt.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/uf68zjd1.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=82036', 'data', 'file=/tmp/tmpygu7z4jl/rc02zqjt.json', 'init=/tmp/tmpygu7z4jl/uf68zjd1.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model1kwl2y2p/prophet_model-20221222121307.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/5arfh7nt.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/8pu6qdvw.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=99799', 'data', 'file=/tmp/tmpygu7z4jl/5arfh7nt.json', 'init=/tmp/tmpygu7z4jl/8pu6qdvw.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model_pyb1lvj/prophet_model-20221222121307.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:07 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/yp65v73z.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/24pzm_qc.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=62171', 'data', 'file=/tmp/tmpygu7z4jl/yp65v73z.json', 'init=/tmp/tmpygu7z4jl/24pzm_qc.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelqwg728oa/prophet_model-20221222121308.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:08 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:08 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/cggclqmk.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/mueyhxbw.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=32988', 'data', 'file=/tmp/tmpygu7z4jl/cggclqmk.json', 'init=/tmp/tmpygu7z4jl/mueyhxbw.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelf_6cde_4/prophet_model-20221222121308.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:08 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:08 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/uznmh7uf.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/v45lk1z3.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=64220', 'data', 'file=/tmp/tmpygu7z4jl/uznmh7uf.json', 'init=/tmp/tmpygu7z4jl/v45lk1z3.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model2dasxdge/prophet_model-20221222121309.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/ljfxbpjs.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/qncou6o9.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=82641', 'data', 'file=/tmp/tmpygu7z4jl/ljfxbpjs.json', 'init=/tmp/tmpygu7z4jl/qncou6o9.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelm45vqm5w/prophet_model-20221222121309.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/t34tly9d.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/rlzpak_k.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=73431', 'data', 'file=/tmp/tmpygu7z4jl/t34tly9d.json', 'init=/tmp/tmpygu7z4jl/rlzpak_k.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model8_r60rxn/prophet_model-20221222121309.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:09 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/t_uo3pvq.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/7_uxupaw.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=11204', 'data', 'file=/tmp/tmpygu7z4jl/t_uo3pvq.json', 'init=/tmp/tmpygu7z4jl/7_uxupaw.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model73zaxpd9/prophet_model-20221222121310.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/k1bop_0p.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/rrmosmdb.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=7755', 'data', 'file=/tmp/tmpygu7z4jl/k1bop_0p.json', 'init=/tmp/tmpygu7z4jl/rrmosmdb.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelkcqeechv/prophet_model-20221222121310.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/8hxrifmp.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/6noe7zja.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=10513', 'data', 'file=/tmp/tmpygu7z4jl/8hxrifmp.json', 'init=/tmp/tmpygu7z4jl/6noe7zja.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_model9wovgma5/prophet_model-20221222121310.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:10 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/r66qiiav.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/mvhy0iqc.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=2844', 'data', 'file=/tmp/tmpygu7z4jl/r66qiiav.json', 'init=/tmp/tmpygu7z4jl/mvhy0iqc.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modeleae4auck/prophet_model-20221222121311.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:11 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:11 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/st8l6xs8.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/6a4su4kl.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=90262', 'data', 'file=/tmp/tmpygu7z4jl/st8l6xs8.json', 'init=/tmp/tmpygu7z4jl/6a4su4kl.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelj8g6b8k_/prophet_model-20221222121311.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:11 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:11 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/vff38_6j.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/62dfxntq.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=37483', 'data', 'file=/tmp/tmpygu7z4jl/vff38_6j.json', 'init=/tmp/tmpygu7z4jl/62dfxntq.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelxx52_4cf/prophet_model-20221222121312.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:12 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:12 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n",
"INFO:prophet:Disabling yearly seasonality. Run prophet with yearly_seasonality=True to override this.\n",
"INFO:prophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/5oofh9xr.json\n",
"DEBUG:cmdstanpy:input tempfile: /tmp/tmpygu7z4jl/nhgg7dt0.json\n",
"DEBUG:cmdstanpy:idx 0\n",
"DEBUG:cmdstanpy:running CmdStan, num_threads: None\n",
"DEBUG:cmdstanpy:CmdStan args: ['/home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages/prophet/stan_model/prophet_model.bin', 'random', 'seed=86598', 'data', 'file=/tmp/tmpygu7z4jl/5oofh9xr.json', 'init=/tmp/tmpygu7z4jl/nhgg7dt0.json', 'output', 'file=/tmp/tmpygu7z4jl/prophet_modelmo3oct1k/prophet_model-20221222121312.csv', 'method=optimize', 'algorithm=lbfgs', 'iter=10000']\n",
"12:13:12 - cmdstanpy - INFO - Chain [1] start processing\n",
"INFO:cmdstanpy:Chain [1] start processing\n",
"12:13:12 - cmdstanpy - INFO - Chain [1] done processing\n",
"INFO:cmdstanpy:Chain [1] done processing\n"
]
}
],
"source": [
"pred_return_rate = []\n",
"pred_return_rate_mean = []\n",
"i = 0 # current index of stock\n",
"N = 20 # number of candidate stocks\n",
"periods = data_test.shape[0] # number of days to forecast, same as data_test length\n",
"\n",
"for i in range(N): # for each stock\n",
" model = Prophet()\n",
" data['ds'] = data.index # index values (dates) as ds \n",
" data['y'] = data.iloc[:,i] # stock close price as y\n",
" stock_data = data[['ds','y']] # training data\n",
" model.fit(stock_data)\n",
" future = model.make_future_dataframe(periods=periods) # predict \n",
" forecast = model.predict(future)['yhat'] # predicted values stored in yhat column\n",
" return_rate = np.zeros(len(forecast)) # return rate to be stored in numpy array\n",
" for j in range(len(forecast)-1): \n",
" return_rate[j+1] = (forecast[j+1] - forecast[j])/forecast[j] # calculate the return rate per day\n",
" pred_return_rate.append(return_rate)\n",
" pred_return_rate_mean.append(np.mean(return_rate)) # predicted mean return rate for the stock\n",
" data.drop(columns=['ds','y'], inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "M7_Ijdiqg8Aj",
"outputId": "338dc52d-17e7-4406-b362-1b6364ead9b3"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-3.295120211693522e-05, 0.0004212352164733634, 0.0006443323333763341, 0.0006180876152479546, 0.0005054489754747284, -0.0011332381597290778, 0.0009594389542370872, -0.0004696525374264205, -0.0002568726047600979, -0.0010414517228963223, 0.00016392271944240828, 0.0006151834219843348, -9.028002715718702e-05, -6.64373367560175e-05, 0.0005966905208675229, -0.00019014641894963337, 0.0011606326569417041, 0.0005490521838582731, 0.0004985535894805963, -0.0005879920495425485]\n"
]
}
],
"source": [
"print(pred_return_rate_mean)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "aPGYYaL-diMU",
"outputId": "af9a60bf-8b6a-414e-fb85-1c1ec20298bb"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.00014317780640250334\n"
]
}
],
"source": [
"print(np.mean(pred_return_rate_mean))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "F2vo0ZanMwXh"
},
"source": [
"We also calculate actual return rates from the historical stock price data. We will use the results for the calculation of covariance that will be used for the second objective function to reduce the risk by diversify the portfolio items to avoid losing it altogether."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"id": "_yAkNF2qM2qP"
},
"outputs": [],
"source": [
"actual_return_rate = []\n",
"actual_return_rate_mean = []\n",
"N = 20 # number of candidate stocks\n",
"for i in range(N): # for each stock\n",
" stock_data = data.iloc[:,i]\n",
" return_rate = np.zeros(len(stock_data)) # return rate to be stored in numpy array\n",
" for j in range(len(stock_data)-1): \n",
" return_rate[j+1] = (stock_data[j+1] - stock_data[j])/stock_data[j] # calculate the return rate per day\n",
" actual_return_rate.append(return_rate)\n",
" actual_return_rate_mean.append(np.mean(return_rate)) # predicted mean return rate for the stock"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "gmgGmzeNWMcr",
"outputId": "75d34b1b-befe-446a-838e-1c527ee07f1d"
},
"outputs": [
{
"data": {
"text/plain": "0.00039797518633514584"
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"actual_return_rate_mean_mean = np.mean(actual_return_rate_mean)\n",
"actual_return_rate_mean_mean"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "TfEthgUNW6IM"
},
"source": [
"Finally, we calculate the actual return rates on the test set for the final evaluation."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"id": "UdmpvABkW_ic"
},
"outputs": [],
"source": [
"actual_return_rate_test = []\n",
"actual_return_rate_mean_test = []\n",
"N = 20 # number of candidate stocks\n",
"for i in range(N): # for each stock\n",
" stock_data = data_test.iloc[:,i]\n",
" return_rate = np.zeros(len(stock_data)) # return rate to be stored in numpy array\n",
" for j in range(len(stock_data)-1): \n",
" return_rate[j+1] = (stock_data[j+1] - stock_data[j])/stock_data[j] # calculate the return rate per day\n",
" actual_return_rate_test.append(return_rate)\n",
" actual_return_rate_mean_test.append(np.mean(return_rate)) # predicted mean return rate for the stock"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "o0vfvXq6Xh7Z",
"outputId": "7b1470e5-cd79-4749-f9f4-029dca6099c3"
},
"outputs": [
{
"data": {
"text/plain": "0.26481830738480705"
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"actual_return_rate_mean_mean_test = np.mean(actual_return_rate_mean_test)\n",
"actual_return_rate_mean_mean_test"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "UvVp2ir5sfZJ"
},
"source": [
"#3. Formulate the cost functions and constraints to obtain a QUBO matrix\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "cqZdX2OCNDs3"
},
"source": [
"## Constraint Terms"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "qCrz-zttNG-z"
},
"source": [
"In our portfolio optimisation example, let's assume 1 constraint: the number of stock to be included in the portfolio.\n",
"\n",
"In this example, we will use PyQubo, a python library that that can conver the cost function to a quadratic unconstraintbinary optimisation matrix that can be sent to D-wave qunatum annealer and JIJ simulated quantum annealer.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"id": "Z7pdl2fZPdFw"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collecting pyqubo\r\n",
" Downloading pyqubo-1.4.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (244 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m244.8/244.8 kB\u001B[0m \u001B[31m1.6 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0ma \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hCollecting Deprecated>=1.2.12\r\n",
" Downloading Deprecated-1.2.13-py2.py3-none-any.whl (9.6 kB)\r\n",
"Collecting dimod<0.13,>=0.9.14\r\n",
" Downloading dimod-0.12.3-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (15.9 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m15.9/15.9 MB\u001B[0m \u001B[31m33.3 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: six>=1.15.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pyqubo) (1.16.0)\r\n",
"Collecting dwave-neal>=0.5.7\r\n",
" Downloading dwave_neal-0.6.0-py3-none-any.whl (8.7 kB)\r\n",
"Requirement already satisfied: numpy>=1.17.3 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pyqubo) (1.23.4)\r\n",
"Requirement already satisfied: wrapt<2,>=1.10 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Deprecated>=1.2.12->pyqubo) (1.14.1)\r\n",
"Collecting dwave-samplers<2.0.0,>=1.0.0\r\n",
" Downloading dwave_samplers-1.0.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (6.1 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m6.1/6.1 MB\u001B[0m \u001B[31m27.1 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: networkx>=2.4.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-samplers<2.0.0,>=1.0.0->dwave-neal>=0.5.7->pyqubo) (2.8.7)\r\n",
"Installing collected packages: dimod, Deprecated, dwave-samplers, dwave-neal, pyqubo\r\n",
"Successfully installed Deprecated-1.2.13 dimod-0.12.3 dwave-neal-0.6.0 dwave-samplers-1.0.0 pyqubo-1.4.0\r\n",
"\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m A new release of pip available: \u001B[0m\u001B[31;49m22.3\u001B[0m\u001B[39;49m -> \u001B[0m\u001B[32;49m22.3.1\u001B[0m\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m To update, run: \u001B[0m\u001B[32;49mpip install --upgrade pip\u001B[0m\r\n"
]
}
],
"source": [
"!pip install pyqubo # https://pyqubo.readthedocs.io/en/latest/index.html"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"id": "9GfTXQ-PPuiK"
},
"outputs": [],
"source": [
"from pyqubo import Array, Constraint, Placeholder"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"id": "Rdza3plZQsQ1"
},
"outputs": [],
"source": [
"x = Array.create('x',shape=N, vartype='BINARY') # array takes binary 0 or 1 indicate to exclude or include a stock in the portfolio, which we need to optimise"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"id": "k13isroNsuNo"
},
"outputs": [],
"source": [
"# number of stocks constraint\n",
"num_stocks = 10 # specify the desired number of stocks\n",
"h_const1 = (num_stocks - np.dot(x,x))**2 # MSE from the budget. This lead to x values that makes the portfolio as close as the budget"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "gt4jvfQsjNE_"
},
"source": [
"## Cost Functions"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ualaqpV7BXFY"
},
"source": [
"We have set two cost functions to optimise our portfolio. One is to maximise the sum of predicted growth rates that we predicted in the feature engineering section. Another is to minimise the covariance between each of the stock items to be selected in the portfolio. We then add up two cost functions for QUBO."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"id": "4efqOsnhT_Gh"
},
"outputs": [],
"source": [
"# cost function 1 to aim the highest return\n",
"h_cost1 = 0\n",
"h_cost1 -= np.dot(x,pred_return_rate_mean) # maximisation problem. negative to make it minimisation problem"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"id": "V5bd5p0HeFAY"
},
"outputs": [],
"source": [
"# cost function 2 to balance the portfolio not to put all eggs in one basket\n",
"# minimising the sum of covariance leads to the safer choice\n",
"h_cost2 = 0\n",
"for i in range(N):\n",
" for j in range(N):\n",
" h_cost2 += x[i]*x[j]*np.sum((actual_return_rate[i]-actual_return_rate_mean[i])*(actual_return_rate[j]-actual_return_rate_mean[j]))/len(actual_return_rate[i])"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"id": "3TCBs2OKRpJ2"
},
"outputs": [],
"source": [
"h_cost = h_cost1 + h_cost2"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "pUh02l4qMA-2"
},
"source": [
"## Prepare QUBO"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "wiwvWGQHCaxr"
},
"source": [
"The problme formulation and representation in QUBO format are habled by PyQubo applications. In the lines below, we compile the model using the objective function that needs to be minimised, then define the constraint coefficient. to_qubo function generates a QUBO matrix in the dictionary type and an offset value which is an costanct value not required for the search for the minimum."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"id": "Oz6z5bXgXIZl"
},
"outputs": [],
"source": [
"h = h_cost + Placeholder('l1')*Constraint(h_const1,label='num_stocks')\n",
"model = h.compile()"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"id": "QtbTdpEja8pM"
},
"outputs": [],
"source": [
"feed_dict = {'l1': 2}\n",
"qubo, offset = model.to_qubo(feed_dict=feed_dict)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "4NgrOS7_v_zB",
"outputId": "b656b259-1a94-4612-bc8c-5c29139b420b"
},
"outputs": [
{
"data": {
"text/plain": "200.0"
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"offset"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "s1HR5k4L_3Rm"
},
"source": [
"# 4. Sample optimisation results from a D-Wave solver\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "m84Qzh9ZDZrU"
},
"source": [
"Install D-Wave's Ocean SDK to request the computation to the quantum annealer. EmbeddingComposite is used to embed the QUBO to the physical quantum anneler in D-Wave.\n"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"id": "Wgkgc61O_t7P"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collecting dwave-ocean-sdk\r\n",
" Downloading dwave_ocean_sdk-6.1.1-py3-none-any.whl (8.5 kB)\r\n",
"Requirement already satisfied: dimod==0.12.3 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-ocean-sdk) (0.12.3)\r\n",
"Requirement already satisfied: dwave-neal==0.6.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-ocean-sdk) (0.6.0)\r\n",
"Collecting dwave-system==1.18.0\r\n",
" Downloading dwave_system-1.18.0-py3-none-any.whl (103 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m103.5/103.5 kB\u001B[0m \u001B[31m1.1 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0ma \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hCollecting dwavebinarycsp==0.2.0\r\n",
" Downloading dwavebinarycsp-0.2.0-py3-none-any.whl (35 kB)\r\n",
"Collecting pyqubo==1.3.1\r\n",
" Downloading pyqubo-1.3.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (236 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m236.5/236.5 kB\u001B[0m \u001B[31m2.8 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hCollecting dwave-greedy==0.3.0\r\n",
" Downloading dwave_greedy-0.3.0-py3-none-any.whl (10 kB)\r\n",
"Collecting dwave-preprocessing==0.5.3\r\n",
" Downloading dwave_preprocessing-0.5.3-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (2.3 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m2.3/2.3 MB\u001B[0m \u001B[31m13.2 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0ma \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hCollecting dwave-hybrid==0.6.9\r\n",
" Downloading dwave_hybrid-0.6.9-py3-none-any.whl (74 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m74.7/74.7 kB\u001B[0m \u001B[31m18.9 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting dwave-cloud-client==0.10.3\r\n",
" Downloading dwave_cloud_client-0.10.3-py3-none-any.whl (111 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m111.4/111.4 kB\u001B[0m \u001B[31m27.5 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting dwave-inspector==0.3.0\r\n",
" Downloading dwave_inspector-0.3.0-py3-none-any.whl (26 kB)\r\n",
"Collecting minorminer==0.2.9\r\n",
" Downloading minorminer-0.2.9-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (5.6 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m5.6/5.6 MB\u001B[0m \u001B[31m27.0 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0ma \u001B[36m0:00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: dwave-samplers==1.0.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-ocean-sdk) (1.0.0)\r\n",
"Collecting dwave-tabu==0.5.0\r\n",
" Downloading dwave_tabu-0.5.0-py3-none-any.whl (9.2 kB)\r\n",
"Collecting dwave-networkx==0.8.12\r\n",
" Downloading dwave_networkx-0.8.12-py3-none-any.whl (97 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m97.5/97.5 kB\u001B[0m \u001B[31m26.8 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting penaltymodel==1.0.2\r\n",
" Downloading penaltymodel-1.0.2-py3-none-any.whl (36 kB)\r\n",
"Requirement already satisfied: numpy<2.0.0,>=1.17.3 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dimod==0.12.3->dwave-ocean-sdk) (1.23.4)\r\n",
"Collecting homebase>=1.0\r\n",
" Downloading homebase-1.0.1-py2.py3-none-any.whl (11 kB)\r\n",
"Collecting plucky>=0.4.3\r\n",
" Downloading plucky-0.4.3-py2.py3-none-any.whl (10 kB)\r\n",
"Collecting diskcache>=5.2.1\r\n",
" Downloading diskcache-5.4.0-py3-none-any.whl (44 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m45.0/45.0 kB\u001B[0m \u001B[31m15.0 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting pydantic>=1.7.3\r\n",
" Downloading pydantic-1.10.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (13.6 MB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m13.6/13.6 MB\u001B[0m \u001B[31m41.6 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m00:01\u001B[0m00:01\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: click>=7.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-cloud-client==0.10.3->dwave-ocean-sdk) (8.1.3)\r\n",
"Requirement already satisfied: requests[socks]>=2.18 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-cloud-client==0.10.3->dwave-ocean-sdk) (2.28.1)\r\n",
"Requirement already satisfied: python-dateutil>=2.7 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-cloud-client==0.10.3->dwave-ocean-sdk) (2.8.2)\r\n",
"Requirement already satisfied: networkx in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-hybrid==0.6.9->dwave-ocean-sdk) (2.8.7)\r\n",
"Collecting Flask>=1.1.1\r\n",
" Downloading Flask-2.2.2-py3-none-any.whl (101 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m101.5/101.5 kB\u001B[0m \u001B[31m28.2 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hRequirement already satisfied: importlib-resources>=3.2.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-inspector==0.3.0->dwave-ocean-sdk) (5.2.0)\r\n",
"Requirement already satisfied: scipy>=1.7.3 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from dwave-system==1.18.0->dwave-ocean-sdk) (1.9.3)\r\n",
"Collecting rectangle-packer>=2.0.1\r\n",
" Downloading rectangle_packer-2.0.1-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (302 kB)\r\n",
"\u001B[2K \u001B[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001B[0m \u001B[32m302.6/302.6 kB\u001B[0m \u001B[31m21.2 MB/s\u001B[0m eta \u001B[36m0:00:00\u001B[0m\r\n",
"\u001B[?25hCollecting fasteners\r\n",
" Downloading fasteners-0.18-py3-none-any.whl (18 kB)\r\n",
"Requirement already satisfied: six>=1.15.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pyqubo==1.3.1->dwave-ocean-sdk) (1.16.0)\r\n",
"Requirement already satisfied: Deprecated>=1.2.12 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pyqubo==1.3.1->dwave-ocean-sdk) (1.2.13)\r\n",
"Requirement already satisfied: wrapt<2,>=1.10 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Deprecated>=1.2.12->pyqubo==1.3.1->dwave-ocean-sdk) (1.14.1)\r\n",
"Requirement already satisfied: Werkzeug>=2.2.2 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Flask>=1.1.1->dwave-inspector==0.3.0->dwave-ocean-sdk) (2.2.2)\r\n",
"Requirement already satisfied: Jinja2>=3.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Flask>=1.1.1->dwave-inspector==0.3.0->dwave-ocean-sdk) (3.0.3)\r\n",
"Requirement already satisfied: importlib-metadata>=3.6.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Flask>=1.1.1->dwave-inspector==0.3.0->dwave-ocean-sdk) (5.0.0)\r\n",
"Collecting itsdangerous>=2.0\r\n",
" Downloading itsdangerous-2.1.2-py3-none-any.whl (15 kB)\r\n",
"Requirement already satisfied: zipp>=3.1.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from importlib-resources>=3.2.0->dwave-inspector==0.3.0->dwave-ocean-sdk) (3.9.0)\r\n",
"Requirement already satisfied: typing-extensions>=4.1.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from pydantic>=1.7.3->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (4.4.0)\r\n",
"Requirement already satisfied: charset-normalizer<3,>=2 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests[socks]>=2.18->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (2.1.1)\r\n",
"Requirement already satisfied: idna<4,>=2.5 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests[socks]>=2.18->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (3.4)\r\n",
"Requirement already satisfied: urllib3<1.27,>=1.21.1 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests[socks]>=2.18->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (1.26.12)\r\n",
"Requirement already satisfied: certifi>=2017.4.17 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests[socks]>=2.18->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (2022.9.24)\r\n",
"Requirement already satisfied: PySocks!=1.5.7,>=1.5.6 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from requests[socks]>=2.18->dwave-cloud-client==0.10.3->dwave-ocean-sdk) (1.7.1)\r\n",
"Requirement already satisfied: MarkupSafe>=2.0 in /home/obm/Prog/miniconda3/envs/qml/lib/python3.8/site-packages (from Jinja2>=3.0->Flask>=1.1.1->dwave-inspector==0.3.0->dwave-ocean-sdk) (2.1.1)\r\n",
"Installing collected packages: rectangle-packer, plucky, homebase, pydantic, itsdangerous, fasteners, diskcache, penaltymodel, Flask, dwave-preprocessing, dwave-networkx, minorminer, dwavebinarycsp, dwave-tabu, dwave-greedy, dwave-cloud-client, pyqubo, dwave-system, dwave-inspector, dwave-hybrid, dwave-ocean-sdk\r\n",
" Attempting uninstall: pyqubo\r\n",
" Found existing installation: pyqubo 1.4.0\r\n",
" Uninstalling pyqubo-1.4.0:\r\n",
" Successfully uninstalled pyqubo-1.4.0\r\n",
"Successfully installed Flask-2.2.2 diskcache-5.4.0 dwave-cloud-client-0.10.3 dwave-greedy-0.3.0 dwave-hybrid-0.6.9 dwave-inspector-0.3.0 dwave-networkx-0.8.12 dwave-ocean-sdk-6.1.1 dwave-preprocessing-0.5.3 dwave-system-1.18.0 dwave-tabu-0.5.0 dwavebinarycsp-0.2.0 fasteners-0.18 homebase-1.0.1 itsdangerous-2.1.2 minorminer-0.2.9 penaltymodel-1.0.2 plucky-0.4.3 pydantic-1.10.2 pyqubo-1.3.1 rectangle-packer-2.0.1\r\n",
"\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m A new release of pip available: \u001B[0m\u001B[31;49m22.3\u001B[0m\u001B[39;49m -> \u001B[0m\u001B[32;49m22.3.1\u001B[0m\r\n",
"\u001B[1m[\u001B[0m\u001B[34;49mnotice\u001B[0m\u001B[1;39;49m]\u001B[0m\u001B[39;49m To update, run: \u001B[0m\u001B[32;49mpip install --upgrade pip\u001B[0m\r\n"
]
}
],
"source": [
"!pip install dwave-ocean-sdk # https://docs.ocean.dwavesys.com/en/stable/"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"id": "9slCue62bZeJ"
},
"outputs": [],
"source": [
"from dwave.system import DWaveSampler, EmbeddingComposite\n",
"endpoint = 'https://cloud.dwavesys.com/sapi'\n",
"token = '***' # Please relace with your own token"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "s0K_qwB2GbeE"
},
"source": [
"Available annelers depend on the geographic location and machine's online status. The list of solvers can be found on the D-Wave Leap dashboard. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "Gv8-8NEZb7Ve",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"sampler_dw = DWaveSampler(solver='Advantage_system4.1',token=token)\n",
"sampler_qa = EmbeddingComposite(sampler_dw)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Y8UAX9eNnKLU",
"outputId": "7baceb91-3cf6-4a9d-ff5b-24b49030157b",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"sampleset_qa = sampler_qa.sample_qubo(qubo,num_reads=10)\n",
"records_qa = sampleset_qa.record\n",
"print(records_qa)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "FcnZ5rnlGx_F"
},
"source": [
"Above is the result of 10 samples from the solver. \n",
"\n",
"First list items are the optimised combination of the stock items, 1 indicates the stock is included in the portfolio. \n",
"\n",
"The second value is the energy state where the lowest number is the best solution. The number shows close to minus two hundred, because the offset value is not included. \n",
"\n",
"Third number is the number of occurance of the solution. We have 10 reads and each read has unique stock selections.\n",
"\n",
"The last number indicate \"chain break\" that the connection between qubits are broken then fixed by the software. Ideal soliution should have 0. chain break ratio.\n",
"\n",
"Using records returned from the D-Wave solver, we can validate if the constraint terms are complied."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "f8CAuYfAoHWY",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_candidates_qa = []\n",
"num_stock_counts_qa = []\n",
"record_len = len(records_qa)\n",
"for i in range(record_len):\n",
" portfolio = []\n",
" num_stock_count = 0\n",
" for j in range(len(records_qa[i][0])):\n",
" if records_qa[i][0][j] == 1:\n",
" portfolio.append(stocks[j])\n",
" num_stock_count += 1\n",
" portfolio_candidates_qa.append(portfolio)\n",
" num_stock_counts_qa.append(num_stock_count)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "fbhM3LzcpfTT",
"outputId": "165ba3db-46bb-428d-f916-dc42364624bd",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"num_stock_counts_qa # number of stocks selected for the portfolio"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "QYi817jh0VkQ"
},
"source": [
"In our example, 4 of the 10 solutions fit the requirement. We may adjust the coeffient value by increasing the penalty."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "-rApXs6WKAWS"
},
"source": [
"Finally, we can select one best solution and display the sotck codes.\n",
"We will use 7th record where there was no chain break."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "UAQWDqxUqwxw",
"outputId": "aceef0e2-1117-435d-d650-0fd040bef612",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"best_portfolio_qa = portfolio_candidates_qa[7]\n",
"print(best_portfolio_qa)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "_O2vKFSj0s5v"
},
"source": [
"D-Wave has hybrid solver that uses classical computer and quantum annealer. \n",
"The solver 'hybrid_binary_quadratic_model_version2' can have up to 1,000,000 variables.\n",
"\n",
"To try the hybrid solver, we need to import LeapHybridSampler library. The difference from the quantum solver are: \n",
"\n",
"- it does not require embed composite\n",
"- only one sample record is returned (no num_reads parameter)\n",
"- it consumes more QPU time (different conversion rate applies QPU: Hybrid 1:20)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Ual7XMplxg0G",
"outputId": "de83b02c-6a00-4f90-8f99-9531d27927f4",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"from dwave.system import LeapHybridSampler\n",
"sampler_dw_hybrid = LeapHybridSampler(solver='hybrid_binary_quadratic_model_version2',token=token)\n",
"sampleset_qa_hybrid = sampler_dw_hybrid.sample_qubo(qubo)\n",
"records_qa_hybrid = sampleset_qa_hybrid.record\n",
"print(records_qa_hybrid)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "fryOP91hMKGd"
},
"source": [
"# 5. Sample optimisation results from a simulated quantum annealing solver\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "-SfZxfL4Klhu"
},
"source": [
"There is publicly available open source QA simulators provided by [OpenJIJ](https://www.openjij.org/). It can be used as an alternative the real quantum annerler. We can use the same QUBO prepared by PyQubo.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "mlbr7BGVk1lG",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"!pip install openjij"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "EdLM1IcYcm9D",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"from openjij import SQASampler\n",
"sampler_jij = SQASampler()\n",
"sampleset_sqa = sampler_jij.sample_qubo(qubo,num_reads=10)\n",
"records_sqa = sampleset_sqa.record"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "NxaMuv38x_A3",
"outputId": "8d0d94db-8a5c-4c90-ef08-4597d25eb1d9",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"print(records_sqa)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "0Z-fhb_1kv9g",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_candidates_sqa = []\n",
"num_stock_counts_sqa = []\n",
"record_len = len(records_sqa)\n",
"for i in range(record_len):\n",
" portfolio = []\n",
" num_stock_count = 0\n",
" for j in range(len(records_sqa[i][0])):\n",
" if records_sqa[i][0][j] == 1:\n",
" portfolio.append(stocks[j])\n",
" num_stock_count += 1\n",
" portfolio_candidates_sqa.append(portfolio)\n",
" num_stock_counts_sqa.append(num_stock_count)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "NjiaDZPpyReI",
"outputId": "a83ee422-0972-4ff1-d396-67725c3f2bb3",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"num_stock_counts_sqa"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Vy2OQFnzxvyJ",
"outputId": "a4ae761d-1870-46b8-9032-a404efff5806",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"best_portfolio_sqa = portfolio_candidates_sqa[-1]\n",
"print(best_portfolio_sqa)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Ai4w_FVITS_K"
},
"source": [
"#6. Evaluate the outcomes\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "QFCaegY63nlS"
},
"source": [
"We can select the best portfolio selection from QA and SQA samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "PJIhTJps3mdB",
"outputId": "64de0934-9d2a-460c-82e4-43a7f8e694e4",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_qa = records_qa[7][0] # best combination from QA hybdrid solver\n",
"print(portfolio_qa)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "MCAoXGr03_2f",
"outputId": "dbaccaf6-90d9-426d-8d45-985d625d6a5a",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_sqa = records_sqa[-1][0] # best combination from SQA hybdrid solver\n",
"print(portfolio_sqa)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "5bE_Q0uI_1cK"
},
"source": [
"We compare the cost 1: increases the return rate.\n",
"The more return is achieved if the number is smaller than the other."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Hp29zsrE4KRc",
"outputId": "d59a9744-eba7-4966-a5d6-2e9a49a85c09",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"h_cost1_qa = 0\n",
"h_cost1_qa -= np.dot(portfolio_qa,actual_return_rate_mean_test) # maximisation problem. negative to make it minimisation problem\n",
"print(h_cost1_qa) # return rate factor by QA"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "3YG-1hsE43S3",
"outputId": "5d33ea5d-e464-4b06-84f9-b70ba58cb73f",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"h_cost1_sqa = 0\n",
"h_cost1_sqa -= np.dot(portfolio_sqa,actual_return_rate_mean_test) # maximisation problem. negative to make it minimisation problem\n",
"print(h_cost1_sqa) # return rate cost by SQA"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "WzwJDZyTAaIO"
},
"source": [
"We also compare the cost 2: decrease the risk.\n",
"The more return is achieved if the number is smaller than the other."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "QKwaUtMq4gLk",
"outputId": "6b2d73d8-03f3-48f1-b7e0-33fba0c97e7b",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"N= 20\n",
"h_cost2_qa = 0\n",
"for i in range(N):\n",
" for j in range(N):\n",
" h_cost2_qa += portfolio_qa[i]*portfolio_qa[j]*np.sum((actual_return_rate_test[i]-actual_return_rate_mean_test[i])*(actual_return_rate_test[j]-actual_return_rate_mean_test[j]))/len(actual_return_rate_test[i])\n",
"print(h_cost2_qa)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Qkyixcj75Ckz",
"outputId": "c0aefee2-7e0f-4338-acd4-5c31bb1de561",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"N= 20\n",
"h_cost2_sqa = 0\n",
"for i in range(N):\n",
" for j in range(N):\n",
" h_cost2_sqa += portfolio_sqa[i]*portfolio_sqa[j]*np.sum((actual_return_rate_test[i]-actual_return_rate_mean_test[i])*(actual_return_rate_test[j]-actual_return_rate_mean_test[j]))/len(actual_return_rate_test[i])\n",
"print(h_cost2_sqa)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "dehQVD1PBGMm"
},
"source": [
"We can check the costs in case we selected all 20 stocks as benchmark.\n",
"The numbers will be divided by two for comparison."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "yqtwWnx_-eRR",
"outputId": "283a93ca-c173-44bd-cc5c-6d2ad12f99d3",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_all = np.ones(20)\n",
"h_cost1_all = 0\n",
"h_cost1_all -= np.dot(portfolio_all,actual_return_rate_mean_test) # maximisation problem. negative to make it minimisation problem\n",
"print(h_cost1_all/2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "9pZVXJmD_Ar7",
"outputId": "52004dac-c26b-43e0-da64-f0ff3cc60f37",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"N= 20\n",
"h_cost2_all = 0\n",
"for i in range(N):\n",
" for j in range(N):\n",
" h_cost2_all += portfolio_all[i]*portfolio_all[j]*np.sum((actual_return_rate_test[i]-actual_return_rate_mean_test[i])*(actual_return_rate_test[j]-actual_return_rate_mean_test[j]))/len(actual_return_rate_test[i])\n",
"print(h_cost2_all/2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "NXOrt37ECBSF"
},
"source": [
"Finally, the costs from Hybrid solver as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Uv_kWSvpBWkB",
"outputId": "d401de8e-8a1c-4194-8568-4024f2451317",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"portfolio_qa_hybrid = records_qa_hybrid[0][0] # best combination from QA hybdrid solver\n",
"print(portfolio_qa_hybrid)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "L0mN6mj4BjjH",
"outputId": "b9acacd9-6282-4c8e-d8b2-8dc9e074746b",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"h_cost1_qa_hybrid = 0\n",
"h_cost1_qa_hybrid -= np.dot(portfolio_qa_hybrid,actual_return_rate_mean_test) # maximisation problem. negative to make it minimisation problem\n",
"print(h_cost1_qa_hybrid) # return rate factor by QA Hybrid Solver"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "owuDyxIGBxa-",
"outputId": "2629f91f-afce-43ba-a68b-ac8a91b0af79",
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"N= 20\n",
"h_cost2_qa_hybrid = 0\n",
"for i in range(N):\n",
" for j in range(N):\n",
" h_cost2_qa_hybrid += portfolio_qa_hybrid[i]*portfolio_qa_hybrid[j]*np.sum((actual_return_rate_test[i]-actual_return_rate_mean_test[i])*(actual_return_rate_test[j]-actual_return_rate_mean_test[j]))/len(actual_return_rate_test[i])\n",
"print(h_cost2_qa_hybrid)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "4vc2E4sVCSMj"
},
"source": [
"With compariso to the benchmark, we can observe sample reocrds from QA/SQA/Hybdrid solvers as: \n",
"\n",
"- QA solver showed high-return but high-risks.\n",
"- SQA solver recommended slightly lower return but minimum risk\n",
"- Hybrid solver showed very similar costs for both return and risk as SQA\n",
"\n",
"Above is just a comparison of 1 sample record. We could analyse more sample records and visualsing the outcome. We could also tune hyper-parameters, like constraint coefficient."
]
}
],
"metadata": {
"colab": {
"authorship_tag": "ABX9TyPipMoNtOXQNeppvssSErAZ",
"collapsed_sections": [],
"include_colab_link": true,
"name": "TDS_QA_Example_Portfolio_Optimisation.ipynb",
"provenance": []
},
"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.9.16"
}
},
"nbformat": 4,
"nbformat_minor": 0
}