<jupyter_start><jupyter_text>The Portfolio Rebalancing Optimization using QAOA with soft constraints In this notebook, we will apply the QAOA algorithm to solve the portfolio optimization problem using a soft constraint$$C_{RR}({\bf z}) = \sum_i (1-\lambda) r_i z_i - \lambda\sum_{ij} \sigma_{ij} z_i z_j$$$$C_{TC}({\bf z}) = \sum_i T (1- \delta(z-y))$$The soft constraint that we impose will be given $$H =C_{RR}({\bf z})+C_{TC}({\bf z})+A\left(D-\sum_i z_i \right)^2$$  Generate the inital ground state wave function of the reference Hamiltonian<jupyter_code># Choose the number of portfolios to optimize
N = 2

# Define the transaction cost
T = sympy.symbols("T")
D = sympy.symbols("D")
A = sympy.symbols("A")

# Define our previous holdings
y = np.zeros(N)

# Define qubits on the grid.
qubits = [cirq.LineQubit(i) for i in range(2*N)]

# Define the indices that we will use
index_grid_pair = [(2*i,2*i+1) for i in range(0,N) ]

print(index_grid_pair)
print()

# Next we take these quibits which are in states |0> and apply a Hadamard gate to each one.
# this generates the eigenvectors of the reference Hamiltonian

print(qubits)

hadamard_operators = [cirq.H(qubits[i]) for i in range(len(qubits))]
measurements = cirq.Moment([cirq.measure(qubits[i],key=str(i)) for i in range(len(qubits))])
#measurements = cirq.Moment( [cirq.measure(qubits[2],qubits[0],qubits[1],qubits[3],key='x') ])

# Create the circuit from this 
circuit = cirq.Circuit((hadamard_operators,measurements))

#print(circuit)
SVGCircuit(circuit)<jupyter_output>[(0, 1), (2, 3)]

[cirq.LineQubit(0), cirq.LineQubit(1), cirq.LineQubit(2), cirq.LineQubit(3)]<jupyter_text>Implementation of the mixing operatorWe have generated the ground state wave function of the reference Hamiltonian. Now, we need to construc the Unitary operator $U({\beta})$ that represents the time evolution operator $$\begin{align}U({\bf \beta}) &= e^{-i\beta B}, \\B &= \sum_{i} X_i.\end{align}$$In the `cirq` api, this operator is implemented as $$\text{cirq.rx} \equiv e^{-i X \beta / 2}.$$<jupyter_code># Generate a list of symbolic arguments
beta_list = ["b_"+str(i) for i in range(1)]

# Convert them to symbolic numpy objects
betas = sympy.symbols(beta_list)


mixer_operators = [cirq.rx(2*betas[0]).on(qubits[i]) for i in range(len(qubits))]

# Create the circuit from this 
circuit = cirq.Circuit((hadamard_operators,mixer_operators,measurements))

#print(circuit)

SVGCircuit(circuit)<jupyter_output><empty_output><jupyter_text>Implementation of the exponential cost-functionFor the portfolio optimization problem we have that the operator $Z_i$ is the difference between the Pauli-$z$ matrices$$\begin{align}Z_i = \frac{1}{2}\left(z^{+}_{i}-z^{-}_{i} \right),\end{align}$$where the $z^+_i$ and $z^-_i$ operate on the $i$-th qubit pair $|q_+,q_- \rangle$ and $q_{\pm}=\lbrace 0,1\rbrace$$$C_{RR}({\bf z}) = \sum_i (1-\lambda) r_i z_i - \lambda\sum_{ij} \sigma_{ij} z_i z_j$$$$\begin{align}e^{-i \gamma C_{RR}(z) } &=   \rm{Exp}\left(  i \gamma (1-\lambda) \sum_{k=1} \mu_k Z_k \right)  \rm{Exp}\left( - i \lambda \gamma  \sum_{ij} \sigma_{ij} Z_i Z_j \right) \\&=   \prod_{k=1} \rm{Exp}\left( i \gamma (1-\lambda) \mu_k Z_k \right) \prod_{ij} \rm{Exp}\left( - i \lambda \gamma  \sigma_{ij} Z_i Z_j \right) \\&=   \prod_{k=1} \rm{Exp}\left( i \gamma (1-\lambda) \frac{1}{2} \mu_k z^{+}_{k} \right)\left( i \gamma (1-\lambda) \frac{1}{2} \mu_k (-z^{-}_{k}) \right) \prod_{ij} \rm{Exp}\left( - i \lambda \frac{1}{4} \gamma  \sigma_{ij} \left(z^{+}_{i}z^{+}_{j}-z^{+}_{i}z^{-}_{j}-z^{-}_{i}z^{+}_{j}+z^{-}_{i}z^{-}_{j} \right) \right)\end{align}$$<jupyter_code># Returns
mu = sympy.symbols(["mu_"+str(k) for k in range(N)])


# Convert them to symbolic numpy objects
lam,gamma = sympy.symbols("lam, gamma")

cost_returns_operators = [cirq.rz((-1)**int(i)*(1-lam)*gamma*mu[int(i/2)]).on(qubits[i]) for i in range(len(qubits))]

# Create the circuit from this 
circuit = cirq.Circuit((hadamard_operators,mixer_operators,cost_returns_operators,measurements))

#print(circuit)

SVGCircuit(circuit)<jupyter_output><empty_output><jupyter_text>Implementation of $e^{ i \gamma z_i z_j}$ This operator only contributes if $i \neq j$. If $i=j$ then the operator is simply an overall constant phase that can be ignored. First, we will implement $e^{ i \gamma \pi z_i z_j}$ by rescaling the $\gamma$ parameter. Furthermore, we note that an overall phase factor does not not change the final result when taking expectation values with wavefunctions, so we will remove an $e^{i\pi\gamma}$ phase from the implementation of this operator to simplify its decomposition.  $$ \begin{align}  \text{Exp}\left(i \pi \gamma z_i z_j \right) &=  \begin{bmatrix}   e^{\left(i \pi \gamma   \right)} & 0 & 0 & 0 \\  0 & e^{-\left(i \pi \gamma   \right)} & 0 & 0 \\  0 & 0 & e^{-\left(i \pi \gamma \right)} & 0 \\  0 & 0 & 0 & e^{\left(i \pi \gamma  \right)}  \end{bmatrix},\\  &=  e^{i \pi \gamma} \begin{bmatrix}   1 & 0 & 0 & 0 \\  0 & e^{-2\left(i \pi \gamma   \right)} & 0 & 0 \\  0 & 0 & e^{-2\left(i \pi \gamma \right)} & 0 \\  0 & 0 & 0 &1  \end{bmatrix}, \end{align} $$  Now, we note that in `cirq`, the above matrix is implemented as  $$ \text{cirq.ZZPowGate}(t) =    \begin{bmatrix}   1 & 0 & 0 & 0 \\  0 & e^{i\pi t} & 0 & 0 \\  0 & 0 &  e^{i\pi t} & 0 \\  0 & 0 & 0 &  1  \end{bmatrix}. $$   For the purposes of implementation also note that  $$ X_2  \equiv \left( I \otimes X \right)  = \begin{bmatrix}  X & {\bf 0}\\  {\bf 0} & X   \end{bmatrix}   = \begin{bmatrix}  0 & 1 & 0& 0\\  1 & 0 & 0& 0\\  0 & 0 & 0& 1\\  0 & 0 & 1& 0\\  \end{bmatrix}, $$   $$ X_1  \equiv \left( X \otimes I \right)  = \begin{bmatrix}  {\bf 0} & I\\  I & {\bf 0}   \end{bmatrix}   = \begin{bmatrix}  0 & 0 & 1& 0\\  0 & 0 & 0& 1\\  1 & 0 & 0& 0\\  0 & 1 & 0& 0\\  \end{bmatrix}, $$    $$ X_1  X_2 = X_2  X_1 = \begin{bmatrix}  {\bf 0} & X\\  X & {\bf 0}   \end{bmatrix}   = \begin{bmatrix}  0 & 0 & 0& 1\\  0 & 0 & 1& 0\\  0 & 1 & 0& 0\\  1 & 0 & 0& 0\\  \end{bmatrix}. $$ Finally, we arrive at the decomposition of our operator of interest in terms of gate operations $$\text{Exp}\left(i \pi \gamma z_i z_j \right) = e^{i\gamma \pi}\text{cirq.ZZPowGate}(-2\gamma)$$<jupyter_code>def exp_zz(gamma,qubit1,qubit2):
    """
    Exponentiates the Exp[]
    """
    
    #D1 = [ cirq.X.on(qubit2), cirq.CZPowGate(exponent= -2*gamma).on(qubit1,qubit2), cirq.X.on(qubit2)]
    #D2 = [ cirq.X.on(qubit1),cirq.X.on(qubit2), cirq.CZPowGate(exponent= -2*gamma).on(qubit1,qubit2),cirq.X.on(qubit1),cirq.X.on(qubit2)]
    
    #decomposition = D2.append(D1)
    
    return [cirq.ZZPowGate(exponent= -2*gamma).on(qubit1,qubit2)]



def remove_digonal_combinations(combination_list):
    """
    Simple utility function to remove the diagonal combinations from a list
    of possible connections
    """
    
    purged_list = []
    
    
    for i in range(len(index_grid_pp)):
        p1,p2 = combination_list[i]
    
        if(p1!=p2):
            purged_list.append((p1,p2))
    
    
    return purged_list

def map_combinations_to_portfolio_index(combination_list):
    
    
    portfolio_index= [ (int(combination_list[i][0]/2), int(combination_list[i][1]/2)) for i in range(len(combination_list))]
    
    return portfolio_index


# Covariance matrix
sigma = sympy.symbols([["sigma_"+str(i)+"_"+str(j) for i in range(N)] for j in range(N)])

print(qubits)
print()
# Define the s^+ , s^- permutations that we will use for summations
index_grid_p = [2*i for i in range(0,N)]
index_grid_m = [2*i+1 for i in range(0,N)]
index_grid_pm = [(2*i,2*j+1) for i in range(0,N) for j in range(0,N) ]
index_grid_pm_ii = [(2*i,2*i+1) for i in range(0,N)]
index_grid_mp = [(2*i+1,2*j) for i in range(0,N) for j in range(0,N) ]
index_grid_pp = [(2*i,2*j) for i in range(0,N) for j in range(0,N)  ]
index_grid_mm = [(2*i+1,2*j+1) for i in range(0,N) for j in range(0,N) ]


index_grid_pp_offdiagonal  =  remove_digonal_combinations(index_grid_pp)
index_grid_mm_offdiagonal  =  remove_digonal_combinations(index_grid_mm)
index_grid_pm_offdiagonal  =  remove_digonal_combinations(index_grid_pm)
index_grid_mp_offdiagonal  =  remove_digonal_combinations(index_grid_mp)


index_grid_pp_portfolio = map_combinations_to_portfolio_index(index_grid_pp_offdiagonal)
index_grid_mm_portfolio = map_combinations_to_portfolio_index(index_grid_mm_offdiagonal)
index_grid_pm_portfolio = map_combinations_to_portfolio_index(index_grid_pm_offdiagonal)
index_grid_mp_portfolio = map_combinations_to_portfolio_index(index_grid_mp_offdiagonal)

print("++ Combinations")
print(index_grid_pp_offdiagonal)
print()
print("-- Combinations")
print(index_grid_mm_offdiagonal)
print()
print("+- Combinations")
print(index_grid_pm_offdiagonal)
print()
print("-+ Combinations")
print(index_grid_mp_offdiagonal)
print()
print("+- Combinations| Same portfolio")
print(index_grid_pm_ii)


k=0
print(sigma[int(index_grid_pp[k][0]/2)][int(index_grid_pp[k][1]/2)])

pp1,pp2 = index_grid_pp_offdiagonal[0]
pp1_portfolio, pp2_portfolio =index_grid_pp_portfolio[0]
circuit = cirq.Circuit(exp_zz(gamma*lam/4*sigma[pp1_portfolio][pp2_portfolio],qubit1=qubits[pp1],qubit2=qubits[pp2]))

print("")
print("Example circut for application of Exp(zi zj) \n")
#print(circuit)
SVGCircuit(circuit)
#print("")<jupyter_output>[cirq.LineQubit(0), cirq.LineQubit(1), cirq.LineQubit(2), cirq.LineQubit(3)]

++ Combinations
[(0, 2), (2, 0)]

-- Combinations
[(1, 3), (3, 1)]

+- Combinations
[(0, 1), (0, 3), (2, 1), (2, 3)]

-+ Combinations
[(1, 0), (1, 2), (3, 0), (3, 2)]

+- Combinations| Same portfolio
[(0, 1), (2, 3)]
sigma_0_0

Example circut for application of Exp(zi zj)<jupyter_text>Implementation of of full volatility term$$\text{Exp}\left( \sum_{ij} \lambda \sigma_{ij}Z_i Z_j \right) = \text{Exp}\left( \sum_{ij} \frac{\lambda}{4}\sigma_{ij}\left(z^{+}_{i} z^{+}_{j}-z^{+}_{i} z^{-}_{j}-z^{-}_{i} z^{+}_{j}+z^{-}_{i} z^{-}_{j} \right) \right)$$<jupyter_code># Now let us build the entire circuit up to now

cost_volatility_operators_pp = []
cost_volatility_operators_mm = []
cost_volatility_operators_pm = []
cost_volatility_operators_mp = []

# Loop over all ++ combinations
for k in range(len(index_grid_pp_offdiagonal)):
    i_qubit,j_qubit = index_grid_pp_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_pp_portfolio[k]
    cost_volatility_operators_pp.append(exp_zz(gamma*lam/4*sigma[i_portfolio][j_portfolio],qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

# Loop over all +- combinations
for k in range(len(index_grid_pm_offdiagonal)):
    i_qubit,j_qubit = index_grid_pm_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_pm_portfolio[k]
    cost_volatility_operators_pm.append(exp_zz((-1)*gamma*lam/4*sigma[i_portfolio][j_portfolio],qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

    # Loop over all -+ combinations
for k in range(len(index_grid_mp_offdiagonal)):
    i_qubit,j_qubit = index_grid_mp_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_mp_portfolio[k]
    cost_volatility_operators_mp.append(exp_zz((-1)*gamma*lam/4*sigma[i_portfolio][j_portfolio],qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

    # Loop over all -- combinations
for k in range(len(index_grid_mm_offdiagonal)):
    i_qubit,j_qubit = index_grid_mm_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_mm_portfolio[k]
    cost_volatility_operators_mm.append(exp_zz(gamma*lam/4*sigma[i_portfolio][j_portfolio],qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))


#cost_returns_operators = [cirq.rz((-1)**int(i)*(1-lam)*gamma*mu[int(i/2)]).on(qubits[i]) for i in range(len(qubits))]

# Create the circuit from this 
circuit = cirq.Circuit((cost_volatility_operators_pp,cost_volatility_operators_pm,cost_volatility_operators_mp,cost_volatility_operators_mm))
SVGCircuit(circuit)<jupyter_output><empty_output><jupyter_text>Implement the exponential operator of the transaction cost$$C_{TC}(z) = \frac{1}{4}T\left( 3+(1-y^2_i-y_i)z^{+}_{i}+(1-y^2_i+y_i)z^{-}_{i}+(2y^2_i-1)z^{+}_{i}z^{-}_{i} \right)$$$$\begin{align}\text{Exp}\left( i\pi \gamma C_{TC}(z) \right) = \prod_{j=1} \text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(1-y^2_j-y_j\right)z^{+}_{j}\right)\text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(1-y^2_j+y_j\right)z^{-}_{j} \right)\text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(2y^2_j-1\right)z^{+}_{j}z^{-}_{j} \right) \\= \prod_{j=1} \text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(1-y^2_j-y_j\right)z^{+}_{j}\right) \prod_{j=1} \text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(1-y^2_j+y_j\right)z^{-}_{j} \right) \prod_{j=1}\text{Exp}\left(  i\pi \gamma \frac{1}{4}T\left(2y^2_j-1\right)z^{+}_{j}z^{-}_{j} \right)\end{align}$$<jupyter_code>transaction_operators_p_i = []
transaction_operators_m_i = []
transaction_operators_pm_ii = []

print(index_grid_p)



for k in range(len(index_grid_p)):
    angle_i = (1/4)*T*(1-y[k]**2-y[k])*gamma
    i_qubit = index_grid_p[k]
    transaction_operators_p_i.append(cirq.rz(angle_i).on(qubits[i_qubit]) )
    
for k in range(len(index_grid_pm_ii)):
    angle_k = (1/4)*T*(2*y[k]**2-1)*gamma
    i_qubit, j_qubit = index_grid_pm_ii[k]
    i_portfolio, j_portfolio =index_grid_pm_ii[k]
    transaction_operators_pm_ii.append(exp_zz( angle_k,qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

for k in range(len(index_grid_m)):
    angle_i = (1/4)*T*(1-y[k]**2+y[k])*gamma
    i_qubit = index_grid_m[k]
    transaction_operators_m_i.append(cirq.rz(angle_i).on(qubits[i_qubit]) )
    

circuit = cirq.Circuit((transaction_operators_p_i,transaction_operators_pm_ii,transaction_operators_m_i))
SVGCircuit(circuit)<jupyter_output>[0, 2]<jupyter_text>Now we implement the soft constraint Term$$P_{\rm INV}(z) = A\left(\sum^N_i z^2_i -D \right)^2$$which can be factored into$$\begin{align}P_{\rm INV}(z) = \frac{A}{4} \left( z^+_i z^+_j-z^+_i z^-_j - z^+_i z^+_j+ z^-_i z^-_j \right) -AD(z^+_j-z^-_j)+AD^2\end{align}$$<jupyter_code># We implement the AD terms of the penalty function just like we did for the returns
penalty_AD_operators = [cirq.rz( (-1)**i * gamma*A*D).on(qubits[i]) for i in range(len(qubits))]

# Create the circuit from this 
circuit = cirq.Circuit((penalty_AD_operators))

#print(circuit)

SVGCircuit(circuit)
# Now let us build the entire circuit up to now

penalty_A_operators_pp = []
penalty_A_operators_mm = []
penalty_A_operators_pm = []
penalty_A_operators_mp = []

# Loop over all ++ combinations
for k in range(len(index_grid_pp_offdiagonal)):
    i_qubit,j_qubit = index_grid_pp_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_pp_portfolio[k]
    penalty_A_operators_pp.append(exp_zz(gamma*A/4,qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

# Loop over all +- combinations
for k in range(len(index_grid_pm_offdiagonal)):
    i_qubit,j_qubit = index_grid_pm_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_pm_portfolio[k]
    penalty_A_operators_pm.append(exp_zz((-1)*A/4*gamma,qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

    # Loop over all -+ combinations
for k in range(len(index_grid_mp_offdiagonal)):
    i_qubit,j_qubit = index_grid_mp_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_mp_portfolio[k]
    penalty_A_operators_mp.append(exp_zz((-1)*A/4*gamma,qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

    # Loop over all -- combinations
for k in range(len(index_grid_mm_offdiagonal)):
    i_qubit,j_qubit = index_grid_mm_offdiagonal[k]
    i_portfolio, j_portfolio =index_grid_mm_portfolio[k]
    penalty_A_operators_mm.append(exp_zz(A/4*gamma,qubit1=qubits[i_qubit],qubit2=qubits[j_qubit]))

# Create the circuit from this 
circuit = cirq.Circuit((penalty_A_operators_pp ,penalty_A_operators_pm,penalty_A_operators_mp,penalty_A_operators_mm))
#circuit = cirq.Circuit((penalty_A_operators_pm))
SVGCircuit(circuit)<jupyter_output><empty_output><jupyter_text>Construct the final circuit and run the simulator for a single instance of the parameters<jupyter_code>penalty = (penalty_AD_operators,penalty_A_operators_pp ,penalty_A_operators_pm,penalty_A_operators_mp,penalty_A_operators_mm)
portfolio_lagrangian = (cost_returns_operators,cost_volatility_operators_pp,cost_volatility_operators_pm,cost_volatility_operators_mp,cost_volatility_operators_mm)
transaction_costs = (transaction_operators_p_i,transaction_operators_pm_ii,transaction_operators_m_i)


circuit = cirq.Circuit(hadamard_operators,mixer_operators,portfolio_lagrangian,transaction_costs,penalty,measurements)
SVGCircuit(circuit)
def convert_to_y(x):
    '''
    Converts the bitstring x to the vector y={-1,0,1}
    for the portfolio
    '''
    
    # Rescales the bit string to [-1 or 1]
    s = 2.0*np.array(x)-1.0
    
    y = np.zeros(N)
    
    for i in range(N):
        y[i] = (s[2*i]-s[2*i+1])/2
    
    return y

def portfolio_penalty_cost(y,A,D):
    penalty = 0.0
    
    
    penalty = A*(np.sum(y)-D)**2
    
    return penalty

def portfolio_cost_function(y):
    '''
    The cost function for a portfolio
    
    y-- a bitstring array giving the results of the portfolio {-1,0,1}
    '''
    cost = 0 
    
    for i in range(N):
        Zi = x[2*i]+x[2*i+1]
    
    
    return cost
def convert_measurements_to_array(results,N,n_repetitions):
    '''
    
    results -- 
    '''
    
    # Instantiate the array
    X = np.zeros((2*N,n_repetitions))
    
    # Retrieve all of the measurements for each Qbit and store the result in the tensor
    for i in range(len(result.measurements.keys())):
        key = list(result.measurements.keys())[i]
        X[i,:] = np.array(result.measurements[key]).flatten()
    
    
    return X


def convert_measurements_to_Z(X,N,n_repetitions):
    '''
    Converts collection of simulations in the computational basis to 
    the Z-values
    
    Z_i = (1/2)(z^+_i - z^-_i)
    
    where z^{+/-} = {-1,1}
    '''
    
    Z = np.zeros((N,n_repetitions))
    
    for i in range(0,n_repetitions):
        x_bitstring = X[:,i]
        
        # Converts the bitstring to portfolio 
        Z[:,i] = convert_to_y(x_bitstring)
    
    return Z


def count_instances(N,results,n_repetitions):
    '''
    
    '''
    
    X = convert_measurements_to_array(results,N,n_repetitions)
    Z = convert_measurements_to_Z(X,N,n_repetitions)
    
    # The dictionary to store the results
    portfolio_results = {}
    
    
    # The values that can be taken by Zi
    portfolio_index = [-1,0,1]
    
    portfolio_keys = itertools.product(portfolio_index,repeat=N)
    
    labels = []
    label_counts = []
    label_penalty = []
    
    
    for state in list(portfolio_keys):
        
        # convert the tuple to numpy integer vector 
        Z_state =  np.array([state[i] for i in range(len(state))]).astype(int)
        
        # convert the Z_state to a string
        state_key = ''.join(Z_state.astype('str'))
        
        state_count = 0
        
        for i in range(n_repetitions):
            
            # retrieve a sample
            Zi = Z[:,i].flatten()
            
            label_penalty.append(portfolio_penalty_cost(Zi,A_val,D_val))
            
            # check if Zi sample matches the Z-label
            if(np.array_equal(Z_state,Zi)==True):
                state_count+=1
        
        # Save all of the values
        labels.append(state_key)
        label_counts.append(state_count)
        
        portfolio_results[state_key] = state_count
        
    portfolio_results['penalty_cost_samples'] = label_penalty 
    portfolio_results['labels']= labels
    portfolio_results['counts']= label_counts
    
    return portfolio_results

print(mu,sigma,lam,gamma,betas)


resolved_params = {}

A_val = 1.0
D_val = 0.0


for k in range(len(mu)):
    
    mu_k = mu[k]
    
    if(k==0):
        resolved_params[mu_k] = 100.0
    else:
        resolved_params[mu_k] = -10.0

for i in range(len(mu)):
    for j in range(len(mu)):
        
        if(i==j):
            resolved_params[sigma[i][j]] = 0.0
        else:
            resolved_params[sigma[i][j]] = 0.0

for beta in betas:
    resolved_params[beta] =  2.3
   
resolved_params[gamma] = 0.0
resolved_params[lam] = 0


resolved_params[A] = A_val
resolved_params[D] = D_val


print(resolved_params)


# Resolve all of the parameters
resolver = cirq.ParamResolver(resolved_params)
import itertools

n_repetitions = 5000

#circuit = cirq.Circuit(hadamard_operators,mixer_operators,portfolio_lagrangian,penalty)
#circuit = cirq.Circuit(hadamard_operators,mixer_operators,penalty)
#penalty_AD_operators
#circuit = cirq.Circuit(hadamard_operators,mixer_operators,penalty_AD_operators)
circuit = cirq.Circuit(hadamard_operators,mixer_operators,penalty_AD_operators)
circuit.append(measurements)

simulator = cirq.Simulator()
result = simulator.run(circuit, resolver,repetitions=n_repetitions)


portfolio_results = count_instances(N,result,n_repetitions)

plt.figure(figsize=(10,10))
plt.title('D-value: '+str(D_val))
plt.bar(np.arange(len(portfolio_results['labels'])), portfolio_results['counts'], tick_label=portfolio_results['labels'])
plt.xticks(rotation=90)
plt.show()
#cirq.plot_state_histogram(result)
# # Generate a grid for the parameters
# sweep_b0 = cirq.Linspace(key=betas[0], start=0.0, stop=1.0, length=10)
# angles_b0 = [angle for angle in sweep_b0]
# #sweep_gamma = cirq.Linspace(key=gamma, start=0.0, stop=1.0, length=10)
# #sweeps = (sweep_b0*sweep_gamma)

# result_sweep = simulator.run_sweep(circuit,params=sweep_b0, repetitions=10)

# #
# X = np.zeros((len(qubits),10))

# for i in range(len(result.measurements.keys())):
#     key = list(result.measurements.keys())[i]
#     X[i,:] = np.array(result.measurements[key]).flatten()

# print(X)


plt.figure(figsize=(10,10))
plt.hist(portfolio_results['penalty_cost_samples'])
plt.show()

plt.clf()
cirq.plot_state_histogram(result)

SVGCircuit(circuit)
print(result_sweep[1].measurements['0'].flatten())
print(result_sweep[1].measurements['1'].flatten())
print(result_sweep[1].measurements['2'].flatten())
print(result_sweep[1].measurements['3'].flatten())
plt.hist(label_penalty)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15604805679439337, 'fraction_numerical': 0.02070628925093292, 'mean_word_length': 4.485521717423865, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '625978', 'license': 'Apache-2.0', 'n_tokens_mistral': 8968, 'n_tokens_neox': 8563, 'n_words': 1665, 'path': '.ipynb_checkpoints/Portfolio_rebalancing_soft_constraints-checkpoint.ipynb', 'repo_name': 'OscarJHernandez/qc_portfolio_optimization', 'size': '21974'}