docs/notebooks/09_example_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization with numpy Arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-12T04:19:15.953542Z",
     "start_time": "2019-02-12T04:19:14.691832Z"
    },
    "attributes": {
     "classes": [],
     "id": "",
     "n": "1"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy            1.17.2\n",
      "scipy            1.3.1\n",
      "matplotlib       3.1.2\n",
      "krotov           1.0.0\n",
      "matplotlib.pylab 1.17.2\n",
      "CPython 3.7.3\n",
      "IPython 7.10.2\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "%load_ext watermark\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib\n",
    "import matplotlib.pylab as plt\n",
    "import krotov\n",
    "# note that qutip is NOT imported\n",
    "%watermark -v --iversions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{tr}[0]{\\operatorname{tr}}\n",
    "\\newcommand{diag}[0]{\\operatorname{diag}}\n",
    "\\newcommand{abs}[0]{\\operatorname{abs}}\n",
    "\\newcommand{pop}[0]{\\operatorname{pop}}\n",
    "\\newcommand{aux}[0]{\\text{aux}}\n",
    "\\newcommand{opt}[0]{\\text{opt}}\n",
    "\\newcommand{tgt}[0]{\\text{tgt}}\n",
    "\\newcommand{init}[0]{\\text{init}}\n",
    "\\newcommand{lab}[0]{\\text{lab}}\n",
    "\\newcommand{rwa}[0]{\\text{rwa}}\n",
    "\\newcommand{bra}[1]{\\langle#1\\vert}\n",
    "\\newcommand{ket}[1]{\\vert#1\\rangle}\n",
    "\\newcommand{Bra}[1]{\\left\\langle#1\\right\\vert}\n",
    "\\newcommand{Ket}[1]{\\left\\vert#1\\right\\rangle}\n",
    "\\newcommand{Braket}[2]{\\left\\langle #1\\vphantom{#2} \\mid\n",
    "#2\\vphantom{#1}\\right\\rangle}\n",
    "\\newcommand{op}[1]{\\hat{#1}}\n",
    "\\newcommand{Op}[1]{\\hat{#1}}\n",
    "\\newcommand{dd}[0]{\\,\\text{d}}\n",
    "\\newcommand{Liouville}[0]{\\mathcal{L}}\n",
    "\\newcommand{DynMap}[0]{\\mathcal{E}}\n",
    "\\newcommand{identity}[0]{\\mathbf{1}}\n",
    "\\newcommand{Norm}[1]{\\lVert#1\\rVert}\n",
    "\\newcommand{Abs}[1]{\\left\\vert#1\\right\\vert}\n",
    "\\newcommand{avg}[1]{\\langle#1\\rangle}\n",
    "\\newcommand{Avg}[1]{\\left\\langle#1\\right\\rangle}\n",
    "\\newcommand{AbsSq}[1]{\\left\\vert#1\\right\\vert^2}\n",
    "\\newcommand{Re}[0]{\\operatorname{Re}}\n",
    "\\newcommand{Im}[0]{\\operatorname{Im}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `krotov` package heavily builds on QuTiP. However, in rare circumstances\n",
    "the overhead of `qutip.Qobj` objects might limit numerical efficiency, in\n",
    "particular when QuTiP's automatic sparse storage is inappropriate. If you know\n",
    "what you are doing, it is possible to replace `Qobj`s with low-level objects\n",
    "such as numpy arrays.  This example revisits the [Optimization of a\n",
    "State-to-State Transfer in a Two-Level-System](01_example_simple_state_to_state.ipynb),\n",
    "but exclusively uses numpy objects for states and operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-level-Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider again the standard Hamiltonian of a two-level system, but now we\n",
    "construct the drift Hamiltonian `H0` and the control Hamiltonian `H1` as numpy\n",
    "matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian(omega=1.0, ampl0=0.2):\n",
    "    \"\"\"Two-level-system Hamiltonian\n",
    "\n",
    "    Args:\n",
    "        omega (float): energy separation of the qubit levels\n",
    "        ampl0 (float): constant amplitude of the driving field\n",
    "    \"\"\"\n",
    "    # .full() converts everything to numpy arrays\n",
    "    H0 = -0.5 * omega * np.array([[-1, 0], [0, 1]], dtype=np.complex128)\n",
    "    H1 = np.array([[0, 1], [1, 0]], dtype=np.complex128)\n",
    "\n",
    "    def guess_control(t, args):\n",
    "        return ampl0 * krotov.shapes.flattop(\n",
    "            t, t_start=0, t_stop=5, t_rise=0.3, func=\"sinsq\"\n",
    "        )\n",
    "\n",
    "    return [H0, [H1, guess_control]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "H = hamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization target\n",
    "\n",
    "By default, the `Objective` initializer checks that the objective is expressed with\n",
    "QuTiP objects. If we want to use low-level objects instead, we have to\n",
    "explicitly disable this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "krotov.Objective.type_checking = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we initialize the initial and target states,\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ket0 = np.array([[1], [0]], dtype=np.complex128)\n",
    "ket1 = np.array([[0], [1]], dtype=np.complex128)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and instantiate the `Objective` for the state-to-state transfer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Objective[a₀[2,1] to a₁[2,1] via [a₂[2,2], [a₃[2,2], u₁(t)]]]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objectives = [krotov.Objective(initial_state=ket0, target=ket1, H=H)]\n",
    "\n",
    "objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how all objects are numpy arrays, as indicated by the symbol `a`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate dynamics under the guess field\n",
    "\n",
    "To simulate the dynamics under the guess pulse, we can use the objective's\n",
    "`propagator` method. However, the propagator we use must take into account the\n",
    "format of the states and operators. We define a simple propagator that solve\n",
    "the dynamics within a single time step my matrix exponentiation of the\n",
    "Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expm(H, state, dt, c_ops=None, backwards=False, initialize=False):\n",
    "    eqm_factor = -1j  # factor in front of H on rhs of the equation of motion\n",
    "    if backwards:\n",
    "        eqm_factor = eqm_factor.conjugate()\n",
    "    A = eqm_factor * H[0]\n",
    "    for part in H[1:]:\n",
    "        A += (eqm_factor * part[1]) * part[0]\n",
    "    return scipy.linalg.expm(A * dt) @ state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will want to analyze the population dynamics, and thus define the projectors\n",
    "on the ground and excited levels, again as numpy matrices:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj0 = np.array([[1, 0],[0, 0]], dtype=np.complex128)\n",
    "proj1 = np.array([[0, 0],[0, 1]], dtype=np.complex128)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will pass these as `e_ops` to the `propagate` method, but since `propagate`\n",
    "assumes that `e_ops` contains `Qobj` instances, we will have to teach it how to\n",
    "calculate expectation values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expect(proj, state):\n",
    "    return complex(state.conj().T @ (proj @ state)).real"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can simulate the dynamics over a time grid from $t=0$ to $T=5$ and plot\n",
    "the resulting dynamics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "tlist = np.linspace(0, 5, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-12T04:19:16.422181Z",
     "start_time": "2019-02-12T04:19:16.265626Z"
    },
    "attributes": {
     "classes": [],
     "id": "",
     "n": "12"
    }
   },
   "outputs": [],
   "source": [
    "guess_dynamics = objectives[0].propagate(tlist, propagator=expm, e_ops=[proj0, proj1], expect=expect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-12T04:19:16.429662Z",
     "start_time": "2019-02-12T04:19:16.424216Z"
    },
    "attributes": {
     "classes": [],
     "id": "",
     "n": "13"
    }
   },
   "outputs": [],
   "source": [
    "def plot_population(result):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(result.times, result.expect[0], label='0')\n",
    "    ax.plot(result.times, result.expect[1], label='1')\n",
    "    ax.legend()\n",
    "    ax.set_xlabel('time')\n",
    "    ax.set_ylabel('population')\n",
    "    plt.show(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-12T04:19:16.606020Z",
     "start_time": "2019-02-12T04:19:16.433034Z"
    },
    "attributes": {
     "classes": [],
     "id": "",
     "n": "14"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXAc533m8e8P5+C+eeEgeJmURFEiCUpyvLG92TiWFVtKJetYyrFRpESpXSvrbDbeUmqzTuStrXKOTSVrO9koiTdrZW2VYycxIyvyqjZSHMc6CB6SSEqUaF4AeAEg7vv47R9vAxiCBDEkMRgA/Xyqunqmu2fmHXLQT/f7vv22uTsiIhJfWZkugIiIZJaCQEQk5hQEIiIxpyAQEYk5BYGISMzlZLoA16u6utobGxszXQwRkWVl//79He5ec7V1yy4IGhsbaW5uznQxRESWFTM7Pdc6VQ2JiMScgkBEJOYUBCIiMacgEBGJOQWBiEjMpS0IzOxLZnbRzA7Psd7M7H+Y2XEze8PMdqWrLCIiMrd0nhH8BXDvNdZ/BNgSTY8Bf5zGsoiIyBzSdh2Bu3/HzBqvsckDwJc9jIP9ipmVm9ladz+XjvLsO3WJf3qnHczIMjDCPCvLAMgyw4zpdeHx5fOcrCwK87IpyMumKC+HgrxsCqPHhfnZlBXkkput2jYRWV4yeUFZLdCS9Lw1WnZFEJjZY4SzBhoaGm7oww6c7uLzLx4n3bdfKMnPobwol4rCPMoKwryiMJeaknxWlyZYXZpgTVmC1SUJSgtyMLP0FkhEZB7L4spid38KeAqgqanphnblv/SBTfzSBzZNvR+TPjOfjNJh0h2PnjvgkzOPJ90Zn3AGR8cZHJ2IpnGGoscDo+N0D47RNThK18AoXYNjdA+OcrpzkK7BUfqGx68oUyI3i9WlCdaWJaivKKShspD66amAmuJ8BYWIpF0mg6ANqE96XhctSzszI9sAFm8nOzw2wcXeEc73DnPhsmmEs91DfOfddi70jlz2moLcbOorC1hfVcSmmmI2rwrTppoiShK5i1Z2EVnZMhkEe4HHzewZ4G6gJ13tA0tBIjebhqpCGqoK59xmeGyC1q5BWi4NcebS4PR0qmOAl45dZGxi5mRodWl+CIYoILatLWXbmhIFhIhct7QFgZl9FfggUG1mrcBvArkA7v4/geeA+4DjwCDw8+kqy3KRyM1m86oSNq8quWLd2MQkZy4N8v2L/Rxv7+f4xX6+3z7ANw600T8yU+3UUFnILWtLuHVtWZivK6W2vEBVTCIyJ1tuN69vampyjT46w9053zvM2+f6OHqul6PnennrbC8nOwemG8ZLEzncsraU22vL2FFfzp115dRXKhxE4sTM9rt709XWLYvGYpmbmbG2rIC1ZQX8y22rppcPjo7z9vk+3jrXy9GzISCefuU0I989CUBFYS476sq5o64szOvLqSnJz9TXEJEMUhCsUIV5OexqqGBXQ8X0srGJSd650MfrLT280drNoZZuvvBiO5PRmcO6sgR31JdzZ305u9dXsL22jERudoa+gYgsFlUNxdzg6DhHzvbyeks3r7f28HpLN2cuDQKQl53F9tpSdq+vYPf6Snavr9BZg8gC6xka41zPEO19I7T3jdDZP8rQ2ASj45OMTkwyMjYxfbD2sTvWcdeGyhv6HFUNyZwK83LY01jJnsaZH1dH/wgHTnexP5r+98un+dN/ClVK66sK2d1Qwe7GCnavr2DLqhKys9TWIDIfd+fMpUH2nerizdZu3r3Yz7sX+2nvG7nq9lkGeTlZ5GVnkRONWLCjruyGg+BadEYg8xoZn+BwWy8HTnfRfPoS+0930dE/CoQrqXeur6BpfQVNjRXcWV9OYZ6OL0TGJyY5eq6X5lPh72bfqa7pnX5xfs70dUGbVxXTUFlIdXE+NSX5VBXnUZibPb3zXyjXOiNQEMh1mzqy2Z901nDsQh/ukJNl3LaulKbGSvY0hiolVSdJHAyMjHPwTDf7Tl2i+fQlDp7pZnB0AoC6igL2NFbS1FjBnsZKNtcUT49ztlgUBJJ2PUNjHDjTRfOpcOTzeks3I+OTADRWFU4HQ1NjJRuri9R1VZa9i73DNJ/uCjv+U10cPdfLxKSTZXDL2tLoLDns/NeWFWS6uAoCWXyj45McPtszHQzNpy7RNTgGQGVRHrvXV0wHw/Z1ZeTlaNRWWbrcnRMdAzSfusRrJ0NVz+nO0KkikZvFzvqZ3/POhvIleYW/gkAyLvkPaSoYTkV/SPk5WdxRXz79h7SroYKygqX3hyTxMVW//9rJS9N1/FPtYlXRgcxdGyppaqzktnWly2L4eQWBLEntfSPsPz0TDEfO9jI+6ZjB1tUl0/WpTY2V1JZn/tRaVq6h0QkOtnSx72So6jlwpmu6fr+hspCmxgruaqxkz4blW7WpIJBlYXB0nEMt3TSfCn+MB890T4+jtK4swe6pdob1lWxdo26rcuO6BkajRt0uXjt5icNtPdMHIdvWlLInOgjZ01jJmrJEpou7IBQEsixNTDpvn++dDobmU12c7x0GQrfVXdPdViu5s76cgjxdBS1Xmph0jl/s58CZLg6c7uJgSzfHL/YD4aLJO+rLpnf6u9av3GpJBYGsCO5OW/fQZcFw7EIfELqtbq8tu6ynRnWxuq3GUffgKAfPdHPgTBcHz4ShVKbOLCsKc9nZUDFdx397jIZRURDIitUzGLqtTgXDodZuRqNuqxuri9i9voJd6yu4vbaM96wuUe+kFWZ0fJJj5/t4vbWbg2e6OXimixMdA0C4MnfbmlJ2rS9nV0MFOxsqaKwqXJb1+wtBQSCxMXUV9FTvpP2nZ7qt5mVnsW1tCdtry7g9mhQOy8fUTv/Nth7ebOvhcFsPx873MToRgr+qKI+dDRXsbAg7/h11ZRTl6yr3KQoCia2pq6DfbOvhzdae6Z3I1D2kZ4fD1jUlbF1doh1IhvWPjPPOhT7ePtd31Z1+aSKH2+vKpv/fdtTqHhvzURCIJHF3TncOTofCm609HD47Ew4A9ZUFbF0dbv+5dU0J29aUsKG6aMHHf4m7sYlJTnYM8Pb5Po6d7+XY+T7ePt9Ha9fQ9DbJO/0dteXcXlumnf4NUBCIzGNy0mnpGuTY+b6wM7oQ5ic7BpiIxgDOy85iY00Rm2qK2VBdRGN1ERuqi9hYXURFUV6Gv8HS1jUwyomOfk60D3CyY2B6frJjYPooPzvL2FhdNB28W9eEIK6r0E5/IWgYapF5ZGUZ66uKWF9VxI/ctmZ6+fDYBN9v7w/VFOf7eOd8uCXo80fOTwcEQFlB7nQo1FcWUltRQG15AevKC1hblljxPVOGxyY42z1EW/cQrV1DtHWFx6c7BzjRMUB31E4DoYdXQ1UhG6uL+OC2mrDTX13KplVF5Oes7H+npUpBIHINidxsbltXxm3ryi5bPjYxSculQU51hqPbU53h6PaVE538zaE2Zp9oVxfnU1ueYF15AatLE1QX51FdnE91cRh2eGoI4qUWGMNjE3QOjNIR3TSlvX9k+gYq7X0jnO8dpq176Iox9bOzjDWlCeoqCvjI9rVsqonOnmqKqasoWBZDMsSJgkDkBuRmZ7GxppiNNcX80LbL142OT3Ih2kG2dQ1xtnuIsz3hSPnYhT6++24HfSPjV33fwrxsShI5lCZyKUnkUJI0L03kkJ+bTX50s5K8nKzwOJqyzZjKn6kg8mjJxKQzMjbJ8PhEmI9NTD8eGpugb3ic3uExeobG6B0ao3d4nJ6hsemuuLNVFOZSU5LPqpIEP7R1FbUVBdRFZ0G1FQWsKU2oPWUZURCILLC8nCzqKwupryycc5vkI+2O/nB7wvb+ES4NjNI3PEbf8Dh9w+N0D47ScmmQ3uFx+obHpof2XghZFs54ErnZlCZyKC3IpTSRy7qyAkoLQhiVFuRSVZRHTUn+9FRVlK8utyuMgkAkAxK52eHo+ToH03N3xiackfHke9qG+WR0GmCEhtWp9lUDzIxEbhaJ6IwikZut6hmZpiAQWUbMjLwc0xG5LCj9mkREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCIiMacgEBGJubQGgZnda2bHzOy4mT1xlfUNZvaimR00szfM7L50lkdERK6UtiAws2zgi8BHgFuBh8zs1lmb/QbwNXffCTwI/FG6yiMiIleXzjOCu4Dj7n7C3UeBZ4AHZm3jQGn0uAw4m8byiIjIVaQzCGqBlqTnrdGyZL8F/IyZtQLPAb98tTcys8fMrNnMmtvb29NRVhGR2Mp0Y/FDwF+4ex1wH/C0mV1RJnd/yt2b3L2ppqZm0QspIrKSpTMI2oD6pOd10bJkjwJfA3D3l4EEUJ3GMomIyCzpDIJ9wBYz22BmeYTG4L2ztjkD/CsAM7uFEASq+xERWURpCwJ3HwceB74NvEXoHXTEzD5rZvdHm/1H4BfN7HXgq8DD7rNv8iciIumU1vsRuPtzhEbg5GWfSXp8FHhfOssgIiLXlunGYhERyTAFgYhIzCkIRERiTkEgIhJzCgIRkZhTEIiIxJyCQEQk5hQEIiIxpyAQEYk5BYGISMwpCEREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCIiMacgEBGJOQWBiEjMKQhERGJOQSAiEnMKAhGRmFMQiIjEnIJARCTmFAQiIjGnIBARiTkFgYhIzCkIRERiTkEgIhJzCgIRkZjLSXVDM8sGVie/xt3PpKNQIiKyeFI6IzCzXwYuAC8A34qmZ1N43b1mdszMjpvZE3Ns85NmdtTMjpjZV66j7CIisgBSPSP4FLDV3TtTfePoDOKLwIeAVmCfme1196NJ22wBfh14n7t3mdmq1IsuIiILIdU2ghag5zrf+y7guLufcPdR4BnggVnb/CLwRXfvAnD3i9f5GSIicpNSPSM4AbxkZt8CRqYWuvvvX+M1tYQAmdIK3D1rm/cAmNk/A9nAb7n787PfyMweAx4DaGhoSLHIIiKSilSD4Ew05UXTQn7+FuCDQB3wHTO73d27kzdy96eApwCampp8AT9fRCT2UgoCd38SwMyKo+f9KbysDahPel4XLUvWCrzq7mPASTN7hxAM+1Ipl4iI3LxUew1tN7ODwBHgiJntN7Pb5nnZPmCLmW0wszzgQWDvrG3+lnA2gJlVE6qKTlxH+UVE5CalWjX0FPCr7v4igJl9EPhT4AfmeoG7j5vZ48C3CfX/X3L3I2b2WaDZ3fdG637EzI4CE8Cnr6dnkojIYhobG6O1tZXh4eFMF2VOiUSCuro6cnNzU36Nuc9f5W5mr7v7HfMtWwxNTU3e3Ny82B8rIsLJkycpKSmhqqoKM8t0ca7g7nR2dtLX18eGDRsuW2dm+9296WqvS7X76Akz+y9m1hhNv4GqcEQkZoaHh5dsCACYGVVVVdd9xpJqEDwC1AB/HU010TIRkVhZqiEw5UbKl1IQuHuXu/97d98VTZ+aughMREQWz/PPP8/WrVvZvHkzn/vc5xbkPa/ZWGxmf+Duv2Jmfwdc0Zjg7vcvSClERGReExMTfPKTn+SFF16grq6OPXv2cP/993Prrbfe1PvO12vo6Wj+ezf1KSIictNee+01Nm/ezMaNGwF48MEH+eY3v5neIHD3/dHDO939D5PXmdmngH+8qU8XEVmmnvy7Ixw927ug73nrulJ+82NzX6LV1tZGff3Mdbp1dXW8+uqrN/25qTYW/9xVlj18058uIiIZN18bwUPATwEbzCz5quAS4FI6CyYispRd68g9XWpra2lpmRnLs7W1ldra2pt+3/naCL4HnAOqgf+etLwPeOOmP11ERFK2Z88e3n33XU6ePEltbS3PPPMMX/nKzd/Pa742gtPAaeC9N/1JIiJyU3JycvjCF77Ahz/8YSYmJnjkkUe47babPzNJaawhM7sH+DxwC2EY6mxgwN1Lb7oEIiKSsvvuu4/77rtvQd8z1cbiLwAPAe8CBcAvEG5DKSIiy1yqQYC7Hwey3X3C3f8XcG/6iiUiIosl1WGoB6N7Chwys98hNCCnHCIiIrJ0pboz/1lCu8DjwADhzmM/ka5CiYjI4kn1VpWno4dDwJPpK46IiCy2+S4oe5OrDDY3xd13LHiJRERkUc13RvDRRSmFiIik5JFHHuHZZ59l1apVHD58eEHe85ptBO5++lrTgpRARERS9vDDD/P8888v6Hum1FhsZn1m1htNw2Y2YWYLO+yeiIjM6/3vfz+VlZUL+p6pNhaXTD22cB+0B4B7FrQkIiLLyd8/AeffXNj3XHM7fGRh7jp2Pa77WgAP/hb4cBrKIyIiiyzVsYZ+POlpFtAEDKelRCIiy0EGjtzTJdUriz+W9HgcOEWoHhIRkWUu1TaCn093QUREZH4PPfQQL730Eh0dHdTV1fHkk0/y6KOP3tR7plo1tBH4Q0IDsQMvA//B3U/c1KeLiMh1+epXv7rg75lqY/FXgK8Ba4F1wF8BC18aERFZdKkGQaG7P+3u49H0l0AinQUTEZHFkWpj8d+b2RPAM4SqoU8Az5lZJYC760b2IiLLVKpB8JPR/JdmLX+QEAwbF6xEIiJLmLsTrqtdmtznHCd0Tqn2Gtpw3e8sIrLCJBIJOjs7qaqqWpJh4O50dnaSSFxfzX2qvYZygX8LvD9a9BLwJ+4+Ns/r7iX0NsoG/szdr3oFhpn9BPB1YI+7N6dWdBGRxVVXV0drayvt7e2ZLsqcEokEdXV11/WaVKuG/hjIBf4oev6z0bJfmOsFZpZNuMH9h4BWYJ+Z7XX3o7O2KwE+Bbx6XSUXEVlkubm5bNiw8ipIUg2CPe5+R9LzfzCz1+d5zV3A8alrDczsGcLVyEdnbfdfgd8GPp1iWUREZAGl2n10wsw2TT2JLjCbmOc1tUBL0vPWaNk0M9sF1Lv7t671Rmb2mJk1m1nzUj4lExFZjlI9I/g08KKZTV1J3Ajc1LATZpYF/D7w8HzbuvtTwFMATU1N198kLiIic0r1jOCfgT8BJoFL0eOX53lNG1Cf9LwuWjalBNgOvGRmpwjDV+w1s6YUyyQiIgsg1SD4MrCBUJ//ecJ1A0/P85p9wBYz22BmeYRrDvZOrXT3HnevdvdGd28EXgHuV68hEZHFlWrV0HZ3vzXp+YtmNrvR9zLuPm5mjwPfJnQf/ZK7HzGzzwLN7r73Wq8XEZHFkWoQHDCze9z9FQAzuxuY98jd3Z8Dnpu17DNzbPvBFMsiIiILKNUg2A18z8zORM8bgGNm9ibh7pU70lI6ERFJu1SD4N60lkJERDIm1bGGTqe7ICIikhmp9hoSEZEVSkEgIhJzCgIRkZhTEIiIxJyCQEQk5hQEIiIxpyAQEYk5BYGISMwpCEREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCIiMacgEBGJOQWBiEjMKQhERGJOQSAiEnMKAhGRmFMQiIjEnIJARCTmFAQiIjGnIBARiTkFgYhIzCkIRERiTkEgIhJzCgIRkZhLaxCY2b1mdszMjpvZE1dZ/6tmdtTM3jCz/2dm69NZHhERuVLagsDMsoEvAh8BbgUeMrNbZ212EGhy9x3A14HfSVd5RETk6tJ5RnAXcNzdT7j7KPAM8EDyBu7+orsPRk9fAerSWB4REbmKdAZBLdCS9Lw1WjaXR4G/v9oKM3vMzJrNrLm9vX0BiygiIkuisdjMfgZoAn73auvd/Sl3b3L3ppqamsUtnIjICpeTxvduA+qTntdFyy5jZj8M/GfgA+4+ksbyiIjIVaTzjGAfsMXMNphZHvAgsDd5AzPbCfwJcL+7X0xjWUREZA5pCwJ3HwceB74NvAV8zd2PmNlnzez+aLPfBYqBvzKzQ2a2d463ExGRNEln1RDu/hzw3Kxln0l6/MPp/HwREZnfkmgsFhGRzFEQiIjEnIJARCTmFAQiIjGnIBARibm09hoSkUUyOQE+GT2xaGYzz82SnotcTkEgkinuMNoP/RdhoAMGLsJgJwz3wHAvjPReOR8fhvGRME2MwPhoWOYT839eTmJmyk1ATgHk5ENuISRKIVEWTeUzjwvKobAaildB8WrIL07/v4ssOgWBSDpMTsJgB3S3QE8L9LTOzPvOwUA79LfD+NDVX29ZkF8adtD5ZWFeug5yCyA7H3Lyonk0ZedDVhb41BtEDzyaT45HIRJNY8Phs8dHYGwIes/CxbdguDuEzswbXS63aCYUildByRooq4fyeihvgLIGKKrW2ccyoyAQuVET49BzBjq/H03Hw9R9GnrawhF7srzisNMsWQNVm6GoJuxMi2qgaBUU10BhVTgSzyvO3M50chJG+2CoOwTDQEcUXBfC2Uv/hTB1vAMn/hFGei5/fU4BlNWFcKhoDN+1agtUbw5Bka3dzlKj/xGR+YwOQvvbcPFoOGqe2ul3nYLJsZnt8sugahOsvRNu+VjY6ZfVRVN92MEvhyPlrKyZqiFSuGngUHc425k6++k+MzO17Q9VXdPvnQuVG0IwVG2Cmq2w6lao2QZ5hWn7SnJtCgKRKRPjcOlEtMM/CheOhPmlk0xXleQkoHITrLoFbvlodLQbTYVVy2NHv9AKysO05vYr17mHdo/O49Dx7sxZU+dxOP4CTIxGGxpUboTVt8Kq22bmlRsgK3tRv04cKQgknibGwlH+2YNw9hCcOxR2/OPDYb1lRTum7bDjE+GodfVtoapDO6bUmYU2g6JqaLjn8nWTEyFkLx6BC0dn5m89y0zwFsCa7eEsa92dsG4nVG9V9dICM/c5GoWWqKamJm9ubs50MWQ5mRgPO/1zh2Z2/OffnKnDzyuBtXeEaXV0NFqzLTTMyuIbHYSOYyEULhyGc6+HabQ/rJ8dDmvvDP9fCodrMrP97t501XUKAllxBi9By6tw5mU480rYiUwd6U/t9KeOLtfeGY78s3Rt5ZI2OQmXvj9z9nb2IJx7IzRqQwiHdTuhfg/U7YG6u6BkdWbLvMQoCGTlcg+NtmdemdnxdxwL67Jyw86hbk+Yr7sz1O9rp78yJIfD2QPQui+E/lS7Q3nDTCjU74HVt4dutzF1rSDQuZQsLxPjcP6Ny4/4+y+EdYkyqL8b7vgENLw37PxVvbNyZWVB9ZYw7fh4WDY2HH4frfug5bXw+zj8jbAuJxHOAOv3hN9J/T2hy67ojECWuJG+8Ed9JtrxtzbD2EBYV94Q/pgb7gk7/pptOtqXK/W0hd/QVDicOzRz1lC1Bda/N/x+Gt4bOgOs0J5fOiOQ5aP3XNjhTx3xn38zjKFjWaEhd+dPhx1//T1QVpvp0spyUFYbptt+LDwfHwnVSWe+F84Yjn4TDnw5rCtZGx1Y/EAIiFW3xqKXmIJAMmdyMtTnn3llpo6/+3RYl1MAdU3wg78W/jDr9oRhFkRuVk4+NNwdJgi/w/a3wu/v9MthfuRvwrr8Mqi/a+asYd2uME7TCqOqIVk8Y8Oht0fLKzM7/+HusK6oZqaKp+EeWLMDsnMzW16JJ/dwhfTpl2fOGtrfDuuy86B298xvtf4uKKjIbHlTpF5DkhnT3Tijnf7ZA5fXzSbv+Cs3rti6WVkBBjqjA5jorOHcoTCQH4Tqo+Tfcln9kvwtKwgk/ZK7cU4d8U8dRWXlhq6b00dRd4crTUWWq9HBMI7SVJVmy2sz1zSU1l4eDEuknUGNxbLwJsbCBT0tr0Y7/leh/3xYN1WvevvHwx9D7S5145SVJa8QNvxgmCAMl3HhyEwwnP7eTLfV/NLw9zAVDut2LbkB9nRGIKkZvBSOelpeDVPbgZmx9MsaLm9Qq7lF3Tgl3tzD6KvJZ8gXj4Z1WblhiIyaW6DmPaHbc/n6MCR5QWXa/nZUNSTXxz2MDnnmlZkdf8c7YV1WTmjIbbgn7Pzr7w43TBGRaxu8FF0T83I4kOp4J9ykKJllh2rT3MLohkN54UK4qaql9z4eRr29Aaoakmsb6AwNuW37ww+0dR8MXQrrEuXR1boPhvkSPK0VWRYKK+E9Hw7TlKHuMDx3b2vSTX8uztySdGLqVqTR/ajT1NagIIibkf4wHkvyjn+q7z4WTlO33hf6WNffHXr3qJpHJD0KysOQF+zJaDEUBCvZUHcYxvf84XCF7tkDoSfP1NFFeUM4wt/zC6FBd+0dkF+S2TKLyKJTEKwEk5PQfSrs8C9EO/3zh8P9dKcU1YRB2G65P1wQs26nBtwSEUBBsLxMjIW++u3HwtAMHe9Gj9+d6cNsWaE6p/4u2PNIuH3g6ts1NruIzElBsNRMjEFPa9jhd50K9fedx6H9nXA/3eSbpZesC93P7vypMCDbmu3h4hX12ReR65DWIDCze4E/BLKBP3P3z81anw98GdgNdAKfcPdT6SxTRk3dyLvvHPSdh96zYd7TEnb4XafCkLk+MfOarFyoWB/u07rtvjCveU846tcgbCKyANIWBGaWDXwR+BDQCuwzs73ufjRps0eBLnffbGYPAr8NfCJdZVoQ7jA2BGODMDoQzQdD1cxQV+grPHgpdL9MnvdfDAGQfEQ/pWhVGAe9/h7YsT48Lo/mpeuWxOXpIrJypfOM4C7guLufADCzZ4AHgOQgeAD4rejx14EvmJl5Oq5yO/A0fO/zUY8ZD3NPms9eNv08adnkeNj5k0Lx8orDqIQFFaH/cNWmMNZ5yVooXTvzuHh1rG+fJyKZl84gqAVakp63AnfPtY27j5tZD1AFdCRvZGaPAY8BNDQ03FhpCith1S1hVEDLAqK5Zc1aZjPPL9vOQjVNXmG46i+veNbjopmdfkFFuCpQRGQZWBaNxe7+FPAUhCEmbuhNtv1omERE5DLpvGS0DahPel4XLbvqNmaWA5QRGo1FRGSRpDMI9gFbzGyDmeUBDwJ7Z22zF/i56PG/Bv4hLe0DIiIyp7RVDUV1/o8D3yZ0H/2Sux8xs88Cze6+F/hz4GkzOw5cIoSFiIgsorS2Ebj7c8Bzs5Z9JunxMPDxdJZBRESuTcNKiojEnIJARCTmFAQiIjGnIBARiblld89iM2sHTs+74dVVM+uq5RjQd44Hfed4uJnvvN7dr3oTkmUXBDfDzJrnunnzSqXvHA/6zvGQru+sqiERkZhTEIiIxFzcguCpTBcgA/Sd40HfOR7S8p1j1UYgIiJXist4Ld0AAAPESURBVNsZgYiIzKIgEBGJudgEgZnda2bHzOy4mT2R6fKkm5l9ycwumtnhTJdlsZhZvZm9aGZHzeyImX0q02VKNzNLmNlrZvZ69J2fzHSZFoOZZZvZQTN7NtNlWQxmdsrM3jSzQ2bWvODvH4c2AjPLBt4BPkS4ZeY+4CF3P3rNFy5jZvZ+oB/4srtvz3R5FoOZrQXWuvsBMysB9gM/tsL/nw0ocvd+M8sFvgt8yt1fyXDR0srMfhVoAkrd/aOZLk+6mdkpoMnd03IBXVzOCO4Cjrv7CXcfBZ4BHshwmdLK3b9DuMdDbLj7OXc/ED3uA94i3Bd7xfKgP3qaG00r+ujOzOqAHwX+LNNlWSniEgS1QEvS81ZW+A4i7sysEdgJvJrZkqRfVE1yCLgIvODuK/07/wHwn4DJTBdkETnwf81sv5k9ttBvHpcgkBgxs2LgG8CvuHtvpsuTbu4+4e53Eu4LfpeZrdiqQDP7KHDR3fdnuiyL7F+4+y7gI8Ano6rfBROXIGgD6pOe10XLZIWJ6sm/Afwfd//rTJdnMbl7N/AicG+my5JG7wPuj+rMnwF+yMz+MrNFSj93b4vmF4G/IVR3L5i4BME+YIuZbTCzPMK9kfdmuEyywKKG0z8H3nL33890eRaDmdWYWXn0uIDQIeLtzJYqfdz91929zt0bCX/H/+DuP5PhYqWVmRVFnR8wsyLgR4AF7Q0YiyBw93HgceDbhAbEr7n7kcyWKr3M7KvAy8BWM2s1s0czXaZF8D7gZwlHiYei6b5MFyrN1gIvmtkbhAOeF9w9Fl0qY2Q18F0zex14DfiWuz+/kB8Qi+6jIiIyt1icEYiIyNwUBCIiMacgEBGJOQWBiEjMKQhERGJOQSByDWZWbmb/Lnq8zsy+nukyiSw0dR8VuYZozKJn4zKCq8RTTqYLILLEfQ7YFA3q9i5wi7tvN7OHgR8DioAtwO8BeYQL2kaA+9z9kpltAr4I1ACDwC+6+4q98leWJ1UNiVzbE8D3o0HdPj1r3Xbgx4E9wH8DBt19J+GK7n8TbfMU8Mvuvhv4NeCPFqXUItdBZwQiN+7F6L4HfWbWA/xdtPxNYEc0CuoPAH8VhkECIH/xiylybQoCkRs3kvR4Mun5JOFvKwvojs4mRJYsVQ2JXFsfUHIjL4zuhXDSzD4OYXRUM7tjIQsnshAUBCLX4O6dwD+b2WHgd2/gLX4aeDQaOfIIK/wWqbI8qfuoiEjM6YxARCTmFAQiIjGnIBARiTkFgYhIzCkIRERiTkEgIhJzCgIRkZj7//ELfbRWuit5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_population(guess_dynamics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is the same as in the original example.\n",
    "\n",
    "## Optimize\n",
    "\n",
    "First, we define the update shape and step width as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def S(t):\n",
    "    \"\"\"Shape function for the field update\"\"\"\n",
    "    return krotov.shapes.flattop(\n",
    "        t, t_start=0, t_stop=5, t_rise=0.3, t_fall=0.3, func='sinsq'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "pulse_options = {\n",
    "    H[1][1]: dict(lambda_a=5, update_shape=S)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the optimization with only small additional adjustments. This is\n",
    "because Krotov's method internally does very little with the states and\n",
    "operators: nearly all of the numerical effort is in the propagator, which we\n",
    "have already defined above for the specific use of numpy arrays.\n",
    "\n",
    "Beyond this, the optimization only needs to know three things: First, it must\n",
    "know how to calculate and apply the operator $\\partial H/\\partial \\epsilon$. We\n",
    "can easily teach it how to do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mu(objectives, i_objective, pulses, pulses_mapping, i_pulse, time_index):\n",
    "    def _mu(state):\n",
    "        return H[1][0] @ state\n",
    "    return _mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, the pulse updates are calculated from an overlap of states, and we\n",
    "define an appropriate function for numpy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def overlap(psi1, psi2):\n",
    "    return complex(psi1.conj().T @ psi2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Third, it must know how to calculate the norm of states, for which we can use `np.linalg.norm`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By passing all these routines to `optimize_pulses`, we get the exact same\n",
    "results as in the original example, except much faster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "15"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter.      J_T    ∫gₐ(t)dt          J       ΔJ_T         ΔJ  secs\n",
      "0     1.00e+00    0.00e+00   1.00e+00        n/a        n/a     0\n",
      "1     7.65e-01    2.33e-02   7.88e-01  -2.35e-01  -2.12e-01     0\n",
      "2     5.56e-01    2.07e-02   5.77e-01  -2.09e-01  -1.88e-01     0\n",
      "3     3.89e-01    1.66e-02   4.05e-01  -1.67e-01  -1.51e-01     0\n",
      "4     2.65e-01    1.23e-02   2.77e-01  -1.24e-01  -1.12e-01     0\n",
      "5     1.78e-01    8.63e-03   1.87e-01  -8.68e-02  -7.82e-02     0\n",
      "6     1.19e-01    5.86e-03   1.25e-01  -5.89e-02  -5.30e-02     0\n",
      "7     8.01e-02    3.91e-03   8.40e-02  -3.92e-02  -3.53e-02     0\n",
      "8     5.42e-02    2.58e-03   5.68e-02  -2.59e-02  -2.33e-02     0\n",
      "9     3.71e-02    1.70e-03   3.88e-02  -1.71e-02  -1.54e-02     0\n",
      "10    2.58e-02    1.13e-03   2.70e-02  -1.13e-02  -1.02e-02     0\n"
     ]
    }
   ],
   "source": [
    "opt_result = krotov.optimize_pulses(\n",
    "    objectives,\n",
    "    pulse_options=pulse_options,\n",
    "    tlist=tlist,\n",
    "    propagator=expm,\n",
    "    chi_constructor=krotov.functionals.chis_re,\n",
    "    info_hook=krotov.info_hooks.print_table(J_T=krotov.functionals.J_T_re),\n",
    "    check_convergence=krotov.convergence.check_monotonic_error,\n",
    "    iter_stop=10,\n",
    "    norm=np.linalg.norm,\n",
    "    mu=mu,\n",
    "    overlap=overlap,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "16"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Krotov Optimization Result\n",
       "--------------------------\n",
       "- Started at 2019-12-15 22:48:47\n",
       "- Number of objectives: 1\n",
       "- Number of iterations: 10\n",
       "- Reason for termination: Reached 10 iterations\n",
       "- Ended at 2019-12-15 22:48:51 (0:00:04)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt_result"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}