from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit import Aer
from qiskit.utils import algorithm_globals, QuantumInstance
from qiskit.algorithms import QAOA, NumPyMinimumEigensolver
import numpy as np

val = [5,6,7,8,9]
wt = [4,5,6,7,8]
W = 18

def dp(W, wt, val, n):
    k = [[0 for x in range(W + 1)] for x in range(n + 1)]
    for i in range(n + 1):
        for w in range(W + 1):
            if i == 0 or w == 0:
                k[i][w] = 0
            elif wt[i-1] <= w:
                k[i][w] = max(val[i-1] + k[i-1][w-wt[i-1]], k[i-1][w])
            else:
                k[i][w] = k[i-1][w]
                
    picks=[0 for x in range(n)]
    volume=W
    for i in range(n,-1,-1):
        if (k[i][volume]>k[i-1][volume]):
            picks[i-1]=1
            volume -= wt[i-1]
    return k[n][W],picks

n = len(val)
print("optimal value:", dp(W, wt, val, n)[0])
print('\n index of the chosen items:')
for i in range(n): 
    if dp(W, wt, val, n)[1][i]: 
        print(i,end=' ')

# import packages necessary for application classes.
from qiskit_optimization.applications import Knapsack

def knapsack_quadratic_program():
    # Put values, weights and max_weight parameter for the Knapsack()
    
    ##############################
    # Provide your code here
    prob = Knapsack(val, wt, W)
    
    #
    ##############################
    
    # to_quadratic_program generates a corresponding QuadraticProgram of the instance of the knapsack problem.
    kqp = prob.to_quadratic_program()
    return prob, kqp

prob,quadratic_program=knapsack_quadratic_program()
quadratic_program

# Numpy Eigensolver
meo = MinimumEigenOptimizer(min_eigen_solver=NumPyMinimumEigensolver())
result = meo.solve(quadratic_program)
print('result:\n', result)
print('\n index of the chosen items:', prob.interpret(result)) 

# QAOA
seed = 123
algorithm_globals.random_seed = seed
qins = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), shots=1000, seed_simulator=seed, seed_transpiler=seed)

meo = MinimumEigenOptimizer(min_eigen_solver=QAOA(reps=1, quantum_instance=qins))
result = meo.solve(quadratic_program)
print('result:\n', result)
print('\n index of the chosen items:', prob.interpret(result)) 

# Check your answer and submit using the following code
from qc_grader import grade_ex4a
grade_ex4a(quadratic_program)

L1 = [5,3,3,6,9,7,1]
L2 = [8,4,5,12,10,11,2]
C1 = [1,1,2,1,1,1,2]
C2 = [3,2,3,2,4,3,3]
C_max = 16

def knapsack_argument(L1, L2, C1, C2, C_max):
      
    ##############################
    # Provide your code here
    values = [j-i for i,j in zip(L1,L2)]
    weights = [j-i for i,j in zip(C1,C2)]
    max_weight = max([i+j for i,j in zip(L1, L2)])
    #
    ##############################
    return values, weights, max_weight
    
values, weights, max_weight = knapsack_argument(L1, L2, C1, C2, C_max)
print(values, weights, max_weight)
prob = Knapsack(values = values, weights = weights, max_weight = max_weight)
qp = prob.to_quadratic_program()
qp

# Check your answer and submit using the following code
from qc_grader import grade_ex4b
grade_ex4b(knapsack_argument)

# QAOA
seed = 123
algorithm_globals.random_seed = seed
qins = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), shots=1000, seed_simulator=seed, seed_transpiler=seed)

meo = MinimumEigenOptimizer(min_eigen_solver=QAOA(reps=1, quantum_instance=qins))
result = meo.solve(qp)
print('result:', result.x)

item = np.array(result.x)
revenue=0
for i in range(len(item)):
    if item[i]==0:
        revenue+=L1[i]
    else:
        revenue+=L2[i]

print('total revenue:', revenue)

instance_examples = [
    {
        'L1': [3, 7, 3, 4, 2, 6, 2, 2, 4, 6, 6],
        'L2': [7, 8, 7, 6, 6, 9, 6, 7, 6, 7, 7],
        'C1': [2, 2, 2, 3, 2, 4, 2, 2, 2, 2, 2],
        'C2': [4, 3, 3, 4, 4, 5, 3, 4, 4, 3, 4],
        'C_max': 33
    },
    {
        'L1': [4, 2, 2, 3, 5, 3, 6, 3, 8, 3, 2],
        'L2': [6, 5, 8, 5, 6, 6, 9, 7, 9, 5, 8],
        'C1': [3, 3, 2, 3, 4, 2, 2, 3, 4, 2, 2],
        'C2': [4, 4, 3, 5, 5, 3, 4, 5, 5, 3, 5],
        'C_max': 38
    },
    {
        'L1': [5, 4, 3, 3, 3, 7, 6, 4, 3, 5, 3],
        'L2': [9, 7, 5, 5, 7, 8, 8, 7, 5, 7, 9],
        'C1': [2, 2, 4, 2, 3, 4, 2, 2, 2, 2, 2],
        'C2': [3, 4, 5, 4, 4, 5, 3, 3, 5, 3, 5],
        'C_max': 35
    }
]

from typing import List, Union
import math
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, assemble
from qiskit.compiler import transpile
from qiskit.circuit import Gate
from qiskit.circuit.library.standard_gates import *
from qiskit.circuit.library import QFT

def phase_return(index_qubits: int, gamma: float, L1: list, L2: list, to_gate=True) -> Union[Gate, QuantumCircuit]:
    qr_index = QuantumRegister(index_qubits, "index")
    qc = QuantumCircuit(qr_index)
    
    ##############################
    ### U_1(gamma * (lambda2 - lambda1)) for each qubit ###
    # Provide your code here
     


    ##############################
    
    return qc.to_gate(label=" phase return ") if to_gate else qc

def subroutine_add_const(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    qc = QuantumCircuit(data_qubits)
    ##############################
    ### Phase Rotation ###
    # Provide your code here
    
    
    
    ##############################
    return qc.to_gate(label=" [+"+str(const)+"] ") if to_gate else qc


def const_adder(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_data)
    
    ##############################
    ### QFT ###
    # Provide your code here
    
    
    
    ##############################

    ##############################
    ### Phase Rotation ###
    # Use `subroutine_add_const`
    
    
    ##############################
    
    ##############################
    ### IQFT ###
    # Provide your code here
    
    
    
    ##############################
    return qc.to_gate(label=" [ +" + str(const) + "] ") if to_gate else qc

def cost_calculation(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_index = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qc = QuantumCircuit(qr_index, qr_data)

    for i, (val1, val2) in enumerate(zip(list1, list2)):
        
        ##############################
        ### Add val2 using const_adder controlled by i-th index register (set to 1) ###
        # Provide your code here
        
        
        
        ##############################
        qc.x(qr_index[i])
        
        ##############################
        ### Add val1 using const_adder controlled by i-th index register (set to 0) ###
        # Provide your code here
        
        
        
        ##############################
        qc.x(qr_index[i])
    
    return qc.to_gate(label=" Cost Calculation ") if to_gate else qc

def constraint_testing(data_qubits: int, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qr_f = QuantumRegister(1, "flag")
    qc = QuantumCircuit(qr_data, qr_f)
    
    ##############################
    ### Set the flag register for indices with costs larger than C_max ###
    # Provide your code here
    
    
    
    ##############################
    
    return qc.to_gate(label=" Constraint Testing ") if to_gate else qc

def penalty_dephasing(data_qubits: int, alpha: float, gamma: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qr_f = QuantumRegister(1, "flag")
    qc = QuantumCircuit(qr_data, qr_f)
    
    ##############################
    ### Phase Rotation ###
    # Provide your code here
    
    
    
    ##############################
    
    return qc.to_gate(label=" Penalty Dephasing ") if to_gate else qc

def reinitialization(index_qubits: int, data_qubits: int, C1: list, C2: list, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_index = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qr_f = QuantumRegister(1, "flag")
    qc = QuantumCircuit(qr_index, qr_data, qr_f)
    
    ##############################
    ### Reinitialization Circuit ###
    # Provide your code here
    
    
    
    ##############################
    
    return qc.to_gate(label=" Reinitialization ") if to_gate else qc

def mixing_operator(index_qubits: int, beta: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_index = QuantumRegister(index_qubits, "index")
    qc = QuantumCircuit(qr_index)
    
    ##############################
    ### Mixing Operator ###
    # Provide your code here
    
    
    
    ##############################
    
    return qc.to_gate(label=" Mixing Operator ") if to_gate else qc

def solver_function(L1: list, L2: list, C1: list, C2: list, C_max: int) -> QuantumCircuit:
    
    # the number of qubits representing answers
    index_qubits = len(L1)
    
    # the maximum possible total cost
    max_c = sum([max(l0, l1) for l0, l1 in zip(C1, C2)])
    
    # the number of qubits representing data values can be defined using the maximum possible total cost as follows:
    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
    
    ### Phase Operator ###
    # return part
    def phase_return():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
    
    # penalty part
    def subroutine_add_const():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################

    # penalty part
    def const_adder():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
    
    # penalty part
    def cost_calculation():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
    
    # penalty part
    def constraint_testing():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
    
    # penalty part
    def penalty_dephasing():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
        
    # penalty part
    def reinitialization():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################

    ### Mixing Operator ###
    def mixing_operator():
        ##############################
        ### TODO ###
        ### Paste your code from above cells here ###
        
        
        
        ##############################
    
    qr_index = QuantumRegister(index_qubits, "index") # index register
    qr_data = QuantumRegister(data_qubits, "data") # data register
    qr_f = QuantumRegister(1, "flag") # flag register
    cr_index = ClassicalRegister(index_qubits, "c_index") # classical register storing the measurement result of index register
    qc = QuantumCircuit(qr_index, qr_data, qr_f, cr_index)
    
    ### initialize the index register with uniform superposition state ###
    qc.h(qr_index)
    
    ### DO NOT CHANGE THE CODE BELOW
    p = 5
    alpha = 1
    for i in range(p):
        
        ### set fixed parameters for each round ###
        beta = 1 - (i + 1) / p
        gamma = (i + 1) / p
        
        ### return part ###
        qc.append(phase_return(index_qubits, gamma, L1, L2), qr_index)
        
        ### step 1: cost calculation ###
        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2), qr_index[:] + qr_data[:])
        
        ### step 2: Constraint testing ###
        qc.append(constraint_testing(data_qubits, C_max), qr_data[:] + qr_f[:])
        
        ### step 3: penalty dephasing ###
        qc.append(penalty_dephasing(data_qubits, alpha, gamma), qr_data[:] + qr_f[:])
        
        ### step 4: reinitialization ###
        qc.append(reinitialization(index_qubits, data_qubits, C1, C2, C_max), qr_index[:] + qr_data[:] + qr_f[:])
        
        ### mixing operator ###
        qc.append(mixing_operator(index_qubits, beta), qr_index)

    ### measure the index ###
    ### since the default measurement outcome is shown in big endian, it is necessary to reverse the classical bits in order to unify the endian ###
    qc.measure(qr_index, cr_index[::-1])
    
    return qc

# Execute your circuit with following prepare_ex4c() function.
# The prepare_ex4c() function works like the execute() function with only QuantumCircuit as an argument.
from qc_grader import prepare_ex4c
job = prepare_ex4c(solver_function)

result = job.result()

# Check your answer and submit using the following code
from qc_grader import grade_ex4c
grade_ex4c(job)