{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resource estimation for simulating a 2D Ising Model Hamiltonian\n",
    "\n",
    "In this Python+Q# notebook we demonstrate how to estimate the resources for quantum dynamics,\n",
    "specifically the simulation of an Ising model Hamiltonian on an $N \\times N$ 2D\n",
    "lattice using a *fourth-order Trotter Suzuki product formula* assuming a 2D\n",
    "qubit architecture with nearest-neighbor connectivity.\n",
    "\n",
    "First, we load the necessary Python packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qsharp\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background: 2D Ising model\n",
    "\n",
    "The Ising model is a mathematical model of ferromagnetism in a lattice (in our case a 2D square lattice) with two kinds of terms in the Hamiltonian: (i) an interaction term between adjacent sites and (ii) an external magnetic field acting at each site. For our purposes, we consider a simplified version of the model where the interaction terms have the same strength and the external field strength is the same at each site.\n",
    "Formally, the Ising model Hamiltonian on an $N \\times N$ lattice we consider is formulated as:\n",
    "\n",
    "$$\n",
    "H = \\underbrace{-J \\sum_{i, j} Z_i Z_j}_{B} + \\underbrace{g \\sum_j X_j}_{A}\n",
    "$$\n",
    "where $J$ is the interaction strength, $g$ is external field strength.\n",
    "\n",
    "The time evolution $e^{-iHt}$ for the Hamiltonian is simulated with the fourth-order product formula so that any errors in simulation are sufficiently small. Essentially, this is done by simulating the evolution for small slices of time $\\Delta$ and repeating this for `nSteps` $= t/\\Delta$ 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} 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.\n",
    "\n",
    "## Implementation\n",
    "\n",
    "### Helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that expanding $U_4(\\Delta)$ to express it in terms of $A, B$ gives:\n",
    "$$\n",
    "U_4(\\Delta) = e^{-iAp\\Delta/2} e^{-iBp\\Delta} e^{-iAp\\Delta} e^{-iBp\\Delta} e^{-iA(1 - 3p)\\Delta/2} e^{-iB(1-4p)\\Delta} e^{-iA(1 - 3p)\\Delta/2} e^{-iBp\\Delta} e^{-iAp\\Delta} e^{-iBp\\Delta} e^{-iAp\\Delta/2}\n",
    "$$\n",
    "\n",
    "The above equation with $11$ exponential terms works for one time step. For `nSteps` $> 1$ time steps, some adjacent terms can be merged to give $10t+1$ exponential terms for $e^{-iHt}$.\n",
    "\n",
    "The function below creates a sequence containing the constant factors that will be applied with $A$ and $B$, respectively, in the exponential sequence of the above formula."
   ]
  },
  {
   "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, dt : Double, J : Double, g : Double) : Double[] {\n",
    "\n",
    "    let len1 = 3;\n",
    "    let len2 = 3;\n",
    "    let valLength = 2*len1+len2+1;\n",
    "    mutable values = [0.0, size=valLength];\n",
    "\n",
    "    let val1 = J*p*dt;\n",
    "    let val2 = -g*p*dt;\n",
    "    let val3 = J*(1.0 - 3.0*p)*dt/2.0;\n",
    "    let val4 = g*(1.0 - 4.0*p)*dt/2.0;\n",
    "\n",
    "    for i in 0..len1 {\n",
    "\n",
    "        if (i % 2 == 0) {\n",
    "            set values w/= i <- val1;\n",
    "        }\n",
    "        else {\n",
    "            set values w/= i <- val2;\n",
    "        }\n",
    "\n",
    "    }\n",
    "\n",
    "    for i in len1+1..len1+len2 {\n",
    "        if (i % 2 == 0) {\n",
    "            set values w/= i <- val3;\n",
    "        }\n",
    "        else {\n",
    "            set values w/= i <- val4;\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",
    "    return values;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantum operations\n",
    "\n",
    "There are two kinds of Pauli exponentials needed for simulating the time evolution of an Ising Model:\n",
    "- The transverse field $e^{-iX\\theta}$ applied to each qubit for an angle $\\theta$;\n",
    "- $e^{-i (Z \\otimes Z)\\theta}$ applied to neighboring pairs of qubits in the lattice.\n",
    "\n",
    "The operation below applies $e^{-iX\\theta}$ on all qubits in the 2D lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "microsoft": {
     "language": "qsharp"
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyAllX(n : Int, qArr : Qubit[][], theta : Double) : Unit {\n",
    "    // This applies `Rx` with an angle of `2.0 * theta` to all qubits in `qs`\n",
    "    // using partial application\n",
    "    for row in 0..n-1 {\n",
    "        ApplyToEach(Rx(2.0 * theta, _), qArr[row]);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next operation below applies $e^{-i(Z \\otimes Z)\\theta}$ on overlapping pairs of neighboring qubits. Observe that unlike the previous case, it is not possible to simultaneously apply all the rotations in one go. For example, while applying the rotation on qubits at $(0, 0)$ and $(0, 1)$, it is not possible to also apply the rotation on qubits $(0, 1)$ and $(0, 2)$. Instead, we try to apply as many rotations as possible. This is broken up as follows:\n",
    "- in the vertical (resp. horizontal) direction of the 2D lattice as chosen by `dir`,\n",
    "- consider pairs starting with an even (resp. odd) index as given by `grp`;\n",
    "- apply the exponential to all such pairs in the lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "microsoft": {
     "language": "qsharp"
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "operation ApplyDoubleZ(n : Int, m : Int, qArr : Qubit[][], theta : Double, dir : Bool, grp : Bool) : Unit {\n",
    "    let start = grp ? 1 | 0;    // Choose either odd or even indices based on group number\n",
    "    let P_op = [PauliZ, PauliZ];\n",
    "    let c_end = dir ? m-1 | m-2;\n",
    "    let r_end = dir ? m-2 | m-1;\n",
    "\n",
    "    for row in 0..r_end {\n",
    "        for col in start..2..c_end {    // Iterate through even or odd columns based on `grp`\n",
    "\n",
    "            let row2 = dir ? row+1 | row;\n",
    "            let col2 = dir ? col | col+1;\n",
    "\n",
    "            Exp(P_op, theta, [qArr[row][col], qArr[row2][col2]]);\n",
    "        }\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next operation puts everything together and calls the operations needed to\n",
    "simulate the Ising model Hamiltonian using a fourth order product formula.\n",
    "Observe that the `ApplyDoubleZ` operation is called four times for different\n",
    "choices of direction and starting index to ensure all possible pairs of qubits\n",
    "are appropriately considered.\n",
    "\n",
    "The various parameters taken in by the operation correspond to:\n",
    "\n",
    "- `N1`, `N2`: row and column size for the lattice.\n",
    "- `J`, `g`: parameters by which the Hamiltonian terms are scaled.\n",
    "- `totTime`: the number of Trotter steps.\n",
    "- `dt` : the step size for the simulation, sometimes denoted as $\\Delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "microsoft": {
     "language": "qsharp"
    },
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "import Std.Math.*;\n",
    "import Std.Arrays.*;\n",
    "\n",
    "operation IsingModel2DSim(N1 : Int, N2 : Int, J : Double, g : Double, totTime : Double, dt : Double) : Unit {\n",
    "\n",
    "    use qs = Qubit[N1*N2];\n",
    "    let qubitArray = Chunks(N2, qs); // qubits are re-arranged to be in an N1 x N2 array\n",
    "\n",
    "    let p = 1.0 / (4.0 - 4.0^(1.0 / 3.0));\n",
    "    let t = Ceiling(totTime / dt);\n",
    "\n",
    "    let seqLen = 10 * t + 1;\n",
    "\n",
    "    let angSeq = SetAngleSequence(p, dt, J, g);\n",
    "\n",
    "    for i in 0..seqLen - 1 {\n",
    "        let theta = (i==0 or i==seqLen-1) ? J*p*dt/2.0 | angSeq[i%10];\n",
    "\n",
    "        // for even indexes\n",
    "        if i % 2 == 0 {\n",
    "            ApplyAllX(N1, qubitArray, theta);\n",
    "        } else {\n",
    "            // iterate through all possible combinations for `dir` and `grp`.\n",
    "            for (dir, grp) in [(true, true), (true, false), (false, true), (false, false)] {\n",
    "                ApplyDoubleZ(N1, N2, qubitArray, theta, dir, grp);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting logical resource counts\n",
    "\n",
    "For the purpose of generating the rQOPS for some target runtime, it suffices to obtain the logical resource estimates to simulate the Heisenberg model Hamiltonian. We consider three problem instances with lattice sizes $\\{10 \\times 10, 20 \\times 20, 30 \\times 30\\}$ with $J = g = 1.0$. These instances are simulated for a total time of $L$ seconds for lattice size $L$, with step size `dt`$ = 0.9$, 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": [
    "N1 = [10, 20, 30]\n",
    "N2 = [10, 20, 30]\n",
    "totTime = [10.0, 20.0, 30.0]\n",
    "J = 1.0\n",
    "g = 1.0\n",
    "dt = 0.9"
   ]
  },
  {
   "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": [
    "results = []\n",
    "for i in range(3):\n",
    "    qsharp_string = f\"IsingModel2DSim({N1[i]}, {N2[i]}, {J}, {g}, {totTime[i]}, {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 $10 \\times 10$ lattice by displaying `results[0]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Displaying estimates for 10x10 lattice size.\n",
    "results[0]\n",
    "# Change index to 1 (resp. 2) for 20x20 (resp. 30x30) lattice size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing and understanding the results\n",
    "\n",
    "### Result summary table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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 = [\"Isi10\", \"Isi20\", \"Isi30\"]\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 2 days i.e., 172800 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 = 172800\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.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": 2
}
