{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to Continue Making a Wormhole\n",
    "\n",
    "## Part 3: A Really Existing Quantum Computer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a quick example using Rigetti's quantum computer. You go on their site: [rigetti.com](https://www.rigetti.com/.com) and make an account. The idea is you have a classical virtual server sitting next to the quantum lattices. You can ssh in, put code to it. Meanwhile, you reserve a lattice at a certain time using their web interface. And when your time arrives, your classical server knows about it, and if you call \"get_qc\" in your python code with the lattice name, then your program runs on the quantum hardware. They have pretty good documentation:\n",
    "\n",
    "https://www.rigetti.com/qcs-docs\n",
    "\n",
    "http://docs.rigetti.com/en/stable/\n",
    "\n",
    "https://grove-docs.readthedocs.io/en/latest/\n",
    "\n",
    "When you run things offline, you have to first make sure you do\n",
    "\n",
    "`qvm -S`\n",
    "\n",
    "`quilc -S`\n",
    "\n",
    "which starts the quantum virtual machine and the compiler. The idea is you can specify your quantum program in terms of a circuit of unitary matrices and measurements. Then this is compiled down into local gates, etc, which are optimized for the specific lattice structure you're running on, wherein only certain qubits can interact. Moreover your python defined program is compiled into an intermediate language.\n",
    "\n",
    "So here we just do a simple wormhole. We just have to translate over our qutip operators into forest operators, which require definitions, etc. We use grove's \"create_arbitrary_state\" to get a unitary that takes us from the fixed initial state to the TFD. In the end, we just measure the output qubit, and repeat the experiment many times.\n",
    "\n",
    "Now I'm sure there are much much more elegant ways of translating the wormhole stuff into something really tuned to Rigetti's system, but that's for the future. As it is, if I do any more qubits that 5 total, the compiler/qvm system times out. You'll find a very weak effect. But that's all we can do for now!\n",
    "\n",
    "It's also recommended to use [pytket](https://github.com/CQCL/pytket) which allows your quantum programs to run on many of the currently available quantum computing platforms.\n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import qutip as qt\n",
    "from qutip.qip.operations import swap\n",
    "\n",
    "from pyquil import Program\n",
    "from pyquil.quil import DefGate\n",
    "from pyquil.gates import MEASURE, SWAP\n",
    "from pyquil import Program\n",
    "from pyquil import get_qc\n",
    "import grove\n",
    "import grove.alpha.arbitrary_state.arbitrary_state\n",
    "\n",
    "####################################################################################################\n",
    "\n",
    "def construct_tfd(E, beta=0):\n",
    "    El, Ev = E.eigenstates()\n",
    "    return (1/np.sqrt((-beta*E).expm().tr()))*\\\n",
    "             sum([np.exp(-(1/2)*beta*l)*\\\n",
    "                 qt.tensor(Ev[i].conj(), Ev[i])\\\n",
    "                    for i, l in enumerate(El)])\n",
    "\n",
    "def construct_coupling(n):\n",
    "    return\\\n",
    "    (1/(n-1))*sum([qt.tensor(*[qt.sigmaz()\\\n",
    "                       if j == i else qt.identity(2) \\\n",
    "                            for j in range(2*n)])*\n",
    "                   qt.tensor(*[qt.sigmaz()\\\n",
    "                        if j == i+n else qt.identity(2) \\\n",
    "                            for j in range(2*n)])\n",
    "                                for i in range(n)])\n",
    "\n",
    "def optimize_g(n, t, state, E, V, ID2, IDn):\n",
    "    def prep_worm(n, t, E, ID2, IDn):\n",
    "        U1 = qt.tensor(qt.identity(2), (1j*E*t).expm(), IDn)\n",
    "        U2 = swap(N=2*n+1, targets=[1,2])\n",
    "        U3 = qt.tensor(qt.identity(2), (-1j*E*t).expm(), IDn)\n",
    "        U5 = qt.tensor(qt.identity(2), IDn, (-1j*E.trans()*t).expm())\n",
    "        return [U5, U3*U2*U1]\n",
    "    def worm(n, g, V, state, ops, ID2):\n",
    "        A, B = ops\n",
    "        U4 = qt.tensor(qt.identity(2), (1j*g*V).expm())\n",
    "        state = A*U4*B*state\n",
    "        return qt.expect(qt.sigmaz(), state.ptrace(n+1))\n",
    "    G = np.linspace(0, 15, 60)\n",
    "    worm_prep = prep_worm(n, t, E, ID2, IDn)\n",
    "    Z = [worm(n, g, V, state, worm_prep, ID2) for g in G]\n",
    "    g = G[np.argmin(np.array(Z))]\n",
    "    return g\n",
    "\n",
    "####################################################################################################\n",
    "\n",
    "n = 2\n",
    "N = 2*n + 1\n",
    "beta = 0\n",
    "t = 10\n",
    "\n",
    "ID2 = qt.identity(2)\n",
    "IDn = qt.identity(2**n)\n",
    "IDn.dims = [[2]*n, [2]*n]\n",
    "E = qt.rand_herm(2**n)\n",
    "E.dims = IDn.dims\n",
    "tfd = construct_tfd(E, beta=beta)\n",
    "tfd.dims = [[2]*(2*n), [1]*(2*n)]\n",
    "\n",
    "msg = qt.basis(2,0)\n",
    "state = qt.tensor(msg, tfd)\n",
    "V = construct_coupling(n)\n",
    "g = optimize_g(n, t, state, E, V, ID2, IDn)\n",
    "\n",
    "####################################################################################################\n",
    "\n",
    "CREATE = grove.\\\n",
    "alpha.\\\n",
    "arbitrary_state.\\\n",
    "arbitrary_state.\\\n",
    "create_arbitrary_state(state.full().T[0])\n",
    "\n",
    "LBACK_ = (1j*E*t).expm()\n",
    "LBACK_DEF = DefGate(\"LBACK\", LBACK_.full())\n",
    "LBACK = LBACK_DEF.get_constructor()\n",
    "\n",
    "LFORWARD_ = (-1j*E*t).expm()\n",
    "LFORWARD_DEF = DefGate(\"LFORWARD\", LFORWARD_.full())\n",
    "LFORWARD = LFORWARD_DEF.get_constructor()\n",
    "\n",
    "def coupling(g, V):\n",
    "    V_ = (1j*g*V).expm().full()\n",
    "    V_DEF = DefGate(\"V\", V_)\n",
    "    V = V_DEF.get_constructor()\n",
    "    return (V_DEF, V)\n",
    "\n",
    "RFORWARD_ = (-1j*E.trans()*t).expm()\n",
    "RFORWARD_DEF = DefGate(\"RFORWARD\", RFORWARD_.full())\n",
    "RFORWARD = RFORWARD_DEF.get_constructor()\n",
    "\n",
    "def wormhole(g):\n",
    "    global CREATE, V, LBACK_DEF, LFORWARD_DEF, RFORWARD_DEF,\\\n",
    "            LBACK, LFORWARD, RFORWARD, N, n\n",
    "    COUPLING_DEF, COUPLING = coupling(g, V)\n",
    "    p = Program()\n",
    "    ro = p.declare('ro', 'BIT', 1)\n",
    "    p += LBACK_DEF\n",
    "    p += LFORWARD_DEF\n",
    "    p += COUPLING_DEF\n",
    "    p += RFORWARD_DEF\n",
    "    p += CREATE\n",
    "    p += LBACK(*range(1, n+1))\n",
    "    p += SWAP(0, 1)\n",
    "    p += LFORWARD(*range(1, n+1))\n",
    "    p += COUPLING(*range(1, 2*n+1))\n",
    "    p += RFORWARD(*range(n+1, 2*n+1))\n",
    "    p += MEASURE(n+1, ro[0])\n",
    "    p.wrap_in_numshots_loop(1000)\n",
    "    return p\n",
    "\n",
    "W = wormhole(g)\n",
    "qc = get_qc('%dq-qvm' % N)  # You can make any 'nq-qvm' this way for any reasonable 'n'\n",
    "executable = qc.compile(W)\n",
    "result = qc.run(executable)\n",
    "nonzero = np.count_nonzero(result)\n",
    "R = -1*(nonzero/len(tresult)) + 1*((len(result)-nonzero)/len(result))\n",
    "print(R)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "VPython",
   "language": "python",
   "name": "vpython"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}