Constant_Depth_w_QSearch_TFIM_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant-depth circuits for simulating the transverse field Ising model using Qsearch\n",
    "\n",
    "In this notebook, we demonstrate how to generate the constant-depth circuits for dynamic simulation of the transverse field Ising model (TFIM) using the circuit synthesis software QSearch, which is part of the **[BQSKit](https://bqskit.lbl.gov/)** tookit.  More details on our constant-depth circuits can be found in our **[paper](https://arxiv.org/abs/2103.07429)**.  We also prepare analgous circuits using IBM's general purpose compiler and show how our circuits perform better on currently available quantum computers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import necessary libraries\n",
    "First we download the necessary libraries.  See instructions for downloading **[numpy](https://numpy.org/install/)**, **[matplotlib](https://matplotlib.org/stable/users/installing.html)**, **[QSearch](https://github.com/BQSKit/qsearch)**, or **[Qiskit](https://qiskit.org/documentation/install.html)** if you don't already have them installed on your compuer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import necessary libraries\n",
    "from scipy.stats import unitary_group\n",
    "import math\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "import qsearch\n",
    "from qsearch.gates import *\n",
    "from qsearch.unitaries import *\n",
    "from qsearch.assemblers import *\n",
    "from qsearch import multistart_solvers, utils, options, leap_compiler, post_processing, assemblers\n",
    "from qsearch.defaults import standard_defaults, standard_smart_defaults\n",
    "import itertools as it\n",
    "import qiskit as qk\n",
    "from qiskit.tools.monitor import job_monitor\n",
    "from qiskit import Aer, IBMQ, execute\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System parameters\n",
    "Next, we define the parameters of the system and the simulation we wish to run. Here, we will be performing a quantum quench of a TFIM.  Set the number of desired qubits; the coupling constant Jx; the strength and frequency of the transverse field hz and freq, respectively; the time-step size; the number of time-steps; and the number of shots for performing experiments on the quantum simulator or the real quantum hardware."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define simulation parameters\n",
    "N = 3 #number of spins/qubits\n",
    "#TFIM Simulation parameters\n",
    "Jx = 0.01183898\n",
    "hz = 2.0*Jx\n",
    "freq = 0.0048\n",
    "delta_t = 3 #time-step size\n",
    "num_steps = 40 #number of time-steps, we start with a small number for speed of demonstration\n",
    "shots = 8192 #number of shots for circuit execution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating the high-level circuits\n",
    "Next, we generate the high-level circuits for simulating the dynamics of the XY model.  These circuits will be used to create our constant-depth circuits, as well as used with the IBM general-purpose circuit compiler for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create the circuit to execute the time-evolution operator for a given time-step\n",
    "def evolution_circuit_TFIM(num_time_steps, Jx, hz, freq, delta_t, N):\n",
    "    hbar = 0.658212    # eV*fs\n",
    "    time_evol_circuit = qk.QuantumCircuit(N)\n",
    "    #define rotation angles for gates in circuit\n",
    "    psiX = -2.0*Jx*delta_t/hbar\n",
    "    for step in range(num_time_steps):\n",
    "        t = (step + 0.5) * delta_t\n",
    "        psiZ = -2.0*hz*np.cos(2*np.pi*freq*t)*delta_t/hbar\n",
    "        #implement XX operator\n",
    "        for q in range(0,N-1):\n",
    "            time_evol_circuit.h(q)\n",
    "            time_evol_circuit.h(q+1)\n",
    "            time_evol_circuit.cx(q,q+1)\n",
    "            time_evol_circuit.rz(psiX,q+1)\n",
    "            time_evol_circuit.cx(q,q+1)\n",
    "            time_evol_circuit.h(q)\n",
    "            time_evol_circuit.h(q+1)\n",
    "        #implement Z operator\n",
    "        for q in range(0,N):\n",
    "            time_evol_circuit.rz(psiZ, q)\n",
    "    return time_evol_circuit\n",
    "\n",
    "#draw circuit for time-step 1\n",
    "#circ = evolution_circuit_TFIM(1,Jx,hz,N)\n",
    "#circ.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create set of final circuits for quantum quench simulations\n",
    "circuits = []\n",
    "for i in range(0, num_steps+1):\n",
    "    #TFIM\n",
    "    circuits.append(evolution_circuit_TFIM(i,Jx,hz,freq,delta_t,N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions for generating the constant-depth circuit structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_matchgate():\n",
    "    rz = ZGate()\n",
    "    rx = XGate()\n",
    "    cnot = CNOTGate()\n",
    "    I = IdentityGate()\n",
    "\n",
    "    RZ_layer = KroneckerGate(rz,rz)\n",
    "    RX_layer = KroneckerGate(rx,rx)\n",
    "    RXRZ_layer = KroneckerGate(rx, rz)\n",
    "    matchgate = ProductGate(RZ_layer, RX_layer, cnot, RXRZ_layer, cnot, RX_layer, RZ_layer)\n",
    "    return matchgate\n",
    "    \n",
    "def make_layertype1(N):\n",
    "    I = IdentityGate()\n",
    "    for i in range(int(N/2)):\n",
    "        if (i==0):\n",
    "            layer = make_matchgate()\n",
    "        else:\n",
    "            layer = KroneckerGate(layer,make_matchgate())\n",
    "    if (N%2 != 0): \n",
    "        layer = KroneckerGate(layer, I)\n",
    "    return layer\n",
    "        \n",
    "def make_layertype2(N):\n",
    "    I = IdentityGate()\n",
    "    layer = I\n",
    "    #N even\n",
    "    if (N%2==0):\n",
    "        for _ in range(int(N/2)-1):\n",
    "            layer = KroneckerGate(layer,make_matchgate())\n",
    "        layer = KroneckerGate(layer,I)\n",
    "    #N odd\n",
    "    else:\n",
    "        for _ in range(int(N/2)):\n",
    "            layer = KroneckerGate(layer,make_matchgate())\n",
    "    return layer\n",
    "\n",
    "def make_MGC(N):\n",
    "    \"\"\"\n",
    "    Make matchgate circuit for N qubits.\n",
    "    Args:\n",
    "        N (int): The number of spins.\n",
    "    Returns\n",
    "        circuit (ProductGate): Circuit of QSearch gates.\n",
    "    \"\"\"\n",
    "    for l in range(N):\n",
    "        #add layer_type1\n",
    "        if (l%2 == 0):\n",
    "            if(l==0):\n",
    "                circuit = make_layertype2(N)\n",
    "            else:\n",
    "                circuit = ProductGate(circuit, make_layertype2(N))\n",
    "        #add layer_type2\n",
    "        else:\n",
    "            circuit = ProductGate(circuit, make_layertype1(N))\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate the constant-depth circuits\n",
    "The constant-depth circuits are generating by computing the target unitary matrix from the high-level circuits and the constant-depth circuit structure.  Given the target unitary and the constant-depth circuit structure, QSearch can optimize the parameters in the constant-depth circuit, which are written to a file for future use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make constant-depth circuits\n",
    "unitary_sim = Aer.get_backend('unitary_simulator')\n",
    "for t in range(num_steps+1):\n",
    "    #get target unitary for the given timestep\n",
    "    job = execute(circuits[t], unitary_sim)\n",
    "    result = job.result()\n",
    "    target_unitary = result.get_unitary(circuits[t], decimals=12)\n",
    "    #get constant-depth circuit structure for N qubits\n",
    "    circ_struct = make_MGC(N)\n",
    "    #optimize parameters of circuit\n",
    "    # use the multistart solver, may want to increase the number of starts for more qubits, but that will also be slower\n",
    "    solv = multistart_solvers.MultiStart_Solver(24)\n",
    "    # set up some options\n",
    "    opts = qsearch.Options()\n",
    "    opts.target = target_unitary\n",
    "    #opts.gateset = gateset\n",
    "    opts.set_defaults(**standard_defaults)\n",
    "    opts.set_smart_defaults(**standard_smart_defaults)\n",
    "    # optimize the circuit structure (circ_struct) for target U\n",
    "    # returns the calculated matrix and the vector of parameters\n",
    "    dist = 1\n",
    "    # run a few times to make sure we find the correct solution\n",
    "    for _ in range(5):\n",
    "        mat, vec = solv.solve_for_unitary(circ_struct, opts)\n",
    "        dist_new = utils.matrix_distance_squared(mat, target_unitary)\n",
    "        print(dist_new)\n",
    "        if dist_new < dist:\n",
    "            dist = dist_new\n",
    "        if dist < 1e-10:\n",
    "            break\n",
    "\n",
    "    print(f'For timestep {t} got distance {dist}')\n",
    "    #get final circuit\n",
    "\n",
    "    result_dict = {}\n",
    "    result_dict[\"structure\"] = circ_struct\n",
    "    result_dict[\"parameters\"] = vec\n",
    "    \n",
    "    opts.assemblydict=assemblydict_ibmopenqasm\n",
    "    out = opts.assembler.assemble(result_dict, opts)\n",
    "    #write results to file\n",
    "    with open(f'{N}q_TFIM_timestep{t}_qsearch.qasm', \"w\") as wfile:\n",
    "        wfile.write(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prepare constant depth circuits to run on quantum backends\n",
    "q_regs = qk.QuantumRegister(N, 'q')\n",
    "c_regs = qk.ClassicalRegister(N, 'c')\n",
    "cd_circuits = []\n",
    "for t in range(num_steps+1):\n",
    "    total_circ = qk.QuantumCircuit(q_regs, c_regs)\n",
    "    #initial state is ground state of XX\n",
    "    total_circ.h(q_regs)\n",
    "    total_circ.barrier()\n",
    "    #get constant-depth evolution circuit\n",
    "    circ = qk.QuantumCircuit.from_qasm_file(f'{N}q_TFIM_timestep{t}_qsearch.qasm')\n",
    "    total_circ.compose(circ, inplace=True)\n",
    "    #flip basis to measure X-mag\n",
    "    total_circ.h(q_regs)\n",
    "    total_circ.measure(q_regs, c_regs)\n",
    "    cd_circuits.append(total_circ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate the IBM-compiled circuits\n",
    "For comparison, we generate circuits using IBM's general-purpose quantum circuit compiler.  These circuits will not be constant in depth with increasing numbers of time-steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For comparison, create circuits using IBM compiler\n",
    "q_regs = qk.QuantumRegister(N, 'q')\n",
    "c_regs = qk.ClassicalRegister(N, 'c')\n",
    "ibm_circuits = []\n",
    "for t in range(num_steps+1):\n",
    "    total_circ = qk.QuantumCircuit(q_regs, c_regs)\n",
    "    #initial state is ground state of XX\n",
    "    total_circ.h(q_regs)\n",
    "    total_circ.barrier()\n",
    "    #add evolution circuit\n",
    "    total_circ.compose(circuits[t], inplace=True)\n",
    "    #flip basis to measure X-mag\n",
    "    total_circ.h(q_regs)\n",
    "    total_circ.measure(q_regs, c_regs)\n",
    "    ibm_circuits.append(total_circ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare sizes of the consant-depth and IBM-compiled circuits\n",
    "Chose any time-step t to see how the volume of our constant-depth circuits compare to those generated by IBM's general-purpose compiler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare lengths of IBM vs constant-depth circuits for time-step t\n",
    "t = 9\n",
    "print(len(circuits[t]))\n",
    "print(len(cd_circuits[t]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connect to IBM quantum backends over the cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make connection to IBM to run constant-depth circuits on quantum backend\n",
    "#!! If this is your first time running this notebook get you IBM accounts's API key and save your account!!\n",
    "#qk.IBMQ.save_account('your_API_key_here')\n",
    "qk.IBMQ.load_account()\n",
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "\n",
    "#check for available quantum processor backend\n",
    "provider = qk.IBMQ.get_provider(hub='ibm-q', group='open', project='main')\n",
    "provider.backends()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select one of the available quantum backend\n",
    "backend = provider.get_backend('ibmq_athens')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transpile all circuits for the target quantum backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compile constant-depth circuits for the simulator\n",
    "cd_circs_sim = qk.transpile(cd_circuits, backend=simulator, optimization_level=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compile constant-depth circuits for real quantum backend\n",
    "cd_circs_qp = qk.transpile(cd_circuits, backend=backend, optimization_level=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compile ibm circuits for the simulator\n",
    "ibm_circs_sim = qk.transpile(ibm_circuits, backend=simulator, optimization_level=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compile ibm circuits for real quantum backend\n",
    "ibm_circs_qp = qk.transpile(ibm_circuits, backend=backend, optimization_level=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execute all circuits on desired qunatum backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Run constant-depth circuits on simulator\n",
    "cd_sim_results = execute(cd_circs_sim, simulator, shots=shots).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Run constant-depth circuits on quantum processor\n",
    "cd_job = qk.execute(cd_circs_qp, backend=backend, shots=shots)\n",
    "job_monitor(cd_job)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For comparison run IBM circuits on simulator\n",
    "ibm_sim_results = execute(ibm_circs_sim, simulator, shots=shots).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For comparison run IBM circuits on real quantum processor\n",
    "ibm_job = qk.execute(ibm_circs_qp, backend=backend, shots=shots)\n",
    "job_monitor(ibm_job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Post-processing\n",
    "To post-process the results from the quantum backends, we define a function to compute the average magnetization for each time-step from the results of the quantum backend.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define post-processing function\n",
    "def average_magnetization(result: dict, shots: int):\n",
    "    \"\"\"Compute average magnetization from results of qk.execution.\n",
    "    Args:\n",
    "    - result (dict): a dictionary with the counts for each qubit, see qk.result.result module\n",
    "    - shots (int): number of trials\n",
    "    Return:\n",
    "    - average_mag (float)\n",
    "      \"\"\"\n",
    "    mag = 0\n",
    "    for spin_str, count in result.items():\n",
    "        spin_int = [1 - 2 * float(s) for s in spin_str]\n",
    "        mag += (sum(spin_int) / len(spin_int)) * count\n",
    "    average_mag = mag / shots\n",
    "    return average_mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Post-processing for constant-depth circuit results from simulator\n",
    "avg_mag_cd_sim = []\n",
    "for c in cd_circs_sim:\n",
    "    result_dict = cd_sim_results.get_counts(c)\n",
    "    avg_mag_cd_sim.append(average_magnetization(result_dict, shots))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Post-processing for constant-depth circuit results from real quantum processor\n",
    "cd_results = cd_job.result()        \n",
    "\n",
    "avg_mag_cd_qp = []\n",
    "for c in cd_circs_qp:\n",
    "    result_dict = cd_results.get_counts(c)\n",
    "    avg_mag_cd_qp.append(average_magnetization(result_dict, shots))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Post-processing for IBM circuit results from simulator\n",
    "avg_mag_ibm_sim = []\n",
    "for c in ibm_circs_sim:\n",
    "    result_dict = ibm_sim_results.get_counts(c)\n",
    "    avg_mag_ibm_sim.append(average_magnetization(result_dict, shots))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Post-processing for IBM circuit results from real quantum processor \n",
    "ibm_results = ibm_job.result()        \n",
    "avg_mag_ibm_qp = []\n",
    "for c in ibm_circs_qp:\n",
    "    result_dict = ibm_results.get_counts(c)\n",
    "    avg_mag_ibm_qp.append(average_magnetization(result_dict, shots))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error Mitigation\n",
    "We only perform minimal error mitigation to combat read-out error noise in the results from the real quantum hardware.  There will be no such noise on the quantum simulator as we chose to use noise-free simulators.  Given that we know the magnetization result for the initial time-step should be equal to 1, we scale all subsequent results by the initial magnetization measured in the first time-step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#error mitigation of readout-noise\n",
    "#scaled all results by result from first time-step\n",
    "cd_scaled = np.asarray(avg_mag_cd_qp)/avg_mag_cd_qp[0]\n",
    "ibm_scaled = np.asarray(avg_mag_ibm_qp)/avg_mag_ibm_qp[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,4))\n",
    "plt.plot(avg_mag_ibm_sim, label=\"Noise Free Simulation\", color=\"blue\")\n",
    "plt.plot(ibm_scaled, label=\"IBM generated HW\", color=\"green\")\n",
    "plt.plot(cd_scaled, label=\"Constant-Depth HW\", color=\"red\")\n",
    "plt.plot(avg_mag_cd_sim, label=\"Constant-Depth Simulation\", color=\"orange\")\n",
    "#plt.legend()\n",
    "plt.xlim(0,39)\n",
    "plt.ylim(-1.0, 1.0)\n",
    "plt.xlabel(\"Simulation Timestep\", fontsize=18, fontname=\"Times\")\n",
    "plt.yticks([-1.0, -0.5, 0, 0.5, 1.0],fontsize=18, fontname=\"Times\")\n",
    "plt.xticks(fontsize=18, fontname=\"Times\")\n",
    "plt.ylabel(\"Average Magnetization\", fontsize=18, fontname=\"Times\")\n",
    "plt.show()\n",
    "#plt.savefig(f'{N}q_TFIM_cd_v_ibm.svg', dpi=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optionally write results from the quantum backend to file for future use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write all results to file for future use\n",
    "with open(f'{N}q_TFIM_ibm_hardware.txt', 'w') as f:\n",
    "    for item in avg_mag_ibm_qp:\n",
    "        f.write(\"%s\\n\" % item)\n",
    "with open(f'{N}q_TFIM_cd_hardware.txt', 'w') as f:\n",
    "    for item in avg_mag_cd_qp:\n",
    "        f.write(\"%s\\n\" % item)\n",
    "with open(f'{N}q_TFIM_ibm_sim.txt', 'w') as f:\n",
    "    for item in avg_mag_ibm_sim:\n",
    "        f.write(\"%s\\n\" % item)\n",
    "with open(f'{N}q_TFIM_cd_sim.txt', 'w') as f:\n",
    "    for item in avg_mag_cd_sim:\n",
    "        f.write(\"%s\\n\" % item)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}