# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + id="RH4U56H_kpyQ"
"""
The following code simulates iVQAEC algorithm on QASM simulator.
iVQAEC algorithm is a variational quantum algorithm for equality constrained semidefinite programs (SDPs) that we proposed in our paper.
Here, we use iVQAEC to solve random instances of an equality constrained SDP.
In the paper, we report simulation results for N = 8, 16, and 32, where N is the dimension of input Hermitian operators of an SDP.
For writing clean, short and understandable code, here we fix N = 8. By simply substituting another value of N and with some other minor changes, 
one can obtain results for that N.
"""

"""
This code is based on qiskit-aqua package which recently got deprecated. We intend to use new packages. However, for now the following code is based on
qiskit-aqua package. To install it, use the following pip command:

!pip install git+https://github.com/Qiskit/qiskit-aqua.git

"""

# + id="rJp_NmzgqDBS"
# import packages.
import pennylane as qml
import numpy as np
import cvxpy as cp
import qiskit
import qiskit.providers.aer.noise as noise
import scipy as sc
import scipy.stats as stats
import scipy.sparse as sparse


# + id="2CfnpyWgrEKG"
def decompose_into_pauli_strings(M):
  """
  This function expresses an Hermitian operator as a linear combination of Pauli strings

  :param M: Hermitian operator
  """
  coeffs, pauli_string_objs = qml.utils.decompose_hamiltonian(M)
  return qml.Hamiltonian(coeffs, pauli_string_objs)


# + id="k3kbsOXzrMRO"
# modeling noise

# probabilities of error for single and two qubit gates
prob_1 = 0.001
prob_2 = 0.01

# Depolarizing errors
error_1 = noise.depolarizing_error(prob_1, 1)
error_2 = noise.depolarizing_error(prob_2, 2)

noise_model = noise.NoiseModel()
noise_model.add_all_qubit_quantum_error(error_2, ['cx'])
noise_model.add_all_qubit_quantum_error(error_1, ['u1', 'u2', 'u3'])


# + id="qAxwamZbUFQt"
def sparse_rand_sym(n, density):
  """
  This function generates a sparse symmetric matrix

  :param n: dimension of the matrix
  :param density: sparsity
  """
  rvs = stats.norm().rvs
  
  # generates a sparse matrix
  X = sparse.random(n, n, density=density, data_rvs=rvs)

  # make it symmetric
  upper_X = sparse.triu(X) 
  result = upper_X + upper_X.T - sparse.diags(X.diagonal())
  return result


# + id="IjcymwGTVSYn"
# constants of an SDP
N = 8
M = 3
R = 10

# + id="KX7HEKbTVbzv"
# generate a weakly constrained random sparse SDP and solve them using cvxpy
C = sparse_rand_sym(N, 0.1)
A = []
b = []
for i in range(M-1):
    A.append(sparse_rand_sym(N, 0.1))
    b.append(np.random.randn())
b.append(R)
A.append(np.eye(N)) 

# define the variable and constraints
X = cp.Variable((N,N), symmetric=True)
# positive semidefinite constraint
constraints = [X >> 0]

constraints += [
    cp.trace(A[i]@X) == b[i] for i in range(M)
]

# solve the problem
prob = cp.Problem(cp.Minimize(cp.trace(C@X)),
                  constraints)
prob.solve()

# result.
print("Optimal value:", prob.value)

# + id="Rpd3sGkfyPhg"
C = C.todense()
A[0] = A[0].todense()
A[1] = A[1].todense()

# + id="pr1Pkftgr8fW"
# Now solving the above SDP using iVQAEC
# quantum circuit settings
num_wires = 3
num_layers = 4

# quantum device
# can also define the number of shots here
device = qml.device("default.qubit", wires=num_wires)

# def the ansatz
def circuit(param):
  """
  This function instantiates quantum circuit based on templates provided by pennylane

  :param param: the parameters of the circuit
  """
  qml.templates.StronglyEntanglingLayers(param, wires=list(range(num_wires)))


# + id="3wbRjJxhtZeJ"
# def the cost function (Just the expectation value of the hamiltonian)
@qml.qnode(device)
def evalute_exp(param, M):
  """
  This function evaluates the expectation value of a Hermitian operator

  :param param: circuit parameters
  :param M: Hermitian operator
  """
  circuit(param)
  M_pauli = decompose_into_pauli_strings(M)
  return qml.expval(M_pauli)

def cost_func_ALM(param, **kwargs):
  """
  This function evaluates the cost function at given parameters

  :param param: circuit parameters
  """
  exp_C = evalute_exp(param, C)
  exp_A1 = evalute_exp(param, A[0])
  exp_A2 = evalute_exp(param, A[1])
  exp_A3 = evalute_exp(param, A[2])
  b_minus_phi_vector = np.subtract(b, [R*exp_A1, R*exp_A2, R*exp_A3])
  return R*exp_C + np.dot(kwargs['y'], b_minus_phi_vector) + (c/2)*(np.linalg.norm(b_minus_phi_vector)**2)

def cost_func_ALM_grad_step_and_cost(param, **kwargs):
  """
  This function evaluates new parameters, as well as returns current cost function value

  :param param: circuit parameters
  """
  gradient = qml.grad(evalute_exp, argnum=0)
  grad_C = gradient(param, C)
  grad_A1 = gradient(param, A[0])
  grad_A2 = gradient(param, A[1])
  grad_A3 = gradient(param, A[2])
  exp_A1 = evalute_exp(param, A[0])
  exp_A2 = evalute_exp(param, A[1])
  exp_A3 = evalute_exp(param, A[2])
  full_gradient = R*grad_C - R*(kwargs['y'][0]*grad_A1 + kwargs['y'][1]*grad_A2 + kwargs['y'][2]*grad_A3) - R*c*(b[0]*grad_A1 - R*exp_A1*grad_A1 + b[1]*grad_A2 - R*exp_A2*grad_A2 + b[2]*grad_A3 - R*exp_A3*grad_A3)
  updated_param = np.subtract(param, 0.01*full_gradient)
  prev_cost_func_value = cost_func_ALM(param, y=y)
  return updated_param, prev_cost_func_value


# + id="MtJ_D7oSuBtn"
# Inner maximizer function
def sup_wrt_theta(y):
  
  # initialize thetas
  params_ALM = np.random.random(qml.templates.StronglyEntanglingLayers.shape(n_layers=4, n_wires=3))
  # cost function storage
  cost_func_ALM_store = [cost_func_ALM(params_ALM, y=y)]

  # store the params
  params_ALM_store = [params_ALM]

  max_iterations = 100
  max_tol = 1e-01

  # iterate
  while True:
    params_ALM, prev_cost_func_ALM_value = cost_func_ALM_grad_step_and_cost(params_ALM, y=y)

    # append to the store
    cost_func_ALM_store.append(cost_func_ALM(params_ALM, y=y))
    params_ALM_store.append(params_ALM)

    # check the tolerance
    tol = np.abs(cost_func_ALM_store[-1] - prev_cost_func_ALM_value)
    #print("Cost function: ", cost_func_ALM_store[-1], " Tolerance: ", tol)

    if tol <= max_tol:
      break
    
  return params_ALM


# + id="3OgPRTa_uKez"
# Outer minimization iteratation
max_outer_iterations = 100
max_outer_tol = 1e-04

# penalty parameter
c = 0

# initialize y
y = [0, 0, 0]

for i in range(max_outer_iterations):
    optimal_thetas = sup_wrt_theta(y)
  # evaluate the expectation values of A1, A2, and A3 wrt optimal thetas
    exp_A1 = evalute_exp(optimal_thetas, A[0])
    exp_A2 = evalute_exp(optimal_thetas, A[1])
    exp_A3 = evalute_exp(optimal_thetas, A[2])

  # update the dual variables
    y[0] += 0.0005*(b[0] - R*exp_A1)
    y[1] += 0.0005*(b[1] - R*exp_A2)
    y[2] += 0.0005*(b[2] - R*exp_A3)

    print("Step:", i , " Cost function value:", cost_func_ALM(optimal_thetas, y=y))

StarcoderdataJupyterScriptsDedupFiltered