samples/azure-quantum/resource-estimation/estimation-time-qubits-constraints.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum resource estimation with time or number of qubits constraints\n",
    "\n",
    "In this Q# notebook we demonstrate resource estimation with constraints \n",
    "on time or the number of physical. We would use a quantum dynamics \n",
    "problem as a sample, more 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. A detailed sample describing this model is available \n",
    "at [Quantum dynamics resource estimation](/~https://github.com/microsoft/Quantum/blob/main/samples/azure-quantum/resource-estimation/estimation-dynamics.ipynb). \n",
    "\n",
    "First, we connect to the Azure quantum service and load the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from azure.quantum import Workspace\n",
    "from azure.quantum.target.microsoft import MicrosoftEstimator, QubitParams, QECScheme\n",
    "import qsharp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace = Workspace(\n",
    "    resource_id=\"\",\n",
    "    location=\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qsharp.packages.add(\"Microsoft.Quantum.Numerics\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Here we just put together the implementation of the Ising model provided at [Quantum dynamics resource estimation](/~https://github.com/microsoft/Quantum/blob/main/samples/azure-quantum/resource-estimation/estimation-dynamics.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Math;\n",
    "\n",
    "function GetQubitIndex(row : Int, col : Int, n : Int) : Int {\n",
    "    return row % 2 == 0             // if row is even,\n",
    "        ? col + n * row             // move from left to right,\n",
    "        | (n - 1 - col) + n * row;  // otherwise from right to left.\n",
    "}\n",
    "\n",
    "function SetSequences(len : Int, p : Double, dt : Double, J : Double, g : Double) : (Double[], Double[]) {\n",
    "    // create two arrays of size `len`\n",
    "    mutable seqA = [0.0, size=len];\n",
    "    mutable seqB = [0.0, size=len];\n",
    "\n",
    "    // pre-compute values according to exponents\n",
    "    let values = [\n",
    "        -J * p * dt,\n",
    "        g * p * dt,\n",
    "        -J * p * dt,\n",
    "        g * p * dt,\n",
    "        -J * (1.0 - 3.0 * p) * dt / 2.0,\n",
    "        g * (1.0 - 4.0 * p) * dt,\n",
    "        -J * (1.0 - 3.0 * p) * dt / 2.0,\n",
    "        g * p * dt,\n",
    "        -J * p * dt,\n",
    "        g * p * dt\n",
    "    ];\n",
    "\n",
    "    // assign first and last value of `seqA`\n",
    "    set seqA w/= 0 <- -J * p * dt / 2.0;\n",
    "    set seqA w/= len - 1 <- -J * p * dt / 2.0;\n",
    "\n",
    "    // assign other values to `seqA` or `seqB`\n",
    "    // in an alternating way\n",
    "    for i in 1..len - 2 {\n",
    "        if i % 2 == 0 {\n",
    "            set seqA w/= i <- values[i % 10];\n",
    "        }\n",
    "        else {\n",
    "            set seqB w/= i <- values[i % 10];\n",
    "        }\n",
    "    }\n",
    "\n",
    "    return (seqA, seqB);\n",
    "}\n",
    "\n",
    "operation ApplyAllX(qs : 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",
    "    ApplyToEach(Rx(2.0 * theta, _), qs);\n",
    "}\n",
    "\n",
    "operation ApplyDoubleZ(n : Int, qs : Qubit[], theta : Double, dir : Bool, grp : Bool) : Unit {\n",
    "    let start = grp ? 0 | 1;    // Choose either odd or even indices based on group number\n",
    "\n",
    "    for i in 0..n - 1 {\n",
    "        for j in start..2..n - 2 {    // Iterate through even or odd `j`s based on `grp`\n",
    "            // rows and cols are interchanged depending on direction\n",
    "            let (row, col) = dir ? (i, j) | (j, i);\n",
    "\n",
    "            // Choose first qubit based on row and col\n",
    "            let ind1 = GetQubitIndex(row, col, n);\n",
    "            // Choose second qubit in column if direction is horizontal and next qubit in row if direction is vertical\n",
    "            let ind2 = dir ? GetQubitIndex(row, col + 1, n) | GetQubitIndex(row + 1, col, n);\n",
    "\n",
    "            within {\n",
    "                CNOT(qs[ind1], qs[ind2]);\n",
    "            } apply {\n",
    "                Rz(2.0 * theta, qs[ind2]);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}\n",
    "\n",
    "operation IsingModel2DSim(N : Int, J : Double, g : Double, totTime : Double, dt : Double, eps : Double) : Unit {\n",
    "    use qs = Qubit[N * N];\n",
    "    let len = Length(qs);\n",
    "\n",
    "    let p = 1.0 / (4.0 - PowD(4.0, 1.0 / 3.0));\n",
    "    let t = Ceiling(totTime / dt);\n",
    "\n",
    "    let seqLen = 10 * t + 1;\n",
    "\n",
    "    let (seqA, seqB) = SetSequences(seqLen, p, dt, J, g);\n",
    "\n",
    "    for i in 0..seqLen - 1 {\n",
    "        // for even indexes\n",
    "        if i % 2 == 0 {\n",
    "            ApplyAllX(qs, seqA[i]);\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(N, qs, seqB[i], dir, grp);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the experiment\n",
    "\n",
    "Next, we are estimating the physical resource estimates to simulate the Ising\n",
    "model Hamiltonian for a $10 \\times 10$ lattice with $J = g = 1.0$, total time\n",
    "$20$, step size $0.25$, and `eps` ${}=0.001$. As configurations for the\n",
    "experiment we use all six pre-defined qubit parameters.  As pre-defined QEC\n",
    "scheme we are using `surface_code` with gate-based qubit parameters (default),\n",
    "and `floquet_code` with Majorana based qubit parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimator = MicrosoftEstimator(workspace)\n",
    "\n",
    "labels = [\"Gate-based µs, 10⁻³\", \"Gate-based µs, 10⁻⁴\", \"Gate-based ns, 10⁻³\", \"Gate-based ns, 10⁻⁴\", \"Majorana ns, 10⁻⁴\", \"Majorana ns, 10⁻⁶\"]\n",
    "\n",
    "params = estimator.make_params(num_items=6)\n",
    "params.arguments[\"N\"] = 10\n",
    "params.arguments[\"J\"] = 1.0\n",
    "params.arguments[\"g\"] = 1.0\n",
    "params.arguments[\"totTime\"] = 20.0\n",
    "params.arguments[\"dt\"] = 0.25\n",
    "params.arguments[\"eps\"] = 0.001\n",
    "params.items[0].qubit_params.name = QubitParams.GATE_US_E3\n",
    "params.items[1].qubit_params.name = QubitParams.GATE_US_E4\n",
    "params.items[2].qubit_params.name = QubitParams.GATE_NS_E3\n",
    "params.items[3].qubit_params.name = QubitParams.GATE_NS_E4\n",
    "params.items[4].qubit_params.name = QubitParams.MAJ_NS_E4\n",
    "params.items[4].qec_scheme.name = QECScheme.FLOQUET_CODE\n",
    "params.items[5].qubit_params.name = QubitParams.MAJ_NS_E6\n",
    "params.items[5].qec_scheme.name = QECScheme.FLOQUET_CODE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are submitting a resource estimation job with all target parameter configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = estimator.submit(IsingModel2DSim, input_params=params)\n",
    "job.wait_until_completed(timeout_secs=300)\n",
    "results = job.get_results()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing and understanding the results\n",
    "\n",
    "### Result summary table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.summary_data_frame(labels=labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we'll define a method for drawing charts to visualize data related to different qubit types. The initial step involves plotting a single point for each qubit type with distinct colors. The data for these points has been obtained earlier.\n",
    "\n",
    "Below, we will take multiple shots for each qubit type. This will introduce variability in the data and lead to more interesting and informative charts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def plot_results(results, labels, num_shots, title):\n",
    "    plt.figure(figsize=(10, 6))\n",
    "\n",
    "    x = []\n",
    "    y = []\n",
    "    label_values = []\n",
    "\n",
    "    for i in range(len(results.data())):\n",
    "        r = results.data()[i]\n",
    "        if 'physicalCounts' in r and 'runtime' in r['physicalCounts'] and 'physicalQubits' in r['physicalCounts']:\n",
    "            x.append(r['physicalCounts']['runtime'])\n",
    "            y.append(r['physicalCounts']['physicalQubits'])\n",
    "            label_values.append(labels[i // num_shots])\n",
    "\n",
    "    for i, label in enumerate(labels):\n",
    "        mask = [l == label for l in label_values]\n",
    "        plt.scatter(np.array(x)[mask], np.array(y)[mask], label=label,s=7)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "\n",
    "    x_min = np.min(x) / 2\n",
    "    x_max = np.max(x) * 2\n",
    "    plt.xlim([x_min, x_max])\n",
    "\n",
    "    y_min = np.min(y) / 2\n",
    "    y_max = np.max(y) * 2\n",
    "    plt.ylim([y_min, y_max])\n",
    "\n",
    "    # for shorter or longer algorithms, it would be worth to modify the xticks\n",
    "    plt.xticks([1e9, 6e10, 3.6e12, 8.64e13, 6.048e14], ['1 sec', '1 min', '1 hour', '1 day', '1 week'])\n",
    "\n",
    "    plt.xlabel('Time')\n",
    "    plt.ylabel('Number of qubits')\n",
    "    plt.title(title)\n",
    "\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing shortest runtimes for each qubit type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_results(results, labels, 1, 'Ising model estimates')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us save params, durations and numbers of qubits for future references."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_shot_params = params\n",
    "single_shot_x = [r['physicalCounts']['runtime'] for r in results.data()]\n",
    "single_shot_y = [r['physicalCounts']['physicalQubits'] for r in results.data()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiments with `max_duration` constraint\n",
    "\n",
    "In this phase, we will conduct experiments by incrementally adjusting the `max_duration` constraint. By default, the Resource Estimator seeks solutions with the shortest runtime possible. However, loosening the `max_duration` constraint provides flexibility in exploring a broader range of possibilities.\n",
    "\n",
    "\n",
    "We will run multiple shots for each type of qubit, systematically increasing the `max_duration` constraint. This approach allows us to observe how the variations in time constraints impact the outcomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiplier = 1.1\n",
    "num_shots = 100\n",
    "\n",
    "params = estimator.make_params(num_items=6 * num_shots)\n",
    "\n",
    "params.arguments = single_shot_params.arguments\n",
    "\n",
    "for param_index in range(6):\n",
    "    max_duration = single_shot_x[param_index]\n",
    "    for shot_index in range(num_shots):   \n",
    "        i = param_index * num_shots + shot_index\n",
    "        params.items[i].constraints.max_duration = f\"{max_duration}ns\"\n",
    "        params.items[i].qubit_params.name = single_shot_params.items[param_index].qubit_params.name\n",
    "        params.items[i].qec_scheme.name = single_shot_params.items[param_index].qec_scheme.name\n",
    "        max_duration *= multiplier\n",
    "\n",
    "job = estimator.submit(IsingModel2DSim, input_params=params)\n",
    "job.wait_until_completed(timeout_secs=300)\n",
    "results_duration = job.get_results()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize resutls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_results(results_duration, labels, num_shots, 'Ising model estimates with duration constraints')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiments with `max_physical_qubits` constraint\n",
    "\n",
    "\n",
    "We will perform experiments by gradually reducing the `max_physical_qubits` constraint. In contrast to the previous `max_duration` constraint experiments, we will take the opposite approach by progressively tightening the constraint from shot to shot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiplier = 1.05\n",
    "num_shots = 100\n",
    "\n",
    "params = estimator.make_params(num_items=6 * num_shots)\n",
    "\n",
    "params.arguments = single_shot_params.arguments\n",
    "\n",
    "for param_index in range(6):\n",
    "    max_physical_qubits = single_shot_y[param_index]\n",
    "    for shot_index in range(num_shots):   \n",
    "        i = param_index * num_shots + shot_index\n",
    "        params.items[i].constraints.max_physical_qubits = int(max_physical_qubits)\n",
    "        params.items[i].qubit_params.name = single_shot_params.items[param_index].qubit_params.name\n",
    "        params.items[i].qec_scheme.name = single_shot_params.items[param_index].qec_scheme.name\n",
    "        max_physical_qubits /= multiplier\n",
    "\n",
    "job = estimator.submit(IsingModel2DSim, input_params=params)\n",
    "job.wait_until_completed(timeout_secs=300)\n",
    "\n",
    "results_qubits = job.get_results()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_results(results_qubits, labels, num_shots, 'Ising model estimates with qubit constraints')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "As we progress, consider exploring alternative algorithms for estimation. The algorithm used for estimation doesn't necessarily have to be an embedded Q# algorithm; you have the flexibility to leverage a variety of examples from the [Quantum GitHub repository](README.md).\n",
    "\n",
    "Feel free to experiment with the multipliers for steps between shots and the number of steps. This customization allows for a more comprehensive exploration of the algorithm's behavior under different configurations.\n",
    "\n",
    "It's important to note that the efficiency of parallelization can vary between algorithms. For parallelizable algorithms, the time-#qubit frontiers may exhibit a wide range, enabling an efficient tradeoff between time and the number of qubits. On the other hand, for algorithms that are not efficiently parallelized, the tradeoff might be less significant.\n",
    "\n",
    "Adjust the parameters and explore different algorithms to gain a deeper understanding of their performance characteristics in the context of resource estimation. Refer to the provided GitHub repository for additional examples and inspiration."
   ]
  }
 ],
 "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}