!python3 -m pip install -q qiskit

!python3 -m pip install -q qiskit_ibm_runtime

from qiskit import IBMQ
from qiskit_ibm_runtime import QiskitRuntimeService, Session, Options, Sampler, Estimator
from qiskit.circuit.quantumcircuit import QuantumCircuit
from qiskit.quantum_info.operators import Operator
from qiskit.extensions.unitary import UnitaryGate
from qiskit.providers.aer import AerSimulator, Aer
from qiskit.execute_function import execute
import numpy as np
from typing import Dict, Optional

# ====================================
# pass IBM API Token
# ====================================

QiskitRuntimeService.save_account(channel='ibm_quantum', token="", overwrite=True)



num_players = 2

# define the payoff matrices. Assume, the 1st matrix is Alice's payoff, the secomd matrix is Bob's payoff
payoff_matrix =[[[3, 0],
                [5, 1]],
                [[3, 5],
                [0, 1]]]

# players' matrices
# |C> strategy
# bob = [
#     [1, 0],
#     [0, 1]]

# |D> strategy
# bob = [
#     [0, 1],
#     [1, 0],
# ]


def alice_plays_quantum(theta:float, phi:float):
  """
    Set up a quantum game for Alice playing quantum, and Bob playing classically. 
    Arg:
      theta: the value of the theta parameter.
      phi: the value of the phi parameter.
    Returns: 
      qc: the class of 'qiskit.circuit.quantumcircuit.QuantumCircuit', 
            a constructed quantum circuit for the quantum game set up.
  """
  alice = np.array([[np.exp(1.0j*phi)*np.cos(theta/2), np.sin(theta/2)],
          [-np.sin(theta/2), np.exp(-1.0j*phi)*np.cos(theta/2)]])

  bob = np.array([
    [1, 0],
    [0, 1]]).astype(complex)

  qc = QuantumCircuit(num_players)
  J = 1/np.sqrt(2)*np.array([[1.0, 0.0, 0.0, 1.0j],
                             [0.0, 1.0, 1.0j, 0.0],
                             [0.0, 1.0j, 1.0, 0.0],
                             [1.0j, 0.0, 0.0, 1.0]]).astype(complex)
  J_unitary = UnitaryGate(Operator(J))
  qc.append(J_unitary, [0,1])
  qc.barrier()


  unitary_alice = UnitaryGate(Operator(alice))
  unitary_bob = UnitaryGate(Operator(bob))

  qc.append(unitary_alice, [0])
  qc.append(unitary_bob, [1])
  qc.barrier()

  Jdagger_unitary = UnitaryGate(Operator(J.conj().T))
  qc.append(Jdagger_unitary, [0,1])
  
  return qc

# test run
# define the simulator
simulator = Aer.get_backend('statevector_simulator')

theta = np.pi
phi = 0.0
qc = alice_plays_quantum(theta, phi)
results = execute(qc, simulator, shots=1024).result().get_counts()
results

# calculating payoffs

def get_payoff(counts):
  """
    Calculate the reward for the players after the game ends.
  """
  payoff_bob = []
  payoff_alice = []

  for strategy, prob in counts.items():
    strategy_bob = int(strategy[1])
    strategy_alice = int(strategy[0])
    payoff_bob.append(prob * payoff_matrix[0][strategy_alice][strategy_bob])
    payoff_alice.append(prob * payoff_matrix[0][strategy_bob][strategy_alice])

  return sum(payoff_alice), sum(payoff_bob)

# change the size on the quantum space to explore more quantum strategies

space_size = 4

for phi in np.linspace(0, np.pi/2, space_size):
  for theta in np.linspace(0, np.pi, space_size):
    qc = alice_plays_quantum(theta, phi)
    results = execute(qc, simulator, shots=1024).result().get_counts()
    payoff_alice, payoff_bob = get_payoff(results)
    print("theta = {}, phi = {}, results = {}, Alice's payoff {}, Bob's payoff {}".format(theta, phi, results, payoff_alice, payoff_bob))
    print("Next Game")

qc = alice_plays_quantum(theta = np.pi/2, phi = np.pi/4)

job = simulator.run(qc, shots=2048)
result = job.result()
outputstate = result.get_statevector(qc, decimals=3)
print(outputstate)

from qiskit.visualization import plot_state_city
plot_state_city(outputstate)

from qiskit.visualization import plot_histogram
counts = result.get_counts()
plot_histogram(counts)

# probabilities of a certain outcome
counts

num_players = 2

def alice_bob_play_quantum(theta_a:float, phi_a:float, theta_b:float, phi_b:float, measure=False):
  """
    Set up a quantum game. Both players have access to quantum strategies.
    Args:
      theta_a: Theta parameter for Alice's unitary.
      phi_a: Phi parameter for Alice's unitary.
      theta_b: Theta parameter for Bob's unitary.
      phi_b: Phi parameter for Bob's unitary.
    Returns:
      qc: the class of 'qiskit.circuit.quantumcircuit.QuantumCircuit', 
            a constructed quantum circuit for the quantum game set up.
  """

  alice = np.array([
      [np.exp(1.0j*phi_a)*np.cos(theta_a/ 2), np.sin(theta_a / 2)],
      [-np.sin(theta_a / 2), np.exp(-1.0j*phi_a)*np.cos(theta_a / 2)]])

  bob = np.array([
      [np.exp(1.0j*phi_b)*np.cos(theta_b / 2), np.sin(theta_b / 2)],
      [-np.sin(theta_b / 2), np.exp(-1.0j*phi_b)*np.cos(theta_b / 2)]])

  qc = QuantumCircuit(num_players)

  J = 1/np.sqrt(2)*np.array([[1.0, 0.0, 0.0, 1.0j],
                             [0.0, 1.0, 1.0j, 0.0],
                             [0.0, 1.0j, 1.0, 0.0],
                             [1.0j, 0.0, 0.0, 1.0]]).astype(complex)
  J_unitary = UnitaryGate(Operator(J))
  qc.append(J_unitary, [0,1])
  qc.barrier()

  unitary_alice = UnitaryGate(Operator(alice))
  unitary_bob = UnitaryGate(Operator(bob))

  qc.append(unitary_alice, [0])
  qc.append(unitary_bob, [1])
  qc.barrier()

  Jdagger_unitary = UnitaryGate(Operator(J.conj().T))
  qc.append(Jdagger_unitary, [0,1])

  if measure:
    qc.measure_all()
    
  return qc

alice_reward = 0
bob_reward = 0
draw = 0

space_size = 2

for phi1 in np.linspace(0, 2*np.pi, space_size):
  for theta1 in np.linspace(0, np.pi, space_size):
    for phi2 in np.linspace(0, 2*np.pi, space_size):
      for theta2 in np.linspace(0, np.pi, space_size):
        qc = alice_bob_play_quantum(theta1, phi1, theta2, phi2)
        results = execute(qc, simulator, shots=1024).result().get_counts()
        payoff_alice, payoff_bob = get_payoff(results)

        # count winning 
        if payoff_alice > payoff_bob:
          alice_reward += 1
        elif payoff_bob > payoff_alice:
          bob_reward += 1
        else:
          draw += 1
        # print results 
        print("theta_alice = {}, phi_alice = {}, theta_bob = {}, phi_bob = {}".format(theta1, phi1, theta2, phi2 ))
        print("results = {}, Alice's raward {}, Bob's reward {}".format(results, payoff_alice, payoff_bob))
        print("Next Game")
        print("===================================================")

print("In {} games Alice gets a higher reward than Bob.". format(alice_reward))
print("In {} games, Bob gets a higher reward than Alice.".format(bob_reward))
print("In {} games Alice and Bob get equal reward.".format(draw))

matrix = np.array([[1.0j, 0], [0, -1.0j]])
qc = QuantumCircuit(2)
gate = UnitaryGate(Operator(matrix))

J = 1/np.sqrt(2)*np.array([[1.0, 0.0, 0.0, 1.0j],
                             [0.0, 1.0, 1.0j, 0.0],
                             [0.0, 1.0j, 1.0, 0.0],
                             [1.0j, 0.0, 0.0, 1.0]]).astype(complex)
J_unitary = UnitaryGate(Operator(J))
qc.append(J_unitary, [0,1])
qc.barrier()

qc.append(gate, [0])
qc.append(gate, [1])
qc.barrier()

Jdagger_unitary = UnitaryGate(Operator(J.conj().T))
qc.append(Jdagger_unitary, [0,1])

results = execute(qc, simulator, shots=1024).result().get_counts()
alice_payoff, bob_payoff = get_payoff(results)
print("Strategy: {}".format(results))
print("Alice's payoff is {}, Bob's payoff is {}".format(alice_payoff, bob_payoff))

import matplotlib.pyplot as plt

space_size = 40

def payoff_plot():
  """
    Plot expected payoff distribution for Alice.
  """
  x = np.linspace(1, -1, space_size)
  y = np.linspace(-1, 1, space_size)
  X, Y = np.meshgrid(x, y)
  Z = np.zeros(X.shape)

  for i in range(0,space_size):
    for inner in range(0, space_size):
      if X[inner][i] < 0 and Y[inner][i] < 0:
        qc = alice_bob_play_quantum(0, X[inner][i]*np.pi/2, 0, Y[inner][i]*np.pi/2)
        payoff_alice, _ = get_payoff(execute(qc, simulator, shots=1024).result().get_counts())
        Z[inner][i] = payoff_alice
      elif X[inner][i] >= 0 and Y[inner][i] >= 0:
        qc = alice_bob_play_quantum(X[inner][i]*np.pi, 0, Y[inner][i]*np.pi, 0)
        payoff_alice, _ = get_payoff(execute(qc, simulator, shots=1024).result().get_counts())
        Z[inner][i] = payoff_alice
      elif X[inner][i] >= 0 and Y[inner][i] < 0:
        qc = alice_bob_play_quantum(X[inner][i]*np.pi, 0, 0, Y[inner][i]*np.pi/2)
        payoff_alice, _ = get_payoff(execute(qc, simulator, shots=1024).result().get_counts())
        Z[inner][i] = payoff_alice
      elif X[inner][i] < 0 and Y[inner][i] >= 0:
        qc = alice_bob_play_quantum(0, X[inner][i]*np.pi/2, Y[inner][i]*np.pi, 0)
        payoff_alice, _ = get_payoff(execute(qc, simulator, shots=1024).result().get_counts())
        Z[inner][i] = payoff_alice

  fig = plt.figure(figsize=(15, 10))
  ax = fig.add_subplot(111, projection='3d')
  ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm, antialiased=True)
  ax.set(xlim=(-1, 1), ylim=(1, -1), xlabel="Bob strategy", ylabel='Alice strategy', zlabel="Alice's expected payoff")
  plt.show()

payoff_plot()

def build_qcir_three_parameters(theta1, phi1, lmbda, theta2, phi2):
  """
    Set up a quantum game. Both players have access to quantum strategies space
  """
  alice = np.array([
      [np.exp(1.0j*phi1)*np.cos(theta1/2), 1.0j*np.exp(1.0j * lmbda) * np.sin(theta1/2)],
      [1.0j*np.exp(-1.0j * lmbda) * np.sin(theta1/2), np.exp(-1.0j*phi1)*np.cos(theta1/2)]])

  bob = np.array([
      [np.exp(1.0j*phi2)*np.cos(theta2/2), np.sin(theta2/2)],
      [-np.sin(theta2/2), np.exp(-1.0j*phi2)*np.cos(theta2/2)]])

  qc = QuantumCircuit(num_players)

  J = 1/np.sqrt(2)*np.array([[1.0, 0.0, 0.0, 1.0j],
                             [0.0, 1.0, 1.0j, 0.0],
                             [0.0, 1.0j, 1.0, 0.0],
                             [1.0j, 0.0, 0.0, 1.0]]).astype(complex)
  J_unitary = UnitaryGate(Operator(J))
  qc.append(J_unitary, [0,1])
  qc.barrier()

  unitary_alice = UnitaryGate(Operator(alice))
  unitary_bob = UnitaryGate(Operator(bob))

  qc.append(unitary_alice, [0])
  qc.append(unitary_bob, [1])
  qc.barrier()
  
  Jdagger_unitary = UnitaryGate(Operator(J.conj().T))
  qc.append(Jdagger_unitary, [0,1])
  
  return qc

alice_reward = 0
bob_reward = 0
draw = 0

space_size = 4

for phi1 in np.linspace(-np.pi, np.pi, space_size):
  for theta1 in np.linspace(0, np.pi, space_size):
    for phi2 in np.linspace(-np.pi, np.pi, space_size):
      for theta2 in np.linspace(0, np.pi, space_size):
        for lmbda in np.linspace(-np.pi, np.pi, 2*space_size):
          qc = build_qcir_three_parameters(theta1, phi1, lmbda, theta2, phi2)
          results = execute(qc, simulator, shots=1024).result().get_counts()
          payoff_alice, payoff_bob = get_payoff(results)

          # count winning 
          if payoff_alice > payoff_bob:
            alice_reward += 1
          elif payoff_bob > payoff_alice:
            bob_reward += 1
          else:
            draw += 1
          # print results 
          print("theta_alice = {}, phi_alice = {}, theta_bob = {}, phi_bob = {}".format(theta1, phi1, theta2, phi2 ))
          print("results = {}, Alice's payoff {}, Bob's payoff {}".format(results, payoff_alice, payoff_bob))
          print("Next Game")
          print("===================================================")

print("In {} games Alice gets a higher reward than Bob.". format(alice_reward))
print("In {} games, Bob gets a higher reward than Alice.".format(bob_reward))
print("In {} games Alice and Bob get equal reward.".format(draw))

!python3 -m pip install -q qiskit-ibm-provider

from qiskit.providers.ibmq import least_busy
from qiskit_ibm_provider import IBMProvider

provider = IBMProvider()
# ====================================
# pass IBM API Token
# ====================================
provider.save_account(token=, overwrite=True)


small_devices = provider.backends(filters=lambda x: x.configuration().n_qubits == 5 and not x.configuration().simulator)
least_busy(small_devices)

backend = least_busy(small_devices)

qc = alice_bob_play_quantum(theta_a = 0, phi_a = 0, theta_b = np.pi/2, phi_b = 0, measure=True)
job = execute(qc, backend, optimization_level=2)
result = job.result()

counts = result.get_counts()
counts

plot_histogram(counts)

job = execute(qc, backend=Aer.get_backend('statevector_simulator'), shots=1024)
result = job.result()
counts = result.get_counts()
print(counts)
plot_histogram(counts)

from functools import reduce
from operator import add
from qiskit import quantum_info as qi

def mixed_strategy(alpha_a:float, alpha_b:float):
  """
    Set up a quantum game. Both players have access to quantum strategies.
    Args:
      theta_a: Theta parameter for Alice's unitary.
      phi_a: Phi parameter for Alice's unitary.
      theta_b: Theta parameter for Bob's unitary.
      phi_b: Phi parameter for Bob's unitary.
    Returns:
      qc: the class of 'qiskit.circuit.quantumcircuit.QuantumCircuit', 
            a constructed quantum circuit for the quantum game set up.
  """

  qc = QuantumCircuit(2)
  # Alice strategies
  unitary1_a = np.eye(2).astype(complex)
  unitary2_a = np.array([[-1.0j, 0],
                        [0, 1.0j]])

  # Bob's strategies
  unitary1_b = np.array([[0, 1],
                        [-1, 0]]).astype(complex)
  unitary2_b = np.array([[0, -1.0j],
                        [-1.0j, 0]]).astype(complex)

  # define probabilities for Alice
  p1_a = np.cos(alpha_a / 2) ** 2
  p2_a = np.sin(alpha_a/ 2) ** 2

  # define probabilities for Bob
  p1_b = np.cos(alpha_b / 2) ** 2
  p2_b = np.sin(alpha_b/ 2) ** 2

  # define the set of actions and their probability distribution.
  mixed_strategy = [[(p1_a, unitary1_a), (p2_a, unitary2_a)], [(p1_b, unitary1_b), (p2_b, unitary2_b)]]

  identity = np.eye(2)

  J = 1/np.sqrt(2)*np.array([[1.0, 0.0, 0.0, 1.0j],
                             [0.0, 1.0, 1.0j, 0.0],
                             [0.0, 1.0j, 1.0, 0.0],
                             [1.0j, 0.0, 0.0, 1.0]]).astype(complex)

  J_unitary = UnitaryGate(Operator(J))
  qc.append(J_unitary, [0, 1])

  rho = qi.DensityMatrix.from_instruction(qc)

  for index, strategy in enumerate(mixed_strategy):
    rho = reduce(add,  (prob * rho.evolve(np.kron(*[strat if index == player else identity for player in range(num_players)])) for prob, strat in strategy))

  rho = rho.evolve(Operator(J.conj().T))
  return rho.probabilities()
  

payoff_matrix =[[[3, 0],
                [5, 1]],
                [[3, 5],
                [0, 1]]]

for alpha_a in np.linspace(0, np.pi, 4):
  for alpha_b in np.linspace(0, np.pi, 2):
    print("alpha_a = ", alpha_a)
    print("alpha_b = ", alpha_b)
    prob = mixed_strategy(alpha_a, alpha_b)
    print("prob = ", prob)
    alice_matrix = np.reshape(payoff_matrix[0], (1,4))
    bob_matrix = np.reshape(payoff_matrix[1], (1,4))
    print("Alice's Payoff = ", np.round(np.dot(prob, alice_matrix[0]), 3))
    print("Bob's Payoff = ", np.round(np.dot(prob, bob_matrix[0]),3))
    print("=======================")


# probability of 1/2 and 1/2 yeilds the payoffs of 2.5 for each player.
prob = mixed_strategy(np.pi/2, np.pi/2)
print("prob = ", prob)
alice_matrix = np.reshape(payoff_matrix[0], (1,4))
bob_matrix = np.reshape(payoff_matrix[1], (1,4))
print("Alice's Payoff = ", np.round(np.dot(prob, alice_matrix[0]), 3))
print("Bob's Payoff = ", np.round(np.dot(prob, bob_matrix[0]),3))