from qiskit.finance.data_providers import RandomDataProvider

import numpy as np
import matplotlib.pyplot as plt
import datetime
import pandas as pd
from nsepy import get_history

from docplex.mp.model import Model
from qiskit.optimization import QuadraticProgram
from qiskit.optimization.converters import LinearEqualityToPenalty

import pennylane as qml

assets = ['INFY', 'TCS', 'DRREDDY']

num_assets = len(assets)

print('The total number of assets {0}'.format(num_assets))

stockStartDate = datetime.date(2019,8,20)
stockEndDate = datetime.date(2020,8,20)

df = pd.DataFrame()

for stock in assets:
    df[stock] = get_history(symbol=stock, start=stockStartDate, end=stockEndDate)['Close']

title = 'Portfolio Adj. Close Price History    '
#Get the stocks
my_stocks = df
#Create and plot the graph
plt.figure(figsize=(12.2,4.5)) #width = 12.2in, height = 4.5
# Loop through each stock and plot the Adj Close for each day
for c in my_stocks.columns.values:
    plt.plot( my_stocks[c],  label=c)#plt.plot( X-Axis , Y-Axis, line_width, alpha_for_blending,  label)
plt.title(title)
plt.xlabel('Date',fontsize=18)
plt.ylabel('Adj. Price Rs.',fontsize=18)
plt.legend(my_stocks.columns.values, loc='upper left')
plt.show()

returns = df.pct_change()
cov_matrix_annual = returns.cov() * 249

mu = returns.mean().to_numpy()
sigma = cov_matrix_annual.to_numpy()

q = 0.5                   # set risk factor
budget = 2       # set budget
penalty = num_assets      # set parameter to scale the budget penalty term

mdl = Model('docplex model')
x = mdl.binary_var_list(num_assets)

# set objective function: 
#
#     maximize { mu^T * x - q * x^T * sigma * x }
#

objective = mdl.sum([mu[i] * x[i] for i in range(num_assets)])  # mu^T * x
objective -= q * mdl.sum([sigma[i][j]*x[i]*x[j] for i in range(num_assets) for j in range(num_assets)])
mdl.maximize(objective)

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

# converting to Quadratic Program
mod = QuadraticProgram()        
mod.from_docplex(mdl)

#removing the constraint to create the QUBO
lineq2penalty = LinearEqualityToPenalty(penalty) 
qubo = lineq2penalty.convert(mod)

#converting QUBO to an Ising Hamiltonian
H, offset = qubo.to_ising()
H.num_qubits

H = H.to_legacy_op()

H.print_details()

def anstanz(theta):
    
    qml.RX(theta[0], wires=0)
    qml.RX(theta[1], wires=1)
    qml.RX(theta[2], wires=2)
    qml.CNOT(wires=[0,1])
    qml.CNOT(wires=[0,2])
    qml.CNOT(wires=[1,2])
    qml.RX(theta[3], wires=0)
    qml.RX(theta[4], wires=1)
    qml.RX(theta[5], wires=2)
    qml.CNOT(wires=[0,1])
    qml.CNOT(wires=[0,2])
    qml.CNOT(wires=[1,2])
    qml.RX(theta[6], wires=0)
    qml.RX(theta[7], wires=1)
    qml.RX(theta[8], wires=2)
    
    return anstanz

dev1 = qml.device("default.qubit", wires=3)

@qml.qnode(dev1)
def circuit_IIZ(params):
    anstanz(params)
    return qml.expval(qml.Identity(0)@qml.Identity(1)@qml.PauliZ(2))

@qml.qnode(dev1)
def circuit_IZI(params):
    anstanz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.Identity(2))

@qml.qnode(dev1)
def circuit_ZII(params):
    anstanz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.Identity(2))

@qml.qnode(dev1)
def circuit_IZZ(params):
    anstanz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.PauliZ(2))

@qml.qnode(dev1)
def circuit_ZIZ(params):
    anstanz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.PauliZ(2))

@qml.qnode(dev1)
def circuit_ZZI(params):
    anstanz(params)
    return qml.expval(qml.PauliZ(0)@qml.PauliZ(1)@qml.Identity(2))


def pauli_operator_to_dict(pauli_operator):
    
    d = pauli_operator.to_dict()
    paulis = d['paulis']
    paulis_dict = {}

    for x in paulis:
        label = x['label']
        coeff = x['coeff']['real']
        paulis_dict[label] = coeff

    return paulis_dict

pauli_dict = pauli_operator_to_dict(H)

def vqe(parameters):
    
    expval_ZII = pauli_dict['ZII'] * circuit_ZII(parameters)  
    expval_IIZ = pauli_dict['IIZ'] * circuit_IIZ(parameters)
    expval_IZI = pauli_dict['IZI'] * circuit_IZI(parameters)
    expval_ZZI = pauli_dict['ZZI'] * circuit_ZZI(parameters)
    expval_ZIZ = pauli_dict['ZIZ'] * circuit_ZIZ(parameters)
    expval_IZZ = pauli_dict['IZZ'] * circuit_IZZ(parameters)
    
    # summing the measurement results
    classical_adder = expval_ZII + expval_IIZ + expval_ZIZ +  expval_ZZI + expval_ZIZ + expval_IZZ
    
    return classical_adder

opt = qml.GradientDescentOptimizer(stepsize=0.5)
value = []
    
params = [0.01, 0.02, 0.01, np.pi, 0.1, np.pi/2, np.pi, np.pi/2, np.pi] ## Random initial parameters

steps = 100
for i in range(steps):
    #print(i)
    params = opt.step(vqe, params)
    value.append(vqe(params))


@qml.qnode(dev1)
def final_circ(params):
    anstanz(params)
    return qml.probs(wires=[0,1,2])

plt.bar(['000', '001', '010', '011', '100', '101', '110', '111'], final_circ(params))

from qiskit.optimization.applications.ising.common import sample_most_likely

exact_eigensolver = NumPyMinimumEigensolver(H)
result = exact_eigensolver.run()
sample_most_likely(result.eigenstate)

 

