{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "959d157b-cd6c-46f5-8762-e30dba930d20",
   "metadata": {},
   "source": [
    "# Variational Basis State Encoder (Time Dependent)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf2da1ef",
   "metadata": {},
   "source": [
    "## 1 Background\n",
    "This tutorial is for the time evolution of variational basis state encoder (VBE). Here, spin-boson model is taken as an example. As shall be presented below, the key is to calculate $\\dot{\\theta}_j$ and $\\dot{C}[l]^*$. The algorithm realization are presented in section 2. For more theoretical derivations, see https://doi.org/10.1103/PhysRevResearch.5.023046."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78152079",
   "metadata": {},
   "source": [
    "## 2 Algorithm Realization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b13a28e",
   "metadata": {},
   "source": [
    "### 2.1 Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3fdceaa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "from opt_einsum import contract\n",
    "import tensorcircuit as tc\n",
    "\n",
    "from tencirchem import set_backend, Op, Mpo, Model, OpSum\n",
    "from tencirchem.dynamic import get_ansatz, get_deriv, get_jacobian_func, qubit_encode_basis, sbm\n",
    "from tencirchem.applications.vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f90cb726",
   "metadata": {},
   "source": [
    "### 2.2 Initialize\n",
    "In this tutorial, we study the time evolution of the following spin-boson model:\n",
    "\n",
    "$$\\hat{H}=\\frac{\\epsilon}{2} \\hat{\\sigma}_z + \\Delta\\hat{\\sigma}_x + \\sum_j g_j \\omega_j \\hat{\\sigma}_z (\\hat{b}^\\dagger_j + \\hat{b}_j)\n",
    "+ \\sum_j \\omega_j \\hat{b}^\\dagger_j \\hat{b}_j$$\n",
    "\n",
    "Here, `epsilon`, `delta`, `omega_list` and `g_list` correspond to $\\epsilon$, $\\Delta$, $\\omega_j$ and $g_j$ in the Hamiltonian, respectively.\n",
    "In this section, the parameters are defined and the circuit is initialized. The schematic diagram of the circuit is plotted in Fig. 1, where blue square corresponds to qubit representing spin and green circles correspond to qubits representing vibrations. \n",
    "![fig1](../statics/vbe_td_Fig1.svg)\n",
    "Fig. 1 Schematic diagram of the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "766e26d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-08-08 16:12:43.927000: E external/org_tensorflow/tensorflow/compiler/xla/stream_executor/cuda/cuda_driver.cc:267] failed call to cuInit: CUDA_ERROR_NO_DEVICE: no CUDA-capable device is detected\n",
      "WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "psi_index_top:  ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'v-2-0-bottom', 'v-2-1-bottom'] \n",
      "psi_index_bottom:  ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'v-2-0-bottom', 'v-2-1-bottom'] \n",
      "b_dof_vidx:  [array([1, 2]), array([3, 4])] \n",
      "psi_shape2:  [2, 2, 2, 2, 2]\n"
     ]
    }
   ],
   "source": [
    "set_backend(\"jax\")\n",
    "\n",
    "epsilon = 0\n",
    "delta = 1\n",
    "omega_list = [0.5, 1]\n",
    "g_list = [0.25, 1]\n",
    "\n",
    "nmode = len(omega_list)  # number of phonon modes\n",
    "# make sure correct input\n",
    "assert nmode == len(g_list)\n",
    "\n",
    "# two qubit for each mode\n",
    "n_qubit_per_mode = 2\n",
    "nbas_v = 1 << n_qubit_per_mode\n",
    "\n",
    "# -1 for electron dof, natural numbers for phonon dof\n",
    "dof_nature = np.array([-1] + [0] * n_qubit_per_mode + [1] * n_qubit_per_mode)\n",
    "b_dof_pidx = np.array([1, 2])  # index of basis that need VBE\n",
    "\n",
    "n_dof = len(dof_nature)\n",
    "psi_shape2 = [2] * n_dof\n",
    "\n",
    "psi_idx_top, psi_idx_bottom, b_dof_vidx = get_psi_indices(dof_nature, b_dof_pidx, n_qubit_per_mode)\n",
    "print(\n",
    "    \"psi_index_top: \",\n",
    "    psi_idx_bottom,\n",
    "    \"\\npsi_index_bottom: \",\n",
    "    psi_idx_bottom,\n",
    "    \"\\nb_dof_vidx: \",\n",
    "    b_dof_vidx,\n",
    "    \"\\npsi_shape2: \",\n",
    "    psi_shape2,\n",
    ")\n",
    "\n",
    "\n",
    "# spin boson model has been define in tencirchem.dynamic.model\n",
    "def get_model(epsilon, delta, nmode, omega_list, g_list, nlevels):\n",
    "    ham_terms = sbm.get_ham_terms(epsilon, delta, nmode, omega_list, g_list)\n",
    "    basis = sbm.get_basis(omega_list, nlevels)\n",
    "    return Model(basis, ham_terms)\n",
    "\n",
    "\n",
    "nbas = 16  # number of phonon levels (basis)\n",
    "\n",
    "b_shape = tuple([2] * n_qubit_per_mode + [nbas])\n",
    "\n",
    "assert len(omega_list) == nmode\n",
    "assert len(g_list) == nmode\n",
    "model = get_model(epsilon, delta, nmode, omega_list, g_list, [nbas] * nmode)\n",
    "\n",
    "h_mpo = Mpo(model)\n",
    "\n",
    "# generate the quantum circuit with defined qubits\n",
    "circuit = tc.Circuit(1 + nmode * n_qubit_per_mode)\n",
    "psi0 = circuit.state()\n",
    "n_layers = 3  # layers of ansatz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41235844",
   "metadata": {},
   "source": [
    "### 2.3 Get variational Hamiltonian ansatz terms\n",
    "In this section, we will generate variational hamiltonian ansatz terms. The following ansatz is adopted:\n",
    "\n",
    "$$\\ket{\\phi} = \\prod_n e^{i \\theta_n \\hat {h}_{a_n}} \\ket{\\phi_0}$$\n",
    "\n",
    "if\n",
    "\n",
    "$$\\hat{H}=\\sum_a J_a \\hat{h}_a$$\n",
    "\n",
    "Note that the phonon operators have been transformed to qubit operators based on Gray Encoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d8bb10fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[BasisHalfSpin(dof: spin, nbas: 2), BasisHalfSpin(dof: ('v0', 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ('v0', 'TCCQUBIT-1'), nbas: 2), BasisHalfSpin(dof: ('v1', 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ('v1', 'TCCQUBIT-1'), nbas: 2)]\n",
      "Op('X', ['spin'], 1.0)\n",
      "Op('X', [('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y', [('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z', [('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X', [('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('X X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y', [('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Y X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z', [('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z', ['spin'], 2.0)\n",
      "Op('Z X', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z X X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z Y X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z Z X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X', [('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y', [('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z', [('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X', [('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('X X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('X Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y', [('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Y X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Y Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z', [('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z X X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z X Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z Y X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Y Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
      "Op('Z Z X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
      "Op('Z Z Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n"
     ]
    }
   ],
   "source": [
    "def get_vha_terms():\n",
    "    basis = sbm.get_basis(omega_list, [nbas_v] * nmode)\n",
    "    spin_basis = qubit_encode_basis(basis, \"gray\")\n",
    "\n",
    "    spin_ham_terms = OpSum([Op(\"X\", [\"spin\"], 1.0)])\n",
    "    for i in range(nmode):\n",
    "        complete_list = []\n",
    "        for j in range(n_qubit_per_mode):\n",
    "            complete = OpSum()\n",
    "            dof = (f\"v{i}\", f\"TCCQUBIT-{j}\")\n",
    "            for symbol in \"IXYZ\":\n",
    "                complete += Op(symbol, dof)\n",
    "            complete_list.append(complete)\n",
    "        complete_real = complete_list[0]\n",
    "        for c in complete_list[1:]:\n",
    "            complete_real = complete_real * c\n",
    "        spin_ham_terms.extend(complete_real)\n",
    "        spin_ham_terms.extend(Op(\"Z\", \"spin\") * complete_real)\n",
    "    spin_ham_terms = OpSum([op.squeeze_identity() for op in spin_ham_terms.simplify() if not op.is_identity]).simplify()\n",
    "    return spin_ham_terms, spin_basis\n",
    "\n",
    "\n",
    "spin_ham_terms, spin_basis = get_vha_terms()\n",
    "\n",
    "print(spin_basis)\n",
    "for i in range(len(spin_ham_terms)):\n",
    "    print(spin_ham_terms[i])\n",
    "\n",
    "\n",
    "theta0 = np.zeros(n_layers * len(spin_ham_terms), dtype=np.float64)\n",
    "ansatz = get_ansatz(spin_ham_terms, spin_basis, n_layers, psi0)\n",
    "jacobian_func = get_jacobian_func(ansatz)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6a646d6",
   "metadata": {},
   "source": [
    "### 2.4 Time Evolution\n",
    "This sections defines the function involved in time evolution. The derivations $\\dot{\\theta}_j$ are calculated by following equation:\n",
    "\n",
    "$$\\sum_j \\textrm{Re} \\frac{\\partial \\bra{\\phi}}{\\partial \\theta_k} \\frac{\\partial \\ket{\\phi}}{\\partial \\theta_j}\\dot{\\theta}_j = \n",
    "\\textrm{Im} \\frac{\\partial \\bra{\\phi}}{\\partial \\theta_k} \\hat {H}_1$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\hat{H}_1=\\prod_{l}\\hat{B}[l] \\hat{H} \\prod_{l}\\hat{B}[l]^\\dagger$$\n",
    "\n",
    ", and $\\frac{\\partial \\ket{\\phi}}{\\partial \\theta_j}$ is calculated via `theta_deriv`.\n",
    "The time derivations of $B[l]$ are calculated by following equations:\n",
    "\n",
    "$$i\\rho[l]\\dot{B}[l]^* = (1-\\hat{P}[l])\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$$\n",
    "$$\\dot{B}[l]^* = -i \\rho[l]^{-1}(1-\\hat{P}[l])\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\tilde{H}'=\\left( \\prod_{k\\neq l}\\hat{B}[l] \\hat{H} \\prod_{k\\neq l}\\hat{B}[k]^\\dagger \\right) \\hat{B}[l]^\\dagger$$\n",
    "\n",
    "the projection operator\n",
    "\n",
    "$$\\hat{P}[l]=\\hat{B}[l]^\\dagger\\hat{B}[l]$$\n",
    "\n",
    "and the reduced density matrix\n",
    "\n",
    "$$\\rho[l]_{nn'}=\\textrm{Tr} \\left[ \\bra{\\phi}\\ket{n}_l \\bra{n'}_l \\ket{\\phi} \\right ]$$\n",
    "\n",
    "An example that illustrates the function in detail is also presented below.\n",
    "![fig2](../statics/vbe_td_Fig2.svg)\n",
    "Fig. 2 Graphic illustration of (a) $\\rho[l]$, (b) $P[l]$, (c) $\\bra\\phi\\tilde{H}'[l]\\ket\\phi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7597768e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv_fun(t, theta_and_b):\n",
    "    # split \\thera and b to independent arrays\n",
    "    theta = theta_and_b[: len(theta0)]\n",
    "    psi = ansatz(theta)\n",
    "    b_array = theta_and_b[len(theta0) :].reshape(nmode, nbas_v, nbas)\n",
    "\n",
    "    # get contracted H\n",
    "    h_contracted = get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom)\n",
    "    # get the derivation of \\theta\n",
    "    theta_deriv = get_deriv(ansatz, jacobian_func, theta, h_contracted)\n",
    "\n",
    "    psi = psi.reshape(psi_shape2)\n",
    "    b_deriv_list = []\n",
    "    for i in range(nmode):\n",
    "        b = b_array[i]\n",
    "        # calculate rho\n",
    "        indices_base = [(\"contract\", ii) for ii in range(n_dof)]\n",
    "        psi_top_indices = indices_base.copy()\n",
    "        psi_bottom_indices = indices_base.copy()\n",
    "        for j in b_dof_vidx[i]:\n",
    "            psi_top_indices[j] = (\"top\", j)\n",
    "            psi_bottom_indices[j] = (\"bottom\", j)\n",
    "        out_indices = [(\"top\", j) for j in b_dof_vidx[i]] + [(\"bottom\", j) for j in b_dof_vidx[i]]\n",
    "        args = [psi.conj(), psi_top_indices, psi, psi_bottom_indices, out_indices]\n",
    "        rho = contract(*args).reshape(1 << n_qubit_per_mode, 1 << n_qubit_per_mode)\n",
    "\n",
    "        from scipy.linalg import pinv\n",
    "\n",
    "        rho += np.eye(len(rho)) * 1e-5\n",
    "        rho_inv = pinv(rho)\n",
    "\n",
    "        b = b.reshape(nbas_v, nbas)\n",
    "        # calculate projector\n",
    "        proj = b.conj().T @ b\n",
    "\n",
    "        # derivative\n",
    "        args = get_contract_args(psi, h_mpo, b_array, i, n_qubit_per_mode, psi_idx_top, psi_idx_bottom, b_dof_pidx)\n",
    "        k = b_dof_pidx[i]\n",
    "        args.append(b_array[i].reshape(b_shape))\n",
    "        args.append([f\"v-{k}-{l}-bottom\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-bottom\"])\n",
    "        # output indices\n",
    "        args.append([f\"v-{k}-{l}-top\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-top\", \"mpo-0\", f\"mpo-{len(h_mpo)}\"])\n",
    "\n",
    "        # take transpose to be compatible with previous code\n",
    "        b_deriv = contract(*args).squeeze().reshape(nbas_v, nbas).T\n",
    "        b_deriv = np.einsum(\"bf, bg -> fg\", b_deriv, np.eye(nbas) - proj)\n",
    "        b_deriv = -1j * np.einsum(\"fg, fh -> hg\", b_deriv, rho_inv.T)\n",
    "        b_deriv_list.append(b_deriv)\n",
    "    return np.concatenate([theta_deriv, np.array(b_deriv_list).ravel()])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c552000f",
   "metadata": {},
   "source": [
    "Let's illustrate the code in detail.\n",
    "Code that calculates $\\rho[l]$ (see Fig. 2(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1477e876",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode  0\n",
      "out_indices:  [('top', 1), ('top', 2), ('bottom', 1), ('bottom', 2)]\n",
      "args: \n",
      " [DeviceArray([[[[[1.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]]],\n",
      "\n",
      "\n",
      "\n",
      "             [[[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]]]], dtype=complex128), [('contract', 0), ('top', 1), ('top', 2), ('contract', 3), ('contract', 4)], DeviceArray([[[[[1.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]]],\n",
      "\n",
      "\n",
      "\n",
      "             [[[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]]]], dtype=complex128), [('contract', 0), ('bottom', 1), ('bottom', 2), ('contract', 3), ('contract', 4)], [('top', 1), ('top', 2), ('bottom', 1), ('bottom', 2)]]\n",
      "rho: \n",
      " [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# split \\thera and b to independent arrays\n",
    "\n",
    "b_list = []\n",
    "for _ in range(nmode):\n",
    "    b = np.eye(nbas)[:nbas_v]  # nbas_v * nbas\n",
    "    b_list.append(b)\n",
    "theta_and_b = np.concatenate([theta0, np.array(b_list).ravel()]).astype(complex)\n",
    "\n",
    "theta = theta_and_b[: len(theta0)]\n",
    "psi = ansatz(theta)\n",
    "b_array = theta_and_b[len(theta0) :].reshape(nmode, nbas_v, nbas)\n",
    "\n",
    "# get contracted H\n",
    "h_contracted = get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom)\n",
    "# get the derivation of \\theta\n",
    "theta_deriv = get_deriv(ansatz, jacobian_func, theta, h_contracted)\n",
    "\n",
    "psi = psi.reshape(psi_shape2)\n",
    "b_deriv_list = []\n",
    "i = 0\n",
    "print(\"mode \", i)\n",
    "b = b_array[i]\n",
    "\n",
    "# calculate rho[l]\n",
    "indices_base = [(\"contract\", ii) for ii in range(n_dof)]\n",
    "psi_top_indices = indices_base.copy()\n",
    "psi_bottom_indices = indices_base.copy()\n",
    "for j in b_dof_vidx[i]:\n",
    "    psi_top_indices[j] = (\"top\", j)\n",
    "    psi_bottom_indices[j] = (\"bottom\", j)\n",
    "out_indices = [(\"top\", j) for j in b_dof_vidx[i]] + [(\"bottom\", j) for j in b_dof_vidx[i]]\n",
    "print(\"out_indices: \", out_indices)\n",
    "args = [psi.conj(), psi_top_indices, psi, psi_bottom_indices, out_indices]\n",
    "print(\"args: \\n\", args)\n",
    "rho = contract(*args).reshape(1 << n_qubit_per_mode, 1 << n_qubit_per_mode)\n",
    "print(\"rho: \\n\", rho)\n",
    "\n",
    "from scipy.linalg import pinv\n",
    "\n",
    "rho += np.eye(len(rho)) * 1e-5\n",
    "rho_inv = pinv(rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e4ebde6",
   "metadata": {},
   "source": [
    "Code that calculates $\\hat{P}[l]$, (see Fig. 2(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6780faa5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proj:  [[1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
      "  0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "b = b.reshape(nbas_v, nbas)\n",
    "# calculate projector P[l]\n",
    "proj = b.conj().T @ b\n",
    "print(\"proj: \", proj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08928a9f",
   "metadata": {},
   "source": [
    "Code that calculates $\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$ (see Fig. 2(c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5115b03a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "args in derivative: \n",
      " [DeviceArray([[[[[1.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]]],\n",
      "\n",
      "\n",
      "\n",
      "             [[[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]],\n",
      "\n",
      "               [[0.-0.j, 0.-0.j],\n",
      "                [0.-0.j, 0.-0.j]]]]], dtype=complex128), ['p-0-top', 'v-1-0-top', 'v-1-1-top', 'v-2-0-top', 'v-2-1-top'], DeviceArray([[[[[1.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]]],\n",
      "\n",
      "\n",
      "\n",
      "             [[[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]],\n",
      "\n",
      "\n",
      "              [[[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]],\n",
      "\n",
      "               [[0.+0.j, 0.+0.j],\n",
      "                [0.+0.j, 0.+0.j]]]]], dtype=complex128), ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'v-2-0-bottom', 'v-2-1-bottom'], <Matrix at 0x7f68b5462690 (1, 2, 2, 3) float64>, ['mpo-0', 'p-0-top', 'p-0-bottom', 'mpo-1'], <Matrix at 0x7f68503148d0 (3, 16, 16, 3) float64>, ['mpo-1', 'p-1-top', 'p-1-bottom', 'mpo-2'], <Matrix at 0x7f68503147d0 (3, 16, 16, 1) float64>, ['mpo-2', 'p-2-top', 'p-2-bottom', 'mpo-3'], array([[[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "        [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]],\n",
      "\n",
      "       [[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "        [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]]]), ['v-2-0-bottom', 'v-2-1-bottom', 'p-2-bottom'], array([[[1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
      "         0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j],\n",
      "        [0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
      "         0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j]],\n",
      "\n",
      "       [[0.-0.j, 0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
      "         0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j],\n",
      "        [0.-0.j, 0.-0.j, 0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
      "         0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j]]]), ['v-2-0-top', 'v-2-1-top', 'p-2-top'], array([[[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "        [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]],\n",
      "\n",
      "       [[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "        [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
      "         0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]]]), ['v-1-0-bottom', 'v-1-1-bottom', 'p-1-bottom'], ['v-1-0-top', 'v-1-1-top', 'p-1-top', 'mpo-0', 'mpo-3']]\n"
     ]
    }
   ],
   "source": [
    "# derivative\n",
    "args = get_contract_args(psi, h_mpo, b_array, i, n_qubit_per_mode, psi_idx_top, psi_idx_bottom, b_dof_pidx)\n",
    "k = b_dof_pidx[i]\n",
    "args.append(b_array[i].reshape(b_shape))\n",
    "args.append([f\"v-{k}-{l}-bottom\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-bottom\"])\n",
    "# output indices\n",
    "args.append([f\"v-{k}-{l}-top\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-top\", \"mpo-0\", f\"mpo-{len(h_mpo)}\"])\n",
    "print(\"args in derivative: \\n\", args)\n",
    "\n",
    "# take transpose to be compatible with previous code\n",
    "b_deriv = contract(*args).squeeze().reshape(nbas_v, nbas).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5197573",
   "metadata": {},
   "source": [
    "Contract tensors and calculate $\\dot{B}[l]^*$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1921f511",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_deriv = np.einsum(\"bf, bg -> fg\", b_deriv, np.eye(nbas) - proj)\n",
    "b_deriv = -1j * np.einsum(\"fg, fh -> hg\", b_deriv, rho_inv.T)\n",
    "b_deriv_list.append(b_deriv)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "075f0e23",
   "metadata": {},
   "source": [
    "## 2.5 Main Structure of the Function\n",
    "This is the main stucture of the code. We set parameters for time evolution and initialize the system. Then, we update `theta_and_b` in time evolution and calcutate $\\langle \\sigma_z \\rangle$ `z` and $\\langle \\sigma_x \\rangle$ `x` step by step as the time evolves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "04e5b336",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1.0\n",
      "1 0.9801368820199942\n",
      "2 0.9221622297057478\n",
      "3 0.8307164281767181\n",
      "4 0.712902625694087\n",
      "5 0.5774908043700016\n",
      "6 0.4340489975856556\n",
      "7 0.29215043357543824\n",
      "8 0.1607516949764221\n",
      "9 0.04771089791814444\n",
      "10 -0.040502642430565255\n",
      "11 -0.09904787191380995\n",
      "12 -0.12493534752914795\n",
      "13 -0.11727456640490394\n",
      "14 -0.07753204265425653\n",
      "15 -0.009732009482324864\n",
      "16 0.0795128404043315\n",
      "17 0.18126828150880098\n",
      "18 0.2849420221890886\n",
      "19 0.37936481084384754\n",
      "20 0.454192768911805\n",
      "21 0.5014051186981502\n",
      "22 0.5165776398039041\n",
      "23 0.49960262655096854\n",
      "24 0.454624729908236\n",
      "25 0.38915769753158064\n",
      "26 0.31259780970942813\n",
      "27 0.2345813095925135\n",
      "28 0.1636633691524335\n",
      "29 0.10653226239863313\n",
      "30 0.06766729976756988\n",
      "31 0.04926876472593326\n",
      "32 0.05139921686178229\n",
      "33 0.07232882932406118\n",
      "34 0.10899831508240451\n",
      "35 0.15745403958628454\n",
      "36 0.21315107699434557\n",
      "37 0.27112011043252915\n",
      "38 0.3261026901074539\n",
      "39 0.3728086768061548\n",
      "40 0.40637666288671287\n",
      "41 0.4229683862535386\n",
      "42 0.42032486089352117\n",
      "43 0.3981218276511652\n",
      "44 0.358045462603464\n",
      "45 0.303614837159146\n",
      "46 0.23991582373561202\n",
      "47 0.17311347674692437\n",
      "48 0.10948894653071074\n",
      "49 0.05440718215103625\n",
      "50 0.011717807210271824\n",
      "51 -0.016427087542680092\n",
      "52 -0.02981545362112685\n",
      "53 -0.030162355384666367\n",
      "54 -0.020379375716308275\n",
      "55 -0.0033909784547891857\n",
      "56 0.018576266600602798\n",
      "57 0.04368717824936481\n",
      "58 0.06971819393413235\n",
      "59 0.09364096597083588\n",
      "60 0.11186039266314592\n",
      "61 0.12076156704155695\n",
      "62 0.11755115511155295\n",
      "63 0.10090752714979544\n",
      "64 0.07138982669111976\n",
      "65 0.03253624871953725\n",
      "66 -0.011552300380283664\n",
      "67 -0.05607764510179323\n",
      "68 -0.09628855507563477\n",
      "69 -0.12836308326138923\n",
      "70 -0.14994039236833612\n",
      "71 -0.1596458713814789\n",
      "72 -0.15576794369843444\n",
      "73 -0.1373071561893098\n",
      "74 -0.10671249014129201\n",
      "75 -0.06743879382842288\n",
      "76 -0.022285878126101024\n",
      "77 0.024280809221682687\n",
      "78 0.06810377510953083\n",
      "79 0.10678673871439198\n",
      "80 0.13914885388461803\n",
      "81 0.1611886576331503\n",
      "82 0.16581373326788112\n",
      "83 0.1542497273174869\n",
      "84 0.13048909287272364\n",
      "85 0.09708896314020657\n",
      "86 0.057310903820287765\n",
      "87 0.015533678436258006\n",
      "88 -0.024192504501293052\n",
      "89 -0.05790366853755518\n",
      "90 -0.0797592370408079\n",
      "91 -0.08582739907891124\n",
      "92 -0.07396770719727741\n",
      "93 -0.048080239560664026\n",
      "94 -0.012411290850183795\n",
      "95 0.031070370091046573\n",
      "96 0.08042891268644264\n",
      "97 0.13037553828085885\n",
      "98 0.17451072357728148\n",
      "99 0.20513840470656897\n"
     ]
    }
   ],
   "source": [
    "theta0 = np.zeros(n_layers * len(spin_ham_terms), dtype=np.float64)\n",
    "ansatz = get_ansatz(spin_ham_terms, spin_basis, n_layers, psi0)\n",
    "jacobian_func = get_jacobian_func(ansatz)\n",
    "\n",
    "b_list = []\n",
    "for _ in range(nmode):\n",
    "    b = np.eye(nbas)[:nbas_v]  # nbas_v * nbas\n",
    "    b_list.append(b)\n",
    "theta_and_b = np.concatenate([theta0, np.array(b_list).ravel()]).astype(complex)\n",
    "z_list = [1]\n",
    "x_list = [0]\n",
    "\n",
    "tau = 0.1\n",
    "steps = 100\n",
    "\n",
    "dummy_model = get_model(epsilon, delta, nmode, omega_list, g_list, [nbas_v] * nmode)\n",
    "z_op = Mpo(dummy_model, Op(\"Z\", \"spin\", factor=1)).todense()\n",
    "x_op = Mpo(dummy_model, Op(\"X\", \"spin\", factor=1)).todense()\n",
    "\n",
    "for n in range(steps):\n",
    "    print(n, float(z_list[-1]))\n",
    "    # update `theta_and_b`\n",
    "    sol = solve_ivp(deriv_fun, [n * tau, (n + 1) * tau], theta_and_b)\n",
    "    theta_and_b = sol.y[:, -1]\n",
    "    theta = theta_and_b[: len(theta0)]\n",
    "    psi = ansatz(theta)\n",
    "    z = psi.conj().T @ (z_op @ psi)\n",
    "    x = psi.conj().T @ (x_op @ psi)\n",
    "    z_list.append(z.real)\n",
    "    x_list.append(x.real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "82642da3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEeCAYAAABrB7XiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBbElEQVR4nO3de1yb9d0//lcSjqWFEHqmJ0LPra0E2qq1altQp7ZOBbRzurmtibrtvncS7m73d7e7t3sV1p3c4deg29Rtdm3SOut0WlLbarW2QHqwZ0joAWgpEC6gnCHX74/kSoEGyOE6JXk/H48+7hvCdeUTB+988v6835+PgmVZFoQQQsKaUuoBEEIIER4Fe0IIiQAU7AkhJAJQsCeEkAgQJfUApLZo0SKkp6f7fV1tbS1SU1NFu06K56TrIvM6KZ6TruP3OqvVipqamsHfZCOcTqdj161bx7755pt+Xbdu3bqAni/Q66R4TrouMq+T4jnpOn6ue/PNN9l169axSUlJNz0W8TP71NRU7N69W7Tn27Bhg+jXBvOcYj4fvT5+rwtGuL/GcH19GzZswIYNG5CZmXnzgwG93YQRKWZNoSLcXyO9vtAX7q+Rz/hEC7QBkmK2JbZwf430+kJfuL9GPl+fgmWl66C12+2wWCwwmUwoLS316Zri4mKo1WoAAMMwKCgo8OvxodavXy9qGocQQoTmLa5JNrO3Wq2wWCxgGAYOh8Ona7hArtfrodfrodPpYDAYfH6cEEIilaQzewAwm83YvHkzKioqRv3Z5ORkVFdXe2buAKBQKMC9hNEe94Zm9oSQcCOrmb2/7HY7GIYZFMg5Fotl1Mf55nTS/nGEkNARUsHeG7VaDYZhRn2cT6Un6rDqx++jsa2L1/sSQohQQr7OXqPRwOFweJ3RD3x8OLW1tVi/fr3na65OdTi9fU4U/t2Kyiut+NofP8VbL9wDlTJk3jMJIWFo27Zt2LZtm+fr2tram34m5IP9aIu7oz3ub1NVdJQSf/v2nVj94gfYd+oqfrrzBF7Mu9Xn6wkhhG9DJ6kDJ7CckJmSarVar99nGAZarXbUx/m0cJoav//6CgDAL985jXcqLvN6f0II4VtIBXu1Wu01N5+dnT3q43zLu30Wnr93HgDAYDyEhlbK3xNC5EvyYD9cmsVut6O4uHjQ9zZt2jSossZsNkOv1/v8ON9+9kQG5qcmoa2rD6Un6gR7HkIICZZkwd5ut6OwsBBGoxFWqxWFhYUoKSnxPG6xWGA0GgddU1BQAIZhUFJSgpKSEpSVlQ36mdEe51t0lBIP6lzbj+47eVWw5yGEkGBJ3lQltWCbqg6cvoqHXvoQk5LiUPnyI1AoFDyOjhBC/BfSTVVydducCYiPUaG+pQtnalukHg4hhHhFwT5IsdEqrJw3EQDwIaVyCCEyFfHBnmuqGtiQ4K/ViycDAPadvMLXsAghxG/btm3D+vXrw7OpKlh8nFS1epEr2B88ew3dvf2IjVbxMTRCCPEL11wV0k1VcrZomhoTEuPQ0dOPI1WNUg+HEEJuQsGeB0qlAqsXTQIA7DtFeXtCiPxQsOfJPYu4vD0Fe0KI/FCw5wmXt7dWO9DR3SfxaAghZDAK9jyZlpKAlHGxcLIsztW1Sj0cQggZhII9jxZOSwIAnK5hpB0IIYQMQcGeRwtT1QBAnbSEENmhYM+jBe6Z/Rma2RNCZCbigz0fHbScBanuYE8ze0KIBKiDdgR8dNByFkxTAwAuN3WgtbMXifHRvNyXEEJ8QR20IklOiMGU5HgAwFma3RNCZISCPc8WplJFDiFEfijY82w+l7evoZk9IUQ+KNjzbKE7b0+LtIQQOaFgz7MF1FhFCJEhCvY8mz/VFezrW7rQ1NYt8WgIIcSFgj3PxsVHY+b4BACUyiGEyEfEB3s+m6o43CItlV8SQsQ0UlNVxAd7rqlqw4YNvN2TW6SlvD0hREwbNmzA7t27kZqaetNjER/shXBjkZZm9oQQeaBgL4B5UxIBALb6NolHQgghLhTsBTBr4lgAwFWmk06tIoTIAgV7ASQnxCBpjGsTtEuN7RKPhhBCKNgLQqFQYNYE1+z+QsN1iUdDCCEU7AUzkwv21yjYE0KkR8FeIDSzJ4TIScQHeyGaqgAgzb1IW00ze0KISOikqhHweVLVQLMmuLZMoJk9IUQsdFKVBLjyywvXroNlWYlHQwiJdBTsBTI9JQEKBdDR04+G1i6ph0MIiXCSp3GKi4uhVqsBAAzDoKCgYMSfz8vLQ05ODrKysjzXcbRaLSwWC4xGI3JycqDValFaWoply5YhNzdXoFfgXWy0CqnJY1Dj6MCFhnZMTIoX9fkJIWQgSYM9F+j1ej0AwGKxwGAwwGg0DnuN1WqF2Wy+6fu5ubkwmUxgGMbzM1qtFoWFhaIHes6siWNdwf7adSyfPV6SMRBCCCBxGmfz5s3Iz8/3fJ2dnY2SkpIRr8nNzQXLsoP+GY1GmEwmz89UVFSAZVnYbDbPG4kUqPySECIXkgV7u90OhmFuSsUArhn+cAwGw00/m5WVxffweEHll4QQuZAsjWO3271+X61Wg2GYYa/TarWe/59hGNjtdmRnZw/6mR07dkCj0QAAysrKUFRUFPyAAzCTyi8JITIh+QLtUBqNBg6Hw6efLSwsvCm/r9VqodPpPG8KDocDeXl5g9I8A3FNVRyuTpUPs2jLBEKICLZt2zaoMTQkmqp8DfRWq9Xr93U63aCv8/PzYTAYhk0ZCdVUBdxI49Q2d6Cnrx8xUSpBnocQEtmGTlJl1VQ1MB0zEMMwwz42kNFoRHp6+k3fH1qpwwX44dJGQpqQGIcxMSqwLG11TAiRlqTBXq1Wew3CQ3Pw3lgslptm6gzDIC8vb9A9ufy/L28gfFMoFLT7JSFEFiQtvdy0adOgyhuz2TyoVNJut6O4uNjrtXa7/aYArlarYTQaB32/pKQEubm5XlM4YvBsm9BAM3tCiHQkzdkXFBSguLjYU1tvs9kGLbhy3bDeumq1Wq2n4mag/Pz8QW8QTU1Nwy7OiiFtApVfEkKkJ/kC7UjbI+j1+mGbomw2m9fvq9XqUbdcEBNXfnmxkYI9IUQ6tBGawKaluIJ9bVOHxCMhhEQyCvYCm+4O9pebKGdPCJFOxAd7oU6q4kxLGQMAuNbahZ6+fkGegxBCADqpakRCNlUBwPhxsYiNVqK714m65k5PVy0hhPCNTqqSkEKhwDSNa3ZfQ6kcQohEKNiLIFXjytvX0CItIUQiFOxFwOXtL1OwJ4RIhIK9CLg0Tq2D0jiEEGlQsBdBagqlcQgh0qJgL4Lp7jROjYOCPSFEGhTsReBJ41A1DiFEIhEf7IVuqgJupHGYjl60dfYK9jyEkMg2UlNVxAd7rqmKr6MIvUmMj0bSmGgAQC2lcgghAtmwYQN2796N1NTUmx6L+GAvlmmeRVpK5RBCxEfBXiSpGqq1J4RIh4K9SG7U2lOwJ4SIj4K9SDxpHAr2hBAJULAXCbdlAuXsCSFSoGAvEk9jFeXsCSESoGAvEm7ny1pHB1iWlXg0hJBIE/HBXoymKgCYmhwPAOjq7UdjW7egz0UIiUzUVDUCMZqqACA2WoVJSXEAqCKHECIMaqqSiRv72tMiLSFEXBTsRcSVX9bSIi0hRGQU7EXkaaxqpmBPCBEXBXsRTUl2Bfs6ytkTQkRGwV5EqRpXRU6to1PikRBCIg0FexFN5Wb2lMYhhIiMgr2IuJ0v65qpsYoQIq6ID/ZiNVUBwBR3Y1V3rxNN16mxihDCL2qqGoFYTVUAEBOlwoREV2NVHeXtCSE8o6YqGbmxSEt5e0KIeCjYi4wWaQkhUoiSegDFxcVQq9UAAIZhUFBQMOLPWywWGI1G5OTkQKvVorS0FMuWLUNubm7A9xSTZ5E2RGf23b39iI1WST0MQoifJJ3Zc0FZr9dDr9dDp9PBYDCMeA3DMLBarTAYDDAYDEhPT/ca6P25p5imerpoQytn393bj6/+4SCmP2fGjk8vSD0cQoifFKyENYDJycmorq72zMIBQKFQjFiWaDabkZ2dPeiaYO65fv167N69O5DhB+TNg3YYSj7D6kWTsbtwjWjPG4yunn58+Xcf44PjdQCAKJUC279zN+5dOlXikRFCvPEW1ySb2dvtdjAM4zVoWywW2dyTb6khdvB4Z08fnvjNAXxwvA7xMSrcs3AS+vpZfPl3H+NwZYPUwyOE+EjSYO+NWq0GwzAjXrtjxw6YzWaYzWYUFhbyck+xcAu0V0JkgfbX/zqNvSevIiE2Cju/fw92/uAe5CyZgs6efuT96gAaWrukHiIhxAeSL9AOpdFo4HA4hn1cq9VCp9NBq9UCABwOB/Ly8mAymQK6J9dUxdmwYYOgNfdczr6tqw+tnb1IjI8W7LmCxbIsdhy6CADY8nQWVi2YBAD467dXYfWLH+BMbQt2fnYRz947T8phEhLxtm3bNqgx1FtTleyC/UiBHgB0Ot2gr/Pz82EwGEacuY90T66pSiwJsVFITohBc3sP6hwdSExNEu25/XXqMgNbfRtio5V4OGu65/sJsVH46j3pKPy7FebDFOwJkdrQSerACSwn4DTOrl27sGvXrkAv98zMh2IYZtjHANcC7UBcft5utwd8T7Fx2ybIPW//z7LLAIDsW6Zi3JBPII8snwGFAjhc2YhLjXTyFiFy51ewP3bsGJ599lksW7YMdrsdNpsNWVlZeO6553Ds2DG/nlir1UKtVnvNs2dnZ3u9hmEY5OXlDbqGm9FrtdqA7imFUFmk/WfZJQDAF5dNv+mxKcljsGq+K62z8/BFUcdFCPHfqMG+tbUVv/jFL5CVlQWj0QiDwYCysjL84Ac/wAsvvIDy8nLo9Xps3boVy5Ytw5YtW9Da2urTk2/atGlQlYzZbIZer/d8bbfbUVxc7PlarVbDaDQOmqWXlJQgNzfXM8Mf7Z5yEAqLtGdqGJyra0VMlBJfyLh5nw0AeOy2mQAA8yEK9oTI3bA5+507d8JoNEKhUCAvLw/l5eXD3iQjIwNbt24FALzyyitYs2YNUlJSYDAY8Oijjw57XUFBAYqLi1FSUgIAsNlsMBqNnse5btmBHbD5+fmD3gCampoGLc6Odk85SA2BxiouhbNm8WQkjYnx+jMPZ03H998ow4lLzThX14J5U+W7/kBIpPPaVJWfn4/09HTo9XqkpaUFdOPq6moYjUZUV1dj+/btQQ9UKGI3VQHA6wds+NafDuPepVOx8/v3iPrcvlrxw3dxuqYFRv1t+NKdw693PPbL/dhzvA6bvrgYP3x0iYgjJIQMx1tc8zqz37FjR9BPlpaWhpdeeino+4Sjqe4FWrnuj3OurgWna1oQrVLigYxpI/5s3m0zsed4HUyfXaRgT4iM8dJUdeHCBT5uEzHkvkC7/9RVAMCqBROhTvCewuE8qJuGKJUCVVfbqCqHEBnjJdgnJyd7qnGOHj2KTZs28XFbUYh5UhWHW6Btbu9BZ0+faM/rqwq7qy9hxezxo/7suPho3DpLAwD49Nw1QcdFCBmZ4CdVVVdXo6CgAB9++CEyMjKQn5+P++67j49bC07Mk6o4SWOikRDryqDVyXCR1lrdBADQaVN8+vk75k4EAHx6nvbKIURKgp9UtXbtWhQXFyMzMxMXLlxAeXk5bDYbH7cOSwqF4sZWx03ySuW0dvbi/BVX6Wymj8H+9rkTANDMnhA54yXYb968GeXl5UhKSoLNZoPVakVFRQUftw5bqdwircxq7Y9VO8CywIzxCZ7zckfDBftzda1obAu9jdFaOnrAtPdIPQxCBMXL3jh6vR4tLS0AXLN8q9WK5uZmJCVR3fVwpsp0kbaCS+GkaXy+JmVcLOanJuFsbQsOnW/AusybO27l6rPKBjz6i31o6+rD1OR4LJymxvfXLcKd8ydKPTRCeMXbFscDA/sLL7yA5uZmvm4dljzHE8psZm+1+5ev59zhSeWETt7+80vNyP3lfrR1uRbJ65o7Yfn8Ch7/9QFcbLgu8egI4ZdfwT4lxfcAkJGR4fdgIsmN8kt5LdBaq12VOJlpfgb7ea5gf+h8aOTtK6+04uHifWjp6MVtcyag8uVHYPl/OVg+ezxaO3uhLzmEfqdT6mESwhu/gn1zc7PP+96QkXHll3JqrGpo7cKlxnYoFMCtfqRxgBsVOccuNKO9W37lpAP19jmR96v9aGjtwtKZyTB//25MVsdjxZwJePXZOzA2LgqfnmvAb987K/VQCeGN32mc6upqbNmyBZs2bcKuXbso+Afoxv448gn2Fe4UztwpiX4fqjJ9fAKmacag38mirKpRiOHxxvTZBdjqr2NiUhzeemH1oL1/0iaORdGTmQCAn+08gRMXKR1JwoPfwV6n02HPnj2oqKhAQUEBkpOT8cMf/lCIsYlCiqYqAJiqcVXjXGvpQk9fv6jPPRxPvt7PFA6HS+XIuQTT6WTxm3fPAACev3ee14qjp+7S4qHMaejtd6Lo7ZNiD5GQgPHaVLV161bs2bMHe/bsQVVVFSorK+FwOEKmiWooKZqqACBlbCxio13/+a/IpLGqgsvXa/1L4XDumCf/5qo9J+pwprYF4+Ki8PU1c7z+jEKhwP97zLXPz3tHa+icXRIyeG2qevzxxwd9rdVqsXXrVuj1erz66quBjzLCKBQKpCbLp/ySZdmAK3E4t81xzeyt9iY4nTdtpioLv373NADgmdVzRtz3Z+E0NbK0KejrZ/HmwWqxhkeIYPwK9tnZ2di7d6/Xxx577DEqt/TTVBmVX15u6kBjWzeiVArcMj05oHvMm5qI+BgV2rr6UFXfxvMIg3e4sgGfnmtAtEqJb943+rm5T9+dDgB444ANXnYCJySk+BXst27dihdeeAHHjx/3+nhycmBBIlJxi7Q1Mtgy4XQNAwCYNzUJcTGqgO4RpVJiyUzX78BRd3OWnPzufVd1zRMrZ3neaEfy2G0zMSZGhfNXWvFZpbwXnQkZjV/BnkvZZGRk4Pnnnx8U9C9cuICqqireBxjOPOWXMpjZn611dUAvSA2u65nrvD3qzv/LRUtHD/591LVo9fy9o8/qASAxPhqPrHAdvfjGAdrriYQ2v3P22dnZcDgcaGhoQEZGBlQqFVQqFfLy8kK6KkcKqe6KHDk0Vp2tc5XQzp+aGNR9MtyVPFaZBfv3jtaip8+JeVMTsXiG759Av+JO5ew6fBGtnb1CDY8QwQW0XYJarYbJZILT6URVVRWqqqpQVlaGxMTgAkWk8eTsZbBAe849s58f5Mw+w723/fELDll1oL515BIA4JHlM/y67rY54zF3SiI6evqxu/yyEEMjRBRB742TlpYW8Dm1kc5TjSNxGodlWZytcwX7YA8NnzNlHMbGRaGjpx/n6+TRcNfS0YO9n18BADzqZ7BXKBSeN4g9x+t4HxshYvEa7Lds2eI5eSoYx44dw5YtW4K+j5CkaqoCbizQXmU60dsn3Sy4pqkD17v6EKVSIH3SuKDupVIqPSdXySWVMzCFs2Ca2u/rc5ZMAQDsO3kFff3y+bRCyFB+N1WtXbsWP//5zzFnzhxs2bLFry0RWltb8Ytf/AJZWVnYvHkzsrOzAx+5CKRqqgKACYlxiFIpwLJAfYt0eXtuVj97ciKio4LfCDVDZou0gaZwOFnpKUhOiAHT0YsjMt8KgkQ2v5uqMjIysGPHDlRWViIpKQlr1qzBfffdh127dg37JDt37sS9996LtWvXIjk5GeXl5di+fTtuvfVW3l5IuFEqFZ6KHCkbq/iqxOFwFTlWGZRfBpPC4aiUSmTf4prdl564wtvYCBHTqNO4jRs3ory8HFu3bsWRI0cwe/ZsPPfcczh27BiOHTuGZ599FnPmzEF5eTmMRiPKysrwjW98Q4yxh4Ub5ZdSzuz5qcThcBU5n19iJE1PAcGncDg5S6cCoLw9CV0+n1SVlpaGl156CS+99BL27t2Ln//851AoFNDr9di6dauQYwxr01LGAJXSzuzPuBuqgq3E4WgnjkXSmGi0dPTibF0LbvGj1JFv77graAJN4XC4mf2JS824ynRisjo+6LERIqaAjiVcu3Yt1q5dy/dYIpLUaRyWZXGOm9nzFOwVCgUyZmmw/3Q9rNUOyYJ9b58T+09dBQDcf+vNOUx/TEiMgy5NA2u1A6Un6vDUXel8DJEQ0fB2LCEJDNdYJVWt/ZXmTrR29kKlDL4SZyBPc5Vdurz94apGtHX1IWVcrKf+Pxj3ulM5lLcnoYiCvcSkPsSEq8TRThqH2OjA9sTxhqvIOXZBuoqc0hOu/Hr2LVOgVCqCvl/OElewpxJMEooo2EtM6i5avitxOFyt/cnLjGSHs1jcM3Au3x6sTK3GU4JZIeEnFkICEfHBXsqmKuBGF+0VplOS7QXOcNsk8FSJw5k1IQHJCTHo6XPiTE0Lr/f2xVWmEycuNUOhANbyFOxVSiVWLZgEADh4Vr6ncZHIxetJVeFGyqYqAJikjoNKqUBfP4t6RvwTkc7W8rs4y1EoFFjKbXcsQSrH4q6tz5il8Xr0YKBWhsDRiyRy8XpSFeGXSqnE1GTXIm2NyKkcVyUOPxugeXOrhHl7iztfz+XZ+XLnfNfM/tD5Bsrbk5BCwV4GpqUkAABqmtpFfd7Gtm40t/dAoQBmT+avEofDVcCIHez7+p348KSr5DJ7CT8pHM6i6UlIGhONtq4+fH6J4fXehAgpoDp7PhUXF0OtVgMAGIZBQUGBT9cAgM3mOlDCaDR6HrNYLDAajcjJyYFWq0VpaSmWLVuG3Nxc/gfPk+kpY3AIwKVGcYN95ZVW9/MnID6G/1+FoYu0MVH8VfuMpMLehOb2HiQnxCArwPN0h6NSKnH73Al4/1gdDp6t91QdESJ3ks7suUCv1+uh1+uh0+lgMBhGvKawsBAFBQUoKCjwBPmcnBzP4wzDwGq1wmAwwGAwID09XdaBHgCmj+dm9uKmcaquus6JFWJWDwBpE8dCPSYa3b3iLtJye+Hcs2gyolT8/4qvnDcRAPDJuQbe702IUCQN9ps3b0Z+fr7n6+zsbJSUlAz781wgZxjG8z2DwQCLxQK73e75XkVFBViWhc1mg16vF2TsfJruTuNcFjmNU3nVNbOfM1mYQ2cUCgWWumf3Yi7Sfujuml2zeLIg979zvivYf3ruGpxOOoichAbJgr3dbgfDMJ4UzkAWi2XY68rLywcFdq1WCwCD3gBCzbQUV/ml2MFe6Jk9cCOVI1bevqWjB+U2Vw38msX85us5S2dqkBAbheb2Hk/pKiFyJ1nOfmDAHkitVg8buNVqNZqbmwd9j3tj4II+AOzYsQMajSvIlJWVoaioiIcRC2dGijRpHC5nP2eKcMdJir1I+/GZa+h3skifNA4z3OkxvkVHKXHbnPHYe/IqDp6tx6LpakGehxA+Sb5AO5RGo4HD4Xtg2Lx5M4xGo+cTglarhU6n8wR/h8OBvLw8mEwmr9dzTVWcDRs2iF5zn+oO9s3tPWjr7MW4+GjBn7Pf6YS9/joAgWf2aTcWaXv7nLwcjjKSfadc+frVi4RJ4XBWzp+EvSev4pNzDTDkzBP0uQgZzbZt2wY1hnprqpJdsPcn0BcWFsJgMAzKy+t0ukE/k5+fD4PBMGzKiGuqklJifDTUY6LBdPSipqk9qH3XfXWxoR29/U7ERis9aSQhDNzu+ExtC5bMFHYHzH2n6gEAqwXK13O45qqDZ6+BZVkoFMHvvUNIoIZOUgdOYDmS5ewHpl0GYhhm2McGMpvNSE9Pv2kB1mw2D/qaC/DDpY3kgqvIuSxSKofL16dPGgeVUrhfA1cnrTiLtDVN7ai80gqlQoG73NsaCCVTm4LYaCUaWrtgq28T9LkI4YOkwV6tVnsNwqOdW8vl6blAzzCMZ8E3Ly9v0D25/L8vbyBSEruxqspdiTNboEqcgXRaV7AXevMwrpFKp9VAnRAj6HPFRquQ6a7hP3SeSjCJ/Elaerlp06ZBlTdms3nQTN1ut3saqDhWqxVWqxU6nQ52ux1WqxWbN2+GRqOBWq2G0WgcFNhLSkqQm5vrNYUjJ9PdqZRLIgX7yiuu2eicKcLl6znL0scDgOCHde/jSi4FztdzbpvjSuVQsCehQNKcfUFBAYqLiz219TabzWs3LNdVyzAM1q5dC4ZhUFhYOOheXMVNfn7+oDeIpqamYRdn5WSayBU5Ys7sl892BfvTNYxgC9BOJ+s5lWq1QCWXQ90+1xXsP6sU9k2MRCarvQnvVNTgf/KW8nI/yRdoR9oegeus5XgrvRxKrVb7tOWC3MwQubGq0p2znyNgJQ5nsjoeM8Yn4FJjO6zVTbh7If8z75OXGTS2dSMhNgrLZ/O7RcJwuDexyiutaGjt4nV3TRK5mtt78L+m4/jTvkqwLLBizvigj9UEaCM02eAqYsSY2bd393nOvBVjZg8Ay9JdAfhIlTB5e25L4zvnTxRtDx7N2FjPoS+fVVIqhwRv38mr0BW8g1c/dAX6x++Y5WlMDBYFe5ngtkyodXQIvnWuzT2rT06IQcq4WEGfi8Pl7ctswqQ8LJ9zWxqLk8LheFI55ymVQ4LT1dOPZ185hMa2bsxPTcJ7m9bi1WfvwGR1PC/3j/hgL/VJVZzJ6nhEq5Tod7K4ynQK+lxcvl7Iztmhls2+sUjLsvzuJ3O9q9cTbPk6lcpXt83lFmnpMBMSnNf2V6GuuROpmjH4+Cf3e05F8wedVDUCqU+q4iiVCqRqXO/gQtfac9skCNk5O9TSmcmIiVKiqa0b1deu83rvj87Uo7ffibSJY5E+SbzXBNyY2R+70IyO7j5Rn5uEj86ePmx55xQA4IX1ixAXE1gqkk6qChFcRc5lgfe15xqqxJzZx0arPMcU8p3K4Q4WX7t4iuidrDPHJ2CyOh69/U5Yq+kQchKYP31YhfqWLswYn4Cn7hKmJ4iCvYyItdXxja2NxZ0Fc9UrZTwv0nL71/N9KpUvFAqFZ3Z/iPL2JADt3X341b9OAwAKHl4sWIEBBXsZmS5CRQ7LsgMaqsSb2QPCLNLa6ttgv3YdUSrht0gYzu2UtydBeG1fFRpau5A2cSy+tDJNsOeRvM6e3MClcYTsor3W0oXWzl4oFIB2ojQz+xOXmtHZ08fLUYjcrP72uRNE2S3UGy7YH65sRL/TKeheQ2R4juvdKD1Rh8tNHXjijlmevye52/ZJNQDg2/fPF3RXWAr2MjJjvPD743ApnJnjEwJeBArUtJQxmKyOx1WmExV2h+fEp2CUnnCVXK4VqWvWm1tmqJEYH43Wzl58fonhrS6a+OZCw3V860+H8fGZa3C6K71efu8MtupvwwMZ0yQe3cjO1rbg+MVmRKkUeHTFTEGfi6YgMsLtfHmpsZ338kTOeYlSOIArv32Hexb80emrQd+vu7cfH59xpU5ylkwN+n6BUimVnhLMg2cplSOm3j4nvvL7gzhwuh5OlsWi6Wosmq5Gc3sPHv/1R/ivv1eg3yls30owTIcuAACyb5kieM8LBXsZmTk+AQoFcL2rD41t3YI8h+d0KpE6Z4da466D33sy+GC///RVtHf3YUpyPBZLfFrUjUPIKdiLafM/P4e12oHkhBiUbX4Qn/3fAzjw4n345n2uA2X+8ME5/P79cxKP0juWZWH67CIAIP/2WYI/X8QHe7k0VQGu8sTUZNcirV2gPdLPi3AU4Ui4HSnLbU1g2nuCutfbZZcBAOsyp0GplPbwEC4l9clZOoRcLJ+eu4ZfvuOqYnn5meWY7966IjZahZeezMSvv7IMAPCznSc8kxw5Kbc3ofradSTERuEBHT/pJmqqGoFcmqo4WndTkJ3nxiNOlfuXfq5EwX76+ATMnZIIJ8viwOn6gO/T1+/Eu1bXL/TDWTP4Gl7AMmZpMCZGRYeQi6StsxcbjYfgZFk8uUqLLy6/+Xfg62tmY83iyejq7ce3/nxYdm/COz69AAB4KHMaEmL5WT6lpqoQkjZxLACgWoCZfXdvPy40uBZ/xdjHfjhr3EcGcufFBuLTcw1wXO+GZmws7nAfESil6CglVrj3t6dUjvBKLOdxqbEdM8cnoPjLmV5/RqFQ4OVnliMhNgqfnmvAqx9WijzK4fX1O7Hz8CUAQN7twi7McijYy4wn2Asws6++dh1OlsXYuCjeNlcKxBp35cyHQeTt3y53/aE8qEtFlEoev8ZcKocWaYXV0d2H379/FgDwo8eWIHGEktuZE8bifx+/FQDw4+3H0NDaJcYQR/XRmXo0tHYhZVws1iwSp5JMHn8lxEPINM75AYuzUh6QvWrBRESrlKi+dj2gtQmnk8U7FTUAgIeXTed7eAFbOSDYC1VNFYhzdS34xyfVOH7Bgd4++Vam+Or1AzY0tnVj1oQE5N02+qz4G2vmIGOWBu3dfdi6Rx6Lte+Uu35/12dNF7S2fiAK9jKjdc/shVig9VTiSJjCAYCxcdGeBqtAZvfl9iZcae7EuLgo3CPAQSiByky7cQg5dziMlBpau/Dd18qwfNN72Gg8hDt//D6mGkx48uWP0dbZK/XwAtLd24/fvncGAPDdBxf69KlOqVTg++sWAgBe2VuJ613Svnank8W7R13Bfl2meH0AFOxlJs09s29s6+b9D7JSgg3QhrP2FleQ/vCk/3l7rgrn/ltTERstbmPYSOJiVJ4tIaRO5bxrrUGG+xAMJ8ti6cxkJI2JRldvP3aXX8YTv/kIXT39ko4xEG8erEatowNTkuPx5CrfNwx7KHMaZk8eh+b2Hry23ybgCEdnrb4xWRFziw8K9jKTGB/taa7gO29/vs5VJSJVJc5AXN7+wOl6v1IL/U4n/lnmytevz5JPCofD5e0/PhN4pVGwymyN+OofPkFLRy+WzEjGe5vW4uBPv4BLf8zFe5vWYmxcFD46U4+v/vGTkErr9Dud+M27rlLL//jCAr/e6FVKJf7zgQUAgN/9+wx6+qR7o/uX1TWrz1kyVdTJCgV7GUoTIJXDsqwkWxsP59ZZyZiQGIfWzl58cLzO5+s+OF6HS43tSE6Iwb1LpeuaHc49i7hPLFcl6dy82HAdj//6I3T19uO+pVNx4Cf3eQ7BUCoVWLVgEnZ8927ERivxrrUG33m9TPQxBqr0xBXYr11HckIMnlk92+/rN6xMw2R1POqaO7Hj0EUBRugbrmT4IRFTOAAFe1k1VXHSBVikbWzrRrO7iUnsAz68USmV+NKdrh3+/rK/yufrSkrPAwCeuisdY3iqTebT8vTxSIyPhuN6N45WO0R97tbOXuT96gAaWrtwyww1/vL8Sq857VULJuGv31oFpUKBNw7YBDsqkm9/3uf6PXlylTaguvTYaJWns/a3752RZBG98korzta2IFqlFGSyQk1VI5BbUxUwsPySv5k9tzg7PWWMbILkV+9JB+DazOySDwe2nL/Sir0nr0KhADZmzxF6eAGJjlJ6ZvfcIehiKfxbBc7UtmCyOh6m790z4i6gX8hIxQb3m+1PTMfFGmLAapra8cEx1yfAQGb1nK+tmYOE2CicrW3BZ5Xiv8m9607hrFowEUljYni/PzVVhRghau3ltDjLmT05EXcvnASWBd44MPqiGTerv//WVMyaMFbo4QWMO/S89IR4wf7Dk1fwt4/tUCiA17+5EqmaMaNe88NHbkFMlBIHTtdjHw97FQnpjQM2OFkWqxZMDGrNKTE+Go+ucHXbvu7D7xzfuHz9Qzxtj+APCvYylObeZ57PnL3UG6AN52vuWdobH9nQ1z98jru1sxd/P2gHADybM1eUsQUq273ZW7mtCY7rwmxoN9D1rl58+89HAACG7Lm4Y55vW0fPGJ+Ar69x/ff/ifmYrHoDBurrd3oC89eCmNVzvnK36xPlW4cvoqUjuP2Z/HGtpRNHqlyfJvjaC8cfFOxlKH2Sa9Za4+hAdy8/VQPcfi3zpsor2D+UOQ3jx8XiSnMn3j82/ELttoN2XO/qw5wpibKqrfdmWkoCFqQmwcmyosyY/9d8HJca2zFjfAL+J2+pX9e+sH4xEmKjUGF3eBrV5Ob9Y3Woa+5EyrhYrMsMvgJr+ezxmJ+ahI6efuz8TLyF2nettWBZQJem8emTF98o2MvQhMQ4JMRGgWWBizwdPn6mhgEALJR4K+ChYqJUnnrpEst5r7PLS43t+NnOEwAAQ/YcyXe49AW3v36pwHn7MlsjtrrTWy8/sxxj4/w7rWtCYpxn0XLL7lO8j48P3AL+l1dpeSlVVCgUntm9mKmcf1W4+kMe4uENKxAU7GVIoVDwukjb0tGDy+5zbRdOUwd9P749s3o2VEoF9p266jl4mdPb58Qzf/wETEcvsrQpQS3OiYnL21tO1AmWHul3OvG918vAsq6ywrW3BLbHynP3zUNMlBJHLzhEryAaTT3TCYt77YML0Hx4YuUsRKuUsFY7cOJiM2/3HU5rZy/2u3d5FbNrdiAK9jLl2SOnPvhF2tM1rhTO1OR4JCfwXwEQrPRJ4/AL986FL5qO42130xQA/HTnCRypakTSmGj85ZsrERMln47Zkdw+dwLGxKhQ39KFzy8xgjzHa/ttOHahGUljovGzJ24N+D7jx8XhYXeDmj9lsGLYefginCyLrPQUXosLxo+L8wTd1w8I/5pLj9ehp8+JOVMSPfvui42CvUzxObPnUjiLZJbCGWhj9lzPwutG4yF860+Hse6lvfi1u2PyD19fIesKnKFio1W4a6GrmenfR/nPhTe2dXlKJn/06BJMTApuF9OvuRdqTYcuyGrfnB3uY/seF+Akp6fdnxRMhy7ytjY2nHe4FI4EC7OciA/2cmyqAm5siGbjYWZ/yh3sF6Sqg76XkDZ/SYd7l0xBZ08/Xj9g83zsff7eeXh4mfQHlPjri+4x/+PTC7yncn5iOo7m9h4snq7GxrXB9xysnDcRc6Yk4npXn+eoPKlVXW1Fhd0BlVKBx3zY3dJf9yyahKnJ8Whu78G/j97chMSX7t5+7HF3ia/LEjbYU1PVCOTYVAUAc91VM2d5OPXo1GXXPRZNl+bjo6+iVEq89s078f11C1H48GJs3Xgb9r94H156Uif10AKyPms64mNUqLrahnJ7E2/3PVLV6FlY3PJ0Fi/7+SsUCjzjbnJ7bZ88Ujkm95YGaxZPxoTEON7vr1IqPY1lXFmvEA6crkdbl+us5My0FMGeB6CmqpDELaReamxHaxAfq1mWxWkujSPDxdmhxsVH48W8W/Hfjy3Bk6u0yNSmSLr3fjDGxUd7Nmv7xyfVvNyzr9+J//zLEbAs8KU70zwHnfNhw51pslmoZVkW293H9gl5GPeX7nRVgpWeuIJrLZ2CPAeXwnkwQ9qzkinYy5RmbCymJLvysFywDsRVphPN7T1QKhSYN1XeM/tw9MQdswAA5s8u8bLT4h/3nMPJywySE2Lwsycygr7fQOPHxeGL7sNgfOloFpK12gFbfRviY1SCbhg2d0oilqWnoN95482FT/3OG2clC53CGY3km6QUFxdDrVYDABiGQUFBQdDXBHJPOVo8XY0rzZ04fZnBbXMCO2f11GUGAJA+eRziYkKjkiWc3LNoMiYlxaG+pQt7jl8JKnBdbmzH/7n7DX76RIYgqY0nV2mx49BF7DpyCcVfzhTtFKWhuMO4H9RN87t3wF9PrtKizNaEvx+sxrfun8/rJ8kDp13HD2rGxmLVfPH2rvdG0pk9F5T1ej30ej10Oh0MBkNQ1wRyT7niUjmngpjZn3bn/BdNo1m9FKJUSuS7Z/f/+DTwVA7LsvjeG2Xo6OnH7XMn4Ck/Du7wx10LJmFiUhwc17uxN4CDZfjgdLJ4q0y8w7gfXTETsdFKnLrM8F5z/49PLgAAHlsxQ7I3To6kz75582bk5+d7vs7OzkZJSUlQ1wRyT7niSiW5BdZAcDN7OZddhrsNK12LgP8+WuvZZtpfr+234f1jdYhWKfHbry4TLPcbpVIid4UrwJrcZY9i+6yyAVeaO5E0JhprFwt/GHdyQgwezHB94vrbx/wt1LZ393ny9dwbvpQkC/Z2ux0Mw3jSLQNZLJaArgnknnLGLaiermECLt07HSJll+HslhnJuGWGGj19Tvzh/bN+X195pRX/9fcKAMCPc5dggcAL7dxs+l8VNZKc1/rWEdes/oGMaaKd5PTUXa5PSts/vcDbcY3vWWtwvasPsyYkYIX7zGUpSRrsvVGr1WAYJqBrArmnnM2bmgiVUoHm9h7UNftfKdDvdOJsrWu3S7mXXYa7wocXAwBe/vcZ1Dk6fL6ut8+Jb2z9FB09/bh74ST8xxcWCDVEj0xtCrQTx6Kjpx/vWYWrP/fG6WTxT/cZw9xWxGJYvXgypqeMQXN7D3aXX+blntyC7+N3pMmiokx21TgajQYOh39lX6NdM9LjXFMV909OzVWx0SrMnuzaNoFLx/jDXn8dXb39iI9ReTpyiTTWZ03Hijnj0dnTj5/tOuHzdT/ecQzWageSE2Jg1N8uSumeQqHwpB22i5zKOVTZgKuMK4WzZrF4u5uqlEo8dRd/m6M1tHZ5Dq8RI4XDNVNx/0KiqcrfQO/LNSM9zjVVcf/k1ly1KIhFWi6FM39qElRK2f1PHVEUCgV+vsHVHPa3j+04eWn0hcCX/30Gv3enfV5+Zrmo2+LmuWvb935+BQ2tXaI971uHXSmcB3XTRN8H6am7tFAogI/O1MMW5FkSuw5fRL+ThS5NE9RhK77imqm4f7JqqtJqvVcTMAwz7GOjXRPIPeVu8Qw1AOB0ADP7oxccg+5BpLV89ng8snwGWBb44bajcDqHX4d586AdP9p2FADwYt5SfHG5uNtFzJ2SiIxZGvQ7WfzzyKXRL+BBv9OJt90plEdEfr2A6xyCHPfOocH0GbAsi79/7Kq8ErIhzF+SBnu1Wu01z56dnR3QNYHcU+6CKb887D5jc7kMFoeIy4t5SxGtUmLfqat4+g8H0dnTN+jxfqcTxtJzeP7VwwCAb90/H997aKEUQ/Us1Iq1V86h842SpHAGevpu14Zwf//Yjt6+4U9OG8nhqkYcveBAbLRSFlU4HEk/22/atGlQlYzZbIZer/d8bbfbUVxc7Nc1oz0eariSyXN1rX798vX2OVHh3o8l0IYswj/tpHEw6m9DTJQSb5ddxgOb9+LUZQaXGttRbmtE9k9L8YO/VqDfyWLDyjT83xMZki3uPbZiJhQK4ND5Blzm6RCdkbx1xPWmIkUKh/OFjKmYkOhqgnsvwN1K//jBOQCuWb0QjW+BkjTYFxQUgGEYlJSUoKSkBGVlZTAajZ7HLRbLoK99uWa0x0PNjJQEjI2LQk+f06884vGLDnT29CM5IUaUnCHxXd7ts/B2wRokJ8Sg3NaE2370HhZ9722s/skelNuaMC4uClueysT/t3GFpHupTNWM8ey9Yz4s7Oy+r9+Jt464Uji5Auxw6auYKJXnkJTfvnfG75Lny43tnmqe5++dx/v4giH5dgkjbWXAdcH6c40vj4cSpVKBBalJKLM14dRlxueDDz5zp3BWzBkfEsf4RZo750/E3h/fC0PJIZyra0VPnxNOlsUDGakoejITUyU4o9SbvNtn4eDZazAfuojvPihcOumjMze2FZD6jOFnc+bid++fQZmtCYfON/h8gDsAlOw9j34ni7sXTsLiGckCjtJ/kgd7MrpF09UoszXhxKVmn/f1PlzZAABYQSkc2ZozJREf/s99Ug9jRA9nTcf33yjDiUvNOFvbItgpS2b3usAXl02XfFuBSep4fOlOLf6yrwq/ee+Mz8G+vbvPsz30czKb1QMyLL0kN1uW7lpgPXj2mk8/z7KsZ2ZP+XoSjJRxsZ4tC8wCLdR29/bjHXfqQ4y9cHzx7fvnQ6FwbXHh65kSf//YDqajF9qJY3H/rVMFHqH/KNiHgNWLXB9rK+xNaOkYfW+Vi43tuMp0IlqlRKZWI/TwSJjjygdNh/g/cQsA9p68AqajF1OS43H7XHlMTuZMScS6TNd2z79978yoP8+092DzW58DAJ6/b54s+1rkNyKRyfVYwoGmj09A+qRx6HeyPs3uuRTO0lnJiI+hTB0JzgO6VMTHqGC/dh1lNv5O3OLsdH9ieHT5DFkFye886NqaYvunF1B1tXXEn/3ZzhNobOvGvKmJeGb1bDGG5xUdSzgCuR5LOBQ3u99/6uqoP/vZeffiLNXXEx6MjYvGw+4Tt/76Eb+HmnR093kO9xDinNlgLEsfj7WLJ6O337U/0XClz8cvOPDK3koAwC+fzpKsbBSgYwnDwj3uYL/Pl2DvntnL5SMxCX1Pu8sRzZ9d5HUnzHetNWjvdu0MmaUV9nzWQPz+6yuQnBCDCrsDP3enaQZyOll8741yOFkWj62YgbslriQaCQX7EHHXwklQKFzNVSPtmtjS0ePptqXFWcKXO+dPhHbiWFzv6vNsQcyH1/a7Pik8sVIeO0MONS0lAb99ZjkA4Jf/OjUojdrd24/vvl6GI1WNSIiN8ux9JFcU7ENEckIMdGmuxdb9p4ef3e8/dRUsC6RNHItJ6nixhkfCnEKhwFPu2f0bB/g54KPqais+OlMPhQJ42r3jpBw9snwGvrxKC5YF8n61H//5lyPYc7wOOT8txZ/dpZY//5JONr0Rw6FgH0Lu8SFvz23gJMVGUiS8PXlnGpQKBT6rbMC5usBPT+NwWwnn3DIF08cnBH0/If3iqUxkajW43tWHP++rwmO/3I+jFxzQjI3Frh/cg69JuCjrKwr2IWS1J29f77UErqapHaXuPbS5k3cI4cuU5DG4d6mr5v6vHwU3u+/p6/fsDPnVEAiUY+Oi8eGP78O7/7UW+be7zqxdMWc8Pvnp/chZIr+aem+oLi+ErJg9AXHRKlxlOnGurvWmbsa/fmQHywKrFkzE7Mm0Hw7h39N3p+P9Y3V482A1/vvRJYiLCazy5D1rLRpauzApKQ73L725ckSOlEoF7lo4CXctnIStfU5EqRSyXGcYDs3sQ0hcjAp3zne1bv9xz7lBj/U7nZ6yuK/eLf+ZEglN9y9NRapmDBpau4I6nPu1/a5c95fv0kq+PUIgoqOUIRXoAQr2IdFUNdAP1i8C4Ppj4bYwBoB9J6/iclMHkhNisN5dE00I36KjlPiuu9noV/86FdCe75VXWrH3pGvd6Ss0MeEVNVWNIFSaqjgr503EhpWzwLLA914vQ7/T9cf2mnux6/E7ZgX80ZoQXzx9dzomJsXhclMH/vFptd/XF719EgDwhYxUOhuZZ9RUFWZ++ngGEuOjYa124Mfbj+PRLfvwdplrIyluL25ChBIfE4X/+IJrdv/Ld055Jhy+OFvbgh3uQ8x/9MgtQgyPDIOCfQiapI7Hfz+2BIDrUOrSE1egUirw3QcXym4PbRKevr5mNjRjY2Grv45dh31vsnrpn5+DZYF1mdOwdBZt0icmCvYhauPaOVg+ezyiVAo8fXc6rEUP4X8fv1XqYZEIMTYuGt+637Vn+093nkBb5+hbKJy6zGCXu/v2hzSrFx2VXoaoKJUS//7hWvT0OTE2Llrq4ZAIZMiZh7/sq0L1tev4wV/LYdTfPuLP/9+uE2BZV8MffQIVH83sQ1hMlIoCPZFMYnw0Xn32DigVCrx5sBomdy7eG2PpObxTUQOlQoFNX1ws3iCJBwV7QkjA7pg3EYUPu8qBv/NaGez1bTf9TOmJOhT8zQoAeDF/KRZMU4s5ROJGwZ4QEpSChxfjtjkT0NrZi7v+5328urcSTieLfqcTB8/W4yu/Pwgny+Kpu7T4zgMLpB5uxFKwQpwzFkIyMzORmpqKDRs2hEytPSFyU+vowJMvf4QKuwMAMHvyODS0dqGlw7Vwe+f8iXi7YLWkB3tEgm3btmHbtm2ora1FRUXFoMciPtivX78eu3fvlnoYhIS8fqcTr+6txE9Mx9HW1QcAGBcXhZwlU/GrryxDyrhYiUcYObzFNarGIYTwQqVUwpAzD+uzpmP/6XrMm5KIJTOTEaWibLEcULAnhPBqSvIYbFiZJvUwyBD0lksIIRGAgj0hhEQACvaEEBIBKNgTQkgEoGAfoFA57CQY4f4a6fWFvnB/jXy+vogP9oGeVBXuv2RA+L9Gen2hL9xfYyBxiU6qGobYJ1UF88sZ6LVi/0GIPU56ffwL99cYrq+PO6nKm4gP9mKjYC+f6wIV7q8vmOcMldcY7q/P28w+4rdLWLRoEdLT/T/Kr7a21us5j0JdJ8Vz0nWReZ0Uz0nX8Xud1WpFTU3NoO9FfLAnhJBIQGkcQgiJABTsCSEkAlCwJ4SQCEC7XvqpuLgYarUaAMAwDAoKCqQdkACKi4sBADabDQBgNBqlHI6gcnJyUFpaKvUweFdYWOgpPNBoNMjNzZV4RPwqKSkBwzBQq9Ww2WzYtGmT5+8y1NjtdlgsFphMJq+/i7zFHJb4rKioiDUajZ6vS0tLWb1eL+GI+FdQUDDoa71ez2ZnZ0s0GmGZTCY23P4EmpubWZ1OxzY3N7Msy7IVFRVh9xqLioo8r49lXa85NzdXugEFoaKigjUajWxRURGr0+luepzPmBNevwUCU6vVg37JWJYNqz+k5uZmNjs7e9Br5IKFzWaTbmACaG5uZo1GY1j978eyrjfroqKiQd8rLS2VaDTC8BYUQ31CYjKZvL4uPmMO5ex9ZLfbPR8bh7JYLOIPSCDl5eWw2+2er7VaLQDXx8dwsmPHDuTn50s9DN4VFxcjOzvbkxoAgOzsbIlHxS+tVoucnBzP76Tdbvf8noYTvmMOBXsfDQyAA6nV6rAJhGq1Gs3NzdDpdJ7vcb9U4fTHZLFYwi4AAjd+R7kgodVqYTAYwmoyAgCvvPIKHA4HkpOTUVhYCIvFEpbrSnzHHAr2QdJoNHA4HFIPQzCbN2+G0WgM2cUvb7hAGG644KBWq6HT6aDValFUVIS8vDyJR8YvtVoNg8GA3NxcFBcXw2Qyhc2EyxeBxhwK9kEK50BfWFgIg8EAvV4v9VB4U1JSEnaVKRyNRgMAyMrK8nyPmwWG0+y+sLAQWq0WJpMJNpsNDocDmZmZUg9LNIHGHAr2PhpuJhius0Sz2Yz09PSwCvRWq3VQIAw33O/h0GCgVquHTQmEGi5FxaXhtFotKioqoNVqYTabJR4dv/iOOVRn7yOtVuv5oxn6Hzrc8r/cLJAL9AzDwOFwhPybmsPhgNVq9bw+ro+guLgYWq025Gf8arUaWq32pt9RhmHC5k3Obrd7TSmGW6oKECDmBFTDE6GG1ryaTKawq7OvqKhgi4qKWJvNxtpsNraiooItKCi4qfwrHIRjDbrJZBrUK2EymUK+LHGooeXBLMuG/N+h0Wj0qc4+mJhDu176aWA3m81mQ1FRkbQD4hHDMEhLS/O62BVuvyZmsxnbt2+H2WyGXq9HXl5e2HxC47pLAaCpqSmsfkcB1+/p5s2bkZKS4lmT0Ov1IVlEYLfbYTQaYbFYYLVaUVBQcFP6lK+YQ8GeEEIiAC3QEkJIBKBgTwghEYCCPSGERAAK9oQQEgEo2BNCSASgYE8IIRGAgj0hhEQACvaEEBIBKNgTEgSz2SybTcYGds4SMhQFe0ICZLFYZHVKkl6vx8aNG6UeBpEpCvYkogU6K2cYBkajEQUFBTyPKDibNm2CwWCQehhEhijYk4hltVphtVoDupY72EVudDod7Ha7bFJLRD4o2JOIFczpTeXl5bLdJdNgMITdTpckeBTsSUSyWq0oLCwM6Fqz2Szrw0Byc3OxY8cOqYdBZIaCPYk4ZrMZRqMRAGA0GmEwGGAwGHyuZNm+fTtycnKGfbykpMTzz2Aw+HRcntVqRU5ODpKTk1FSUuL5fl5eHpKTk5Genu7T2DgajSbgFBUJUwEdeUJIGADAmkwmv6/T6XSszWbz+pjJZGJzc3MHfU+tVvv8PAAGnUzEsixbUFDAarVav8aYm5t7031IZKOZPSF+stvt0Gg0wz4+dEadnZ2N0tJSn+7t7bQlf2f1gGtmTzX3ZCAK9oT4iWGYYY/Ay83N9RxkzjAMrFar58B2X4z0JuIPtVqNpqYmXu5FwgMFe0J4ZjabkZmZiY0bN8LhcITk2agk/FCwJ8Rt4MLoSLhDroe7x8aNG2EymWAymZCdnR30bD2QdAzDMEhJSQnqeUl4oWBPItbQGbevQVWj0QzbtMTVuA/cQmFgCseXN5Sh4+DSQv5wOByy2caByAMFexKxsrKyUFZWBsC1qKrT6Xy6jutS9QUXuL3938zMzJvKMrOzswcFd4ZhvHbEDnc9R0579hB5oGBPIpbRaITVakVxcTEsFovPHbGPP/74sNU1FRUVKC0tRXFxMUpKSmCxWGAymaBWq1FYWIjc3FzPz9rt9psWbouKiuBwOAZdn5eXBwDIzMwcVOnj7XoOwzA+v3mRyKBgWZaVehCEhJrMzExUVFRIPQyvLBYLSktLacsEMgjN7AkJQFZWVlB76wipqKhIlpu0EWlRsCckAEVFRZ4tF+SEW3ugfD0ZitI4hATIbDbD4XBAr9dLPRSPvLw8mEwmqYdBZIhm9oQEKDc3d8QyTLGVlJTglVdekXoYRKZoZk8IIRGAZvaEEBIBKNgTQkgEoGBPCCERgII9IYREgP8fNeTUNlh7B3cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the outcome\n",
    "from matplotlib import pyplot as plt\n",
    "import mpl_config\n",
    "\n",
    "plt.plot(tau * np.arange(101), np.array(z_list))\n",
    "plt.xlabel(\"t (a.u.)\")\n",
    "plt.ylabel(r\"$\\langle \\sigma_z \\rangle$\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}