Workshop Part 4 - Time Evolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals of this section\n",
    "1. Even *more* physics! Time evolution and \"qunatum quenches\" as a way to...\n",
    "2. Learn more about Julia's parallel computing features like `pmap` and module importation\n",
    "3. Exploit various linear algebra tricks to make things faster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far all we've looked at have been static problems. Now we will do something *time dependent*, as a motivator to explore more of Julia's linear algebra and parallel features.\n",
    "\n",
    "We will start in a situation with non-zero transverse magnetic field and turn it off, seeing what happens. This is going to require doing *a lot* of independent matrix-vector operations - a good use case for parallelism!\n",
    "\n",
    "First, we'll make our Hamiltonian from the previous part, with a little time-dependent spice:\n",
    "\n",
    "$$ \\hat{H}(t) = -\\sum_{\\langle i, j \\rangle} \\hat{\\sigma}_i^z \\hat{\\sigma}_j^z - h(t)\\sum_i \\hat{\\sigma}_i^x $$\n",
    "\n",
    "Now $h(t)$ is some time dependent function. This means that the lowest energy state (the groundstate) will change with time as well. We need a way to simulate this. In a closed quantum system, a wavefunction $|\\Psi\\rangle$ will undergo *unitary time evolution* so that:\n",
    "\n",
    "$$ \\left| \\Psi(t) \\right\\rangle = \\hat{U}(t)\\left| \\Psi(t = 0)\\right\\rangle $$\n",
    "\n",
    "so that\n",
    "\n",
    "$$ \\left| \\left\\langle \\Psi(t) \\left|\\right. \\Psi(t) \\right\\rangle\\right|_2 = 1 \\forall t $$\n",
    "\n",
    "$\\hat{U}(t)$ is a *unitary* operator - it's norm preserving. The particular form it takes is:\n",
    "\n",
    "$$ \\hat{U}(t) = \\exp\\left\\{- i t\\hat{H}(t) \\right\\} $$\n",
    "\n",
    "Since $\\hat{H}(t)$ is always Hermitian, $\\hat{U}(t)$ is always unitary. For an explanation of where all this comes from you can consult a textbook on quantum mechanics. For now, if it's confusing, we're just going to\n",
    "  1. Calculate $\\hat{U}(t)$ for various times (and perhaps use a few shortcuts)\n",
    "  2. Multiply it by $|\\Psi\\rangle$ to find the groundstate at various times\n",
    "  3. Make pretty pictures, learn some things\n",
    "  \n",
    "Remember that $|\\Psi\\rangle$ is \"just\" some vector, and $\\hat{U}$ and $\\hat{H}$ are matrices - underneath all the jargon, we're still just doing linear algebra! We'll be able to reuse our types from the previous part. We can do this by hoovering up all our code into a file called `timeevolution.jl` (or whatever you feel like calling it, it's a free country). If you're not sure how to do this quickly, take a look at [`nbconvert`](https://nbconvert.readthedocs.io/en/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_groundstate (generic function with 2 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "include(\"timeevolve.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write some single-node code to generate the wavefunction at various times, which we'll parallelize in a moment. It's good to have a working single-node version as a proof of concept and something to test against. I wrote a cheap helper function `get_groundstate` to automate finding the lowest energy eigenvector for the `TransverseFieldIsing` Hamiltonian. The first argument to my function is the length $L$, and the second argument is the strength of the transverse field $h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "timeSeries (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function timeEvolve(ψ, H, t)\n",
    "    U = expm(-im*t*Hermitian(H)) # want to use the optimized method\n",
    "    return U*ψ\n",
    "end\n",
    "\n",
    "function timeSeries(ψ, H, start, stop, n)\n",
    "    times = linspace(start, stop, n)\n",
    "    map((t,)->timeEvolve(t, ψ, H), times)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Here's your parallelism\n",
    "\n",
    "Now we can use Julia's [parallel map](https://docs.julialang.org/en/latest/stdlib/parallel.html#Base.Distributed.pmap) function to make this faster! For now, all we have to do is add the \"`p`\": "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mOn worker 2:\n\u001b[91mUndefVarError: #timeEvolve not defined\u001b[39m\ndeserialize_datatype at ./serialize.jl:968\nhandle_deserialize at ./serialize.jl:674\ndeserialize at ./serialize.jl:634\nhandle_deserialize at ./serialize.jl:681\ndeserialize_global_from_main at ./distributed/clusterserialize.jl:154\nforeach at ./abstractarray.jl:1724\ndeserialize at ./distributed/clusterserialize.jl:56\nhandle_deserialize at ./serialize.jl:722\ndeserialize at ./serialize.jl:634\ndeserialize_datatype at ./serialize.jl:966\nhandle_deserialize at ./serialize.jl:676\ndeserialize at ./serialize.jl:634\nhandle_deserialize at ./serialize.jl:681\ndeserialize_msg at ./distributed/messages.jl:98\nmessage_handler_loop at ./distributed/process_messages.jl:161\nprocess_tcp_streams at ./distributed/process_messages.jl:118\n#99 at ./event.jl:73\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mOn worker 2:\n\u001b[91mUndefVarError: #timeEvolve not defined\u001b[39m\ndeserialize_datatype at ./serialize.jl:968\nhandle_deserialize at ./serialize.jl:674\ndeserialize at ./serialize.jl:634\nhandle_deserialize at ./serialize.jl:681\ndeserialize_global_from_main at ./distributed/clusterserialize.jl:154\nforeach at ./abstractarray.jl:1724\ndeserialize at ./distributed/clusterserialize.jl:56\nhandle_deserialize at ./serialize.jl:722\ndeserialize at ./serialize.jl:634\ndeserialize_datatype at ./serialize.jl:966\nhandle_deserialize at ./serialize.jl:676\ndeserialize at ./serialize.jl:634\nhandle_deserialize at ./serialize.jl:681\ndeserialize_msg at ./distributed/messages.jl:98\nmessage_handler_loop at ./distributed/process_messages.jl:161\nprocess_tcp_streams at ./distributed/process_messages.jl:118\n#99 at ./event.jl:73\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1m#568\u001b[22m\u001b[22m at \u001b[1m./asyncmap.jl:178\u001b[22m\u001b[22m [inlined]",
      " [2] \u001b[1mforeach\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Base.##568#570, ::Array{Any,1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./abstractarray.jl:1724\u001b[22m\u001b[22m",
      " [3] \u001b[1mmaptwice\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Function, ::Channel{Any}, ::Array{Any,1}, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Vararg{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},N} where N\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./asyncmap.jl:178\u001b[22m\u001b[22m",
      " [4] \u001b[1mwrap_n_exec_twice\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Channel{Any}, ::Array{Any,1}, ::Base.Distributed.##202#205{WorkerPool}, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Vararg{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},N} where N\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./asyncmap.jl:154\u001b[22m\u001b[22m",
      " [5] \u001b[1m#async_usemap#553\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Function, ::Void, ::Function, ::Base.Distributed.##186#188, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Vararg{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},N} where N\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./asyncmap.jl:103\u001b[22m\u001b[22m",
      " [6] \u001b[1m(::Base.#kw##async_usemap)\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Any,1}, ::Base.#async_usemap, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Vararg{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},N} where N\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./<missing>:0\u001b[22m\u001b[22m",
      " [7] \u001b[1m(::Base.#kw##asyncmap)\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Any,1}, ::Base.#asyncmap, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./<missing>:0\u001b[22m\u001b[22m",
      " [8] \u001b[1m#pmap#201\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Bool, ::Int64, ::Void, ::Array{Any,1}, ::Void, ::Function, ::WorkerPool, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./distributed/pmap.jl:126\u001b[22m\u001b[22m",
      " [9] \u001b[1mpmap\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::WorkerPool, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./distributed/pmap.jl:101\u001b[22m\u001b[22m",
      " [10] \u001b[1m#pmap#211\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Any,1}, ::Function, ::Function, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./distributed/pmap.jl:156\u001b[22m\u001b[22m",
      " [11] \u001b[1mtimeSeriesParallel\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Float64,1}, ::TransverseFieldIsing{Float64,Array{T,2} where T}, ::Float64, ::Float64, ::Int64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./In[3]:3\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "function timeSeriesParallel(ψ, H, start, stop, n)\n",
    "    times = linspace(start, stop, n)\n",
    "    pmap((t,) -> timeEvolve(ψ, H.Mat, t), times)\n",
    "end\n",
    "\n",
    "addprocs(6)\n",
    "ψ, H = get_groundstate(10, 1.0)\n",
    "ψ_quench, H_quench = get_groundstate(10, 0.0)\n",
    "timeSeriesParallel(ψ, H_quench, 0.0, 10.0, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oops! The workers don't know about our type and other functions. We'll need to load them onto each worker to be able to use `pmap`. For a more detailed discussion of why this is, consult the [docs](https://docs.julialang.org/en/latest/manual/parallel-computing.html#Code-Availability-and-Loading-Packages-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Array{Complex{Float64},1},1}:\n",
       " Complex{Float64}[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im  …  0.149587+0.0im, 0.1981+0.0im, 0.253379+0.0im, 0.194948+0.0im, 0.163508+0.0im, 0.206921+0.0im, 0.263587+0.0im, 0.216095+0.0im, 0.278388+0.0im, 0.372904+0.0im]                                                               \n",
       " Complex{Float64}[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im  …  0.0424323-0.143443im, 0.0686685+0.185818im, -0.212603-0.137843im, 0.067576+0.182862im, 0.0463809-0.156791im, 0.0717262+0.194092im, -0.221168-0.143397im, 0.0749063+0.202698im, -0.233588-0.151449im, 0.372083-0.0247317im]  \n",
       " Complex{Float64}[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im  …  -0.125514-0.0813787im, -0.150494+0.128822im, 0.103399+0.231321im, -0.1481+0.126773im, -0.137195-0.0889516im, -0.157196+0.134558im, 0.107565+0.24064im, -0.164165+0.140524im, 0.113605+0.254153im, 0.369623-0.0493545im]     \n",
       " Complex{Float64}[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im  …  -0.11364+0.0972748im, -0.173002-0.0965093im, 0.039084-0.250346im, -0.170249-0.0949739im, -0.124215+0.106327im, -0.180705-0.100807im, 0.0406586-0.260432im, -0.188717-0.105276im, 0.0429418-0.275056im, 0.365536-0.0737599im]\n",
       " Complex{Float64}[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im  …  0.0610439+0.136565im, 0.0305572-0.195729im, -0.168988+0.188796im, 0.0300711-0.192615im, 0.0667245+0.149274im, 0.0319179-0.204445im, -0.175796+0.196402im, 0.033333-0.213509im, -0.185668+0.207431im, 0.359839-0.0978405im]  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "addprocs(6)\n",
    "@everywhere include(\"timeevolve.jl\")\n",
    "\n",
    "@everywhere function timeEvolve(ψ, H, t)\n",
    "    U = expm(UniformScaling(-im*t)*Hermitian(H)) # want to use the optimized method\n",
    "    return U*ψ\n",
    "end\n",
    "ψ, H = get_groundstate(8, 1.0)\n",
    "ψ_quench, H_quench = get_groundstate(8, 0.0)\n",
    "\n",
    "#now we can try running with the workers we added\n",
    "ψs_t = timeSeriesParallel(ψ, H_quench, 0.0, 10.0, )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an *embarrassingly parallel* problem. How convenient that the workers don't need to coordinate with each other, only with the driver node. \n",
    "\n",
    "This might, depending on your choice of $L$ and how many time samples you want, take *ages*. $L = 6$ and 5 time samples is an OK choice to make sure that your code works at all. Now we can make an initial plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mMethodError: no method matching magnetization(::Array{Complex{Float64},1})\u001b[0m\nClosest candidates are:\n  magnetization(::Array{T,1} where T, \u001b[91m::Any\u001b[39m) at /Users/kshyatt/Projects/juliacon2017/timeevolve.jl:26\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mMethodError: no method matching magnetization(::Array{Complex{Float64},1})\u001b[0m\nClosest candidates are:\n  magnetization(::Array{T,1} where T, \u001b[91m::Any\u001b[39m) at /Users/kshyatt/Projects/juliacon2017/timeevolve.jl:26\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1m_collect\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Array{Complex{Float64},1},1}, ::Base.Generator{Array{Array{Complex{Float64},1},1},#magnetization}, ::Base.EltypeUnknown, ::Base.HasShape\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./array.jl:442\u001b[22m\u001b[22m",
      " [2] \u001b[1mmap\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Function, ::Array{Array{Complex{Float64},1},1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./abstractarray.jl:1858\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "using GR\n",
    "\n",
    "mags_t = map(magnetization, ψs_t)\n",
    "plot(linspace(0.0, 10.0, 5), mags_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are clearly a lot of inefficiencies in this code. Let's enumerate some of them:\n",
    "  1. We are constructing each Hamiltonian matrix on the head node and sending the whole thing to the workers. All the workers need is the value of $h(t)$.\n",
    "  2. We are sending the entire groundstate back when, for now, all we need is the magnetization (a single `float`).\n",
    "  3. The Hamiltonian is actually diagonal in the basis we have picked. We can use a linear algebra trick to speed up the hard work of computing $\\hat{U}(t)$.\n",
    "\n",
    "$$ \\exp\\left\\{ \\hat{A} \\right\\} = \\hat{P}^\\dagger \\exp\\left\\{ \\hat{D} \\right\\} \\hat{P} $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\hat{A} = \\hat{P}^\\dagger \\hat{D} \\hat{P} $$\n",
    "\n",
    "and $\\hat{D}$ is a diagonal matrix whose entries are the eigenvalues of $\\hat{A}$, and $\\hat{P}$ is the matrix which diagonalizes $\\hat{A}$ (its columns are the eigenvectors of $\\hat{A}$). So:\n",
    "\n",
    "$$ \\hat{U}(t) = \\exp\\left\\{- i t\\hat{H}(t) \\right\\} = \\hat{P}^\\dagger \\exp\\left\\{ -it\\hat{D} \\right\\} \\hat{P} $$\n",
    "\n",
    "and we have access to $\\hat{P}$ and $\\hat{D}$ from our previous work...\n",
    "\n",
    "### Exercises:\n",
    "\n",
    "1. Try different batch sizes for `pmap` - does that speed anything up?\n",
    "2. Use the linear algebra trick! For extra points, let's play with `Diagonal` and see if that speeds anything up.\n",
    "3. Implement the speedup in 2 above.\n",
    "4. (Harder) have the plot be generated/updated in real time\n",
    "5. How does this scale with the number of workers?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### If you have extra time:\n",
    "\n",
    "It's worth trying the quench with the $XXZ$ model, quenching to and from the $XY$ model to the Heisenberg model. What happens?\n",
    "\n",
    "We've implemented an *instantaneous* quench. What happens if you change your time-evolution function to support *slow* ([adiabatic](https://en.wikipedia.org/wiki/Adiabatic_theorem)) quenches?\n",
    "\n",
    "Quantum quenches are a good example of a process that is a) relatively easy to simulate on a classical computer and b) illuminating of lots of [interesting](https://arxiv.org/abs/1404.6848) [quantum](https://arxiv.org/abs/1704.01974) [properties](https://arxiv.org/abs/1706.01917)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0-rc1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}