# useful additional packages 
import matplotlib.pyplot as plt
import matplotlib.axes as axes
import numpy as np
import networkx as nx
from qiskit import Aer, execute, QuantumCircuit
from qiskit.quantum_info import Statevector

# auxilliary function to plot graphs
def plot_result(G, x):
    colors = ['r' if x[i] == 0 else 'b' for i in range(n)]
    pos, default_axes = nx.spring_layout(G), plt.axes(frameon=True)
    nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, pos=pos)

# create graph
G = nx.Graph()

# add nodes
n = 5
G.add_nodes_from(range(n))

# add edges: tuple is (i,j,weight) where (i,j) is the edge
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0), (2, 4, 1.0), (3, 4, 1.0)]
G.add_weighted_edges_from(edges)

# plot graph
plot_result(G, [0]*n)

from docplex.mp.model import Model

mdl = Model('MaxCut')
x = mdl.binary_var_list('x{}'.format(i) for i in range(n))

objective = mdl.sum([ w * (x[i] + x[j] - 2*x[i]*x[j]) for (i, j, w) in edges])
mdl.maximize(objective)

mdl.prettyprint()

b = 2
mdl.add_constraint(mdl.sum(x) == b)
mdl.prettyprint()

conda install -c ibmdecisionoptimization cplex

from qiskit.optimization import QuadraticProgram
from qiskit.optimization.algorithms import CplexOptimizer

# convert from DOcplex model to Qiskit Quadratic Program
qp = QuadraticProgram()
qp.from_docplex(mdl)

# Solve Quadratic Program using CPLEX
cplex = CplexOptimizer()
result = cplex.solve(qp)
print(result)
plot_result(G, result.x)

from qiskit.optimization.converters import (
    InequalityToEquality,     # converts inequality constraints to equality constraints by adding slack variables
    LinearEqualityToPenalty,  # converts linear equality constraints to quadratic penalty terms 
    IntegerToBinary,          # converts integer variables to binary variables
    QuadraticProgramToQubo    # combines the previous three converters
)

lineq2penalty = LinearEqualityToPenalty(penalty=1)
qp_wo_constr = lineq2penalty.convert(qp)
qp_wo_constr

# Solve converted Quadratic Program using CPLEX
result = cplex.solve(qp_wo_constr)
print(result)
plot_result(G, result.x)

H, offset = qp_wo_constr.to_ising()
print('offset =', offset)
print()
print('H =', H)

# print Ising Hamiltonian as matrix
H_matrix = np.real(H.to_matrix())
print('dim(H):', H_matrix.shape)
print(H_matrix)

# plot diagonal of matrix
opt_indices = list(np.where(H_matrix.diagonal() == min(H_matrix.diagonal())))[0]
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), H_matrix.diagonal())
plt.bar(opt_indices, H_matrix.diagonal()[opt_indices], color='g')
plt.xticks(range(2**n), ['('+str(i)+') {0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14)
plt.show()

from qiskit.circuit.library import RealAmplitudes
qc = RealAmplitudes(5, reps=1)
qc.draw(fold=120)

# run VQE
from qiskit.aqua.algorithms import VQE
vqe = VQE(H, qc, quantum_instance=Aer.get_backend('statevector_simulator'))
result = vqe.run()
print('optimal value:', np.round(result.eigenvalue, decimals=4))

# plot probabilities
probabilities = np.abs(result.eigenstate)**2
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), probabilities)
plt.bar(opt_indices, probabilities[opt_indices], color='g')
plt.xticks(range(2**n), ['('+str(i)+') {0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14)
plt.show()

from qiskit.circuit import Parameter
gamma, beta = Parameter('gamma'), Parameter('beta')

qc = QuantumCircuit(1)
qc.h(0)
qc.barrier()
qc.rz(gamma, 0)
qc.rx(beta, 0)
qc.barrier()
qc.draw()

def construct_schedule(T, N):
    delta_t = T/N
    gammas, betas = [], []  # H_C, H_X parameters
    for i in range(N+1):
        t = i * delta_t
        gammas += [ 2 * delta_t * t/T ]  # H_C
        betas += [ -2 * delta_t * (1 - t/T) ]  # H_X
    return gammas, betas

T = 5
N = 10
gammas, betas = construct_schedule(T, N)

plt.figure(figsize=(10, 4))
plt.plot(np.linspace(0, T, N+1), gammas, label='gammas')
plt.plot(np.linspace(0, T, N+1), betas, label='betas')
plt.legend(fontsize=14)
plt.xticks(fontsize=14)
plt.xlabel('time', fontsize=14)
plt.yticks(fontsize=14)
plt.ylabel('parameters', fontsize=14);

# track probabilities during trotterized annealing
probabilities = np.zeros((2, N+1))

# construct circuit
qc = QuantumCircuit(1)
qc.h(0)
qc.barrier()
for i, (gamma, beta) in enumerate(zip(gammas, betas)):
    qc.rz(gamma, 0)
    qc.rx(beta, 0)
    qc.barrier()
    probabilities[:, i] = Statevector.from_instruction(qc).probabilities()
qc.draw()

plt.figure(figsize=(12, 7))
plt.plot(np.linspace(0, T, N+1), probabilities[1, :], 'gd-', label=r'$|1\rangle$')
plt.plot(np.linspace(0, T, N+1), probabilities[0, :], 'bo-', label=r'$|0\rangle$')
plt.legend(fontsize=14)
plt.xticks(fontsize=14)
plt.xlabel('time $t$', fontsize=14)
plt.yticks(fontsize=14)
plt.ylabel('probabilities', fontsize=14);

plt.figure(figsize=(12, 7))
plt.plot(np.linspace(0, T, N+1), probabilities[0, :] - probabilities[1, :], 'gd-')
plt.xticks(fontsize=14)
plt.xlabel('time $t$', fontsize=14)
plt.yticks(fontsize=14)
plt.ylabel('objective value', fontsize=14);

gamma = Parameter('gamma')
qc1 = QuantumCircuit(2)
qc1.cx(0, 1)
qc1.rz(gamma, 1)
qc1.cx(0, 1)
qc1.draw()

qc2 = QuantumCircuit(2)
qc2.rzz(gamma, 0, 1)
qc2.draw()

from qiskit.aqua.algorithms.minimum_eigen_solvers.qaoa.var_form import QAOAVarForm

# construct parameters from annealing schedule
T = 10
N = 20
gammas, betas = construct_schedule(T, N)

# construct variational form
var_form = QAOAVarForm(H, N+1)

# evaluate circuit
qc = var_form.construct_circuit(gammas + betas)
sv = Statevector.from_instruction(qc)

# plot probabilities
probabilities = sv.probabilities()
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), probabilities)
plt.bar(opt_indices, probabilities[opt_indices], color='g')
plt.xticks(range(2**n), ['('+str(i)+') {0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14);

# draw trotterized annealing circuit for MaxCut
qc.draw()

from qiskit.aqua.components.optimizers import COBYLA

# QAOA circuit for p = 1
gamma, beta = Parameter('gamma'), Parameter('beta')
qc = QuantumCircuit(1)
qc.h(0)
qc.rz(gamma, 0)
qc.rx(beta, 0)

def objective(params):
    qc_ = qc.bind_parameters({gamma: params[0], beta: params[1]})
    probs = Statevector.from_instruction(qc_).probabilities()
    return probs @ [1, -1]

optimizer = COBYLA()
result = optimizer.optimize(num_vars=2, objective_function=objective, initial_point=[0.5, 0.5])
print('optimal params:', np.round(result[0], decimals=4))
print('optimal value: ', np.round(result[1], decimals=4))
print('optimal state: ', np.round(Statevector.from_instruction(qc.bind_parameters(
    {gamma: result[0][0], beta: result[0][1]})).data, decimals=4))

# illustrating the QAOA var form
p = 1
var_form = QAOAVarForm(H, p)
var_form.construct_circuit([gamma, beta]).draw()

from qiskit.aqua.algorithms import QAOA

qaoa_mes = QAOA(H, p=1, optimizer=optimizer, quantum_instance=Aer.get_backend('statevector_simulator'))
result = qaoa_mes.run()

print('optimal params:      ', result.optimal_parameters)
print('optimal value:       ', result.optimal_value)
print('optimal probability: ', sum(np.abs(result.eigenstate[opt_indices])**2))

# plot probabilities
probabilities = np.abs(result['eigenstate'])**2
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), probabilities)
plt.bar(opt_indices, probabilities[opt_indices], color='g')
plt.xticks(range(2**n), ['('+str(i)+') {0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14)
plt.show()

from qiskit.optimization.algorithms import MinimumEigenOptimizer

# construct QAOA as Minimum Eigensolver
qaoa_mes = QAOA(p=1, optimizer=optimizer, quantum_instance=Aer.get_backend('statevector_simulator'))

# construct Minimum Eigen Optimizer based on QAOA
qaoa = MinimumEigenOptimizer(qaoa_mes)

# solve Quadratic Program
result = qaoa.solve(qp)

print(result)
plot_result(G, result.x)

# converts two's complement bit string to corresponding integer
def twos_complement(val, num_bits):
    val = int(val, 2)
    if (val & (1 << (num_bits - 1))) != 0:
        val = val - (1 << num_bits)     
    return val   

print(twos_complement('0011', 4))
print(twos_complement('1011', 4))

from qiskit.circuit.library import QFT

def encode(num_qubits, k):
    qc = QuantumCircuit(num_qubits, name='enc({})'.format(k))
    for i in range(num_qubits):
        theta = 2*np.pi * 2**i / 2**num_qubits * k
        qc.rz(theta, i)
    return qc

num_value_qubits = 4
qc = QuantumCircuit(num_value_qubits, num_value_qubits)
qc.h(range(num_value_qubits))
qc.barrier()
qc.append(encode(num_value_qubits, 2), range(num_value_qubits))
qc.barrier()
qc.append(QFT(num_value_qubits, do_swaps=False).inverse(), qc.qubits)
qc.measure(qc.qregs[0], qc.cregs[0][::-1])
qc.draw(fold=120)

num_value_qubits = 4
qc = QuantumCircuit(num_value_qubits, num_value_qubits)
qc.h(range(num_value_qubits))
qc.barrier()
qc.extend(encode(num_value_qubits, 2))
# qc.extend(encode(num_value_qubits, -3))
# qc.extend(encode(num_value_qubits, 4))
qc.barrier()
qc.append(QFT(num_value_qubits, do_swaps=False, inverse=True), qc.qubits)
qc.measure(qc.qregs[0], qc.cregs[0][::-1])
qc.draw(fold=120)

counts = execute(qc, Aer.get_backend('qasm_simulator')).result().get_counts()
for key in counts:
    print(key, ' -->', twos_complement(key, num_value_qubits))

num_input_qubits = 3
num_total_qubits = num_input_qubits + num_value_qubits
qc = QuantumCircuit(num_total_qubits)
qc.h([0, 1, 2])
qc.barrier()
qc.h(range(num_input_qubits, num_total_qubits))
qc.append(encode(num_value_qubits,  2).control(2), [0, 2] + list(range(num_input_qubits, num_total_qubits)))
qc.append(encode(num_value_qubits, -3).control(), [1] + list(range(num_input_qubits, num_total_qubits)))
qc.append(encode(num_value_qubits,  4).control(), [2] + list(range(num_input_qubits, num_total_qubits)))
qc.append(QFT(num_value_qubits, do_swaps=False).inverse(), range(num_input_qubits, num_total_qubits))
qc.measure_all()
qc.draw(fold=120)

counts = execute(qc, Aer.get_backend('qasm_simulator')).result().get_counts()
for key, value in counts.items():
    x = key[num_value_qubits:]
    y_bin = key[:num_value_qubits][::-1]
    y_int = twos_complement(y_bin, num_value_qubits)
    print(x, '-->', y_bin, '-->', y_int, '\t(counts: {})'.format(value))

from qiskit.circuit.library import QuadraticForm

# get quadratic / linear / constant part of quadratic program
A = qp_wo_constr.objective.quadratic.to_array()
b = qp_wo_constr.objective.linear.to_array()
c = qp_wo_constr.objective.constant

# set number of results qubits
num_value_qubits = 5

# construct circuit to evaluate quadratic form
qf = QuadraticForm(num_value_qubits, A, b, c)
qf.draw(fold=115)

qc = QuantumCircuit(n + num_value_qubits)
qc.h(range(n))
qc.append(qf, range(n + num_value_qubits))
qc.measure_all()
qc.draw()

counts = execute(qc, backend=Aer.get_backend('qasm_simulator')).result().get_counts()
for key, value in counts.items():
    x_ = key[num_value_qubits:]
    x = [0 if x__ == '0' else 1 for x__ in x_][::-1]
    y_bin = key[:num_value_qubits]
    y_int = twos_complement(y_bin, num_value_qubits)
    qx = qp_wo_constr.objective.evaluate(x)    
    print('x =', x_, '\ty_bin =', y_bin, '\ty_int =', y_int, '\tQ(x) =', qx, '\t(counts: {})'.format(value))

qc = QuantumCircuit(n + num_value_qubits, name='U_f')
qc.append(qf, range(n + num_value_qubits))            # 1. compute Q(x)
qc.z(qc.qubits[-1])                                   # 2. multiply all |x>|Q(x)> by -1 where Q(x) < 0.
qc.append(qf.inverse(), range(n + num_value_qubits))  # 3. uncompute Q(x).
qc.draw()

qc_grover = QuantumCircuit(n + num_value_qubits)
qc_grover.h(range(n))
qc_grover.append(qc, range(qc_grover.num_qubits))
qc_grover.draw()

data = Statevector.from_instruction(qc_grover).data
x = ['{0:05b}'.format(i) for i in range(2**n)]
y = [qp_wo_constr.objective.evaluate([0 if x__ == '0' else 1 for x__ in reversed(x_)]) for x_ in x]
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), np.real(data[:2**n]))
plt.xticks(range(2**n), ['{} $\leftarrow$ '.format(y[i]) + '{0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14)
plt.show()

reflection = QuantumCircuit(n, name='reflection')
reflection.h(range(reflection.num_qubits))
reflection.barrier()
reflection.x(range(reflection.num_qubits))
reflection.barrier()
reflection.h(-1)
reflection.mct(list(range(reflection.num_qubits - 1)), -1)
reflection.h(-1)
reflection.barrier()
reflection.x(range(reflection.num_qubits))
reflection.barrier()
reflection.h(range(reflection.num_qubits))
reflection.draw()

qc_grover.append(reflection, range(n))
qc_grover.draw()

data = Statevector.from_instruction(qc_grover).data
x = ['{0:05b}'.format(i) for i in range(2**n)]
y = [qp_wo_constr.objective.evaluate([0 if x__ == '0' else 1 for x__ in reversed(x_)]) for x_ in x]
plt.figure(figsize=(12, 5))
plt.bar(range(2**n), -np.real(data[:2**n]))  # multiply by -1, since reflection is implemented up to global phase -1
plt.xticks(range(2**n), ['{} $\leftarrow$ '.format(y[i]) + '{0:05b}'.format(i) for i in range(2**n)], rotation=90, fontsize=14)
plt.yticks(fontsize=14)
plt.show()

from qiskit.optimization.algorithms import GroverOptimizer

# set up Grover Optimizer
grover = GroverOptimizer(num_value_qubits=5, quantum_instance=Aer.get_backend('statevector_simulator'))
grover._qubo_converter.penalty = 1  # set to small value to reduce required number of value qubits

# solver problem
result = grover.solve(qp)

print(result)
plot_result(G, result.x)

from qiskit.finance.applications.ising import portfolio

# define tickers and number of assets
n = 6

# set budget
budget = n // 2

# sample expected returns and covariance matrix
mu, sigma = portfolio.random_model(n, seed=4321)

# plot mu/sigma
plt.figure(figsize=(12, 5))
plt.subplot(1,2,1)
plt.title('Expected Returns', fontsize=14)
plt.bar(range(n), mu)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.subplot(1,2,2)
plt.title('Covariance Matrix', fontsize=14)
plt.imshow(sigma, interpolation='nearest')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.show()

def get_qp(n, mu, sigma, budget, q):
    
    # create model with docplex
    mdl = Model('portfolio_optimization')

    # define variables: 
    #
    #     x_i in {0, 1}
    #
    x = mdl.binary_var_list(['x{}'.format(i) for i in range(n)])

    # set objective function: 
    #
    #     maximize { mu^T * x - q * x^T * sigma * x }
    #
    objective = mdl.sum([mu[i] * x[i] for i in range(n)])  # mu^T * x
    objective -= q * mdl.sum([sigma[i][j]*x[i]*x[j] for i in range(n) for j in range(n)])
    mdl.maximize(objective)

    # add budget constraint: 
    #
    #     1^T * x == budget
    #
    cost = mdl.sum([x[i] for i in range(n)])
    mdl.add_constraint(cost == budget, ctname='budget')

    # convert to quadratic program
    qp = QuadraticProgram()
    qp.from_docplex(mdl)
    
    return qp

# track resulting risk and return
mus = []
sigmas = []

# screen different risk factors
for q in [0, 1, 10]:
    
    # solve the portfolio optimization problem using QAOA
    qp = get_qp(n, mu, sigma, budget, q)
    result = qaoa.solve(qp)
    
    # track risk and return
    mus += [portfolio.portfolio_expected_value(result.x, mu)]
    sigmas += [portfolio.portfolio_variance(result.x, sigma)]
    
    # print resulting portfolio
    print(result.x)

# plot risk/return efficient frontier
plt.figure(figsize=(8, 5))
plt.plot(sigmas, mus, '--', linewidth=2)
plt.plot(sigmas, mus, 'rD', markersize=7)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.xlabel('Risk (Variance)', fontsize=14)
plt.ylabel('Return', fontsize=14)
plt.title('Efficient Frontier', fontsize=14)
plt.show()

