{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8d428726-b143-45b7-9c23-d7aef8f3690b",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# How to use the approximate model\n",
    "\n",
    "You have already seen a first example of the [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) function in the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial.\n",
    "In this guide, we are going to revisit that function in more detail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9e6b25b-6f6f-49c6-84ac-9493c151a968",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Setting up a simple model problem\n",
    "\n",
    "For this guide, we will reuse the simple model problem from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial; the Ising model on a line of 10 sites:\n",
    "\n",
    "$$\n",
    "\\hat{\\mathcal{H}}_{\\text{Ising}} = \\sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \\sum_{i=1}^{10} h_i X_i \\, ,\n",
    "$$\n",
    "\n",
    "where $J$ is the coupling strength between two sites and $h$ is the external magnetic field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3c8d6feb-d0fb-45ac-8bf6-213b29bcb8a4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],\n",
      "              coeffs=[1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j,\n",
      " 1. +0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,\n",
      " 0.4+0.j, 0.4+0.j, 0.4+0.j])\n"
     ]
    }
   ],
   "source": [
    "from qiskit.transpiler import CouplingMap\n",
    "from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian\n",
    "\n",
    "coupling_map = CouplingMap.from_line(10, bidirectional=False)\n",
    "\n",
    "hamiltonian = generate_xyz_hamiltonian(\n",
    "    coupling_map,\n",
    "    coupling_constants=(0.0, 0.0, 1.0),\n",
    "    ext_magnetic_field=(0.4, 0.0, 0.0),\n",
    ")\n",
    "print(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46ab2bb9-fd98-4c59-8d9c-66993652ff13",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Choosing $k_j$\n",
    "\n",
    "In this guide, we are going to focus on the tweaks that we can apply to the [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) and its solver that are returned by the [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html). Therefore, we are not concerned with choosing any specific $k_j$ values and can simply reuse the ones from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2ec34c2f-1f6f-4346-9ba6-f7f1c0e3620e",
   "metadata": {},
   "outputs": [],
   "source": [
    "trotter_steps = (8, 12, 19)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c7ad3c-7a5a-4dfe-81a5-9457d68fbd9c",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Computing $x$ analytically\n",
    "\n",
    "First, we have to construct our LSE, $Ax=b$. Once again, we reuse the settings from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "62d48e9c-8b00-43ac-9ee6-90f6418cfb05",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit_addon_mpf.static import setup_lse\n",
    "\n",
    "lse = setup_lse(trotter_steps, order=2, symmetric=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b5c4989-a75d-40e2-bb8a-9070188baef4",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Next, we compute the analytical solution of $x$ and its L1-norm as our reference values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e279b0a3-9639-4e68-8138-90ff5ae6640d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.17239057 -1.19447005  2.02207947]\n"
     ]
    }
   ],
   "source": [
    "coeff_analytical = lse.solve()\n",
    "print(coeff_analytical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a668e3bc-0cd1-4995-bc81-9ac8ca053f02",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.388940092165899\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "print(np.linalg.norm(coeff_analytical, ord=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abda0fe8-c188-49e5-b121-e7d9b64b5313",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Setting up the approximate model\n",
    "\n",
    "Now, we will construct the approximate model.\n",
    "The idea of using approximate solutions $\\tilde{x}$ in order to ensure a smaller L1-norm was proposed by [Zhuk et al., 2023].\n",
    "The model is already explained in detail in the API documentation of [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) so we refrain from repeating it here.\n",
    "Suffice to say, that the model consists of two constraints which ensure that the coefficients, $\\tilde{x}$, sum to $1$ and it enforces their L1-norm to be smaller than the user-provided `max_l1_norm` value.\n",
    "The optimization objective is to minimize the deviation of $A\\tilde{x}=b$.\n",
    "\n",
    "We can see all of this in the output below.\n",
    "Here, we are setting `max_l1_norm=3` because we already know that the L1-norm of the exact solution is approximately $3.4$ and we want to find a $\\tilde{x}$ that has a smaller L1-norm (since otherwise this exercise would be rather pointless).\n",
    "\n",
    "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "de60da28-ce4f-427d-b8c6-11a91cae3472",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "minimize quad_over_lin(Vstack([1. 1. 1.] @ x + -1.0, [0.015625   0.00694444 0.00277008] @ x + -0.0, [2.44140625e-04 4.82253086e-05 7.67336039e-06] @ x + -0.0), 1.0)\n",
      "subject to Sum(x, None, False) == 1.0\n",
      "           norm1(x) <= 3.0\n"
     ]
    }
   ],
   "source": [
    "from qiskit_addon_mpf.static import setup_approximate_model\n",
    "\n",
    "model, coeff_var = setup_approximate_model(lse, max_l1_norm=3)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc5da15f-7771-446b-b0cf-c9320836102a",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Solving the model with default settings\n",
    "\n",
    "First, let us see what happens when we solve the `model` with only default settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "153755d2-4c04-410c-a251-b939f460643a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                     v1.5.3                                    \n",
      "===============================================================================\n",
      "(CVXPY) Sep 05 09:28:06 AM: Your problem has 3 variables, 2 constraints, and 0 parameters.\n",
      "(CVXPY) Sep 05 09:28:06 AM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Sep 05 09:28:06 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Sep 05 09:28:06 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "(CVXPY) Sep 05 09:28:06 AM: Your problem is compiled with the CPP canonicalization backend.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:06 AM: Compiling problem (target solver=OSQP).\n",
      "(CVXPY) Sep 05 09:28:06 AM: Reduction chain: CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP\n",
      "(CVXPY) Sep 05 09:28:06 AM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Sep 05 09:28:06 AM: Applying reduction Qp2SymbolicQp\n",
      "(CVXPY) Sep 05 09:28:06 AM: Applying reduction QpMatrixStuffing\n",
      "(CVXPY) Sep 05 09:28:06 AM: Applying reduction OSQP\n",
      "(CVXPY) Sep 05 09:28:06 AM: Finished problem compilation (took 1.053e-02 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:06 AM: Invoking solver OSQP  to obtain a solution.\n",
      "-----------------------------------------------------------------\n",
      "           OSQP v0.6.3  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2021\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 9, constraints m = 11\n",
      "          nnz(P) + nnz(A) = 33\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-05, eps_rel = 1.0e-05,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on, time_limit: off\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1   0.0000e+00   1.00e+00   2.00e+02   1.00e-01   6.26e-05s\n",
      " 200   4.0661e-09   7.52e-06   1.44e-08   1.22e-05   1.28e-04s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 200\n",
      "optimal objective:    0.0000\n",
      "run time:             1.53e-04s\n",
      "optimal rho estimate: 2.09e-06\n",
      "\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Problem status: optimal\n",
      "(CVXPY) Sep 05 09:28:07 AM: Optimal value: 4.092e-09\n",
      "(CVXPY) Sep 05 09:28:07 AM: Compilation took 1.053e-02 seconds\n",
      "(CVXPY) Sep 05 09:28:07 AM: Solver (including time spent in interface) took 1.028e-03 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4.092215822096244e-09"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.solve(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aafe4c0-c2cc-4de2-824f-ea3594960ca3",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can see from the final value of the cost function that the sum of squares of $A\\tilde{x}-b$ is almost zero.\n",
    "\n",
    "In the output above, we also get a lot more information from [cvxpy.Problem.solve](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem.solve), some of which we will get back to later.\n",
    "\n",
    "For now, let us inspect the final values of $\\tilde{x}$ and its L1-norm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "662b35c2-cd7a-42ff-96a5-792184723f7c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.39646076  0.55556835  0.84089365]\n",
      "1.7929227539822414\n"
     ]
    }
   ],
   "source": [
    "coeff_approx_plain = coeff_var.value\n",
    "print(coeff_approx_plain)\n",
    "print(np.linalg.norm(coeff_approx_plain, ord=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed22d40-140e-4a06-9f6d-90be062fd0cb",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can clearly see that the L1-norm has decreased compared to the analytical solution that we have computed before.\n",
    "\n",
    "Nonetheless, let us see if we can tweak our solution a bit.\n",
    "In the end, we will compare the results obtained for this set of coefficients with all the other ones."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3484a3-ce6e-460a-a144-2bc94ab0ace4",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Tweaking the approximate model\n",
    "\n",
    "Rather than tweaking the input parameters like the chosen `trotter_steps` and `max_l1_norm` values, we will focus on the solver settings in this guide."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a078b0b7-93d0-418f-aaf7-f8584d6889db",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Using a different solver\n",
    "\n",
    "Notice, that we used the default solver method: `OSQP`. We can change that solver to see whether this has any effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "318369c1-7c54-4035-87b9-1e729aa0c048",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                     v1.5.3                                    \n",
      "===============================================================================\n",
      "(CVXPY) Sep 05 09:28:07 AM: Your problem has 3 variables, 2 constraints, and 0 parameters.\n",
      "(CVXPY) Sep 05 09:28:07 AM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Sep 05 09:28:07 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Sep 05 09:28:07 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "(CVXPY) Sep 05 09:28:07 AM: Your problem is compiled with the CPP canonicalization backend.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Compiling problem (target solver=CLARABEL).\n",
      "(CVXPY) Sep 05 09:28:07 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CLARABEL\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction Dcp2Cone\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction ConeMatrixStuffing\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction CLARABEL\n",
      "(CVXPY) Sep 05 09:28:07 AM: Finished problem compilation (took 8.657e-03 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Invoking solver CLARABEL  to obtain a solution.\n",
      "-------------------------------------------------------------\n",
      "           Clarabel.rs v0.8.1  -  Clever Acronym                \n",
      "\n",
      "                   (c) Paul Goulart                          \n",
      "                University of Oxford, 2022                   \n",
      "-------------------------------------------------------------\n",
      "\n",
      "problem:\n",
      "  variables     = 9\n",
      "  constraints   = 11\n",
      "  nnz(P)        = 3\n",
      "  nnz(A)        = 30\n",
      "  cones (total) = 2\n",
      "    :        Zero = 1,  numel = 4\n",
      "    : Nonnegative = 1,  numel = 7\n",
      "\n",
      "settings:\n",
      "  linear algebra: direct / qdldl, precision: 64 bit\n",
      "  max iter = 200, time limit = Inf,  max step = 0.990\n",
      "  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,\n",
      "  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32\n",
      "  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7\n",
      "  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,\n",
      "               max iter = 10, stop ratio = 5.0\n",
      "  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4\n",
      "               max iter = 10\n",
      "\n",
      "iter    pcost        dcost       gap       pres      dres      k/t        μ       step      \n",
      "---------------------------------------------------------------------------------------------\n",
      "  0  +7.1341e-05  -2.3333e+00  2.33e+00  2.40e-01  7.10e-01  1.00e+00  1.88e+00   ------   \n",
      "  1  +7.1330e-05  -1.5432e-01  1.54e-01  1.64e-02  1.12e-01  3.28e-02  1.80e-01  9.61e-01  \n",
      "  2  +6.9440e-05  -1.7266e-03  1.80e-03  2.10e-04  1.89e-03  5.65e-04  2.43e-03  9.87e-01  \n",
      "  3  +2.2053e-05  -5.0002e-05  7.21e-05  5.56e-06  4.83e-05  1.98e-05  7.80e-05  9.87e-01  \n",
      "  4  +1.1119e-06  -2.5543e-05  2.67e-05  1.80e-07  1.41e-06  3.21e-06  1.04e-05  9.90e-01  \n",
      "  5  +7.9250e-09  -2.0484e-06  2.06e-06  2.06e-09  1.59e-08  2.11e-07  6.43e-07  9.90e-01  \n",
      "  6  +2.2543e-09  -2.8390e-08  3.06e-08  2.07e-11  1.59e-10  3.12e-09  9.42e-09  9.90e-01  \n",
      "  7  +1.8480e-09  -8.5733e-10  2.71e-09  1.79e-12  2.36e-10  3.19e-10  9.40e-10  9.15e-01  \n",
      "---------------------------------------------------------------------------------------------\n",
      "Terminated with status = Solved\n",
      "solve time = 113.289µs\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Problem status: optimal\n",
      "(CVXPY) Sep 05 09:28:07 AM: Optimal value: 1.848e-09\n",
      "(CVXPY) Sep 05 09:28:07 AM: Compilation took 8.657e-03 seconds\n",
      "(CVXPY) Sep 05 09:28:07 AM: Solver (including time spent in interface) took 8.504e-04 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.8479854472132682e-09"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.solve(verbose=True, solver=\"CLARABEL\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "207a1f1d-00e1-40ef-bc2a-f904a58fd1ba",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.21284352 -0.00792683  1.22077035]\n",
      "1.4415406931636157\n"
     ]
    }
   ],
   "source": [
    "coeff_approx_clarabel = coeff_var.value\n",
    "print(coeff_approx_clarabel)\n",
    "print(np.linalg.norm(coeff_approx_clarabel, ord=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b672be-9f72-4a11-a9ee-0d4701a8fc86",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Indeed, the outcome has changed. But now we have one major contribution from the largest $k_j$ value which is also not ideal.\n",
    "\n",
    "An exhaustive search of the various available solvers is beyond the scope of this guide. Instead we encourage you to try this yourself."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ba957c-f083-4975-b793-5a96d77a8a37",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Changing the convergence thresholds\n",
    "\n",
    "Rather than changing the solver altogether, we can also tweak the convergence parameters.\n",
    "For the `OSQP` solver the relevant ones are called `eps_abs` and `eps_rel` and they are both set to `1e-5` by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "74528a7b-2bba-445d-b324-a61a261d07c3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                     v1.5.3                                    \n",
      "===============================================================================\n",
      "(CVXPY) Sep 05 09:28:07 AM: Your problem has 3 variables, 2 constraints, and 0 parameters.\n",
      "(CVXPY) Sep 05 09:28:07 AM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Sep 05 09:28:07 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Sep 05 09:28:07 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "(CVXPY) Sep 05 09:28:07 AM: Your problem is compiled with the CPP canonicalization backend.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Compiling problem (target solver=OSQP).\n",
      "(CVXPY) Sep 05 09:28:07 AM: Reduction chain: CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction Qp2SymbolicQp\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction QpMatrixStuffing\n",
      "(CVXPY) Sep 05 09:28:07 AM: Applying reduction OSQP\n",
      "(CVXPY) Sep 05 09:28:07 AM: Finished problem compilation (took 8.835e-03 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Invoking solver OSQP  to obtain a solution.\n",
      "-----------------------------------------------------------------\n",
      "           OSQP v0.6.3  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2021\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 9, constraints m = 11\n",
      "          nnz(P) + nnz(A) = 33\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-08, eps_rel = 1.0e-08,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on, time_limit: off\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1   0.0000e+00   1.00e+00   2.00e+02   1.00e-01   5.02e-05s\n",
      " 200   4.0661e-09   7.52e-06   1.44e-08   1.22e-05   2.10e-04s\n",
      " 400   3.8209e-09   4.95e-06   2.88e-09   2.09e-06   3.72e-04s\n",
      " 600   3.3607e-09   6.09e-06   1.96e-09   2.09e-06   4.33e-04s\n",
      " 800   2.8975e-09   6.24e-06   1.84e-09   2.09e-06   4.87e-04s\n",
      "1000   2.4704e-09   6.08e-06   1.70e-09   2.09e-06   5.40e-04s\n",
      "1200   2.0921e-09   5.78e-06   1.56e-09   2.09e-06   5.92e-04s\n",
      "1400   1.7645e-09   5.40e-06   1.44e-09   2.09e-06   6.44e-04s\n",
      "1600   1.4845e-09   5.01e-06   1.32e-09   2.09e-06   7.01e-04s\n",
      "1800   1.2471e-09   4.63e-06   1.21e-09   2.09e-06   7.53e-04s\n",
      "2000   1.0467e-09   4.26e-06   1.11e-09   2.09e-06   8.05e-04s\n",
      "2200   8.7801e-10   3.91e-06   1.01e-09   2.09e-06   8.57e-04s\n",
      "2400   7.3890e-10   3.25e-06   7.74e-10   2.09e-06   9.10e-04s\n",
      "2600   6.4003e-10   2.47e-06   7.11e-10   2.09e-06   9.62e-04s\n",
      "2800   5.6692e-10   2.03e-06   6.68e-10   2.09e-06   1.01e-03s\n",
      "3000   5.0872e-10   1.76e-06   6.33e-10   2.09e-06   1.07e-03s\n",
      "3200   4.5991e-10   1.58e-06   6.02e-10   2.09e-06   1.12e-03s\n",
      "3400   4.1754e-10   1.46e-06   5.74e-10   2.09e-06   1.17e-03s\n",
      "3600   3.7999e-10   1.37e-06   5.47e-10   2.09e-06   1.22e-03s\n",
      "3800   3.4628e-10   1.29e-06   5.22e-10   2.09e-06   1.32e-03s\n",
      "4000   3.1580e-10   1.22e-06   4.99e-10   2.09e-06   1.37e-03s\n",
      "4200   2.8813e-10   1.16e-06   4.76e-10   2.09e-06   1.43e-03s\n",
      "4400   2.6295e-10   1.11e-06   4.55e-10   2.09e-06   1.48e-03s\n",
      "4600   2.4000e-10   1.06e-06   4.35e-10   2.09e-06   1.53e-03s\n",
      "4800   2.1907e-10   1.01e-06   4.15e-10   2.09e-06   1.58e-03s\n",
      "5000   1.9997e-10   9.67e-07   3.97e-10   2.09e-06   1.63e-03s\n",
      "5200   1.8255e-10   9.23e-07   3.79e-10   2.09e-06   1.69e-03s\n",
      "5400   1.6664e-10   8.82e-07   3.62e-10   2.09e-06   1.74e-03s\n",
      "5600   1.5219e-10   1.01e-06   2.60e-10   2.09e-06   1.79e-03s\n",
      "5800   1.4049e-10   6.59e-07   2.47e-10   2.09e-06   1.84e-03s\n",
      "6000   1.3096e-10   5.75e-07   2.39e-10   2.09e-06   1.90e-03s\n",
      "6200   1.2274e-10   5.24e-07   2.31e-10   2.09e-06   1.95e-03s\n",
      "6400   1.1538e-10   4.90e-07   2.24e-10   2.09e-06   2.00e-03s\n",
      "6600   1.0863e-10   4.67e-07   2.17e-10   2.09e-06   2.05e-03s\n",
      "6800   1.0237e-10   4.48e-07   2.11e-10   2.09e-06   2.10e-03s\n",
      "7000   9.6519e-11   4.32e-07   2.05e-10   2.09e-06   2.16e-03s\n",
      "7200   9.1024e-11   4.19e-07   1.99e-10   2.09e-06   2.21e-03s\n",
      "7400   8.5854e-11   4.06e-07   1.93e-10   2.09e-06   2.29e-03s\n",
      "7600   8.0984e-11   3.94e-07   1.88e-10   2.09e-06   2.36e-03s\n",
      "7800   7.6393e-11   3.82e-07   1.82e-10   2.09e-06   2.41e-03s\n",
      "8000   7.2065e-11   3.71e-07   1.77e-10   2.09e-06   2.46e-03s\n",
      "8200   6.7982e-11   3.60e-07   1.72e-10   2.09e-06   2.51e-03s\n",
      "8400   6.4131e-11   3.50e-07   1.67e-10   2.09e-06   2.56e-03s\n",
      "8600   6.0499e-11   3.40e-07   1.62e-10   2.09e-06   2.62e-03s\n",
      "8800   5.7072e-11   3.30e-07   1.57e-10   2.09e-06   2.67e-03s\n",
      "9000   5.3859e-11   8.77e-06   8.55e-12   2.09e-06   2.72e-03s\n",
      "9200   5.1785e-11   5.78e-07   2.18e-13   2.09e-06   2.79e-03s\n",
      "9400   5.0736e-11   7.72e-08   8.32e-14   2.09e-06   2.87e-03s\n",
      "9550   5.0303e-11   4.69e-08   8.73e-14   2.09e-06   2.96e-03s\n",
      "plsh   4.9640e-11   1.94e-14   1.50e-12   --------   3.00e-03s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      successful\n",
      "number of iterations: 9550\n",
      "optimal objective:    0.0000\n",
      "run time:             3.00e-03s\n",
      "optimal rho estimate: 2.96e-06\n",
      "\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Sep 05 09:28:07 AM: Problem status: optimal\n",
      "(CVXPY) Sep 05 09:28:07 AM: Optimal value: 4.964e-11\n",
      "(CVXPY) Sep 05 09:28:07 AM: Compilation took 8.835e-03 seconds\n",
      "(CVXPY) Sep 05 09:28:07 AM: Solver (including time spent in interface) took 3.988e-03 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4.9640211694779584e-11"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.solve(verbose=True, eps_abs=1e-8, eps_rel=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a5c9fe95-1752-4717-bb74-30766be22f98",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.10925065 -1.          1.89074935]\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "coeff_approx_tight = coeff_var.value\n",
    "print(coeff_approx_tight)\n",
    "print(np.linalg.norm(coeff_approx_tight, ord=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c841fd73-87fa-49d5-b1a8-f0de4b4c53ff",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "As you can see, we have now obtained a comparatively higher L1-norm (which is still constrained by our choice of `max_l1_norm=3`). While the coefficient of the smallest $k_j$ value is still fairly small in magnitude, at least the other two provide a reasonable amount of mixing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f852761-5ef5-4cc6-ab62-d7660f41e82a",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Comparing our results\n",
    "\n",
    "We will now briefly repeat the computations from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial in order to assess the quality of the various coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "71e34633-b044-4fcc-aa63-cebc155a848f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit.synthesis import SuzukiTrotter\n",
    "from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit\n",
    "\n",
    "time = 8.0\n",
    "circuits = []\n",
    "for k in trotter_steps:\n",
    "    circ = generate_time_evolution_circuit(\n",
    "        hamiltonian,\n",
    "        synthesis=SuzukiTrotter(order=2, reps=k),\n",
    "        time=time,\n",
    "    )\n",
    "    circuits.append(circ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "acfc17f7-9817-4576-9709-1f73aa26cf05",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],\n",
      "              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,\n",
      " 0.05+0.j, 0.05+0.j, 0.05+0.j])\n"
     ]
    }
   ],
   "source": [
    "from qiskit.quantum_info import SparsePauliOp\n",
    "\n",
    "L = 10\n",
    "observable = SparsePauliOp.from_sparse_list([(\"Z\", [i], 1 / L / 2) for i in range(L)], num_qubits=L)\n",
    "print(observable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0ea7cac7-0af8-4909-ad23-7b1acb16b607",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit.primitives import StatevectorEstimator\n",
    "\n",
    "estimator = StatevectorEstimator()\n",
    "job = estimator.run([(circ, observable) for circ in circuits])\n",
    "result = job.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "2c2b6c4c-6a27-4601-846a-3e6757ad5156",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array(0.23799162), array(0.35754312), array(0.38649906)]\n"
     ]
    }
   ],
   "source": [
    "evs = [res.data.evs for res in result]\n",
    "print(evs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "655f4663-d19e-4617-89c3-f0355563d915",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import expm\n",
    "\n",
    "exp_H = expm(-1j * time * hamiltonian.to_matrix())\n",
    "\n",
    "initial_state = np.zeros(exp_H.shape[0])\n",
    "initial_state[0] = 1.0\n",
    "\n",
    "time_evolved_state = exp_H @ initial_state\n",
    "\n",
    "exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "35bb0451-091b-4ed2-8776-a1aa7ccb2760",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytical coeff solution: 0.39548478559794153\n",
      "Plain approx.    solution: 0.42928990901998165\n",
      "CLARABEL approx. solution: 0.41833743748333746\n",
      "Tight approx.    solution: 0.3992304700751579\n",
      "Exact            solution: 0.40060242487900255\n"
     ]
    }
   ],
   "source": [
    "print(\"Analytical coeff solution:\", evs @ coeff_analytical)\n",
    "print(\"Plain approx.    solution:\", evs @ coeff_approx_plain)\n",
    "print(\"CLARABEL approx. solution:\", evs @ coeff_approx_clarabel)\n",
    "print(\"Tight approx.    solution:\", evs @ coeff_approx_tight)\n",
    "print(\"Exact            solution:\", exact_obs.real)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adaea843-98fb-48ed-ad1a-0f6ad57f9a7b",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can clearly see now, that the various approximate solutions can yield quite different results.\n",
    "Luckily, testing out different expansion coefficients for a fixed choice of $k_j$ values does not require us to rerun our quantum experiments.\n",
    "Therefore, we encourage you to play around and tweak your expansion coefficients when using the approximate solutions."
   ]
  }
 ],
 "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
}