{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "# Resource estimation for simulating a 2D Hubbard model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "In this Python + Q# notebook we demonstrate how to *efficiently* estimate the resources required to simulate a Hubbard model Hamiltonian on an $N \\times N$ 2D\n",
    "lattice using a *fourth-order Trotter Suzuki product formula* assuming a 2D planar qubit architecture with nearest-neighbor connectivity.\n",
    "\n",
    "First, load the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "import qsharp\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Background: 2D Hubbard Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "The Hubbard Model is a simple model of interacting particles in a lattice (in our case a 2D square lattice) with two kinds of terms in the Hamiltonian: (i) a hopping term between adjacent sites and (ii) a potential term for onsite interactions.\n",
    "Formally, a Hubbard model Hamiltonian on an $N \\times N$ lattice where each site has two spins {up, dn} is given as:\n",
    "\n",
    "$$\n",
    "H = \\underbrace{U \\sum_{i,j} c^{\\dagger}_{(i, j), \\text{up}} c_{(i, j), \\text{up}}  c^{\\dagger}_{(i, j), \\text{dn}} c_{(i, j), \\text{dn}}}_{B} - \\underbrace{t \\sum_{\\sigma \\in \\{ \\text{up}, \\text{dn}\\} , i, j} ( c^\\dagger_{(i,j),\\sigma} c_{(i+1,j),\\sigma} + c^\\dagger_{(i+1,j),\\sigma} c_{(i,j),\\sigma})}_{A} - \\underbrace{t \\sum_{\\sigma \\in \\{ \\text{up}, \\text{dn}\\} , i, j} ( c^\\dagger_{(i,j),\\sigma} c_{(i,j+1),\\sigma} + c^\\dagger_{(i,j+1),\\sigma} c_{(i,j),\\sigma})}_{C}\n",
    "$$\n",
    "where $c$ is the annihilation operator, $c^\\dagger$ is the creation operator, the first set of terms are the on-site repulsion terms, the second set are the vertical row hops and last set is the horizontal column hops."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "### Converting to Pauli terms\n",
    "\n",
    "The first step needed to simulate $H$ is to map the fermionic terms into *local* Pauli operations involving only a few qubits. When the Hamiltoinan acts on a 1D lattice, using the Jordan-Wigner would suffice. However, for the 2D lattice, either row or column hops will end up being non-local terms involving $O(N)$ qubits.\n",
    "To overcome this, we use the [Verstraete-Cirac](http://arxiv.org/abs/cond-mat/0508353v3) encoding to convert the repulsion and hopping terms into local Pauli operators involving at most 4 qubits at a time. This is done by adding one auxiliary qubit for each spin type in each lattice site such that the $N^2$ lattice sites are mapped to $4 N^2$ qubits. Specifically, each lattice site $m$ maps to the qubits $\\{u_m, d_m, \\tilde{u}_m \\tilde{d}_m\\}$ where the first two qubits correspond to the up and down spin types and the latter two to the corresponding auxiliary qubits.\n",
    "\n",
    "This converts the terms as follows:\n",
    "$$\n",
    "\\text{Repulsion Terms:} \\qquad   c^{\\dagger}_{(i, j), \\text{up}} c_{(i, j), \\text{up}}  c^{\\dagger}_{(i, j), \\text{dn}} c_{(i, j), \\text{dn}} \\mapsto \\frac{1}{4} \\left(\\mathbb{I} - Z_{(i,j), \\text{up}} - Z_{(i,j), \\text{dn}} + Z_{(i,j), \\text{up}} \\otimes Z_{(i,j), \\text{dn}}\\right)\n",
    "$$\n",
    "$$\n",
    "\\text{Horizontal Hopping Terms:} \\qquad c^\\dagger_{(i,j),\\sigma} c_{(i,j+1),\\sigma} + c^\\dagger_{(i,j+1),\\sigma} c_{(i,j),\\sigma} \\mapsto \\frac{1}{2} \\left(X_{(i,j),\\sigma} X_{(i,j+1),\\sigma} + Y_{(i,j),\\sigma} Y_{(i,j+1),\\sigma} \\right) \\tilde{Z}_{(i,j),\\sigma}\n",
    "$$\n",
    "$$\n",
    "\\text{Vertical Hopping Terms:} \\qquad c^\\dagger_{(i,j),\\sigma} c_{(i+1,j),\\sigma} + c^\\dagger_{(i+1,j),\\sigma} c_{(i,j),\\sigma} \\mapsto \\frac{(-1)^{i+1}}{2}  \\begin{cases} \\left(X_{(i,j),\\sigma} Y_{(i+1,j), \\sigma} - Y_{(i,j),\\sigma} X_{(i+1,j),\\sigma} \\right) \\tilde{X}_{(i,j),\\sigma} \\tilde{Y}_{(i+1,j), \\sigma} & (j \\text{ is even}); \\\\\n",
    "\\left(X_{(i,j),\\sigma} Y_{(i+1,j), \\sigma} -  Y_{(i,j),\\sigma} X_{(i+1,j),\\sigma} \\right) \\tilde{Y}_{(i,j),\\sigma} \\tilde{X}_{(i+1,j), \\sigma}  & (j \\text{ is odd}). \\end{cases}\n",
    "$$\n",
    "where the $\\tilde{(\\cdot)}$ terms are Paulis acting on the auxiliary system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "### Simulation\n",
    "\n",
    "Now that the Hamiltonian has been converted into $H = \\sum_k h_k$ where each $h_k$ is a Local Hamiltonian term, we can consider how to implement $e^{-iHT}$ using the fourth order Trotter-Suzuki product formula. For this, we simulate the evolution for small slices of time $\\Delta$ and repeating this for `nSteps` $= \\lceil T/\\Delta \\rceil$ to obtain the full time evolution. The Trotter-Suzuki formula for higher orders can be recursively defined using a *fractal decomposition* as discussed in Section 3 of [Hatanao and Suziki's survey](https://link.springer.com/chapter/10.1007/11526216_2). Then the fourth order formula $U_4(\\Delta)$ can be constructed using the second-order one $U_2(\\Delta)$ as follows.\n",
    "$$\n",
    "\\begin{aligned}\n",
    "U_2(\\Delta) & = e^{-iA\\Delta/2} e^{-iB\\Delta/2} e^{-iC\\Delta} e^{-iB\\Delta/2} e^{-iA\\Delta/2}; \\\\\n",
    "U_4(\\Delta) & = U_2(p\\Delta)U_2(p\\Delta)U_2((1 - 4p)\\Delta)U_2(p\\Delta)U_2(p\\Delta); \\\\\n",
    "p & = (4 - 4^{1/3})^{-1}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "For the rest of the notebook, we will present the code that computes the time evolution in a step by step fashion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Implementation\n",
    "\n",
    "### Helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "We will allocate all qubits in the 2D lattice in a one-dimensional array where every lattice point has 4 qubits ordered as $(u, u', d, d')$. Further we assume a snake-like order on the 2D lattice i.e., the ordering goes left-to-right on even rows and right-to-left on odd rows. The function `GetQubitIndex` converts a qubit identified by its 2D lattice point coordinates, `(row, col)`, its spin flag `f_spin` and auxiliary flag `f_aux` to an index in the qubit array. Note that `f_spin = True` for the down spin and `f_aux = True` for auxiliary qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "function GetQubitIndex(k : Int, l : Int, f_spin : Bool, f_aux : Bool, n : Int, m : Int) : Int {\n",
    "\n",
    "    let sshift = f_spin ? 2 | 0;\n",
    "    let ashift = f_aux ? 1 | 0;\n",
    "\n",
    "    return k % 2 == 0    //Check if row is even and return expression accordingly\n",
    "        ? 4*m*k + 4*l + ashift + sshift\n",
    "        | 4*m*k + 4*(m-l-1) + (ashift) + (sshift);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Expanding $U_4(\\Delta)$ to express it in terms of $A, B, C$ gives:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "U_4(\\Delta) & = e^{-iAp\\Delta/2} \\ T_1 \\ T_2 \\ T_1 \\ e^{-iAp\\Delta/2}, \\\\\n",
    "\\text{where } T_1 & = e^{-iBp\\Delta/2} e^{-iCp\\Delta} e^{-iBp\\Delta/2} e^{-iAp\\Delta} e^{-iBp\\Delta/2} e^{-iCp\\Delta} e^{-iBp\\Delta/2}, \\\\\n",
    "\\text{and } T_2 & = e^{-iA(1 - 3p)\\Delta/2} e^{-iB(1-4p)\\Delta/2} e^{-iC(1-4p)\\Delta} e^{-iB(1-4p)\\Delta/2} e^{-iA(1 - 3p)\\Delta/2}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The above equation has $21$ exponential terms and works for one time step. For `nSteps` $> 1$ time steps, some adjacent terms can be merged to give $20\\lceil T/\\Delta \\rceil+1$ exponential terms for time evolution $e^{-iHT}$.\n",
    "\n",
    "The function below sets the rotation angles used to apply the exponentials involving $A, B$ and $C$ in the formula above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "function SetAngleSequence(p : Double, hop_angle : Double, rep_angle : Double, dt : Double) : Double[] {\n",
    "    let len1 = 8;\n",
    "    let len2 = 4;\n",
    "    let valLength = 2*len1 + len2;\n",
    "\n",
    "    let val1 = -p*hop_angle*dt;\n",
    "    let val2 = -p*rep_angle*dt/2.0;\n",
    "\n",
    "    let valB = -(1.0-4.0*p)*rep_angle*dt/2.0;\n",
    "    let valA = -(1.0-3.0*p)*hop_angle*dt/2.0;\n",
    "    let valC = -(1.0-4.0*p)*hop_angle*dt;\n",
    "\n",
    "    mutable values = [0.0, size=valLength];\n",
    "\n",
    "    // Angles bookending Term2\n",
    "    set values w/= len1 <- valA;\n",
    "    set values w/= len1+len2 <- valA;\n",
    "\n",
    "    for i in 0..len1-1 {\n",
    "        if (i%2 == 0) {\n",
    "            set values w/= i <- val1;\n",
    "        }\n",
    "        else {\n",
    "            set values w/= i <- val2;\n",
    "        }\n",
    "    }\n",
    "\n",
    "    for i in len1+1..len1+len2-1 {\n",
    "        if (i%2 == 0) {\n",
    "            set values w/= i <- valC;\n",
    "        }\n",
    "        else {\n",
    "            set values w/= i <- valB;\n",
    "        }\n",
    "    }\n",
    "\n",
    "    for i in len1+len2+1..valLength-1 {\n",
    "        if (i%2 == 0) {\n",
    "            set values w/= i <- val1;\n",
    "        }\n",
    "        else {\n",
    "            set values w/= i <- val2;\n",
    "        }\n",
    "    }\n",
    "\n",
    "    return values;\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "### Quantum operations\n",
    "\n",
    "The first term we want to apply are the horizontal hops. For a given set of qubits and angle $\\theta$, the following operation applies the horizontal hop's evolution operator $e^{i(XZX+YZY)\\theta}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyHorizontalTerm (qslice : Qubit[], theta : Double) : Unit {\n",
    "\n",
    "    let P_x = [PauliX, PauliZ, PauliX];\n",
    "    let P_y = [PauliY, PauliZ, PauliY];\n",
    "\n",
    "    Exp(P_x, theta, qslice);\n",
    "    Exp(P_y, theta, qslice);\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Observe that it is not possible to simultaneously apply all the rotations in one go. For example, while applying the rotation on the up spin qubits at $(0, 0)$ and $(0, 1)$, it is not possible to also apply the rotation on the up spin qubits $(0, 1)$ and $(0, 2)$. Instead, we try to apply as many rotations as possible while also giving consideration to some uniformity with respect to the spins. The operation is broken up as follows:\n",
    "- consider pairs where the first qubit is in an odd (resp. even) column as given by `grp`;\n",
    "- choose the spin up (resp. down) qubits in the array as chosen by `f_spin`,\n",
    "- apply the exponential to all such sets of qubits in the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyHorizontalHops(n : Int, m : Int, g_num : Bool, f_spin : Bool, qs : Qubit[], theta : Double) : Unit {\n",
    "\n",
    "    let start = g_num ? 1 | 0; //Gives the group number for horizontal hop terms\n",
    "\n",
    "    for row in 0..n-1 { // for all rows\n",
    "        for col in start..2..m-2 { // for alternate columns\n",
    "            //Message($\"Creating Pauli expr for Hop terms ({row}, {col}) <--> ({row}, {col+1})\");\n",
    "\n",
    "            let aux_col = (row%2 == 0)? col | col+1;\n",
    "            let i1 = GetQubitIndex(row, col, f_spin, false, n, m);\n",
    "            let aux = GetQubitIndex(row, aux_col, f_spin, true, n, m);\n",
    "            let i2 = GetQubitIndex(row, col+1, f_spin, false, n, m);\n",
    "\n",
    "            ApplyHorizontalTerm([qs[i1], qs[aux], qs[i2]], theta);\n",
    "\n",
    "        } //end for col\n",
    "    } //end for row\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Next, we consider the vertical hopping terms. For a given set of qubits and angle $\\theta$, the following operation applies $e^{i (XXYY - YXXY) \\theta}$ for even columns and $e^{i (XYYX - YYXX) \\theta}$ for odd columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyVerticalTerm(qslice : Qubit[], col : Int, angle : Double) : Unit {\n",
    "\n",
    "    let P_xxyy = [PauliX, PauliX, PauliY, PauliY];\n",
    "    let P_yxxy = [PauliY, PauliX, PauliX, PauliY];\n",
    "    let P_xyyx = [PauliX, PauliY, PauliY, PauliX];\n",
    "    let P_yyxx = [PauliY, PauliY, PauliX, PauliX];\n",
    "\n",
    "    if col%2 == 0 {\n",
    "        Exp(P_xxyy, angle, qslice);\n",
    "        Exp(P_yxxy, -angle, qslice);\n",
    "    }\n",
    "    else {\n",
    "        Exp(P_xyyx, angle, qslice);\n",
    "        Exp(P_yyxx, -angle, qslice);\n",
    "    }\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Like with the horizontal hops, it is not possible to simultaneously apply all the rotations in one go. For example, while applying the rotation on the up spin qubits at $(0, 0)$ and $(1, 0)$, it is not possible to also apply the rotation on the up spin qubits $(1, 0)$ and $(2, 0)$. The operation is accordingly broken up as follows:\n",
    "- consider pairs where the first qubit is in an odd (resp. even) *row* as given by `grp`;\n",
    "- choose the spin up (resp. down) qubits in the array as chosen by `f_spin`,\n",
    "- apply the exponential to all such sets of qubits in the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyVerticalHops(row : Int, n : Int, m : Int, f_spin : Bool, qs : Qubit[], theta : Double, grp : Bool) : Unit { //\n",
    "\n",
    "    let start = grp ? 1 | 0;\n",
    "    for i in start..2..m-1 { // for each column // 0..m-1 {\n",
    "        let angle = i%2 == 0 ? theta | -theta;\n",
    "\n",
    "        let i1_u = GetQubitIndex(row, i, f_spin, false, n, m);\n",
    "        let i1_ua = GetQubitIndex(row, i, f_spin, true, n, m);\n",
    "        let i2_u = GetQubitIndex(row+1, i, f_spin, false, n, m);\n",
    "        let i2_ua = GetQubitIndex(row+1, i, f_spin, true, n, m);\n",
    "\n",
    "        ApplyVerticalTerm([qs[i1_u], qs[i1_ua], qs[i2_u], qs[i2_ua]], i, angle);\n",
    "\n",
    "    } // end for i loop\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Finally, there are two kinds of Pauli exponentials needed for applying the repulsion term for an angle $\\theta$:\n",
    "- an $e^{iZ\\theta}$ applied to each up and down spin qubit;\n",
    "- $e^{-i (Z \\otimes Z)\\theta}$ applied to pairs of $(u, d)$ qubits for each point in the lattice.\n",
    "\n",
    "The operation below applies $e^{iZ\\theta}$ on the up and down spin qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyZRot(n : Int, m : Int, qs : Qubit[], angle : Double) : Unit {\n",
    "\n",
    "    for row in 0..n-1 {\n",
    "        for col in 0..m-1 {\n",
    "            let up = GetQubitIndex(row, col, false, false, n, m);\n",
    "            let down = GetQubitIndex(row, col, true, false, n, m);\n",
    "\n",
    "            // Rz(θ, q) applies e^{-iZθ/2} to qubit q\n",
    "            Rz(-2.0*angle, qs[up]);\n",
    "            Rz(-2.0*angle, qs[down]);\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "The next operation applies $e^{-i(ZZ)\\theta}$ for every $(u, d)$ pair of qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyZoZRot(n : Int, m : Int, qs : Qubit[], angle : Double) : Unit {\n",
    "\n",
    "    let P_op = [PauliZ, PauliZ];\n",
    "\n",
    "    for row in 0..n-1 {\n",
    "        for col in 0..m-1 {\n",
    "            let up = GetQubitIndex(row, col, false, false, n, m);\n",
    "            let down = GetQubitIndex(row, col, true, false, n, m);\n",
    "\n",
    "            Exp(P_op, angle, [qs[up], qs[down]]);\n",
    "\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "The next operation puts everything together and calls the operations needed to\n",
    "simulate the Hubbard model Hamiltonian using a fourth order product formula.\n",
    "Observe that the operations `ApplyHorizontalHops` and `ApplyVerticalHops` are called four times for for different\n",
    "choices of direction and starting index to ensure all possible pairs of qubits are appropriately considered.\n",
    "\n",
    "The various parameters taken in by the operation correspond to:\n",
    "\n",
    "- `N1`, `N2`: row and column size of the lattice.\n",
    "- `t, u`: parameters by which the Hamiltonian terms are scaled.\n",
    "- `totTime`: the total time for the Trotter simulation.\n",
    "- `dt` : the step size for the simulation, sometimes denoted as $\\Delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "microsoft": {
     "language": "qsharp"
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "import Std.Math.*;\n",
    "\n",
    "operation HubModel2DSim(N1 : Int, N2 : Int, t: Double, u : Double, totTime : Double, dt : Double) : Unit {\n",
    "\n",
    "    //Hubbard Hamiltonian: -t*Hopping terms + u*Repulsion Terms\n",
    "\n",
    "    //Represent a 2D array of qubits with 4 qubits per lattice point laid out as a 4N1 x N2 array\n",
    "    use qs = Qubit[4*N1*N2];\n",
    "\n",
    "    let p = 1.0/(4.0 - 4.0^(1.0/3.0));\n",
    "\n",
    "    let nsteps = Ceiling(totTime/dt);\n",
    "\n",
    "    let seqLen = 20*nsteps+1;\n",
    "\n",
    "    // Calculate rotation angles\n",
    "    let hop_angle = -t/2.0;\n",
    "    let rep_angle = u/4.0;\n",
    "\n",
    "    let angSeq = SetAngleSequence(p, hop_angle, rep_angle, dt);\n",
    "\n",
    "    for i in 0..seqLen-1 {\n",
    "        let theta = (i==0 or i==seqLen-1) ? -p*hop_angle*dt/2.0 | angSeq[i%20];\n",
    "\n",
    "        if (i%4 == 0) { // A = Vertical Hops\n",
    "            for row in 0..N1-2 {\n",
    "                for(f_spin, grp) in [(false, false), (false, true), (true, false), (true, true)] {\n",
    "                    ApplyVerticalHops(row, N1, N2, f_spin, qs, theta, grp);\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "        elif (i%2 == 1) { // B = Repulsion Term\n",
    "            ApplyZRot(N1, N2, qs, theta);\n",
    "            ApplyZoZRot(N1, N2, qs, theta);\n",
    "        }\n",
    "        else { // C = Horizontal Hops\n",
    "            for(f_spin, grp) in [(false, false), (false, true), (true, false), (true, true)] {\n",
    "                ApplyHorizontalHops(N1, N2, grp, f_spin, qs, theta);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    //Resetting qubits before release;\n",
    "    ResetAll(qs);\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for typical choices of `totTime` and `dt` used in real-world applications, the `seqLen` can range in the hundreds or thousands of terms each of whose resources will be individually estimated. This results in a lot of time consumed (> 15 min) to compute the resource estimates. However, from the expansion of $U_4(\\Delta)$, it is clear that there are a few unique terms that repeat `nsteps` number of times. Hence, to efficiently estimate resources, we compute the resources for running one step of $U_4(\\Delta)$ and multiply the cost by `nsteps` repetitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "import Std.ResourceEstimation.*;\n",
    "\n",
    "operation EstimateHubModel2DSim(N1 : Int, N2 : Int, t: Double, u : Double, totTime : Double, dt : Double) : Unit {\n",
    "\n",
    "    //Represent a 2D array of qubits with 4 qubits per lattice point laid out as a 4N1 x N2 array\n",
    "    use qs = Qubit[4*N1*N2];\n",
    "\n",
    "    let p = 1.0/(4.0 - 4.0^(1.0/3.0));\n",
    "\n",
    "    let nsteps = Ceiling(totTime/dt);\n",
    "\n",
    "    // Calculate rotation angles\n",
    "    let hop_angle = -t/2.0;\n",
    "    let rep_angle = u/4.0;\n",
    "\n",
    "    let angSeq = SetAngleSequence(p, hop_angle, rep_angle, dt);\n",
    "\n",
    "    within {\n",
    "        RepeatEstimates(nsteps);\n",
    "    }\n",
    "    apply {\n",
    "        for i in 0..19 {\n",
    "            let theta = angSeq[i%20];\n",
    "\n",
    "            if (i%4 == 0) { // A = Vertical Hops\n",
    "                for row in 0..N1-2 {\n",
    "                    for(f_spin, grp) in [(false, false), (false, true), (true, false), (true, true)] {\n",
    "                        ApplyVerticalHops(row, N1, N2, f_spin, qs, theta, grp);\n",
    "                    }\n",
    "                }\n",
    "            }\n",
    "            elif (i%2 == 1) { // B = Repulsion Term\n",
    "                ApplyZRot(N1, N2, qs, theta);\n",
    "                ApplyZoZRot(N1, N2, qs, theta);\n",
    "            }\n",
    "            else { // C = Horizontal Hops\n",
    "                for(f_spin, grp) in [(false, false), (false, true), (true, false), (true, true)] {\n",
    "                    ApplyHorizontalHops(N1, N2, grp, f_spin, qs, theta);\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    let angle = -p*hop_angle*dt/2.0;\n",
    "    for row in 0..N1-2 {\n",
    "        for(f_spin, grp) in [(false, false), (false, true), (true, false), (true, true)] {\n",
    "            ApplyVerticalHops(row, N1, N2, f_spin, qs, angle, grp);\n",
    "        }\n",
    "    }\n",
    "\n",
    "    //Resetting qubits before release;\n",
    "    ResetAll(qs);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Getting logical resource counts\n",
    "\n",
    "For the purpose of generating the rQOPS for some target runtime, it suffices that we obtain the logical resource estimates to simulate the Heisenberg model Hamiltonian. We consider three problem instances with lattice sizes $\\{20 \\times 20, 30 \\times 30, 40 \\times 40\\}$ with $t = u = 1.0$. These instances are simulated for a total time of $1000$ s, with step size `dt`${}= 0.1$, and overall probability of failure $\\varepsilon = 0.01$. Any one of the six pre-defined qubit parameters will do to obtain the logical coounts and in this notebook we choose a Majorana based qubit with the `floquet code`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up problem parameters\n",
    "N1 = [20, 30, 40]\n",
    "N2 = [20, 30, 40]\n",
    "t = 1.0\n",
    "u = 1.0\n",
    "totTime = 1000.0\n",
    "dt = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We submit a resource estimation job with all the problem instances sequentially and collect the estimates in `results`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submit jobs\n",
    "results = []\n",
    "\n",
    "for i in range(3):\n",
    "    qsharp_string = f\"EstimateHubModel2DSim({N1[i]}, {N2[i]}, {t}, {u}, {totTime}, {dt})\"\n",
    "\n",
    "    result = qsharp.estimate(qsharp_string, params={\"errorBudget\": 0.01, \"qubitParams\": {\"name\": \"qubit_maj_ns_e6\"}, \"qecScheme\": {\"name\": \"floquet_code\"}, \"constraints\": {\"logicalDepthFactor\": 4}})\n",
    "    results.append(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the complete information provided when invoking the resource estimator, we output the result for the $20 \\times 20$ lattice by displaying `results[0]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Displaying estimates for 20x20 lattice size.\n",
    "results[0]\n",
    "# Change index to 1 (resp. 2) for 30x30 (resp. 40x40) lattice size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Visualizing and understanding the results\n",
    "\n",
    "### Result summary table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "# Define function to display information in summary format\n",
    "def get_summary_table(results, labels):\n",
    "    logical_qubits = []\n",
    "    logical_depth = []\n",
    "    t_states = []\n",
    "    code_distance = []\n",
    "    t_factories = []\n",
    "    t_factory_fraction = []\n",
    "    physical_qubits = []\n",
    "    rqops = []\n",
    "    runtime = []\n",
    "    logical_error = []\n",
    "\n",
    "    for i in range(3):\n",
    "        logical_qubits.append(results[i]['physicalCounts']['breakdown']['algorithmicLogicalQubits'])\n",
    "        logical_depth.append(results[i]['physicalCountsFormatted']['logicalDepth'])\n",
    "        t_states.append(results[i]['physicalCountsFormatted']['numTstates'])\n",
    "        t_factories.append(results[i]['physicalCounts']['breakdown']['numTfactories'])\n",
    "        logical_error.append(results[i]['physicalCountsFormatted']['requiredLogicalQubitErrorRate'])\n",
    "        physical_qubits.append(results[i]['physicalCountsFormatted']['physicalQubits'])\n",
    "        rqops.append(results[i]['physicalCountsFormatted']['rqops'])\n",
    "        runtime.append(results[i]['physicalCountsFormatted']['runtime'])\n",
    "        code_distance.append(results[i]['logicalQubit']['codeDistance'])\n",
    "        t_factory_fraction.append(results[i]['physicalCountsFormatted']['physicalQubitsForTfactoriesPercentage'])\n",
    "\n",
    "    data = pd.DataFrame()\n",
    "    pd.options.display.float_format = '{:.2E}'.format\n",
    "    data['Logical qubits'] = logical_qubits\n",
    "    data['Logical depth'] = logical_depth\n",
    "    data['Logical error'] = logical_error\n",
    "    data['T states'] = t_states\n",
    "    # data['T states'] = data['T states'].astype('float64')\n",
    "    data['Code Distance'] = code_distance\n",
    "    data['T factories'] = t_factories\n",
    "    data['T factory fraction'] = t_factory_fraction\n",
    "    data['Physical qubits'] = physical_qubits\n",
    "    data['rQOPS'] = rqops\n",
    "    data['Physical runtime'] = runtime\n",
    "    data.index = labels\n",
    "\n",
    "    return data\n",
    "\n",
    "# Display summarized information for all problem instances\n",
    "labels = [\"Hub20\", \"Hub30\", \"Hub40\"]\n",
    "table = get_summary_table(results, labels)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that in general, there is a trade-off between the logical depth and number of T factories used. \n",
    "\n",
    "> To ensure that T factories do not dominate the resource requirements, we set the `logical_depth_factor`${}=4$ adding some number of `noops` to increase the logical depth."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting the target rQOPS\n",
    "\n",
    "While the resource estimator generates a runtime for the given hardware profile, we want to set a target runtime of 1 week i.e., 604800 seconds to obtain a practical quantum advantage. We collect the previous results to compute the corresponding target rQOPS as \n",
    "$$ \\text{Target rQOPS} = \\frac{\\text{Logical qubits}\\cdot\\text{Logical Depth}}{\\text{Target runtime}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_target_rqops(results, labels):\n",
    "\n",
    "    target_runtime = 604800\n",
    "    logical_qubits = []\n",
    "    logical_depth = []\n",
    "    target_rqops = []\n",
    "    logical_error = []\n",
    "\n",
    "    for i in range(3):\n",
    "        logical_qubits.append(results[i]['physicalCounts']['breakdown']['algorithmicLogicalQubits'])\n",
    "        logical_depth.append(results[i]['physicalCountsFormatted']['logicalDepth'])\n",
    "        logical_error.append(results[i]['physicalCountsFormatted']['requiredLogicalQubitErrorRate'])\n",
    "        target_rqops.append(round(results[i]['physicalCounts']['breakdown']['algorithmicLogicalQubits'] * results[i]['physicalCounts']['breakdown']['logicalDepth'] / target_runtime))\n",
    "\n",
    "    data = pd.DataFrame()\n",
    "    pd.options.display.float_format = '{:.2e}'.format\n",
    "    data['Logical qubits'] = logical_qubits\n",
    "    data['Logical depth'] = logical_depth\n",
    "    data['Logical error'] = logical_error\n",
    "    data['Target rQOPS'] = target_rqops\n",
    "    data['Target rQOPS'] = data['Target rQOPS'].astype('float64')\n",
    "    data.index = labels\n",
    "\n",
    "    return data\n",
    "\n",
    "rQOPS_table = get_target_rqops(results, labels)\n",
    "rQOPS_table\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "Feel free to use this notebook as a starting point for your own experiments.  For\n",
    "example, you can\n",
    "\n",
    "* explore how the results change considering other problem instances of the Heisenberg model\n",
    "* explore space- and time-trade-offs by changing the value for\n",
    "  `logical_depth_factor` or `max_t_factories`\n",
    "* visualize these trade-offs with the space and time diagrams\n",
    "* use other or customized qubit parameters"
   ]
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "nteract": {
   "version": "nteract-front-end@1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
