import warnings
#warnings.filterwarnings('ignore', 'DeprecationWarning')
warnings.filterwarnings('once')

#%pip uninstall qiskit

#%pip install qiskit #==0.16

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

# useful additional packages 
import matplotlib.pyplot as plt
import matplotlib.axes as axes
%matplotlib inline
import numpy as np
import networkx as nx

from qiskit import Aer
from qiskit.tools.visualization import plot_histogram

#VGG todo 1: the equivalent run_algorithm and EnergyInput versions updates
#from qiskit.aqua.translators.ising import max_cut, tsp
#from qiskit.aqua import run_algorithm
#from qiskit.aqua.input import EnergyInput

#old v0.16# from qiskit.optimization.ising import max_cut, tsp
#old v0.16# from qiskit.optimization.ising.common import sample_most_likely
#older# from qiskit.optimization.ising import docplex

from qiskit.optimization.applications.ising import max_cut, tsp 
from qiskit.optimization.applications.ising.common import sample_most_likely
from qiskit.optimization.applications.ising import docplex

from qiskit.aqua.algorithms import VQE 
from qiskit.aqua.algorithms import NumPyEigensolver as ExactEigensolver  
from qiskit.aqua.components.optimizers import SPSA
#from qiskit.aqua.components.variational_forms import  RY #RealAmplitudes
from qiskit.circuit.library import RealAmplitudes as RY
from qiskit.aqua import QuantumInstance

# setup aqua logging
import logging
from qiskit.aqua import set_qiskit_aqua_logging
# set_qiskit_aqua_logging(logging.DEBUG)  # choose INFO, DEBUG to see the log

from qiskit import IBMQ
provider = IBMQ.load_account()

#VGG select the backend for coupling_map
try: 
    backend = provider.get_backend('ibmq_16_melbourne')#'ibmq_16_melbourne')#'ibmq_essex')
    #backend = provider.get_backend('ibmq_london')#'ibmq_16_melbourne')#'ibmq_essex')
    #backend = provider.get_backend('ibmq_5_yorktown')#'ibmq_london')#'ibmq_16_melbourne')#'ibmq_essex')
    
except:
    backend = Aer.get_backend('qasm_simulator') #VGG it was 'BasicAer.get_backend' ibmq_16_melbourne
    
coupling_map = backend.configuration().coupling_map
print(coupling_map)

#VGG we will generate a diferent coupling_map for the graph
#coupling_map=None

from typing import List, Tuple

seed = 19120623
np.random.seed(seed)

#VGG: function adopted from the Rigetti's MaxCutQAOA.ipynb
def generate_ising_graph(edges: List[Tuple[int, int]]) -> nx.Graph:
    graph = nx.from_edgelist(edges)
    weights:  np.ndarray =  np.random.rand(graph.number_of_edges()) #VGG the old [-1,1] range into [0,1]
    weights /= np.linalg.norm(weights)
    nx.set_edge_attributes(graph, {e: {'weight': w} for e, w in zip(graph.edges, weights)})
    return graph

if coupling_map != None:
    G=generate_ising_graph(coupling_map)
    n=G.number_of_nodes()
    print(n)

# Generating a graph if there were no coupling_map
if coupling_map== None:
    #define the edges / coupling_map
    
    #'ibmq_16_melbourne'
    elist=[[0, 1], [0, 14], [1, 0], [1, 2], [1, 13], [2, 1], [2, 3], [2, 12], 
           [3, 2], [3, 4], [3, 11], [4, 3], [4, 5], [4, 10], [5, 4], [5, 6], 
           [5, 9], [6, 5], [6, 8], [7, 8], [8, 6], [8, 7], [8, 9], [9, 5], [9, 8], 
           [9, 10], [10, 4], [10, 9], [10, 11], [11, 3], [11, 10], [11, 12], [12, 2], 
           [12, 11], [12, 13], [13, 1], [13, 12], [13, 14], [14, 0], [14, 13]]
   
    
    #elist=[[0,1],[0,2],[0,3],[1,2],[2,3],[0,4],[0,2], [4, 3],[1,5],[3,5]]
    elist=[[0,1],[0,2],[0,3],[1,2],[2,3],[0,4],[0,2], [4, 3]]

    G=generate_ising_graph(elist)
    n=G.number_of_nodes()
    
#other ways to define the graph   
    #n=5 # Number of nodes in graph
    #G=nx.Graph()
    #G.add_nodes_from(np.arange(0,n,1))
    
    #ewlist=[(0,1,1.),(0,2,.5),(0,3,0),(1,2,1.0),(0,3,1.0)]
    #G1 = nx.from_edgelist(elist)
    #G1.add_weighted_edges_from(ewlist)
    
    
coupling_map = backend.configuration().coupling_map


#Visulaize
print(G.number_of_nodes(),G.number_of_edges())
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)

nx.drawing.nx_pylab.draw(G)

elist=G.edges()
print("elist=",elist)

ewlist=[(i,j,G.get_edge_data(i,j,default=0)['weight']) for i,j in G.edges()]
print('ewlist=',ewlist)

def issymmetric(Matrix):
  dim=Matrix.shape[0]
  if Matrix.shape[1] != dim:
    print("Shape Error!")
    return False 
  for i in range(dim):
    for j in range(i,dim):
      if Matrix[i,j]!=Matrix[j,i]:
        print("Shape Error:",(i,j),Matrix[i,j],Matrix[j,i],"difference:",Matrix[i,j]-Matrix[j,i])
        return False
  return True
  

# Computing the weight matrix from the random graph
w = np.zeros([n,n])
w = np.eye(n)
for i in range(n):
    for j in range(n):
        temp = G.get_edge_data(i,j,default=0)
        if temp != 0:
            w[i,j] = temp['weight'] 

w/=np.linalg.det(w)**(1/n)
print("Symmetric:",issymmetric(w),"Norm:",np.linalg.norm(w))
print("Eignvlues:",np.linalg.eigvals(w),"det:",np.linalg.det(w))
print(w)

np.sum(w)/4 #the offset value

def Max_Cut_BF(W,*x0):
  best_cost_brute = 0
  xbest_brute=np.array([1]*n)

  for b in range(2**n):
      x = [int(t) for t in reversed(list(bin(b)[2:].zfill(n)))]
      cost = 0
      for h in x0: 
        cost -= np.dot(h,x)/n #VGG don't give free samples to those with h==1
      for i in range(n):
          cost +=(2-np.dot(x,x))/n/2 #VGG try to favor fewer free samples 
          for j in range(n):
              cost += W[i,j]*x[i]*(1-x[j])
      
              if np.isclose(cost,best_cost_brute):
                if sum(x)<sum(xbest_brute):
                  best_cost_brute = cost
                  xbest_brute = x
              else:
                if best_cost_brute < cost:
                  best_cost_brute = cost
                  xbest_brute = x
      if 1==2:
          print('case = ' + str(x)+ ' cost = ' + str(cost))
  return best_cost_brute, xbest_brute
   

%%time 
if n < 10: 
  best_cost_brute, xbest_brute = Max_Cut_BF(w)

  colors = ['r' if xbest_brute[i] == 0 else 'b' for i in range(n)]
  nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, pos=pos)
  print('\nBest solution = ' + str(xbest_brute) + ' cost = ' + str(best_cost_brute)) 

np.set_printoptions(precision=3)
print(w)

def market_simulations(m,*opt):

  free_samples=0
  boughten=0
  Mw2 = np.zeros([n,n])+w
  relations=np.zeros([n,n])
  x_free_samples=np.int_(np.zeros(n))
  np.set_printoptions(precision=2)

  if 'q' in opt: 
    print("Using Max_Cut option:",'q')
    print("submiting for results using:",backend)
  elif 'dcplx' in opt:
    print("Using Max_Cut option:",'Docplex')
  else:
    print("Using Max_Cut_BF")
    if n > 10 : 
      print("It may take too long to do Brute Force Calulations - skiping!")
      return
            
  print("day"," [free samples]"," [buyers distribution]"," [to be used as constrain]"," the two totals"," w-det")
  for i in range(m):
    
    if sum(x_free_samples)>n*2/3:
      #VGG In the future one can consider giving out to these who have not recived yet!

      x_free_samples = np.array([1]*n)-np.int_(x_free_samples!=0)
      
      if (x_free_samples==0).all:
        x=np.array([0]*n)
        xbest_brute=np.array([0]*n)
        tmp1=Mw2/np.linalg.norm(Mw2) #select only those rows that recived free samples
        tmp=sum(tmp1[:,]) #sum probaility contributions  
      else:
        x=x_free_samples
        free_samples+=sum(x)
        tmp1=Mw2[x==1]/np.linalg.norm(Mw2) #select only those rows that recived free samples
        tmp=sum(tmp1[:,],(np.array([1]*n)-x)) #sum probaility contributions  
        tmp-=np.array([1]*n) #push to negative those with free samples 

    else:
      
      if 'q' in opt:
        best_cost_brute, xbest_brute = Max_Cut_IBMQ(Mw2,x_free_samples)
      elif 'dcplx' in opt:
        best_cost_brute, xbest_brute = Max_Cut_Dcplx(Mw2,x_free_samples)
      else:
        best_cost_brute, xbest_brute = Max_Cut_BF(Mw2,x_free_samples)
      
      x=np.array(xbest_brute)
      x_free_samples+=x
      free_samples+=sum(x)
      tmp1=Mw2[x==1]/np.linalg.norm(Mw2) #select only those rows that recived free samples
      tmp=sum(tmp1[:,],(np.array([1]*n)-x)) #sum probaility contributions  
      tmp-=np.array([1]*n) #push to negative those with free samples 
    
    #print(tmp)
    ab=sum(tmp[tmp > 0])

    for j in range(n):
      test=np.random.uniform()*ab/2
      if tmp[j] > test: #buy the product
        x[j]+=1
        boughten+=1
    
    x0=np.array(xbest_brute)
    x-=x0
    relation_today=x.reshape(n,1) @ x0.reshape(1,n)
    relation_today+=relation_today.T
    relations+=relation_today
    
    #print(x0,x,"\n",relation_today)
    #print(x0,x,x_free_samples,tmp)
    print(i,x0,x,x_free_samples,free_samples, boughten,'{:6.4f}'.format(np.linalg.det(Mw2)))

    if i%4==0 : #weekely updates of the w matrix 
      Mw2+=(np.eye(n)+relations)/n/100 #update the w matrix
      relations=np.zeros([n,n])

      if issymmetric(Mw2) and np.linalg.det(Mw2)>0.:
        Mw2/=np.linalg.det(Mw2)**(1/n)
      else:
        Mw2/=np.linalg.norm(Mw2,ord='fro')

  
  print("\nlast day configuration record:\n")
  print(x0,tmp)
  print()
  print(x,free_samples, boughten, '{:6.4f}'.format(np.linalg.norm(Mw2)),'{:6.4f}'.format(np.linalg.det(Mw2)))
  print()
  print(Mw2)
  return

%time market_simulations(10)

#VGG qubitOp, offset = max_cut.get_max_cut_qubitops(w)
qubitOp, offset = max_cut.get_operator(w)
#algo_input = EnergyInput(qubitOp)
offset

#Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = ExactEigensolver(qubitOp, k=3)
result = ee.run()

print("energys:",result['eigenvalues'].real)

#VGG# x = max_cut.sample_most_likely(result['eigvecs'][0])
x = sample_most_likely(result['eigenstates'][0]) 

print('energy:', result['eigenvalues'][0],', offset:',offset)
print('max-cut objective:', result['eigenvalues'][0] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

#VGG note that the other runs had implemeneted soft constrains!

from docplex.mp.model import Model

#VGG from qiskit.aqua.translators.ising import docplex
#older# from qiskit.optimization.ising import docplex

# Create an instance of a model and variables.
mdl = Model(name='max_cut')
y = {i: mdl.binary_var(name='x_{0}'.format(i)) for i in range(n)}

# Object function
#VGG added y[i]/100 term to split the degenerate 1<->0 states in favor of less free samples

max_cut_func = mdl.sum(y[i]/50+w[i,j]* y[i] * ( 1 - y[j] ) for i in range(n) for j in range(n))
mdl.maximize(max_cut_func)

# No constraints for Max-Cut problems.

qubitOp_docplex, offset_docplex = docplex.get_operator(mdl)
offset_docplex

#VGG define the above as a function 
def set_up_Dcplx_model(W,*c_x0):
  mdl = Model(name='max_cut') 
  y = {i: mdl.binary_var(name='y_{0}'.format(i)) for i in range(n)}

  #VGG try to favor fewer free samples using (2-np.dot(x,x))/n/2 
  #VGG split the degenerate 1<->0 states in favor of less free samples using x[i]/n**2
  max_cut_func=mdl.sum((-1)*(2-y[i])*0.5+(-0)*y[i]/100 for i in range(n)) 
  
  #VGG don't give free samples to those with h==1
  max_cut_func+=mdl.sum(h[i]*y[i]*0.55 for i in range(n) for h in c_x0) 
  
  max_cut_func+=mdl.sum(W[i,j]* y[i] * ( 1 - y[j] ) for i in range(n) for j in range(n))
  
  mdl.maximize(max_cut_func)
  qubitOp, offset = docplex.get_operator(mdl)
  return qubitOp, offset
 
qubitOp_docplex, offset_docplex = set_up_Dcplx_model(w)
print(offset_docplex,x)

#Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = ExactEigensolver(qubitOp_docplex, k=3)
result = ee.run()

print("energys:",result['eigenvalues'].real)

x = sample_most_likely(result['eigenstates'][0])

print('energy:', result['eigenvalues'][0].real)
print('max-cut objective:', result['eigenvalues'][0].real + offset_docplex)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)
x=max_cut.get_graph_solution(x).tolist()

qubitOp_docplex, offset_docplex = set_up_Dcplx_model(w,x)
print(x,offset_docplex)

#VGG note if you keep executing this cell you can see diferent configurations 

#VGG define the above as a function 
def Max_Cut_Dcplx(W,*c_x0): 

  qubitOp_docplex, offset_docplex = set_up_Dcplx_model(W)
  
  for h in c_x0:
    qubitOp_docplex, offset_docplex = set_up_Dcplx_model(W,h)

  ee = ExactEigensolver(qubitOp_docplex,k=3)
  result = ee.run()
  x=sample_most_likely(result['eigenstates'][0])
  x_dcplx=max_cut.get_graph_solution(x).tolist()
  cost_dcplx=result['eigenvalues'][0].real
  #cost_dcplx=max_cut.max_cut_value(x, W)

  return cost_dcplx, x_dcplx

%time market_simulations(10,'dcplx')

warnings.filterwarnings('once')
model=qubitOp_docplex
#model=qubitOp

backend1 = Aer.get_backend('statevector_simulator')
backend2 = Aer.get_backend('qasm_simulator') #VGG it was 'BasicAer.get_backend' ibmq_16_melbourne


seed = 10598

spsa = SPSA(max_trials=10)
ry = RY(model.num_qubits, entanglement='linear') #depth=5, 
vqe = VQE(model, ry, spsa)

#VGG backend1 = Aer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend2, seed_simulator=seed, seed_transpiler=seed)
print(backend2)
result = vqe.run(quantum_instance)

#VGG# x = max_cut.sample_most_likely(result['eigvecs'][0])
x = sample_most_likely(result['eigenstate'])

print('energy:', result['eigenvalue'])
print('time:', result['optimizer_time'])
print('max-cut objective:', result['eigenvalue'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

x=max_cut.get_graph_solution(x).tolist()

qubitOp_docplex, offset_docplex = set_up_Dcplx_model(w,x)
print(x,offset_docplex)


#warnings.filterwarnings('ignore', 'DeprecationWarning')
warnings.filterwarnings('once')

# run quantum algorithm with shots
seed = 10598

spsa = SPSA(max_trials=30)
ry = RY(model.num_qubits,  entanglement='linear') #depth=5,
vqe = VQE(model, ry, spsa)

backend2 = Aer.get_backend('qasm_simulator') #VGG it was 'BasicAer.get_backend' ibmq_16_melbourne
quantum_instance = QuantumInstance(backend2, shots=1024, seed_simulator=seed, seed_transpiler=seed)

print(backend2)
result = vqe.run(quantum_instance)

"""declarative approach, update the param from the previous cell.
params['backend']['provider'] = 'qiskit.BasicAer'
params['backend']['name'] = 'qasm_simulator'
params['backend']['shots'] = 1024
result = run_algorithm(params, algo_input)
"""

#VGG# x = max_cut.sample_most_likely(result['eigvecs'][0])
x = sample_most_likely(result['eigenstate'])

print('energy:', result['eigenvalue'])
print('time:', result['optimizer_time'])
print('max-cut objective:', result['eigenvalue'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))
plot_histogram(result['eigenstate'])

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

x=max_cut.get_graph_solution(x).tolist()
qubitOp_docplex, offset_docplex = set_up_Dcplx_model(w,x)
print(x,offset_docplex)



plot_histogram(result['eigenstate'])

backend_old=backend
#backend=backend2

#warnings.filterwarnings('ignore', 'DeprecationWarning')
# run quantum algorithm with shots
seed = 10598

spsa = SPSA(max_trials=10) #VGG 300
ry = RY(model.num_qubits,  entanglement='linear') #depth=5,
vqe = VQE(model, ry, spsa)

#backend = provider.get_backend('ibmq_16_melbourne')#'ibmq_16_melbourne')#'ibmq_essex' # ibmq_london
#backend = provider.get_backend('ibmq_qasm_simulator')
quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed, seed_transpiler=seed)

print("submiting for results using:",backend)
result = vqe.run(quantum_instance)


#VGG# x = max_cut.sample_most_likely(result['eigvecs'][0])
x = sample_most_likely(result['eigenstate'])

print('energy:', result['eigenvalue'])
print('time:', result['optimizer_time'])
print('max-cut objective:', result['eigenvalue'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

plot_histogram(result['eigenstate'])

#VGG define the above as a function 

def Max_Cut_IBMQ(W,*c_x0): 

  qubitOp_docplex, offset_docplex = set_up_Dcplx_model(W)
  
  for h in c_x0:
    qubitOp_docplex, offset_docplex = set_up_Dcplx_model(W,h)

  model=qubitOp_docplex

  spsa = SPSA(max_trials=10) #VGG 300
  ry = RY(model.num_qubits,  entanglement='linear') #depth=5,
  vqe = VQE(model, ry, spsa)

  quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed, seed_transpiler=seed)

  result = vqe.run(quantum_instance)
  x = sample_most_likely(result['eigenstate'])
  cost_vqe=max_cut.max_cut_value(x, W)
  x_vqe =np.int_(max_cut.get_graph_solution(x)).tolist()

  return cost_vqe, x_vqe

%time market_simulations(10)

%time market_simulations(10,'dcplx')

print(backend)
backend = provider.get_backend('ibmq_16_melbourne')#'ibmq_16_melbourne')#'ibmq_essex' # ibmq_london
#backend = provider.get_backend('ibmq_qasm_simulator')
#backend = Aer.get_backend('qasm_simulator')
#backend = Aer.get_backend('statevector_simulator')
print(backend)

backend = provider.get_backend('ibmq_16_melbourne')#'ibmq_16_melbourne')#'ibmq_essex' # ibmq_london
#backend = provider.get_backend('ibmq_qasm_simulator')
#backend = Aer.get_backend('qasm_simulator')
#backend = Aer.get_backend('statevector_simulator')
%time market_simulations(10,'q')

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright