challenge-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dc9e6d3c",
   "metadata": {},
   "source": [
    "## IBM Quantum Challenge Fall 2021\n",
    "# Challenge 4: Battery revenue optimization\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "We recommend that you switch to **light** workspace theme under the Account menu in the upper right corner for optimal experience."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "092340ed",
   "metadata": {},
   "source": [
    "# Introduction to QAOA\n",
    "\n",
    "When it comes to optimization problems, a well-known algorithm for finding approximate solutions to combinatorial-optimization problems is **QAOA (Quantum approximate optimization algorithm)**. You may have already used it once in the finance exercise of Challenge-1, but still don't know what it is. In this challlenge we will further learn about QAOA----how does it work? Why we need it?\n",
    "\n",
    "First off, what is QAOA?  Simply put, QAOA is a classical-quantum  hybrid algorithm that combines a parametrized quantum circuit known as ansatz, and a classical part to optimize those circuits proposed by Farhi, Goldstone, and Gutmann (2014)[**[1]**](https://arxiv.org/abs/1411.4028). \n",
    "\n",
    "It is a variational algorithm that uses a unitary $U(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ characterized by the parameters $(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ to prepare a quantum state $|\\psi(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})\\rangle$. The goal of the algorithm is to find optimal parameters $(\\boldsymbol{\\beta}_{opt}, \\boldsymbol{\\gamma}_{opt})$ such that the quantum state $|\\psi(\\boldsymbol{\\beta}_{opt}, \\boldsymbol{\\gamma}_{opt})\\rangle$ encodes the solution to the problem. \n",
    "\n",
    "The unitary $U(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ has a specific form and is composed of two unitaries $U(\\boldsymbol{\\beta}) = e^{-i \\boldsymbol{\\beta} H_B}$ and $U(\\boldsymbol{\\gamma}) = e^{-i \\boldsymbol{\\gamma} H_P}$ where $H_{B}$ is the mixing Hamiltonian and $H_{P}$ is the problem Hamiltonian. Such a choice of unitary drives its inspiration from a related scheme called quantum annealing.\n",
    "\n",
    "The state is prepared by applying these unitaries as alternating blocks of the two unitaries applied $p$ times such that \n",
    "\n",
    "$$\\lvert \\psi(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma}) \\rangle = \\underbrace{U(\\boldsymbol{\\beta}) U(\\boldsymbol{\\gamma}) \n",
    "                                            \\cdots U(\\boldsymbol{\\beta}) U(\\boldsymbol{\\gamma})}_{p \\; \\text{times}} \n",
    "\\lvert \\psi_0 \\rangle$$\n",
    "\n",
    "where $|\\psi_{0}\\rangle$ is a suitable initial state.\n",
    "\n",
    "<center><img src=\"resources/qaoa_circuit.png\" width=\"600\"></center>\n",
    "\n",
    "The QAOA implementation of Qiskit directly extends VQE and inherits VQE’s general hybrid optimization structure.\n",
    "To learn more about QAOA, please refer to the [**QAOA chapter**](https://qiskit.org/textbook/ch-applications/qaoa.html) of Qiskit Textbook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b58deef",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Goal**\n",
    "\n",
    "Implement the quantum optimization code for the battery revenue problem. <br>\n",
    "<br>\n",
    "**Plan**\n",
    "\n",
    "First, you will learn about QAOA and knapsack problem.\n",
    "<br><br>\n",
    "**Challenge 4a** - Simple knapsack problem with QAOA: familiarize yourself with a typical knapsack problem and  find the optimized solution with QAOA.\n",
    "<br><br>\n",
    "**Final Challenge 4b** - Battery revenue optimization with Qiskit knapsack class: learn the battery revenue optimization problem and find the optimized solution with QAOA. You can receive a badge for solving all the challenge exercises up to 4b.\n",
    "<br><br>\n",
    "**Final Challenge 4c** - Battery revenue optimization with your own quantum circuit: implement the battery revenue optimization problem to find the lowest circuit cost and circuit depth. Achieve better accuracy with smaller circuits. you can obtain a score with ranking by solving this exercise.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "Before you begin, we recommend watching the [**Qiskit Optimization Demo Session with Atsushi Matsuo**](https://youtu.be/claoY57eVIc?t=104) and check out the corresponding [**demo notebook**](https://github.com/qiskit-community/qiskit-application-modules-demo-sessions/tree/main/qiskit-optimization) to learn how to do classifications using QSVM.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe37bc4",
   "metadata": {},
   "source": [
    "As we just mentioned, QAOA is an algorithm which can be used to find approximate solutions to combinatorial optimization problems, which includes many specific problems, such as:\n",
    "\n",
    "- TSP (Traveling Salesman Problem) problem\n",
    "- Vehicle routing problem\n",
    "- Set cover problem\n",
    "- Knapsack problem\n",
    "- Scheduling problems,etc. \n",
    "\n",
    "Some of them are hard to solve (or in another word, they are NP-hard problems), and it is impractical to find their exact solutions in a reasonable amount of time, and that is why we need the approximate algorithm. Next, we will introduce an instance of using QAOA to solve one of the combinatorial optimization problems----**knapsack problem**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d22c66c",
   "metadata": {},
   "source": [
    "# Knapsack Problem #\n",
    "\n",
    "[**Knapsack Problem**](https://en.wikipedia.org/wiki/Knapsack_problem) is an optimization problem that goes like this: given a list of items that each has a weight and a value and a knapsack that can hold a maximum weight. Determine which items to take in the knapsack so as to maximize the total value taken without exceeding the maiximum weight the knapsack can hold. The most efficient approach would be a greedy approach, but that is not guaranteed to give the best result.\n",
    "\n",
    "\n",
    "<center><img src=\"resources/Knapsack.png\" width=\"400\"></center>\n",
    "    \n",
    "Image source: [Knapsack.svg.](https://commons.wikimedia.org/w/index.php?title=File:Knapsack.svg&oldid=457280382)\n",
    "\n",
    "<div id='problem'></div>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Note: Knapsack problem have many variations, here we will only discuss the 0-1 Knapsack problem: either take an item or not (0-1 property), which is a NP-hard problem. We can not divide one item, or take multiple same items.    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d504c427",
   "metadata": {},
   "source": [
    "## Challenge 4a: Simple knapsack problem with QAOA \n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    \n",
    "**Challenge 4a** <br>\n",
    "You are given a knapsack with a capacity of 18 and 5 pieces of luggage. When the weights of each piece of luggage $W$ is $w_i = [4,5,6,7,8]$ and the value $V$ is $v_i = [5,6,7,8,9]$, find the packing method that maximizes the sum of the values of the luggage within the capacity limit of 18.\n",
    "\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f0c1e977",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_optimization.algorithms import MinimumEigenOptimizer\n",
    "from qiskit import Aer\n",
    "from qiskit.utils import algorithm_globals, QuantumInstance\n",
    "from qiskit.algorithms import QAOA, NumPyMinimumEigensolver\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399c5a92",
   "metadata": {},
   "source": [
    "## Dynamic Programming Approach ##\n",
    "\n",
    "A typical classical method for finding an exact solution, the Dynamic Programming approach is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e7395687",
   "metadata": {},
   "outputs": [],
   "source": [
    "val = [5,6,7,8,9]\n",
    "wt = [4,5,6,7,8]\n",
    "W = 18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49787a13",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal value: 21\n",
      "\n",
      " index of the chosen items:\n",
      "1 2 3 "
     ]
    }
   ],
   "source": [
    "def dp(W, wt, val, n):\n",
    "    k = [[0 for x in range(W + 1)] for x in range(n + 1)]\n",
    "    for i in range(n + 1):\n",
    "        for w in range(W + 1):\n",
    "            if i == 0 or w == 0:\n",
    "                k[i][w] = 0\n",
    "            elif wt[i-1] <= w:\n",
    "                k[i][w] = max(val[i-1] + k[i-1][w-wt[i-1]], k[i-1][w])\n",
    "            else:\n",
    "                k[i][w] = k[i-1][w]\n",
    "                \n",
    "    picks=[0 for x in range(n)]\n",
    "    volume=W\n",
    "    for i in range(n,-1,-1):\n",
    "        if (k[i][volume]>k[i-1][volume]):\n",
    "            picks[i-1]=1\n",
    "            volume -= wt[i-1]\n",
    "    return k[n][W],picks\n",
    "\n",
    "n = len(val)\n",
    "print(\"optimal value:\", dp(W, wt, val, n)[0])\n",
    "print('\\n index of the chosen items:')\n",
    "for i in range(n): \n",
    "    if dp(W, wt, val, n)[1][i]: \n",
    "        print(i,end=' ')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58cd150f",
   "metadata": {},
   "source": [
    "The time complexity of this method $O(N*W)$, where $N$ is the number of items and $W$ is the maximum weight of the knapsack. We can solve this problem using an exact solution approach within a reasonable time since the number of combinations is limited, but when the number of items becomes huge, it will be impractical to deal with by using a exact solution approach. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ec9bdb",
   "metadata": {},
   "source": [
    "## QAOA approach ##\n",
    "\n",
    "Qiskit provides application classes for various optimization problems, including the knapsack problem so that users can easily try various optimization problems on quantum computers. In this exercise, we are going to use the application classes for the `Knapsack` problem here. \n",
    "\n",
    "There are application classes for other optimization problems available as well. See [**Application Classes for Optimization Problems**](https://qiskit.org/documentation/optimization/tutorials/09_application_classes.html#Knapsack-problem) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "83af8b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages necessary for application classes.\n",
    "from qiskit_optimization.applications import Knapsack"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb966cf",
   "metadata": {},
   "source": [
    "To represent Knapsack problem as an optimization problem that can be solved by QAOA, we need to formulate the cost function for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6b2b3b3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\\ This file has been generated by DOcplex\n",
       "\\ ENCODING=ISO-8859-1\n",
       "\\Problem name: Knapsack\n",
       "\n",
       "Maximize\n",
       " obj: 5 x_0 + 6 x_1 + 7 x_2 + 8 x_3 + 9 x_4\n",
       "Subject To\n",
       " c0: 4 x_0 + 5 x_1 + 6 x_2 + 7 x_3 + 8 x_4 <= 18\n",
       "\n",
       "Bounds\n",
       " 0 <= x_0 <= 1\n",
       " 0 <= x_1 <= 1\n",
       " 0 <= x_2 <= 1\n",
       " 0 <= x_3 <= 1\n",
       " 0 <= x_4 <= 1\n",
       "\n",
       "Binaries\n",
       " x_0 x_1 x_2 x_3 x_4\n",
       "End"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def knapsack_quadratic_program():\n",
    "    # Put values, weights and max_weight parameter for the Knapsack()\n",
    "    \n",
    "    ##############################\n",
    "    # Provide your code here\n",
    "    \n",
    "    prob = Knapsack(values=val, weights=wt, max_weight=W)\n",
    "    \n",
    "    #\n",
    "    ##############################\n",
    "    \n",
    "    # to_quadratic_program generates a corresponding QuadraticProgram of the instance of the knapsack problem.\n",
    "    kqp = prob.to_quadratic_program()\n",
    "    return prob, kqp\n",
    "\n",
    "prob,quadratic_program=knapsack_quadratic_program()\n",
    "quadratic_program"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be543eeb",
   "metadata": {},
   "source": [
    "We can solve the problem using the classical `NumPyMinimumEigensolver` to find the minimum eigenvector, which may be useful as a reference without doing things by Dynamic Programming; we can also apply QAOA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "095fd7cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result:\n",
      " optimal function value: 21.0\n",
      "optimal value: [0. 1. 1. 1. 0.]\n",
      "status: SUCCESS\n",
      "\n",
      " index of the chosen items: [1, 2, 3]\n"
     ]
    }
   ],
   "source": [
    "# Numpy Eigensolver\n",
    "meo = MinimumEigenOptimizer(min_eigen_solver=NumPyMinimumEigensolver())\n",
    "result = meo.solve(quadratic_program)\n",
    "print('result:\\n', result)\n",
    "print('\\n index of the chosen items:', prob.interpret(result)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9c8e6e93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result:\n",
      " optimal function value: 21.0\n",
      "optimal value: [0. 1. 1. 1. 0.]\n",
      "status: SUCCESS\n",
      "\n",
      " index of the chosen items: [1, 2, 3]\n"
     ]
    }
   ],
   "source": [
    "# QAOA\n",
    "seed = 123\n",
    "algorithm_globals.random_seed = seed\n",
    "qins = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), shots=1000, seed_simulator=seed, seed_transpiler=seed)\n",
    "\n",
    "meo = MinimumEigenOptimizer(min_eigen_solver=QAOA(reps=1, quantum_instance=qins))\n",
    "result = meo.solve(quadratic_program)\n",
    "print('result:\\n', result)\n",
    "print('\\n index of the chosen items:', prob.interpret(result)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f05fa585",
   "metadata": {},
   "source": [
    "You will submit the quadratic program created by your `knapsack_quadratic_program` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "28b68aac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitting your answer for 4a. Please wait...\n",
      "Congratulations 🎉! Your answer is correct and has been submitted.\n"
     ]
    }
   ],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex4a\n",
    "grade_ex4a(quadratic_program)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b84d7f78",
   "metadata": {},
   "source": [
    "<div id='problem'></div>\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Note: QAOA finds the approximate solutions, so the solution by QAOA is not always optimal. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e13be264",
   "metadata": {},
   "source": [
    "# Battery Revenue Optimization Problem #\n",
    "\n",
    "<div id='problem'></div>\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "In this exercise we will use a quantum algorithm to solve a real-world instance of a combinatorial optimization problem: Battery revenue optimization problem.   \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e3295b",
   "metadata": {},
   "source": [
    "Battery storage systems have provided a solution to flexibly integrate large-scale renewable energy (such as wind and solar) in a power system. The revenues from batteries come from different types of services sold to the grid. The process of energy trading of battery storage assets is as follows: A regulator asks each battery supplier to choose a market in advance for each time window. Then, the batteries operator will charge the battery with renewable energy and release the energy to the grid depending on pre-agreed contracts. The supplier makes therefore forecasts on the return and the number of charge/discharge cycles for each time window to optimize its overall return. \n",
    "\n",
    "How to maximize the revenue of battery-based energy storage is a concern of all battery storage investors. Choose to let the battery always supply power to the market which pays the most for every time window might be a simple guess, but in reality, we have to consider many other factors. \n",
    "\n",
    "What we can not ignore is the aging of batteries, also known as **degradation**. As the battery charge/discharge cycle progresses, the battery capacity will gradually degrade (the amount of energy a battery can store, or the amount of power it can deliver will permanently reduce). After a number of cycles, the battery will reach the end of its usefulness. Since the performance of a battery decreases while it is used, choosing the best cash return for every time window one after the other, without considering the degradation, does not lead to an optimal return over the lifetime of the battery, i.e. before the number of charge/discharge cycles reached.\n",
    "\n",
    "Therefore, in order to optimize the revenue of the battery, what we have to do is to select the market for the battery in each time window taking both **the returns on these markets (value)**, based on price forecast, as well as expected battery **degradation over time (cost)** into account ——It sounds like solving a common optimization problem, right?\n",
    "\n",
    "We will investigate how quantum optimization algorithms could be adapted to tackle this problem.\n",
    "\n",
    "\n",
    "<div>\n",
    "<p></p>\n",
    "<center><img src=\"resources/renewable-g7ac5bd48e_640.jpg\" width=\"600\"></center>\n",
    "\n",
    "</div>\n",
    "\n",
    "Image source: [pixabay](https://pixabay.com/photos/renewable-energy-environment-wind-1989416/)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22a19f05",
   "metadata": {},
   "source": [
    "## Problem Setting\n",
    "\n",
    "Here, we have referred to the problem setting in de la Grand'rive and Hullo's paper [**[2]**](https://arxiv.org/abs/1908.02210).\n",
    "\n",
    "Considering two markets $M_{1}$ , $M_{2}$, during every time window (typically a day), the battery operates on one or the other market, for a maximum of $n$ time windows. Every day is considered independent and the intraday optimization is a standalone problem: every morning the battery starts with the same level of power so that we don’t consider charging problems. Forecasts on both markets being available for the $n$ time windows, we assume known for each time window $t$ (day) and for each market:\n",
    "\n",
    "- the daily returns $\\lambda_{1}^{t}$ , $\\lambda_{2}^{t}$\n",
    "\n",
    "- the daily degradation, or health cost (number of cycles), for the battery $c_{1}^{t}$, $c_{2}^{t}$ \n",
    "\n",
    "We want to find the optimal schedule, i.e. optimize the life time return with a cost less than $C_{max}$ cycles. We introduce $d = max_{t}\\left\\{c_{1}^{t}, c_{2}^{t}\\right\\} $.\n",
    "\n",
    "We introduce the decision variable $z_{t}, \\forall t \\in [1, n]$ such that $z_{t} = 0$ if the supplier chooses $M_{1}$ , $z_{t} = 1$ if choose $M_{2}$, with every possible vector $z = [z_{1}, ..., z_{n}]$ being a possible schedule. The previously formulated problem can then be expressed as:\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\underset{z \\in \\left\\{0,1\\right\\}^{n}}{max} \\displaystyle\\sum_{t=1}^{n}(1-z_{t})\\lambda_{1}^{t}+z_{t}\\lambda_{2}^{t}\n",
    "\\end{equation}\n",
    "<br>\n",
    "\\begin{equation}\n",
    "    s.t. \\sum_{t=1}^{n}[(1-z_{t})c_{1}^{t}+z_{t}c_{2}^{t}]\\leq C_{max}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85efc6d2",
   "metadata": {},
   "source": [
    "This does not look like one of the well-known combinatorial optimization problems, but no worries! we are going to give hints on how to solve this problem with quantum computing step by step."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1938c6e1",
   "metadata": {},
   "source": [
    "# Challenge 4b: Battery revenue optimization with Qiskit knapsack class #\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4b** <br>\n",
    "We will optimize the battery schedule using Qiskit optimization knapsack class with QAOA to maximize the total return with a cost within $C_{max}$ under the following conditions;\n",
    "    <br>\n",
    "- the time window $t = 7$<br>\n",
    "- the daily return $\\lambda_{1} = [5, 3, 3, 6, 9, 7, 1]$<br>\n",
    "- the daily return $\\lambda_{2} = [8, 4, 5, 12, 10, 11, 2]$<br>\n",
    "- the daily degradation for the battery $c_{1} = [1, 1, 2, 1, 1, 1, 2]$<br>\n",
    "- the daily degradation for the battery $c_{2} = [3, 2, 3, 2, 4, 3, 3]$<br>\n",
    "- $C_{max} = 16$<br>\n",
    "<br>\n",
    "    \n",
    "Your task is to find the argument, `values`, `weights`, and `max_weight` used for the Qiskit optimization knapsack class, to get a solution which \"0\" denote the choice of market $M_{1}$, and \"1\" denote the choice of market $M_{2}$. We will check your answer with another data set of $\\lambda_{1}, \\lambda_{2}, c_{1}, c_{2}, C_{max}$.\n",
    "    <br>\n",
    "You can receive a badge for solving all the challenge exercises up to 4b.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6ba82106",
   "metadata": {},
   "outputs": [],
   "source": [
    "L1 = [5,3,3,6,9,7,1]\n",
    "L2 = [8,4,5,12,10,11,2]\n",
    "C1 = [1,1,2,1,1,1,2]\n",
    "C2 = [3,2,3,2,4,3,3]\n",
    "C_max = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c4ed4006",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 1, 2, 6, 1, 4, 1] [2, 1, 1, 1, 3, 2, 1] 7\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\\ This file has been generated by DOcplex\n",
       "\\ ENCODING=ISO-8859-1\n",
       "\\Problem name: Knapsack\n",
       "\n",
       "Maximize\n",
       " obj: 3 x_0 + x_1 + 2 x_2 + 6 x_3 + x_4 + 4 x_5 + x_6\n",
       "Subject To\n",
       " c0: 2 x_0 + x_1 + x_2 + x_3 + 3 x_4 + 2 x_5 + x_6 <= 7\n",
       "\n",
       "Bounds\n",
       " 0 <= x_0 <= 1\n",
       " 0 <= x_1 <= 1\n",
       " 0 <= x_2 <= 1\n",
       " 0 <= x_3 <= 1\n",
       " 0 <= x_4 <= 1\n",
       " 0 <= x_5 <= 1\n",
       " 0 <= x_6 <= 1\n",
       "\n",
       "Binaries\n",
       " x_0 x_1 x_2 x_3 x_4 x_5 x_6\n",
       "End"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def knapsack_argument(L1, L2, C1, C2, C_max):\n",
    "      \n",
    "    ##############################\n",
    "    # Provide your code here\n",
    "    length = len(L1)\n",
    "    \n",
    "    values = [L2[i]-L1[i] for i in range(length)]\n",
    "    weights = [C2[i]-C1[i] for i in range(length)]\n",
    "    max_weight = C_max - sum(C1)\n",
    "    #\n",
    "    ##############################\n",
    "    return values, weights, max_weight\n",
    "    \n",
    "values, weights, max_weight = knapsack_argument(L1, L2, C1, C2, C_max)\n",
    "print(values, weights, max_weight)\n",
    "prob = Knapsack(values = values, weights = weights, max_weight = max_weight)\n",
    "qp = prob.to_quadratic_program()\n",
    "qp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "179661e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running \"knapsack_argument\" (1/3)... \n",
      "Running \"knapsack_argument\" (2/3)... \n",
      "Running \"knapsack_argument\" (3/3)... \n",
      "Submitting your answer for 4b. Please wait...\n",
      "Congratulations 🎉! Your answer is correct and has been submitted.\n"
     ]
    }
   ],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex4b\n",
    "grade_ex4b(knapsack_argument)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356027b6",
   "metadata": {},
   "source": [
    "We can solve the problem using QAOA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "eed75fda",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result: [1. 1. 1. 1. 0. 1. 0.]\n",
      "total revenue: 50\n"
     ]
    }
   ],
   "source": [
    "# QAOA\n",
    "seed = 123\n",
    "algorithm_globals.random_seed = seed\n",
    "qins = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), shots=1000, seed_simulator=seed, seed_transpiler=seed)\n",
    "\n",
    "meo = MinimumEigenOptimizer(min_eigen_solver=QAOA(reps=1, quantum_instance=qins))\n",
    "result = meo.solve(qp)\n",
    "print('result:', result.x)\n",
    "\n",
    "item = np.array(result.x)\n",
    "revenue=0\n",
    "for i in range(len(item)):\n",
    "    if item[i]==0:\n",
    "        revenue+=L1[i]\n",
    "    else:\n",
    "        revenue+=L2[i]\n",
    "\n",
    "print('total revenue:', revenue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99ca6f63",
   "metadata": {},
   "source": [
    "## Challenge 4c:  Battery revenue optimization with adiabatic quantum computation\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "\n",
    "Here we come to the final exercise! The final challenge is for people to compete in ranking.\n",
    "\n",
    "</div>\n",
    "\n",
    "## Background\n",
    "\n",
    "QAOA was developed with inspiration from adiabatic quantum computation. In adiabatic quantum computation, based on the quantum adiabatic theorem, the ground state of a given Hamiltonian can ideally be obtained. Therefore, by mapping the optimization problem to this Hamiltonian, it is possible to solve the optimization problem with adiabatic quantum computation.\n",
    "\n",
    "Although the computational equivalence of adiabatic quantum computation and quantum circuits has been shown, simulating adiabatic quantum computation on quantum circuits involves a large number of gate operations, which is difficult to achieve with current noisy devices. QAOA solves this problem by using a quantum-classical hybrid approach.\n",
    "\n",
    "In this extra challenge, you will be asked to implement a quantum circuit that solves an optimization problem without classical optimization, based on this adiabatic quantum computation framework. In other words, the circuit you build is expected to give a good approximate solution in a single run.\n",
    "\n",
    "Instead of using the Qiskit Optimization Module and Knapsack class, let's try to implement a quantum circuit with as few gate operations as possible, that is, as small as possible. By relaxing the constraints of the optimization problem, it is possible to find the optimum solution with a smaller circuit. We recommend that you follow the solution tips.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c**<br>\n",
    "We will optimize the battery schedule using the adiabatic quantum computation to maximize the total return with a cost within $C_{max}$ under the following conditions;\n",
    "<br>\n",
    "- the time window $t = 11$<br>\n",
    "- the daily return $\\lambda_{1} = [3, 7, 3, 4, 2, 6, 2, 2, 4, 6, 6]$<br>\n",
    "- the daily return $\\lambda_{2} = [7, 8, 7, 6, 6, 9, 6, 7, 6, 7, 7]$<br>\n",
    "- the daily degradation for the battery $c_{1} = [2, 2, 2, 3, 2, 4, 2, 2, 2, 2, 2]$<br>\n",
    "- the daily degradation for the battery $c_{2} = [4, 3, 3, 4, 4, 5, 3, 4, 4, 3, 4]$<br>\n",
    "- $C_{max} = 33$<br>\n",
    "    - **Note:** $\\lambda_{1}[i] < \\lambda_{2}[i]$ and $c_{1}[i] < c_{2}[i]$ holds for $i \\in \\{1,2,...,t\\}$\n",
    "<br>\n",
    "    \n",
    "\n",
    "Let \"0\" denote the choice of market $M_{1}$ and \"1\" denote the choice of market $M_{2}$, the optimal solutions are \"00111111000\", and \"10110111000\" with return value $67$ and cost $33$.\n",
    "Your task is to implement adiabatic quantum computation circuit to meet the accuracy below. \n",
    "We will check your answer with other data set of $\\lambda_{1}, \\lambda_{2}, c_{1}, c_{2}, C_{max}$. We show examples of inputs for checking below. We will use similar inputs with these examples.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "8808362f",
   "metadata": {},
   "outputs": [],
   "source": [
    "instance_examples = [\n",
    "    {\n",
    "        'L1': [3, 7, 3, 4, 2, 6, 2, 2, 4, 6, 6],\n",
    "        'L2': [7, 8, 7, 6, 6, 9, 6, 7, 6, 7, 7],\n",
    "        'C1': [2, 2, 2, 3, 2, 4, 2, 2, 2, 2, 2],\n",
    "        'C2': [4, 3, 3, 4, 4, 5, 3, 4, 4, 3, 4],\n",
    "        'C_max': 33\n",
    "    },\n",
    "    {\n",
    "        'L1': [4, 2, 2, 3, 5, 3, 6, 3, 8, 3, 2],\n",
    "        'L2': [6, 5, 8, 5, 6, 6, 9, 7, 9, 5, 8],\n",
    "        'C1': [3, 3, 2, 3, 4, 2, 2, 3, 4, 2, 2],\n",
    "        'C2': [4, 4, 3, 5, 5, 3, 4, 5, 5, 3, 5],\n",
    "        'C_max': 38\n",
    "    },\n",
    "    {\n",
    "        'L1': [5, 4, 3, 3, 3, 7, 6, 4, 3, 5, 3],\n",
    "        'L2': [9, 7, 5, 5, 7, 8, 8, 7, 5, 7, 9],\n",
    "        'C1': [2, 2, 4, 2, 3, 4, 2, 2, 2, 2, 2],\n",
    "        'C2': [3, 4, 5, 4, 4, 5, 3, 3, 5, 3, 5],\n",
    "        'C_max': 35\n",
    "    }\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0fdda3d",
   "metadata": {},
   "source": [
    "### IMPORTANT: Final exercise submission rules\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "<p style=\"font-size:20px\">\n",
    "For solving this problem:\n",
    "</p>\n",
    "\n",
    "- Do not optimize with classical methods.\n",
    "- Create a quantum circuit by filling source code in the functions along the following steps.\n",
    "- As for the parameters $p$ and $\\alpha$, please **do not change the values from $p=5$ and $\\alpha=1$.**\n",
    "- Please implement the quantum circuit within 28 qubits.\n",
    "- You should submit a function that takes (L1, L2, C1, C2, C_max) as inputs and returns a QuantumCircuit. (You can change the name of the function in your way.)\n",
    "- Your circuit should be able to solve different input values. We will validate your circuit with several inputs.   \n",
    "- Create a circuit that gives precision of 0.8 or better with lower cost. The precision is explained below. The lower the cost, the better.\n",
    "\n",
    "- Please **do not run jobs in succession** even if you are concerned that your job is not running properly. This can create a long queue and clog the backend. You can check whether your job is running properly at:\n",
    "[**https://quantum-computing.ibm.com/jobs**](https://quantum-computing.ibm.com/jobs)\n",
    "        \n",
    "- Judges will check top 10 solutions manually to see if their solutions adhere to the rules. **Please note that your ranking is subject to change after the challenge period as a result of the judging process.**\n",
    "\n",
    "- Top 10 participants will be recognized and asked to submit a write up on how they solved the exercise.\n",
    "    \n",
    "**Note: In this challenge, please be aware that you should solve the problem with a quantum circuit, otherwise you will not have a rank in the final ranking.**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cebe1ca",
   "metadata": {},
   "source": [
    "### Scoring Rule\n",
    "\n",
    "\n",
    "The score of submitted function is computed by two steps.\n",
    "\n",
    "1. In the first step, the precision of output of your quantum circuit is checked.\n",
    "To pass this step, your circuit should output a probability distribution whose **average precision is more than 0.80** for eight instances; four of them are fixed data, while the remaining four are randomly selected data from multiple datasets.\n",
    "If your circuit cannot satisfy this threshold **0.8**, you will not obtain a score.\n",
    "We will explain how the precision of a probability distribution will be calculated when the submitted quantum circuit solves one instance.\n",
    "    \n",
    "    1. This precision evaluates how the values of measured feasible solutions are close to the value of optimal solutions.\n",
    "    2. Firstly the number of measured feasible solutions is very low, the precision will be 0 (Please check **\"The number of feasible solutions\"** below). <br>Before calculating precision, the values of solutions will be normalized so that the precision of the solution whose value is the lowest would be always 0 by subtracting the lowest value.\n",
    "    \n",
    "    Let $N_s$, $N_f$, and $\\lambda_{opt}$ be the total shots (the number of execution), the shots of measured feasible solutions, the optimial solution value. Also let $R(x)$ and $C(x)$ be value and cost of a solution $x\\in\\{0,1\\}^n$ respectively. We normalize the values by subtracting the lowest value of instance, which can be calculated by the summation of $\\lambda_{1}$.\n",
    "    Given a probability distribution, the precision is computed with the following formula:\n",
    "    \n",
    "    \\begin{equation*}\n",
    "    \\text{precision} = \\frac 1 {N_f\\cdot (\\lambda_{opt}-\\mathrm{sum}(\\lambda_{1}) )} \\sum_{x, \\text{$\\mathrm{shots}_x$}\\in \\text{ prob.dist.}} (R(x)-\\mathrm{sum}(\\lambda_{1})) \\cdot \\text{$\\mathrm{shots}_x$} \\cdot 1_{C(x) \\leq C_{max}}\n",
    "    \\end{equation*}\n",
    "                                                                         \n",
    "    Here, $\\mathrm{shots}_x$ is the counts of measuring the solution $x$. For example, given a probability distribution {\"1000101\": 26, \"1000110\": 35, \"1000111\": 12, \"1001000\": 16, \"1001001\": 11} with shots $N_s = 100$,\n",
    "    the value and the cost of each solution are listed below.\n",
    "    | Solution | Value | Cost | Feasible or not | Shot counts |\n",
    "    |:-------:|:-------:|:-------:|:-------:|:--------------:|\n",
    "    | 1000101 | 46 | 16 | 1 | 26 |\n",
    "    | 1000110 | 48 | 17 | 0 | 35 |\n",
    "    | 1000111 | 45 | 15 | 1 | 12 |\n",
    "    | 1001000 | 45 | 18 | 0 | 16 |\n",
    "    | 1001001 | 42 | 16 | 1 | 11 |\n",
    "\n",
    "    Since $C_{max}= 16$, the solutions \"1000101\", \"1000111\", and \"1001001\" are feasible, but the solutions \"1000110\" and \"1001000\" are infeasible. So, the shots of measured feasbile solutions $N_f$ is calculated as $N_f = 26+12+11=49$. And the lowest value is $ \\mathrm{sum}(\\lambda_{1}) = 5+3+3+6+9+7+1=34$. \n",
    "    Therefore, the precision becomes \n",
    "    \n",
    "    $$((46-34) \\cdot 26 \\cdot 1 + (48-34) \\cdot 35 \\cdot 0 + (45-34) \\cdot 12 \\cdot 1 + (45-34) \\cdot 16 \\cdot 0 + (42-34) \\cdot 11 \\cdot 1) / (49\\cdot (50-34)) = 0.68$$\n",
    "    \n",
    "    **The number of feasible solutions**: If $N_f$ is less than 20 ($ N_f < 20$), the precision will be calculated as 0.\n",
    "\n",
    "2. In the second step, the score of your quantum circuit will be evaluated only if your solution passes the first step.\n",
    "The score is the sum of circuit costs of four instances, where the circuit cost is calculated as below.\n",
    "\n",
    "    1. Transpile the quantum circuit without gate optimization and decompose the gates into the basis gates of \"rz\", \"sx\", \"cx\".\n",
    "    2. Then the score is calculated by \n",
    "    \n",
    "    \\begin{equation*}\n",
    "    \\text{score} = 50 \\cdot depth + 10 \\cdot \\#(cx) + \\#(rz) + \\#(sx)\n",
    "    \\end{equation*}\n",
    "    \n",
    "    where $\\#(gate)$ denotes the number of $gate$ in the circuit.\n",
    "    \n",
    "Your circuit will be executed <span style=\"color: deepskyblue; \">512 times</span>, which means <span style=\"color: deepskyblue; \">$N_s = 512$</span> here.<br>\n",
    "The smaller the score become, the higher you will be ranked."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80281422",
   "metadata": {},
   "source": [
    "## General Approach\n",
    "\n",
    "Here we are making the answer according to the way shown in [**[2]**](https://arxiv.org/abs/1908.02210), which is solving the \"relaxed\" formulation of knapsack problem.\n",
    "The relaxed problem can be defined as follows:\n",
    "\\begin{equation*}\n",
    "\\text{maximize } f(z)=return(z)+penalty(z)\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "\\text{where} \\quad return(z)=\\sum_{t=1}^{n} return_{t}(z) \\quad \\text{with} \\quad return_{t}(z) \\equiv\\left(1-z_{t}\\right) \\lambda_{1}^{t}+z_{t} \\lambda_{2}^{t}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "\\quad \\quad \\quad \\quad \\quad \\quad penalty(z)=\\left\\{\\begin{array}{ll}\n",
    "0 & \\text{if}\\quad  cost(z)<C_{\\max } \\\\\n",
    "-\\alpha\\left(cost(z)-C_{\\max }\\right) & \\text{if} \\quad cost(z) \\geq C_{\\max }, \\alpha>0 \\quad \\text{constant}\n",
    "\\end{array}\\right.\n",
    "\\end{equation*}\n",
    "\n",
    "A non-Ising target function to compute a linear penalty is used here.\n",
    "This may reduce the depth of the circuit while still achieving high accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e21ae75",
   "metadata": {},
   "source": [
    "The basic unit of relaxed approach consisits of the following items.\n",
    "1. Phase Operator $U(C, \\gamma_i)$\n",
    "    1. return part\n",
    "    2. penalty part\n",
    "        1. Cost calculation (data encoding)\n",
    "        2. Constraint testing (marking the indices whose data exceed $C_{max}$)\n",
    "        3. Penalty dephasing (adding penalty to the marked indices)\n",
    "        4. Reinitialization of constraint testing and cost calculation (clean the data register and flag register)\n",
    "2. Mixing Operator $U(B, \\beta_i)$\n",
    "\n",
    "This procedure unit $U(B, \\beta_i)U(C, \\gamma_i)$ will be totally repeated $p$ times in the whole relaxed QAOA procedure.\n",
    "<br>\n",
    "Let's take a look at each function one by one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2694a2eb",
   "metadata": {},
   "source": [
    "The quantum circuit we are going to make consists of three types of registers: index register, data register, and flag register.\n",
    "Index register and data register are used for QRAM which contain the cost data for every possible choice of battery.\n",
    "Here these registers appear in the function templates named as follows:\n",
    "- `qr_index`: a quantum register representing the index (the choice of 0 or 1 in each time window)\n",
    "- `qr_data`: a quantum register representing the total cost associated with each index\n",
    "- `qr_f`: a quantum register that store the flag for penalty dephasing\n",
    "<br>\n",
    "\n",
    "We also use the following variables to represent the number of qubits in each register.\n",
    "- `index_qubits`: the number of qubits in `qr_index`\n",
    "- `data_qubits`: the number of qubits in `qr_data`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bc0afab",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 1**\n",
    "    \n",
    "</div>\n",
    "\n",
    "## Phase Operator $U(C, \\gamma_i)$\n",
    "### Return Part\n",
    "The return part $return (z)$ can be transformed as follows:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "e^{-i \\gamma_i . return(z)}\\left|z\\right\\rangle \n",
    "&=\\prod_{t=1}^{n} e^{-i \\gamma_i return_{t}(z)}\\left|z\\right\\rangle \\\\\n",
    "&=e^{i \\theta} \\bigotimes_{t=1}^{n} e^{-i \\gamma_i z_{t}\\left(\\lambda_{2}^{t}-\\lambda_{1}^{t}\\right)}\\left|z_{t}\\right\\rangle \\\\\n",
    "\\text{with}\\quad \\theta &=\\sum_{t=1}^{n} \\lambda_{1}^{t}\\quad \\text{constant}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "Since we can ignore the constant phase rotation, the return part $return (z)$ can be realized by rotation gate $U_1\\left(\\gamma_i \\left(\\lambda_{2}^{t}-\\lambda_{1}^{t}\\right)\\right)= e^{-i \\frac{\\gamma_i \\left(\\lambda_{2}^{t}-\\lambda_{1}^{t}\\right)}  2}$ for each qubit.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `phase_return` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "89ed0535",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List, Union\n",
    "import math\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, assemble\n",
    "from qiskit.compiler import transpile\n",
    "from qiskit.circuit import Gate\n",
    "from qiskit.circuit.library.standard_gates import *\n",
    "from qiskit.circuit.library import QFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "52fa0f7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def phase_return(index_qubits: int, gamma: float, L1: list, L2: list, to_gate=True): #  -> Union[Gate, QuantumCircuit]\n",
    "    qr_index = QuantumRegister(index_qubits, \"index\")\n",
    "    qc = QuantumCircuit(qr_index)\n",
    "    \n",
    "    ##############################\n",
    "    ### U_1(gamma * (lambda2 - lambda1)) for each qubit ###\n",
    "    # Provide your code here\n",
    "    length = len(qr_index)\n",
    "    assert length == len(L1) == len(L2), \"The arrays' lengths are different!\"\n",
    "    for i in range(length):\n",
    "        qc.u1(gamma * (L1[i]-L2[i]),qr_index[i])\n",
    "\n",
    "    ##############################\n",
    "    \n",
    "    return qc.to_gate(label=\" phase return \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "eec09377-365a-4e77-a221-47d2930339ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.gate.Gate at 0x7f6067c620a0>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For debugging purposes\n",
    "phase_return(11, 0.1, instance_examples[0]['L1'], instance_examples[0]['L2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e380f679",
   "metadata": {},
   "source": [
    "## Phase Operator $U(C, \\gamma_i)$\n",
    "### Penalty Part\n",
    "\n",
    "In this part, we are considering how to add penalty to the quantum states in index register whose total cost exceed the constraint $C_{max}$.\n",
    "\n",
    "As shown above, this can be realized by the following four steps.\n",
    "\n",
    "1. Cost calculation (data encoding)\n",
    "2. Constraint testing (marking the indices whose data value exceed $C_{max}$)\n",
    "3. Penalty dephasing (adding penalty to the marked indices)\n",
    "4. Reinitialization of constraint testing and cost calculation (clean the data register and flag register)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bba3d122",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 2**\n",
    "    \n",
    "</div>\n",
    "\n",
    "#### Cost calculation (data encoding)\n",
    "\n",
    "To represent the sum of cost for every choice of answer, we can use QRAM structure.\n",
    "In order to implement QRAM by quantum circuit, the addition function would be helpful.\n",
    "Here we will first prepare a function for constant value addition.\n",
    "<br>\n",
    "<br>\n",
    "To add a constant value to data we can use\n",
    "- Series of full adders\n",
    "- Plain adder network [**[3]**](https://arxiv.org/abs/quant-ph/9511018)\n",
    "- Ripple carry adder [**[4]**](https://arxiv.org/abs/quant-ph/0410184)\n",
    "- QFT adder **[[5](https://arxiv.org/abs/quant-ph/0008033), [6](https://arxiv.org/abs/1411.5949)]**\n",
    "- etc...\n",
    "<br>\n",
    "\n",
    "Each adder has its own characteristics. Here, for example, we will briefly explain how to implement QFT adder, which is less likely increase circuits cost when the number of additions increases.\n",
    "1. QFT on the target quantum register\n",
    "2. Local phase rotation on the target quantum register controlled by quantum register for the constant\n",
    "3. IQFT on the target quantum register\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `const_adder` and `subroutine_add_const` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "693f9fbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def subroutine_add_const(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:\n",
    "    qc = QuantumCircuit(data_qubits)\n",
    "    ##############################\n",
    "    ### Phase Rotation ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    return qc.to_gate(label=\" [+\"+str(const)+\"] \") if to_gate else qc\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8b8652e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def const_adder(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_data = QuantumRegister(data_qubits, \"data\")\n",
    "    qc = QuantumCircuit(qr_data)\n",
    "    \n",
    "    ##############################\n",
    "    ### QFT ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "\n",
    "    ##############################\n",
    "    ### Phase Rotation ###\n",
    "    # Use `subroutine_add_const`\n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    \n",
    "    ##############################\n",
    "    ### IQFT ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    return qc.to_gate(label=\" [ +\" + str(const) + \"] \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46a63366",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 3**\n",
    "    \n",
    "</div>\n",
    "\n",
    "Here we want to store the cost in a QRAM form: \n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_{x\\in\\{0,1\\}^t} \\left|x\\right\\rangle\\left|cost(x)\\right\\rangle\n",
    "\\end{equation*}\n",
    "\n",
    "where $t$ is the number of time window (=size of index register), and $x$ is the pattern of battery choice through all the time window.\n",
    "\n",
    "Given two lists $C^1 = \\left[c_0^1, c_1^1, \\cdots\\right]$ and $C^2 = \\left[c_0^2, c_1^2, \\cdots\\right]$, \n",
    "we can encode the \"total sum of each choice\" of these data using controlled gates by each index qubit.\n",
    "If we want to add $c_i^1$ to the data whose $i$-th index qubit is $0$ and $c_i^2$ to the data whose $i$-th index qubit is $1$, \n",
    "then we can add $C_i^1$ to data register when the $i$-th qubit in index register is $0$,\n",
    "and $C_i^2$ to data register when the $i$-th qubit in index register is $1$.\n",
    "These operation can be realized by controlled gates.\n",
    "If you want to create controlled gate from gate with type `qiskit.circuit.Gate`, the `control()` method might be useful.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `cost_calculation` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b54d86a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost_calculation(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_index = QuantumRegister(index_qubits, \"index\")\n",
    "    qr_data = QuantumRegister(data_qubits, \"data\")\n",
    "    qc = QuantumCircuit(qr_index, qr_data)\n",
    "\n",
    "    for i, (val1, val2) in enumerate(zip(list1, list2)):\n",
    "        \n",
    "        ##############################\n",
    "        ### Add val2 using const_adder controlled by i-th index register (set to 1) ###\n",
    "        # Provide your code here\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "        qc.x(qr_index[i])\n",
    "        \n",
    "        ##############################\n",
    "        ### Add val1 using const_adder controlled by i-th index register (set to 0) ###\n",
    "        # Provide your code here\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "        qc.x(qr_index[i])\n",
    "    \n",
    "    return qc.to_gate(label=\" Cost Calculation \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a1f7be5",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 4**\n",
    "    \n",
    "</div>\n",
    "\n",
    "#### Constraint Testing\n",
    "\n",
    "After the cost calculation process, we have gained the entangled QRAM with flag qubits set to zero for all indices:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_{x\\in\\{0,1\\}^t} \\left|x\\right\\rangle\\left|cost(x)\\right\\rangle\\left|0\\right\\rangle\n",
    "\\end{equation*}\n",
    "\n",
    "In order to selectively add penalty to those indices with cost values larger than $C_{max}$, we have to prepare the following state:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sum_{x\\in\\{0,1\\}^t} \\left|x\\right\\rangle\\left|cost(x)\\right\\rangle\\left|cost(x)\\geq C_{max}\\right\\rangle\n",
    "\\end{equation*}\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `constraint_testing` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c75b888",
   "metadata": {},
   "outputs": [],
   "source": [
    "def constraint_testing(data_qubits: int, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_data = QuantumRegister(data_qubits, \"data\")\n",
    "    qr_f = QuantumRegister(1, \"flag\")\n",
    "    qc = QuantumCircuit(qr_data, qr_f)\n",
    "    \n",
    "    ##############################\n",
    "    ### Set the flag register for indices with costs larger than C_max ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    \n",
    "    return qc.to_gate(label=\" Constraint Testing \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20032242",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 5**\n",
    "    \n",
    "</div>\n",
    "\n",
    "#### Penalty Dephasing\n",
    "\n",
    "We also have to add penalty to the indices with total costs larger than $C_{max}$ in the following way.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\quad \\quad \\quad \\quad \\quad \\quad penalty(z)=\\left\\{\\begin{array}{ll}\n",
    "0 & \\text{if}\\quad  cost(z)<C_{\\max } \\\\\n",
    "-\\alpha\\left(cost(z)-C_{\\max }\\right) & \\text{if} \\quad cost(z) \\geq C_{\\max }, \\alpha>0 \\quad \\text{constant}\n",
    "\\end{array}\\right.\n",
    "\\end{equation*}\n",
    "\n",
    "This penalty can be described as quantum operator $e^{i \\gamma \\alpha\\left(cost(z)-C_{\\max }\\right)}$.\n",
    "<br>\n",
    "To realize this unitary operator as quantum circuit, we focus on the following property.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\alpha\\left(cost(z)-C_{m a x}\\right)=\\sum_{j=0}^{k-1} 2^{j} \\alpha A_{1}[j]-2^{c} \\alpha\n",
    "\\end{equation*}\n",
    "\n",
    "where $A_1$ is the quantum register for qram data, $A_1[j]$ is the $j$-th qubit of $A_1$, and $k$ and $c$ are appropriate constants.\n",
    "\n",
    "Using this property, the penalty rotation part can be realized as rotation gates on each digit of data register of QRAM controlled by the flag register.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `penalty_dephasing` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6b96d4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def penalty_dephasing(data_qubits: int, alpha: float, gamma: float, to_gate = True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_data = QuantumRegister(data_qubits, \"data\")\n",
    "    qr_f = QuantumRegister(1, \"flag\")\n",
    "    qc = QuantumCircuit(qr_data, qr_f)\n",
    "    \n",
    "    ##############################\n",
    "    ### Phase Rotation ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    \n",
    "    return qc.to_gate(label=\" Penalty Dephasing \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8380cb5f",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 6**\n",
    "    \n",
    "</div>\n",
    "\n",
    "#### Reinitialization\n",
    "\n",
    "The ancillary qubits such as the data register and the flag register should be reinitialized to zero states when the operator $U(C, \\gamma_i)$ finishes.\n",
    "<br>\n",
    "If you want to apply inverse unitary of a `qiskit.circuit.Gate`, the `inverse()` method might be useful.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "Fill in the blank in the following cell to complete the `reinitialization` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9ab61f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def reinitialization(index_qubits: int, data_qubits: int, C1: list, C2: list, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_index = QuantumRegister(index_qubits, \"index\")\n",
    "    qr_data = QuantumRegister(data_qubits, \"data\")\n",
    "    qr_f = QuantumRegister(1, \"flag\")\n",
    "    qc = QuantumCircuit(qr_index, qr_data, qr_f)\n",
    "    \n",
    "    ##############################\n",
    "    ### Reinitialization Circuit ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    \n",
    "    return qc.to_gate(label=\" Reinitialization \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dbd7d41",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 7**\n",
    "    \n",
    "</div>\n",
    "\n",
    "### Mixing Operator $U(B, \\beta_i)$\n",
    "\n",
    "Finally, we have to add the mixing operator $U(B,\\beta_i)$ after phase operator $U(C,\\gamma_i)$.\n",
    "The mixing operator can be represented as follows.\n",
    "\\begin{equation*}\n",
    "U(B, \\beta_i)=\\exp (-i \\beta_i B)=\\prod_{i=j}^{n} \\exp \\left(-i \\beta_i \\sigma_{j}^{x}\\right)\n",
    "\\end{equation*}\n",
    "\n",
    "This operator can be realized by $R_x(2\\beta_i)$ gate on each qubits in index register.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "Fill in the blank in the following cell to complete the `mixing_operator` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d18fc934",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixing_operator(index_qubits: int, beta: float, to_gate = True) -> Union[Gate, QuantumCircuit]:\n",
    "    \n",
    "    qr_index = QuantumRegister(index_qubits, \"index\")\n",
    "    qc = QuantumCircuit(qr_index)\n",
    "    \n",
    "    ##############################\n",
    "    ### Mixing Operator ###\n",
    "    # Provide your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##############################\n",
    "    \n",
    "    return qc.to_gate(label=\" Mixing Operator \") if to_gate else qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d57cc41",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 4c - Step 8**\n",
    "    \n",
    "</div>\n",
    "\n",
    "Finally, using the functions we have created above, we will make the submit function `solver_function` for whole relaxed QAOA process.\n",
    "\n",
    "\n",
    "Fill the TODO blank in the following cell to complete the answer function.\n",
    "- You can copy and paste the function you have made above.\n",
    "- you may also adjust the number of qubits and its arrangement if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2a921f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solver_function(L1: list, L2: list, C1: list, C2: list, C_max: int) -> QuantumCircuit:\n",
    "    \n",
    "    # the number of qubits representing answers\n",
    "    index_qubits = len(L1)\n",
    "    \n",
    "    # the maximum possible total cost\n",
    "    max_c = sum([max(l0, l1) for l0, l1 in zip(C1, C2)])\n",
    "    \n",
    "    # the number of qubits representing data values can be defined using the maximum possible total cost as follows:\n",
    "    data_qubits = math.ceil(math.log(max_c, 2)) + 1 if not max_c & (max_c - 1) == 0 else math.ceil(math.log(max_c, 2)) + 2\n",
    "    \n",
    "    ### Phase Operator ###\n",
    "    # return part\n",
    "    def phase_return():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "    \n",
    "    # penalty part\n",
    "    def subroutine_add_const():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "\n",
    "    # penalty part\n",
    "    def const_adder():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "    \n",
    "    # penalty part\n",
    "    def cost_calculation():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "    \n",
    "    # penalty part\n",
    "    def constraint_testing():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "    \n",
    "    # penalty part\n",
    "    def penalty_dephasing():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "        \n",
    "    # penalty part\n",
    "    def reinitialization():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "\n",
    "    ### Mixing Operator ###\n",
    "    def mixing_operator():\n",
    "        ##############################\n",
    "        ### TODO ###\n",
    "        ### Paste your code from above cells here ###\n",
    "        \n",
    "        \n",
    "        \n",
    "        ##############################\n",
    "    \n",
    "    qr_index = QuantumRegister(index_qubits, \"index\") # index register\n",
    "    qr_data = QuantumRegister(data_qubits, \"data\") # data register\n",
    "    qr_f = QuantumRegister(1, \"flag\") # flag register\n",
    "    cr_index = ClassicalRegister(index_qubits, \"c_index\") # classical register storing the measurement result of index register\n",
    "    qc = QuantumCircuit(qr_index, qr_data, qr_f, cr_index)\n",
    "    \n",
    "    ### initialize the index register with uniform superposition state ###\n",
    "    qc.h(qr_index)\n",
    "    \n",
    "    ### DO NOT CHANGE THE CODE BELOW\n",
    "    p = 5\n",
    "    alpha = 1\n",
    "    for i in range(p):\n",
    "        \n",
    "        ### set fixed parameters for each round ###\n",
    "        beta = 1 - (i + 1) / p\n",
    "        gamma = (i + 1) / p\n",
    "        \n",
    "        ### return part ###\n",
    "        qc.append(phase_return(index_qubits, gamma, L1, L2), qr_index)\n",
    "        \n",
    "        ### step 1: cost calculation ###\n",
    "        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2), qr_index[:] + qr_data[:])\n",
    "        \n",
    "        ### step 2: Constraint testing ###\n",
    "        qc.append(constraint_testing(data_qubits, C_max), qr_data[:] + qr_f[:])\n",
    "        \n",
    "        ### step 3: penalty dephasing ###\n",
    "        qc.append(penalty_dephasing(data_qubits, alpha, gamma), qr_data[:] + qr_f[:])\n",
    "        \n",
    "        ### step 4: reinitialization ###\n",
    "        qc.append(reinitialization(index_qubits, data_qubits, C1, C2, C_max), qr_index[:] + qr_data[:] + qr_f[:])\n",
    "        \n",
    "        ### mixing operator ###\n",
    "        qc.append(mixing_operator(index_qubits, beta), qr_index)\n",
    "\n",
    "    ### measure the index ###\n",
    "    ### since the default measurement outcome is shown in big endian, it is necessary to reverse the classical bits in order to unify the endian ###\n",
    "    qc.measure(qr_index, cr_index[::-1])\n",
    "    \n",
    "    return qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b04bdb5",
   "metadata": {},
   "source": [
    "Validation function contains four input instances.\n",
    "The output should pass the precision threshold 0.80 for the eight inputs before scored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "317d7ef1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute your circuit with following prepare_ex4c() function.\n",
    "# The prepare_ex4c() function works like the execute() function with only QuantumCircuit as an argument.\n",
    "from qc_grader import prepare_ex4c\n",
    "job = prepare_ex4c(solver_function)\n",
    "\n",
    "result = job.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d18349c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex4c\n",
    "grade_ex4c(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c64d2dee",
   "metadata": {},
   "source": [
    "### References\n",
    "1. Edward Farhi and Jeffrey Goldstone and Sam Gutmann (2014). A Quantum Approximate Optimization Algorithm. (https://arxiv.org/abs/1411.4028)\n",
    "2. Grand'rive, Pierre & Hullo, Jean-Francois (2019). Knapsack Problem variants of QAOA for battery revenue optimisation.  (https://arxiv.org/abs/1908.02210)\n",
    "3. V. Vedral, A. Barenco, A. Ekert (1995). Quantum Networks for Elementary Arithmetic Operations. (https://arxiv.org/abs/quant-ph/9511018)\n",
    "4. Steven A. Cuccaro, Thomas G. Draper, Samuel A. Kutin, David Petrie Moulton (2004). A new quantum ripple-carry addition circuit. (https://arxiv.org/abs/quant-ph/0410184)\n",
    "5. Thomas G. Draper (2000). Addition on a Quantum Computer (https://arxiv.org/abs/quant-ph/0008033)\n",
    "6. Lidia Ruiz-Perez, Juan Carlos Garcia-Escartin (2014). Quantum arithmetic with the Quantum Fourier Transform. (https://arxiv.org/abs/1411.5949)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2480f364",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "**Created by:** Bo Yang, Hyungseok Chang, Sitong Liu, Kifumi Numata\n",
    "\n",
    "**Version:** 1.0.1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Qiskit v0.31.0 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}