{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dafa4d64-0da8-4e18-aea8-73f13cac3d29",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# How to choose the Trotter steps for an MPF\n",
    "\n",
    "This guide teaches you how to find a good set of Trotter steps for an MPF.\n",
    "We assume that you have already read and understood the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb022d51-3187-4822-8536-7778df8620f1",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Understanding some heuristics\n",
    "\n",
    "In principle, any $k_j$ values can be chosen but some choices of Trotter steps will lead to a larger noise amplification on real devices than others.\n",
    "Thus, it is important that one tries to find _\"good\"_ values of $k_j$.\n",
    "In the following we will discuss some heuristics to help you guide in doing so.\n",
    "\n",
    "1. In practice, our **largest** $k_j$ ($k_{\\text{max}}$) value is simply bound by what we can afford to execute.\n",
    "\n",
    "2. Since the Trotter error is only well behaved for $\\text{d}t = t/k \\lt 1$ (see also our explanations page on [Understanding the stability of MPFs](https://qiskit.github.io/qiskit-addon-mpf/explanations/mpf_stability.html), our **smallest** $k_j$ ($k_{\\text{min}}$) should satisfy: $t/k_{\\text{min}} \\lt 1$ (Note: empirically $\\text{d}t \\leq 1$ seems to work fine, too).\n",
    "\n",
    "3. None of our $x_j$ values should be close to $0$, as this would imply that we might as well exclude the corresponding $k_j$ from the expansion.\n",
    "\n",
    "4. Similarly, we do not want the $x_j$ for $k_{\\text{max}}$ to be the only dominant one, since that implies that we are essentially just recovering a single product formula with $k=k_{\\text{max}}$.\n",
    "\n",
    "5. Finally, we want the L1-norm of our coefficients, $x_j$, to be small, as this indicates a well-conditioned MPF (see [Carrera Vazquez et al., 2023]). Furthermore, small coefficients minimize the variance amplification of the MPF.\n",
    "\n",
    "[Carrera Vazquez et al., 2023]: https://quantum-journal.org/papers/q-2023-07-25-1067/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0f2096f-5fb6-4861-a05e-3ded9ac65a82",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Applying these heuristics\n",
    "\n",
    "Reusing the example from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial, we want to simulate up to time $t=8.0$.\n",
    "Thus, $k_{\\text{min}}=8$.\n",
    "\n",
    "For the sake of this example, let us assume that our deepest circuit is bound by $k_{\\text{max}}=20$.\n",
    "\n",
    "In addition to the constraints above, we will follow the heuristic presented by [Zhuk et al., 2023], where they minimize a re-weighted L1-norm of the coefficients:\n",
    "\n",
    "$$\n",
    "\\text{min}_{k_j} \\sum_j \\frac{|c_j|}{k_j^{2\\chi}} \\, ,\n",
    "$$\n",
    "\n",
    "where $\\chi$ is the order of our individual product formulas.\n",
    "\n",
    "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ff9ff06-47bf-451d-815f-11d331cafaa3",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Preparing our LSE\n",
    "\n",
    "In a first step, we prepare the LSE, $Ax=b$, leveraging a feature of the [setup_lse](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_lse.html) function that allows us to use a [cvxpy.Parameter](https://www.cvxpy.org/api_reference/cvxpy.expressions.html#parameter) object for our $k_j$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bda8c5d3-b48a-4c9e-bef3-fb9f26846ad6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "order = 2\n",
    "symmetric = True\n",
    "\n",
    "min_k = 8\n",
    "max_k = 20\n",
    "\n",
    "max_l1_norm = 5.0\n",
    "min_coeff = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "3b68ec15-b8ed-45a1-bd46-ced7b29b159e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "from qiskit_addon_mpf.static import setup_lse\n",
    "\n",
    "ks = cp.Parameter(3, integer=True)\n",
    "lse = setup_lse(ks, order=order, symmetric=symmetric)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92356e74-97b6-4231-bd35-43d5cd274dca",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Next, we loop over all possible 3-tuples of $k_j$ values within our specified range.\n",
    "For each one, we compute the analytical coefficients, $x_j$, and then apply our filter criteria.\n",
    "We keep track of those $k_j$ tuples which satisfy our criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "da4c0610-8d1c-4596-b082-8cfd0098bc7c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "possible_ks = {}\n",
    "\n",
    "for k1 in range(min_k, max_k):\n",
    "    for k2 in range(k1 + 1, max_k):\n",
    "        for k3 in range(k2 + 1, max_k):\n",
    "            ks.value = [k1, k2, k3]\n",
    "            coeffs = lse.solve()\n",
    "\n",
    "            if np.any(np.isclose(coeffs, 0.0, atol=min_coeff)):\n",
    "                continue\n",
    "\n",
    "            norm1 = np.linalg.norm(coeffs, ord=1)\n",
    "            if norm1 > max_l1_norm:\n",
    "                continue\n",
    "\n",
    "            weighted = np.sum([np.abs(c) / k ** (2 * order) for c, k in zip(coeffs, ks.value)])\n",
    "            possible_ks[tuple(int(k) for k in ks.value)] = {\n",
    "                \"coeffs\": tuple(coeffs),\n",
    "                \"weighted\": float(weighted),\n",
    "                \"norm1\": float(norm1),\n",
    "            }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f37b658-0e73-409b-add2-4cbf83708730",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Finally, we sort our $k_j$ by the re-weighted L1-norm presented by [Zhuk et al., 2023]:\n",
    "\n",
    "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "16d05ed1-c1ed-445d-afce-543f55275e5e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8, 14, 19) {'coeffs': (0.10447913478216586, -1.7638200183654775, 2.6593408835833117), 'weighted': 9.182736455463762e-05, 'norm1': 4.527640036730955}\n",
      "(8, 13, 19) {'coeffs': (0.13134519801186473, -1.4167162698412699, 2.285371071829405), 'weighted': 9.920634920634922e-05, 'norm1': 3.8334325396825397}\n",
      "(8, 12, 19) {'coeffs': (0.17239057239057254, -1.1944700460829496, 2.022079473692377), 'weighted': 0.00011520737327188946, 'norm1': 3.388940092165899}\n",
      "(9, 13, 19) {'coeffs': (0.2662743506493499, -1.6904000946969673, 2.4241257440476174), 'weighted': 0.00011837121212121191, 'norm1': 4.380800189393934}\n",
      "(8, 13, 18) {'coeffs': (0.15003663003663004, -1.7549001536098316, 2.6048635235732016), 'weighted': 0.00012288786482334872, 'norm1': 4.509800307219663}\n",
      "(8, 12, 18) {'coeffs': (0.19692307692307653, -1.4399999999999984, 2.243076923076922), 'weighted': 0.0001388888888888887, 'norm1': 3.879999999999997}\n",
      "(8, 11, 19) {'coeffs': (0.24195168054817162, -1.070248538011695, 1.8282968574635234), 'weighted': 0.00014619883040935662, 'norm1': 3.14049707602339}\n",
      "(9, 12, 19) {'coeffs': (0.37193877551020393, -1.5167873601053319, 2.144848584595128), 'weighted': 0.00014629507717065315, 'norm1': 4.033574720210664}\n",
      "(8, 12, 17) {'coeffs': (0.22755555555555573, -1.7875862068965531, 2.5600306513409974), 'weighted': 0.00017241379310344843, 'norm1': 4.575172413793107}\n",
      "(8, 11, 18) {'coeffs': (0.27638326585694917, -1.265318468585254, 1.988935202728305), 'weighted': 0.00017284590787313073, 'norm1': 3.530636937170508}\n",
      "(9, 12, 18) {'coeffs': (0.42857142857142816, -1.8285714285714276, 2.3999999999999995), 'weighted': 0.00017636684303350958, 'norm1': 4.657142857142855}\n",
      "(9, 11, 19) {'coeffs': (0.5858035714285708, -1.5251041666666654, 1.9393005952380946), 'weighted': 0.00020833333333333313, 'norm1': 4.05020833333333}\n",
      "(8, 11, 17) {'coeffs': (0.3193762183235871, -1.5289264828738525, 2.2095502645502654), 'weighted': 0.00020885547201336693, 'norm1': 4.0578529657477045}\n",
      "(8, 10, 19) {'coeffs': (0.383090160867938, -1.0642826734780744, 1.6811925126101364), 'weighted': 0.00021285653469561486, 'norm1': 3.1285653469561487}\n",
      "(9, 11, 18) {'coeffs': (0.6750000000000009, -1.8030788177339916, 2.1280788177339907), 'weighted': 0.0002463054187192121, 'norm1': 4.606157635467984}\n",
      "(8, 10, 18) {'coeffs': (0.43760683760683694, -1.2400793650793638, 1.8024725274725268), 'weighted': 0.0002480158730158727, 'norm1': 3.4801587301587276}\n",
      "(8, 11, 16) {'coeffs': (0.3742690058479534, -1.9026640675763484, 2.528395061728395), 'weighted': 0.0002599090318388565, 'norm1': 4.805328135152697}\n",
      "(8, 10, 17) {'coeffs': (0.5056790123456789, -1.4697236919459142, 1.9640446796002353), 'weighted': 0.00029394473838918284, 'norm1': 3.9394473838918285}\n",
      "(8, 10, 16) {'coeffs': (0.592592592592593, -1.7806267806267817, 2.1880341880341887), 'weighted': 0.0003561253561253563, 'norm1': 4.561253561253563}\n",
      "(8, 9, 19) {'coeffs': (0.8112497524262232, -1.3783613445378153, 1.5671115921115921), 'weighted': 0.00042016806722689084, 'norm1': 3.756722689075631}\n",
      "(8, 9, 18) {'coeffs': (0.9266968325791858, -1.5882352941176474, 1.6615384615384616), 'weighted': 0.00048414427499394834, 'norm1': 4.176470588235295}\n",
      "(8, 9, 17) {'coeffs': (1.0708496732026151, -1.855486425339368, 1.7846367521367528), 'weighted': 0.0005656108597285072, 'norm1': 4.710972850678736}\n"
     ]
    }
   ],
   "source": [
    "for ks_val in sorted(possible_ks, key=lambda ks_val: possible_ks[ks_val].get(\"weighted\")):\n",
    "    print(ks_val, possible_ks[ks_val])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3b709d3-ed27-4e6f-a2ba-f3f17326099e",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Recalling the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial, from these possible values we picked `ks=(8, 12, 19)` since it strikes a good balance of minimizing the re-weighted and plain L1-norms."
   ]
  }
 ],
 "metadata": {
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}