2018-04-10_mg_construct_bit_representation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projecting operators to the single-excitation subspace\n",
    "\n",
    "In this notebook, we develop a routine that project and/or construct quantum\n",
    "operators in the\n",
    "single-excitation subspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "1"
    }
   },
   "outputs": [],
   "source": [
    "% load_ext autoreload\n",
    "% autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import time\n",
    "from collections import OrderedDict\n",
    "from src.dicke_half_model_v2 import num_vals as num_vals_func, dicke_guess_controls\n",
    "from src.algebra_v1 import split_hamiltonian\n",
    "import qutip\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T=5;theta=0; E0_cycles=2; mcwf=False; non_herm=False;\n",
    "lambda_a=1.0; J_T_conv=1e-4; iter_stop=5000; nt=None;\n",
    "max_ram_mb=100000; kappa=0.01; seed=None; observables='all';\n",
    "keep_pulses='prev'; config='config'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qobj_extract_single_excitation(qobj):\n",
    "    \"\"\"Given a Qutip quantum object (operator or state), that is defined in the\n",
    "    full N-node chain Hilbert space, return a new quantum object in the\n",
    "    zero/single-excitation subspace.\n",
    "\n",
    "    The input quantum object must have the tensor structure (TLS_q ⊗ TLS_c)^N,\n",
    "    as reflected in the `dims` property.\n",
    "    \"\"\"\n",
    "    assert qobj.dims[0] == [2 for _ in qobj.dims[0]]\n",
    "    SINGLE_EXCITATION = 1\n",
    "    subspace_levels = [\n",
    "        i for (i, nums)  # nums = tuple of quantum numbers in the full basis\n",
    "        in enumerate(itertools.product(*[range(n) for n in qobj.dims[0]]))\n",
    "        if sum(nums) <= SINGLE_EXCITATION]\n",
    "    return qobj.extract_states(subspace_levels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def symbolic_hamiltonians_v1(n_nodes):\n",
    "    from src.single_sided_network_v1 import network_slh as network_slh_func_v1\n",
    "\n",
    "    slh = network_slh_func_v1(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    controls = dicke_guess_controls(\n",
    "        slh=slh, theta=theta, T=T, E0_cycles=E0_cycles, nt=nt, kappa=kappa)\n",
    "    control_syms = list(controls.keys())\n",
    "    ham_parts = split_hamiltonian(slh.H, use_cc=False, controls=control_syms)\n",
    "    return ham_parts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def symbolic_hamiltonians_v2(n_nodes):\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    controls = dicke_guess_controls(\n",
    "        slh=slh, theta=theta, T=T, E0_cycles=E0_cycles, nt=nt, kappa=kappa)\n",
    "    control_syms = list(controls.keys())\n",
    "    ham_parts = split_hamiltonian(slh.H, use_cc=False, controls=control_syms)\n",
    "    return ham_parts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def symbolic_lindblads(n_nodes):\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    return slh.Ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_hamiltonians(n_nodes):\n",
    "\n",
    "    from src.single_sided_network_v1 import network_slh as network_slh_func_v2\n",
    "    # assuming symbolic_hamiltonians_v2 is tested, there's no different between\n",
    "    # a potential \"v1\" and \"v2\"\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    controls = dicke_guess_controls(\n",
    "        slh=slh, theta=theta, T=T, E0_cycles=E0_cycles, nt=nt, kappa=kappa)\n",
    "    H_num = slh.H.substitute(num_vals)\n",
    "    control_syms = list(controls.keys())\n",
    "    ham_parts = split_hamiltonian(H_num, use_cc=False, controls=control_syms)\n",
    "    num_hams = OrderedDict()\n",
    "    for (key, H) in ham_parts.items():\n",
    "        num_hams[key] = H.substitute({sym: 1 for sym in control_syms})\n",
    "    return num_hams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_lindblads(n_nodes):\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    Ls = symbolic_lindblads(n_nodes)\n",
    "    return [L.substitute(num_vals) for L in Ls]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_hamiltonians_v1(n_nodes):\n",
    "    \"\"\"Construct the Hamiltonians by constructing them in the full space, and\n",
    "    then converting them to the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v1 import network_slh as network_slh_func_v1\n",
    "    from qnet.convert.to_qutip import convert_to_qutip\n",
    "\n",
    "    slh = network_slh_func_v1(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    controls = dicke_guess_controls(\n",
    "        slh=slh, theta=theta, T=T, E0_cycles=E0_cycles, nt=nt, kappa=kappa)\n",
    "    H_num = slh.H.substitute(num_vals)\n",
    "    control_syms = list(controls.keys())\n",
    "    ham_parts = split_hamiltonian(H_num, use_cc=False, controls=control_syms)\n",
    "    hs = H_num.space\n",
    "    for (key, H) in ham_parts.items():\n",
    "        H = H.substitute({sym: 1 for sym in control_syms})\n",
    "        ham_parts[key] = qobj_extract_single_excitation(\n",
    "                convert_to_qutip(H, full_space=hs))\n",
    "    return ham_parts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_lindblads_v1(n_nodes):\n",
    "    \"\"\"Construct the Lindblad ops by constructing them in the full space, and\n",
    "    then converting them to the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v1 import network_slh as network_slh_func_v1\n",
    "    from qnet.convert.to_qutip import convert_to_qutip\n",
    "\n",
    "    slh = network_slh_func_v1(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    hs = slh.H.space\n",
    "    res = []\n",
    "    for L in slh.Ls:\n",
    "        L_num = L.substitute(num_vals)\n",
    "        res.append(qobj_extract_single_excitation(\n",
    "            convert_to_qutip(L_num, full_space=hs)))\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_spin_index(n_nodes):\n",
    "    \"\"\"Return the spin_index for the given number of nodes (for\n",
    "    debugging/inspection)\"\"\"\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "    from src.to_single_excitation_qutip import construct_bit_index\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    return construct_bit_index(slh.H.space)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_hamiltonians_v2(n_nodes):\n",
    "    \"\"\"Construct the Hamiltonians directly in the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "    from src.to_single_excitation_qutip import (\n",
    "            convert_to_single_excitation_qutip, construct_bit_index)\n",
    "    from qnet.algebra import pattern, wc, LocalSigma, Create, Destroy\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    controls = dicke_guess_controls(\n",
    "        slh=slh, theta=theta, T=T, E0_cycles=E0_cycles, nt=nt, kappa=kappa)\n",
    "    H_num = slh.H.substitute(num_vals)\n",
    "    control_syms = list(controls.keys())\n",
    "    ham_parts = split_hamiltonian(H_num, use_cc=False, controls=control_syms)\n",
    "    hs = H_num.space\n",
    "    bit_index = construct_bit_index(hs)\n",
    "\n",
    "    pat_localsigma = pattern(LocalSigma, wc('i'), wc('j'), hs=wc('hs'))\n",
    "\n",
    "    def localsigma_to_create_destroy(i, j, hs):\n",
    "        if (i, j) == ('g', 'e'):\n",
    "            return Destroy(hs=hs)\n",
    "        elif (i, j) == ('e', 'g'):\n",
    "            return Create(hs=hs)\n",
    "        else:\n",
    "            raise ValueError(\"Invalid (i, j) = %s\" % (str((i, j))))\n",
    "\n",
    "    for (key, H) in ham_parts.items():\n",
    "        H = H.substitute({sym: 1 for sym in control_syms})\n",
    "        H = H.simplify([(pat_localsigma, localsigma_to_create_destroy)])\n",
    "        ham_parts[key] = convert_to_single_excitation_qutip(\n",
    "            H, bit_index, full_space=hs)\n",
    "    return ham_parts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_lindblads_v2(n_nodes):\n",
    "    \"\"\"Construct the Lindblad ops by constructing them in the full space, and\n",
    "    then converting them to the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "    from src.to_single_excitation_qutip import (\n",
    "            convert_to_single_excitation_qutip, construct_bit_index)\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    num_vals = num_vals_func(theta=theta, n_nodes=n_nodes, kappa=kappa)\n",
    "    hs = slh.H.space\n",
    "    bit_index = construct_bit_index(hs)\n",
    "    res = []\n",
    "    for L in slh.Ls:\n",
    "        L_num = L.substitute(num_vals)\n",
    "        res.append(convert_to_single_excitation_qutip(\n",
    "            L_num, bit_index, full_space=hs))\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_states_v1(n_nodes):\n",
    "    \"\"\"Construct the states in the full space, and\n",
    "    then converting them to the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v1 import network_slh as network_slh_func_v1\n",
    "    from src.qdyn_model_v2 import dicke_state, logical_2q_state\n",
    "\n",
    "    slh = network_slh_func_v1(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    hs = slh.H.space\n",
    "    dicke1 = qobj_extract_single_excitation(dicke_state(hs, excitations=1))\n",
    "    psi10 = qobj_extract_single_excitation(logical_2q_state(hs, 1, 0))\n",
    "    return {'dicke1': dicke1, 'psi10': psi10}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _state_to_fmt(qutip_state, fmt):\n",
    "    import QDYN.linalg\n",
    "    if fmt == 'qutip':\n",
    "        return qutip_state\n",
    "    elif fmt == 'numpy':\n",
    "        return QDYN.linalg.vectorize(qutip_state.data.todense())\n",
    "    else:\n",
    "        raise ValueError(\"Unknown fmt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def single_excitation_dicke_state(bit_index, fmt='qutip'):\n",
    "    \"\"\"Return the zero/single-excitation dicke state for the given system,\n",
    "    directly using the single-exctiation subspace encoding. Assumes\n",
    "    that the Hilbert spaces for the qubits at the different nodes are labeled\n",
    "    'q1', 'q2', ...\n",
    "\n",
    "    For example for a Hilbert space ``('q1', 'c1', 'q2', 'c2', 'q3', 'c3')``,\n",
    "    the Dicke state is ``(|000010> + |001000> + |100000>) / sqrt(3)``. The\n",
    "    corresponding bit-index is\n",
    "\n",
    "        {'000000': 0,\n",
    "         '000001': 1,\n",
    "         '000010': 2,\n",
    "         '000100': 3,\n",
    "         '001000': 4,\n",
    "         '010000': 5,\n",
    "         '100000': 6}\n",
    "\n",
    "    and only the eigenstates 2, 4, and 6 contribute to the Dicke state.\n",
    "    \"\"\"\n",
    "    N = len(bit_index)\n",
    "    terms = [qutip.basis(N, i) for i in range(2, N, 2)]\n",
    "    res = sum(terms) / np.sqrt(len(terms))\n",
    "    return _state_to_fmt(res, fmt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def single_excitation_initial_state(bit_index, fmt='qutip'):\n",
    "    \"\"\"Return the initial state \"e0....0\" directly using the single-excitation\n",
    "    subspace enconding\"\"\"\n",
    "    N = len(bit_index)\n",
    "    i = bit_index[\"1\" + \"0\" * (N-2)]\n",
    "    assert i == N - 1\n",
    "    return _state_to_fmt(qutip.basis(N, i), fmt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qutip_states_v2(n_nodes):\n",
    "    \"\"\"Construct the states in the full space, and\n",
    "    then converting them to the single-excitation subspace\"\"\"\n",
    "\n",
    "    from src.single_sided_network_v2 import network_slh as network_slh_func_v2\n",
    "    from src.to_single_excitation_qutip import construct_bit_index\n",
    "\n",
    "    slh = network_slh_func_v2(n_cavity=2, n_nodes=n_nodes, topology='open')\n",
    "    hs = slh.H.space\n",
    "    bit_index = construct_bit_index(hs)\n",
    "    dicke1 = single_excitation_dicke_state(bit_index)\n",
    "    psi10 = single_excitation_initial_state(bit_index)\n",
    "    return {'dicke1': dicke1, 'psi10': psi10}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_symbolic_hamiltonians(n_nodes):\n",
    "    t1 = time.time()\n",
    "    hams_v1 = symbolic_hamiltonians_v1(n_nodes)\n",
    "    t2 = time.time()\n",
    "    hams_v2 = symbolic_hamiltonians_v2(n_nodes)\n",
    "    t3 = time.time()\n",
    "    print(\"v1: %.2f seconds\" % (t2-t1))\n",
    "    print(\"v2: %.2f seconds\" % (t3-t2))\n",
    "    assert hams_v1.keys() == hams_v2.keys()\n",
    "    for key in hams_v1:\n",
    "        assert hams_v1[key] == hams_v2[key]\n",
    "    print(\"OK\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qutip_hamiltonians(n_nodes):\n",
    "    t1 = time.time()\n",
    "    hams_v1 = qutip_hamiltonians_v1(n_nodes)\n",
    "    t2 = time.time()\n",
    "    hams_v2 = qutip_hamiltonians_v2(n_nodes)\n",
    "    t3 = time.time()\n",
    "    print(\"v1: %.2f seconds\" % (t2-t1))\n",
    "    print(\"v2: %.2f seconds\" % (t3-t2))\n",
    "    assert hams_v1.keys() == hams_v2.keys()\n",
    "    for key in hams_v1:\n",
    "        assert hams_v1[key] == hams_v2[key], \"Ham %s does not match\" % key\n",
    "    print(\"OK\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qutip_states(n_nodes):\n",
    "    t1 = time.time()\n",
    "    states_v1 = qutip_states_v1(n_nodes)\n",
    "    t2 = time.time()\n",
    "    states_v2 = qutip_states_v2(n_nodes)\n",
    "    t3 = time.time()\n",
    "    print(\"v1: %.2f seconds\" % (t2-t1))\n",
    "    print(\"v2: %.2f seconds\" % (t3-t2))\n",
    "    assert states_v1.keys() == states_v2.keys()\n",
    "    for key in states_v1:\n",
    "        assert states_v1[key] == states_v2[key], \\\n",
    "            \"state %s does not match\" % key\n",
    "    print(\"OK\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qutip_lindblads(n_nodes):\n",
    "    t1 = time.time()\n",
    "    Ls_v1 = qutip_lindblads_v1(n_nodes)\n",
    "    t2 = time.time()\n",
    "    Ls_v2 = qutip_lindblads_v2(n_nodes)\n",
    "    t3 = time.time()\n",
    "    print(\"v1: %.2f seconds\" % (t2-t1))\n",
    "    print(\"v2: %.2f seconds\" % (t3-t2))\n",
    "    for i, L in enumerate(Ls_v1):\n",
    "        assert L.norm('max') > 0.14\n",
    "        assert L == Ls_v2[i], \"L does not match\"\n",
    "    print(\"OK\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_symbolic_hamiltonians(4)\n",
    "test_qutip_hamiltonians(4)\n",
    "test_qutip_states(4)\n",
    "test_qutip_lindblads(4)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}