2020-09-01_quantum-sys-ID-interactive_01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Interactive Quantum System Identification Applet\n",
    "## 1.1 Background\n",
    "The following notebook demonstrates the current version of quantum system identification I've been working on. This notebook demonstrates the algorithm's ability to identify a random quantum system from simulated data.\n",
    "- Select the order of the system\n",
    "- Generate an appropriate basis\n",
    "- Generate a random quantum state transition matrix\n",
    "- Simulate data\n",
    "- Add noise\n",
    "- Regress data onto basis as $\\dot{x}_{data}=\\Theta w$\n",
    "- Estimate original matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Import Functions\n",
    "Mostly ipywidgets, bqplot, and numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "\n",
    "from bqplot import pyplot as plt\n",
    "import random\n",
    "from IPython.display import display\n",
    "from ipywidgets.widgets import IntSlider, HBox, Button, interact, VBox\n",
    "import ipywidgets.widgets as widgets\n",
    "\n",
    "from scipy import signal\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Select Sytem Order\n",
    "Should work for any order, as long as you have RAM for the basis. Mostly tested with 2 through 5.\n",
    "## n="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6f23a2fa318476aa0fb5337cd2c0bc9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "IntText(value=3)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n=widgets.IntText(\n",
    "    value=3,\n",
    "    disabled=False\n",
    ")\n",
    "display(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Generate a Basis for the Space\n",
    "Not pauli like, but spans the space.\n",
    "\n",
    "First, generates matricies with a single element in one of the diagonals.\n",
    "\n",
    "Then, generate diagonal element for the real component, ie (1,3 & 3,1) with the same sign, for Hamiltonian structure.\n",
    "\n",
    "Finally, generate diagonal element pairs for imaginary component with opposite sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis_button = widgets.Button(description=\"Generate Basis\")\n",
    "output = widgets.Output()\n",
    "#display(basis_button, output)\n",
    "def on_basis_button_clicked(b):\n",
    "    #diagonals\n",
    "    output.clear_output()\n",
    "    diag=np.zeros((n.value,n.value,n.value),dtype=complex);\n",
    "    for i in range(n.value):\n",
    "        diag[i,i,i]=0+1j\n",
    "        \n",
    "    idx_start=1;\n",
    "    counter=-1;\n",
    "    count=int(((n.value*n.value)-n.value)/2)\n",
    "    off_complex=np.zeros((count,n.value,n.value),dtype=complex);\n",
    "    for i in range(0,n.value-1):\n",
    "        for ii in range(idx_start,n.value):\n",
    "            if i != ii and i<ii:\n",
    "                off_complex[counter,i,ii]=0+1j\n",
    "                off_complex[counter,ii,i]=0+1j\n",
    "                counter=counter+1\n",
    "                \n",
    "    idx_start=1;\n",
    "    counter=-1;\n",
    "    count=int(((n.value*n.value)-n.value)/2)\n",
    "    off_real=np.zeros((count,n.value,n.value),dtype=complex);\n",
    "    for i in range(0,n.value-1):\n",
    "        for ii in range(idx_start,n.value):\n",
    "            if i != ii and i<ii:\n",
    "                off_real[counter,i,ii]=1+0j\n",
    "                off_real[counter,ii,i]=-1+0j\n",
    "                counter=counter+1\n",
    "\n",
    "    basis=np.concatenate((diag,off_complex,off_real));\n",
    "    basis_button.value=basis\n",
    "    with output:\n",
    "        display(basis)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "basis_button.on_click(on_basis_button_clicked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Generate a Random Matrix\n",
    "## 5.1 Matrix Options\n",
    "First, select from a set of options. Default will generate a matrix with random entries in every possible element of the state transition matrix. Off diagonal only leaves zeros on the diagonal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be9ffd15e109428db130eb90f1fbaa69",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "RadioButtons(description='Settings:', options=('default', 'Off Diagonal Matrix Only?', 'Upper Diagonal Only?',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "option=widgets.RadioButtons(\n",
    "    options=['default','Off Diagonal Matrix Only?', 'Upper Diagonal Only?', 'Lower Diagonal Only?'],\n",
    "#     value='pineapple',\n",
    "    description='Settings:',\n",
    "    disabled=False\n",
    ")\n",
    "display(option)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Generate Random State Matrix\n",
    "This button can be clicked repeatedly. Adds a random component of each basis to form a state matrix for identification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_button = widgets.Button(description=\"Generate Matrix\")\n",
    "output = widgets.Output()\n",
    "#display(matrix_button, output)\n",
    "def on_matrix_button_clicked(b):\n",
    "    output.clear_output()\n",
    "    \n",
    "    if option.value == 'default':\n",
    "        A=np.zeros((n.value,n.value),dtype=complex)\n",
    "        for i in range(len(basis_button.value)):\n",
    "            A=A+np.squeeze((basis_button.value[i]*random.random()))\n",
    "            \n",
    "    if option.value == 'Upper Diagonal Only?':\n",
    "        A=np.zeros((n.value,n.value),dtype=complex)\n",
    "        for i in range(0,1):\n",
    "            A=A+np.squeeze((basis_button.value[i]*random.random()))\n",
    "            \n",
    "    if option.value == 'Lower Diagonal Only?':\n",
    "        A=np.zeros((n.value,n.value),dtype=complex)\n",
    "        A=A+np.squeeze((basis_button.value[n.value-1]*random.random()))\n",
    "        \n",
    "    if option.value == 'Off Diagonal Matrix Only?':\n",
    "        A=np.zeros((n.value,n.value),dtype=complex)\n",
    "        for i in range((len(basis_button.value)-n.value),len(basis_button.value)):\n",
    "            A=A+np.squeeze((basis_button.value[i]*random.random()))\n",
    "            \n",
    "    matrix_button.value=A\n",
    "    \n",
    "    \n",
    "\n",
    "    with output:\n",
    "        display(A)\n",
    "\n",
    "\n",
    "\n",
    "matrix_button.on_click(on_matrix_button_clicked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Generate Data\n",
    "Form the system in state space terms, and use ODE45 to generate 1000 points of data over 10 seconds. Plot the timeseries data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix_button.value=np.array([[0-1j, 0],[0,0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "simdata_button = widgets.Button(description=\"Simulate Data\")\n",
    "\n",
    "t = np.linspace(0, 10, 1000, endpoint=False)\n",
    "output = widgets.Output()\n",
    "#display(simdata_button, output)\n",
    "def on_simdata_button_clicked(b):\n",
    "    output.clear_output()\n",
    "    Xinit=np.ones(n.value)-(1j*np.ones(n.value))\n",
    "\n",
    "    sys = signal.StateSpace(matrix_button.value, np.zeros((n.value,1)), np.zeros((1,n.value)),0)\n",
    "    width = '500px'\n",
    "\n",
    "    tout, yout, xout = signal.lsim(sys,U=0,T=t,X0=Xinit)\n",
    "    plt.clear()\n",
    "    \n",
    "    fig_circle = plt.figure(title=\"Sim Real\")\n",
    "    plt.plot(x=tout, y=np.real(np.transpose(xout)))\n",
    "    fig_circle.layout.height = width\n",
    "    fig_circle.layout.width = width\n",
    "    \n",
    "    fig_ratios = plt.figure(title=\"Sim Imag\")\n",
    "    plt.plot(x=tout, y=np.imag(np.transpose(xout)))\n",
    "    fig_ratios.layout.height = width\n",
    "    fig_ratios.layout.width = width\n",
    "\n",
    "    simdata_button.value=xout\n",
    "    box = HBox([fig_circle, fig_ratios])\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "    with output:\n",
    "        display(box)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "simdata_button.on_click(on_simdata_button_clicked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.5 Add Noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adds random, normally distributed noise to the signal as $x_{noised}=x_{clean}+\\epsilon*random(-1,1)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b6aaa4e4f8134a639e526127d2681f69",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatSlider(value=0.1, continuous_update=False, description='Noise Coef:', max=1.0, readout_format='.1f', step…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eps=widgets.FloatSlider(\n",
    "    value=0.1,\n",
    "    min=0,\n",
    "    max=1,\n",
    "    step=0.05,\n",
    "    description='Noise Coef:',\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.1f',\n",
    ")\n",
    "display(eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "addnoise_button = widgets.Button(description=\"Add Noise\")\n",
    "\n",
    "t = np.linspace(0, 10, 1000, endpoint=False)\n",
    "output = widgets.Output()\n",
    "def on_addnoise_button_clicked(b):\n",
    "    output.clear_output()\n",
    "    x_clean=simdata_button.value\n",
    "    x_noise=x_clean+(eps.value*np.random.normal(-1,1,[len(x_clean),n.value]))+(1j*eps.value*np.random.normal(-1,1,[len(x_clean),n.value]))\n",
    "\n",
    "    plt.clear()\n",
    "    width = '500px'\n",
    "    \n",
    "    fig = plt.figure(title=\"Real Noise\")\n",
    "    fig_circle = plt.figure()\n",
    "    plt.plot(x=t, y=np.real(np.transpose(x_noise)))\n",
    "    fig_circle.layout.height = width\n",
    "    fig_circle.layout.width = width\n",
    "    \n",
    "    fig_ratios = plt.figure(title=\"Imag Noise\")\n",
    "    plt.plot(x=t, y=np.imag(np.transpose(x_noise)))\n",
    "    fig_ratios.layout.height = width\n",
    "    fig_ratios.layout.width = width\n",
    "\n",
    "    simdata_button.value=x_noise\n",
    "    box = HBox([fig_circle, fig_ratios])\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "    with output:\n",
    "        display(box)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "addnoise_button.on_click(on_addnoise_button_clicked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Regress Data Onto Basis Library Using Only Output Data\n",
    "In SINDY like fashion, pass data through each basis function to generate $\\Theta$. Regress onto the numerical derivative of the data, $\\dot{x}$, as $\\dot{x}=\\Theta w$. Result, $w$ is the coefficient in front of each basis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "regress_button = widgets.Button(description=\"Estimate State Matrix\")\n",
    "output = widgets.Output()\n",
    "#display(regress_button, output)\n",
    "def on_regress_button_clicked(b):\n",
    "    output.clear_output()\n",
    "    Theta_cust=np.empty((len(simdata_button.value), 0))\n",
    "    for i in range(len(basis_button.value)):\n",
    "        x_sigma=np.matmul(simdata_button.value,basis_button.value[i])\n",
    "        Theta_cust=np.append(Theta_cust,x_sigma,axis=1)\n",
    "    Theta_cust=np.reshape(Theta_cust,[-1,len(basis_button.value)],'F')\n",
    "    Theta=np.concatenate((np.real(Theta_cust),np.imag(Theta_cust)),axis=0) \n",
    "\n",
    "    dx=np.gradient(simdata_button.value,0.01,axis=0)\n",
    "    dx2=np.squeeze(np.concatenate((np.real(dx.flatten('F')),np.imag(dx.flatten('F')))))\n",
    "    Xi=np.linalg.lstsq(Theta,dx2,rcond=None)\n",
    "    \n",
    "    lambda_knob=0.1\n",
    "    smallinds=np.argwhere(abs(Xi[0])<lambda_knob)\n",
    "    Theta_small=np.delete(Theta,smallinds,axis=1);\n",
    "    basis_small=np.delete(basis_button.value,smallinds,axis=0)\n",
    "    Xi=np.linalg.lstsq(Theta_small,dx2,rcond=None)\n",
    "    A_solve=np.zeros([n.value,n.value])\n",
    "    for i in range(len(basis_small)):\n",
    "        if np.all(np.isreal(basis_small[i])):\n",
    "            Xi[0][i]=-Xi[0][i]\n",
    "        A_solve=A_solve+(Xi[0][i]*basis_small[i])\n",
    "        \n",
    "    regress_button.value=A_solve\n",
    "        \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "    with output:\n",
    "        display(A_solve)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "regress_button.on_click(on_regress_button_clicked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Xi=regress_button.value; \n",
    "# basis=basis_button.value; basis\n",
    "# for i in range(len(basis)):\n",
    "#     if np.all(np.isreal(basis[i])):\n",
    "#         Xi[0][i]=-Xi[0][i]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Compare Results & Demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b7b3242515bd4f33aa78b4ee02937053",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Generate Basis', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f5a8ef3af3d4de7817aa7796145eab0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Generate Matrix', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "06b79842c8a64519af306e323a8365c4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Simulate Data', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bc3401c07cbe49a79c23e948f688cabb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Add Noise', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9e96a484bc764a3385c92bdd675b57f3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Estimate State Matrix', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87a89edb9b9f4b4da47daa37bf886489",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Compare Results', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56a439af76cc47bd9c1f144c175f04dc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "results_button = widgets.Button(description=\"Compare Results\")\n",
    "output = widgets.Output()\n",
    "\n",
    "display(basis_button,matrix_button,simdata_button,addnoise_button,regress_button,results_button, output)\n",
    "# display(matrix_button.value)\n",
    "# display(regress_button.value)\n",
    "\n",
    "\n",
    "      \n",
    "def on_results_button_clicked(b):\n",
    "    output.clear_output()\n",
    "    A_solve=regress_button.value\n",
    "    Xinit=np.ones(n.value)-(1j*np.ones(n.value))\n",
    "    norm_error=np.linalg.norm((matrix_button.value-regress_button.value))\n",
    "    sys_org = signal.StateSpace(matrix_button.value, np.zeros((n.value,1)), np.zeros((1,n.value)),0)\n",
    "    sys_id = signal.StateSpace(regress_button.value, np.zeros((n.value,1)), np.zeros((1,n.value)),0)\n",
    "    tout, yout, xout = signal.lsim(sys_org,U=0,T=t,X0=Xinit)\n",
    "    tout_id, yout_id, xout_id = signal.lsim(sys_id,U=0,T=t,X0=Xinit)\n",
    "    plt.clear()\n",
    "    result=[]\n",
    "    result_sim_real=np.real(np.transpose(xout))\n",
    "    result_sim_imag=np.imag(np.transpose(xout))\n",
    "    result_id_real=np.real(np.transpose(xout_id))\n",
    "    result_id_imag=np.imag(np.transpose(xout_id))  \n",
    "    \n",
    "    result_sim_real_noise=np.real(np.transpose(simdata_button.value))\n",
    "    result_sim_imag_noise=np.imag(np.transpose(simdata_button.value))\n",
    "    \n",
    "    width = '500px'\n",
    "    \n",
    "    \n",
    "    fig_circle = plt.figure(title=\"Real Comp.\")\n",
    "    plt.plot(tout, result_id_real, 'ks--') \n",
    "    plt.plot(tout, result_sim_real)\n",
    "    fig_circle.layout.height = width\n",
    "    fig_circle.layout.width = width\n",
    "    #plt.plot(x=tout, y=result_) \n",
    "    \n",
    "    fig_ratios = plt.figure(title=\"Imag. Comp\")\n",
    "    plt.plot(tout, result_id_imag, 'ks--') \n",
    "    plt.plot(tout, result_sim_imag)\n",
    "    box = HBox([fig_circle, fig_ratios])\n",
    "    fig_ratios.layout.height = width\n",
    "    fig_ratios.layout.width = width\n",
    "\n",
    "    \n",
    "    fig_noise_real=plt.figure(title='Noise Real')\n",
    "    plt.plot(tout, result_id_real, 'ks--') \n",
    "    plt.plot(tout, result_sim_real_noise)\n",
    "    fig_noise_real.layout.height = width\n",
    "    fig_noise_real.layout.width = width\n",
    "    \n",
    "    fig_noise_imag=plt.figure(title='Noise Imag')\n",
    "    plt.plot(tout, result_id_imag, 'ks--') \n",
    "    plt.plot(tout, result_sim_imag_noise)\n",
    "    fig_noise_imag.layout.height = width\n",
    "    fig_noise_imag.layout.width = width\n",
    "\n",
    "    box2=HBox([fig_noise_real, fig_noise_imag])\n",
    "    box_disp=VBox([box,box2])\n",
    "    \n",
    "    with output:\n",
    "        display(box_disp,A_solve,matrix_button.value,['Norm Error= '+ str(norm_error)])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "results_button.on_click(on_results_button_clicked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}