4_quantum/VQE_action.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving A transverse field Ising model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Yao\n",
    "using Statistics: mean\n",
    "using LinearAlgebra\n",
    "using Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Hamiltonians"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transverse field Ising Model\n",
    "\n",
    "$H = \\sum\\limits_{i=1}^{N-1} s_i^z s^z_{i+1} + h\\sum\\limits^{N}_{i=1} s_i^x$\n",
    "\n",
    "where $s^\\alpha = \\frac 1 2 \\sigma^\\alpha$ is the spin operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nspin (generic function with 1 method)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "for simplicity, we require an AbstractModel contains `size` and `periodic` members.\n",
    "\"\"\"\n",
    "abstract type AbstractModel{D} end\n",
    "\n",
    "nspin(model::AbstractModel) = prod(model.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "transverse field ising model, `h` is the strength of transverse field.\n",
    "\"\"\"\n",
    "struct TFI{D} <:AbstractModel{1}\n",
    "    size::NTuple{D, Int}\n",
    "    h::Float64\n",
    "    periodic::Bool\n",
    "    TFI(size::Int...; h::Real, periodic::Bool) = new{length(size)}(size, Float64(h), periodic)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "get the bonds of a chain model.\n",
    "\"\"\"\n",
    "function get_bonds(model::AbstractModel{1})\n",
    "    nbit, = model.size\n",
    "    [(i, i%nbit+1) for i in 1:(model.periodic ? nbit : nbit-1)]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get the hamiltonian of a TFI model\n",
    "\"\"\"\n",
    "function hamiltonian(model::TFI)\n",
    "    nbit = nspin(model)\n",
    "    sum(repeat(nbit, Z, (i,j)) for (i,j) in get_bonds(model))*0.25 +\n",
    "    sum(put(nbit, i=>X) for i=1:nbit)*0.5model.h\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tfi_model = TFI(4; h=0.5, periodic=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tfi_h = hamiltonian(tfi_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A quantum circuit as an ansatz\n",
    "\n",
    "Construct the following circuit as our ansatz\n",
    "\n",
    "<div style=\"display:inline-block\"><img style=\"float:left\", src=\"images/fourqubit.png\" width=\"500\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dump(Measure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rxz() = chain(Rx(0.0), Rz(0.0))\n",
    "MeasureOp(nbit::Int, i::Int) = Measure{nbit, 1, AbstractBlock}(Z, (i,), nothing, false)\n",
    "\n",
    "function ansatz_circuit(nbit::Int, d::Int)\n",
    "    circuit = chain(nbit)\n",
    "\n",
    "    for i=1:nbit-1\n",
    "        unit = chain(nbit)\n",
    "        for j=1:d\n",
    "            push!(unit, put(nbit, i=>Rxz()))\n",
    "            push!(unit, put(nbit, nbit=>Rxz()))\n",
    "            push!(unit, control(nbit, 1, 2=>shift(0.0)))\n",
    "        end\n",
    "        push!(circuit, unit)\n",
    "    end\n",
    "    push!(circuit, chain([MeasureOp(nbit, i) for i=1:nbit]))\n",
    "    circuit\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "example_circuit = ansatz_circuit(4, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nparameters(example_circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    gensample(circuit, operator; nbatch=1024) -> Vector of Measure\n",
    "\n",
    "Generate samples from MPS-inspired circuit. Here, `nbatch` means nshot.\n",
    "`operator` is the operator to measure.\n",
    "This function returns a vector of `Measure` gates, results are stored in `m.results`.\n",
    "\"\"\"\n",
    "function gensample(circuit, operator; nbatch=1024)\n",
    "    mblocks = collect_blocks(Measure, circuit)\n",
    "    for m in mblocks\n",
    "        m.operator = operator\n",
    "    end\n",
    "    reg = zero_state(nqubits(circuit); nbatch=nbatch)\n",
    "    reg |> circuit\n",
    "    mblocks\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = gensample(example_circuit, X; nbatch=1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res[1].results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the energy through sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    VQE{MT, BT<:AbstractBlock}\n",
    "\n",
    "Variational quantum eigensolver setup\n",
    "* `MT` is the type of hamiltonian model,\n",
    "* `BT` is the type of an input circuit.\n",
    "\"\"\"\n",
    "struct VQE{MT, BT<:AbstractBlock}\n",
    "    model::MT\n",
    "    circuit::BT\n",
    "    nshots::Int\n",
    "end\n",
    "\n",
    "VQE(model, circuit::AbstractBlock; nshots::Int=1024) = VQE(model, circuit, nshots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loss Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Obtain the energy of a VQE application.\n",
    "\"\"\"\n",
    "function energy(vqe::VQE{<:TFI})\n",
    "    # measuring Z\n",
    "    eng = ising_energy(vqe.circuit, get_bonds(vqe.model), Z; nshots=vqe.nshots)\n",
    "    # measuring transverse field term on basis X\n",
    "    mblocks = gensample(vqe.circuit, X; nbatch=vqe.nshots)\n",
    "    engx = sum(mean.([m.results for m in mblocks]))\n",
    "    eng + vqe.model.h*engx/2\n",
    "end\n",
    "\n",
    "function ising_energy(circuit, bonds, basis; nshots::Int)\n",
    "    mblocks = gensample(circuit, basis; nbatch=nshots)\n",
    "    nspin = length(mblocks)\n",
    "    local eng = 0.0\n",
    "    for (a, b) in bonds\n",
    "        eng += mean(mblocks[a].results .* mblocks[b].results)\n",
    "    end\n",
    "    eng/=4\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# obtain the energy through sampling\n",
    "vqe = VQE(tfi_model, example_circuit; nshots=1024)\n",
    "dispatch!(vqe.circuit, :random)\n",
    "energy(vqe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training\n",
    "\n",
    "Sequential optimization\n",
    "1. obtain the gradient of one parameter a time $\\frac{\\partial E_\\theta}{\\partial \\theta} = \\frac 1 2 (E_{\\theta+\\pi/2}-E_{\\theta-\\pi/2})$\n",
    "2. perform gradient descent update of this parameter $\\theta\\rightarrow \\theta-\\alpha \\frac{\\partial E_\\theta}{\\partial \\theta}$\n",
    "3. sweep all parameters until convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "obtain the gradient of a parameter in a rotation gate/cphase gate.\n",
    "\"\"\"\n",
    "function opgrad(vqe::VQE, r)\n",
    "    dispatch!(+, r, π/2)\n",
    "    E₊ = energy(vqe)\n",
    "    dispatch!(-, r, π)\n",
    "    E₋ = energy(vqe)\n",
    "    dispatch!(+, r, π/2)\n",
    "    0.5*(E₊ - E₋)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Flux: ADAM, Optimise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "train! (generic function with 1 method)"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function train!(vqe::VQE; maxiter::Int=200, α::Real=0.3)\n",
    "    circuit, model = vqe.circuit, vqe.model\n",
    "    rots = collect_blocks(Union{RotationGate, ControlBlock{<:Any, <:ShiftGate}}, circuit)\n",
    "    loss_history = Float64[]\n",
    "    params = vcat(parameters.(rots)...)\n",
    "    optimizer = ADAM(α)\n",
    "    for i in 0:maxiter\n",
    "        grad = opgrad.(Ref(vqe), rots)\n",
    "        Optimise.update!(optimizer, params, grad)\n",
    "        dispatch!.(rots, params)\n",
    "        push!(loss_history, energy(vqe)/nspin(model))\n",
    "        \n",
    "        i%10 == 0 && println(\"Iter $i, E/N = $(loss_history[end])\")\n",
    "    end\n",
    "    loss_history, circuit\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice_size = 4\n",
    "vqe = VQE(TFI(lattice_size; h=0.5, periodic=false),\n",
    "    ansatz_circuit(lattice_size, 2);\n",
    "    nshots=1024);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EG = -0.29742315519647605\n"
     ]
    }
   ],
   "source": [
    "# obtain the exact ground state energy\n",
    "res = eigen(mat(hamiltonian(vqe.model)) |> Matrix)\n",
    "EG = res.values[1]/nspin(vqe.model)\n",
    "@show EG\n",
    "VG = res.vectors[:,1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nparameters(vqe.circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 0, E/N = -0.0887451171875\n",
      "Iter 10, E/N = -0.2269287109375\n",
      "Iter 20, E/N = -0.2728271484375\n",
      "Iter 30, E/N = -0.282470703125\n",
      "Iter 40, E/N = -0.2672119140625\n",
      "Iter 50, E/N = -0.2684326171875\n",
      "Iter 60, E/N = -0.2769775390625\n",
      "Iter 70, E/N = -0.277099609375\n",
      "Iter 80, E/N = -0.272705078125\n",
      "Iter 90, E/N = -0.278564453125\n",
      "Iter 100, E/N = -0.2783203125\n"
     ]
    }
   ],
   "source": [
    "dispatch!(vqe.circuit, :random)\n",
    "loss_history, circuit = train!(vqe; maxiter=100, α=0.3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip5800\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip5800)\" points=\"\n",
       "0,1600 2400,1600 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip5801\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip5800)\" points=\"\n",
       "279.259,1503.47 2321.26,1503.47 2321.26,47.2441 279.259,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip5802\">\n",
       "    <rect x=\"279\" y=\"47\" width=\"2043\" height=\"1457\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  337.051,1503.47 337.051,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  818.655,1503.47 818.655,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1300.26,1503.47 1300.26,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1781.86,1503.47 1781.86,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2263.47,1503.47 2263.47,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  279.259,1479.22 2321.26,1479.22 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  279.259,1150.06 2321.26,1150.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  279.259,820.889 2321.26,820.889 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  279.259,491.721 2321.26,491.721 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  279.259,162.553 2321.26,162.553 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,1503.47 2321.26,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,1503.47 279.259,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  337.051,1503.47 337.051,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  818.655,1503.47 818.655,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1300.26,1503.47 1300.26,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1781.86,1503.47 1781.86,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2263.47,1503.47 2263.47,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,1479.22 309.889,1479.22 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,1150.06 309.889,1150.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,820.889 309.889,820.889 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,491.721 309.889,491.721 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  279.259,162.553 309.889,162.553 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 337.051, 1557.47)\" x=\"337.051\" y=\"1557.47\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 818.655, 1557.47)\" x=\"818.655\" y=\"1557.47\">25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1300.26, 1557.47)\" x=\"1300.26\" y=\"1557.47\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1781.86, 1557.47)\" x=\"1781.86\" y=\"1557.47\">75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2263.47, 1557.47)\" x=\"2263.47\" y=\"1557.47\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 255.259, 1496.72)\" x=\"255.259\" y=\"1496.72\">-0.30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 255.259, 1167.56)\" x=\"255.259\" y=\"1167.56\">-0.25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 255.259, 838.389)\" x=\"255.259\" y=\"838.389\">-0.20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 255.259, 509.221)\" x=\"255.259\" y=\"509.221\">-0.15</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 255.259, 180.053)\" x=\"255.259\" y=\"180.053\">-0.10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 156.859, 775.359)\" x=\"156.859\" y=\"775.359\">Energy</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#009af9; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  337.051,88.4582 356.315,406.697 375.579,652.608 394.844,930.665 414.108,944.327 433.372,944.327 452.636,1165.33 471.9,1207.11 491.164,1156.49 510.429,1090.59 \n",
       "  529.693,998.17 548.957,1108.27 568.221,1169.34 587.485,1117.91 606.749,1166.13 626.014,1320.43 645.278,1268.19 664.542,1267.39 683.806,1238.46 703.07,1262.57 \n",
       "  722.334,1300.34 741.598,1243.28 760.863,1260.96 780.127,1288.28 799.391,1271.41 818.655,1277.03 837.919,1293.91 857.183,1300.34 876.448,1255.33 895.712,1240.06 \n",
       "  914.976,1363.82 934.24,1335.7 953.504,1312.39 972.768,1275.42 992.033,1286.67 1011.3,1338.91 1030.56,1331.68 1049.83,1301.94 1069.09,1359.81 1088.35,1359 \n",
       "  1107.62,1263.37 1126.88,1336.5 1146.15,1323.64 1165.41,1346.95 1184.67,1326.05 1203.94,1293.1 1223.2,1407.22 1242.47,1339.71 1261.73,1352.57 1281,1351.77 \n",
       "  1300.26,1271.41 1319.52,1345.34 1338.79,1274.62 1358.05,1241.67 1377.32,1297.12 1396.58,1335.7 1415.84,1362.22 1435.11,1364.63 1454.37,1376.68 1473.64,1343.73 \n",
       "  1492.9,1327.66 1512.17,1335.7 1531.43,1329.27 1550.69,1367.04 1569.96,1316.41 1589.22,1385.52 1608.49,1362.22 1627.75,1327.66 1647.01,1297.12 1666.28,1359 \n",
       "  1685.54,1328.46 1704.81,1366.23 1724.07,1334.89 1743.33,1317.21 1762.6,1309.18 1781.86,1335.7 1801.13,1352.57 1820.39,1281.05 1839.66,1322.03 1858.92,1338.91 \n",
       "  1878.18,1299.53 1897.45,1306.77 1916.71,1373.47 1935.98,1358.2 1955.24,1326.86 1974.5,1289.89 1993.77,1356.59 2013.03,1297.93 2032.3,1328.46 2051.56,1345.34 \n",
       "  2070.83,1338.11 2090.09,1359 2109.35,1341.32 2128.62,1357.39 2147.88,1302.75 2167.15,1350.97 2186.41,1322.84 2205.67,1389.54 2224.94,1273.82 2244.2,1356.59 \n",
       "  2263.47,1336.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5802)\" style=\"stroke:#e26f46; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  337.051,1462.26 356.315,1462.26 375.579,1462.26 394.844,1462.26 414.108,1462.26 433.372,1462.26 452.636,1462.26 471.9,1462.26 491.164,1462.26 510.429,1462.26 \n",
       "  529.693,1462.26 548.957,1462.26 568.221,1462.26 587.485,1462.26 606.749,1462.26 626.014,1462.26 645.278,1462.26 664.542,1462.26 683.806,1462.26 703.07,1462.26 \n",
       "  722.334,1462.26 741.598,1462.26 760.863,1462.26 780.127,1462.26 799.391,1462.26 818.655,1462.26 837.919,1462.26 857.183,1462.26 876.448,1462.26 895.712,1462.26 \n",
       "  914.976,1462.26 934.24,1462.26 953.504,1462.26 972.768,1462.26 992.033,1462.26 1011.3,1462.26 1030.56,1462.26 1049.83,1462.26 1069.09,1462.26 1088.35,1462.26 \n",
       "  1107.62,1462.26 1126.88,1462.26 1146.15,1462.26 1165.41,1462.26 1184.67,1462.26 1203.94,1462.26 1223.2,1462.26 1242.47,1462.26 1261.73,1462.26 1281,1462.26 \n",
       "  1300.26,1462.26 1319.52,1462.26 1338.79,1462.26 1358.05,1462.26 1377.32,1462.26 1396.58,1462.26 1415.84,1462.26 1435.11,1462.26 1454.37,1462.26 1473.64,1462.26 \n",
       "  1492.9,1462.26 1512.17,1462.26 1531.43,1462.26 1550.69,1462.26 1569.96,1462.26 1589.22,1462.26 1608.49,1462.26 1627.75,1462.26 1647.01,1462.26 1666.28,1462.26 \n",
       "  1685.54,1462.26 1704.81,1462.26 1724.07,1462.26 1743.33,1462.26 1762.6,1462.26 1781.86,1462.26 1801.13,1462.26 1820.39,1462.26 1839.66,1462.26 1858.92,1462.26 \n",
       "  1878.18,1462.26 1897.45,1462.26 1916.71,1462.26 1935.98,1462.26 1955.24,1462.26 1974.5,1462.26 1993.77,1462.26 2013.03,1462.26 2032.3,1462.26 2051.56,1462.26 \n",
       "  2070.83,1462.26 2090.09,1462.26 2109.35,1462.26 2128.62,1462.26 2147.88,1462.26 2167.15,1462.26 2186.41,1462.26 2205.67,1462.26 2224.94,1462.26 2244.2,1462.26 \n",
       "  2263.47,1462.26 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip5800)\" points=\"\n",
       "1676.17,312.204 2249.26,312.204 2249.26,130.764 1676.17,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1676.17,312.204 2249.26,312.204 2249.26,130.764 1676.17,130.764 1676.17,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#009af9; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  1700.17,191.244 1844.17,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1868.17, 208.744)\" x=\"1868.17\" y=\"208.744\">QMPS-Expand</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip5800)\" style=\"stroke:#e26f46; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  1700.17,251.724 1844.17,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1868.17, 269.224)\" x=\"1868.17\" y=\"269.224\">Exact</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = length(loss_history)\n",
    "plot(0:M-1, [loss_history, fill(EG, M)], label=[\"QMPS-Expand\", \"Exact\"], lw=3, ylabel=\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: The problem of vanishing gradient\n",
    "Change the circuit size (as well as model lattice size), show the variance of gradient items.\n",
    "\n",
    "#### Hint: useful functions\n",
    "* using Statistics: var\n",
    "* Use `energy_exact` instead of `energy` to avoid sampling error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "wave_function (generic function with 1 method)"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wave_function(circuit) = zero_state(nqubits(circuit)) |> circuit[1:end-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "energy_exact (generic function with 1 method)"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# obtain the exact <H> through wave function\n",
    "energy_exact(vqe) = expect(hamiltonian(vqe.model), wave_function(vqe.circuit)) |> real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "opgrad_exact (generic function with 1 method)"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function opgrad_exact(vqe::VQE, r)\n",
    "    dispatch!(+, r, π/2)\n",
    "    E₊ = energy_exact(vqe)\n",
    "    dispatch!(-, r, π)\n",
    "    E₋ = energy_exact(vqe)\n",
    "    dispatch!(+, r, π/2)\n",
    "    0.5*(E₊ - E₋)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8-element Array{Float64,1}:\n",
       " 0.07549153741128896\n",
       " 0.09583226381779932\n",
       " 0.11916702798334385\n",
       " 0.10823850341118278\n",
       " 0.09526693782285131\n",
       " 0.09561843871114195\n",
       " 0.10251840223865968\n",
       " 0.10697407003096228"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Statistics\n",
    "\n",
    "function rand_gradient(nbit::Int)\n",
    "    d = 3\n",
    "    circuit = ansatz_circuit(nbit, d)\n",
    "    model = TFI(nbit; h=0.5, periodic=false)\n",
    "    vqe = VQE(model, circuit, nshots=1)\n",
    "    dispatch!(circuit, :random)\n",
    "    r = rand(collect_blocks(RotationGate, circuit))\n",
    "    opgrad_exact(vqe, r)\n",
    "end\n",
    "\n",
    "maxN = 9\n",
    "grad_vars = [(println(N); mean([abs(rand_gradient(N)) for i=1:300])) for N=2:maxN]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VPW9//H3OWeSQAgQIhCiIewJm4ogSJFABAQRQUVcSxUVoa2I9vYi1nu73N7WH13sNS7VgAoqirigtCoCAmGxSGVTARMgEPZNIISwBOac7++PVAqyGMhMZns9H488HmRmkvmGfHO+7/mcz3yPZYwxAgAAQMDYoR4AAABAtCFgAQAABFhEBqzDhw9r+fLlOnz4cKiHAgAAcJqIDFgFBQXq1KmTCgoKLvh7HDhwIIAjQqxjPiGQmE8IJOZTaERkwAoE13VDPQREEeYTAon5hEBiPoVGzAYsAACAYCFgAQAABJgv1AMAAADRq7y8XIWFhfL7/aEeSpU4jqMmTZooOTm5Uo8nYAEAgKAoKipShw4dVFZWFuqhBMwDDzygF154QbZ97pOABCwAABBwnufp/vvvV/369TVjxgwlJiaGekhVcuzYMS1YsECPPfaYJGn8+PHnfDwBCwAABNyOHTs0f/58vfHGG+revXuohxMQXbt2lSSNHTtWf/zjH895upAmdwAAEHB79uyRJLVo0SLEIwmsHj16SJI2bdp0zscRsIAq4nrpAHA6z/MkST5fdJ0si4+Pl/T9+4sRsIAqGrPEU795idp3lKAFAJGouLhYOTk5qlu3rjp06BCQ70nAAqpgzxGjZ9d4WrbP1rUz/NpfTsgCgEhTp04d/e53v9Mbb7wRsO9JwAKqYHyBJ0vStOwjKj4o9Zvh6sAxQhYAhKM///nPGjFixInPS0pKVL9+fUlS9+7dVatWrYA9V3SdGAWq0XHP6K9rPA1taalHQ1ezr/ep94d+XTfD1cz+jurEW6EeIgCEjcN+o4KS4D5H62Qp0Xf2Y+/w4cOVmZl54h2AEydO1I033qiUlJSAj4WABVygdzYYbT8sjW7vSJI61rc0+3pHfT5ydf3Hrj7u7ygpjpAFAJJUUCJ1ei+4u7kvu9mnjvXPfn9ycrKGDBmil19+WT/72c/0/PPPa+rUqUEZCwELuEC5qz31utjSpSmW9u2ruO3KBrZm9peu/cjVgI9dfXSdo1qELABQ6+SKABTs5/g+o0eP1qBBg9SmTRs1aNBAV1xxRVDGQsACLsCS3Z6W7Daa3tc57b6rGtr6uH9FP9bAma4+uM45Z8kaAGJBos86Z3WpurRu3VrNmzfXiBEj9Mc//jFoz0OTO3ABcld5al5bGtD4zMGpW6qtj/o5WrLH6MZZro74aXwHgHDxwAMPyO/3a8iQIZKkw4cPKz09XbfeeqvWrFmj9PR0/eIXv6jSc1DBAs7TtkNGb28w+tNVthz77JWp7DRbH/aTrv/Y1eDZrt671lENKlkAEHLz5s3TT3/6U8XFxUmSEhMTtXXr1oA+BxUs4Dw9v8ZTTZ90X9b3//nkXGzr7/0c5e8wGvKJq3KXShYAhMr27dvVunVrLV++XI888khQn4sKFnAejvqN8go83ZtpV3obht6X2Hr/WmnQLFe3z3H1Vm9H8Q6VLACobhdffLEKCgqq5bmoYAHn4Y0io71HpYfand+fTr/Gtt671tFHW4zunOvquEclCwCiGQELqCRjjHJXuRqQYall3fOvQF2fYevdPo7+vtlo6DxXfkIWgCjmOBXvsj527FiIRxJYhw8flqQT/VtnwylCoJLm7zD6cp/0ZNcLf10ysImtqb2kW+e48lmuXs1xztkoDwCRqmnTpqpRo4Z++9vf6le/+pXi4+NDPaQq8fv9Kioq0i9+8QvVrl1bmZmZ53w8AQuopNxVntrVk3pfXLVAdHMzW1N6SXfOdeVYrib2JGQBiD5169bV9OnTdeONN2rGjBmhHk7A5OTkaO7cuUpISDjn4whYQCVsLDWavsnohe6OLKvqYejW5rZcI/1wniuf7erFHo7sAHxfAAgnffv21c6dO1VcXCzXdUM9nCqxbVsNGzZUo0aNZNvffyaDgAVUwrNrPNVLkIa2ClwIuqOFLb8n3Z1fEbJe6E7IAhB96tatq8svvzzUw6h2BCzge5QdN3qp0NOP29gBv+TN0Fa2/Ea6b74rn+XpuavtgFTIAAChRcACvscraz2VHZcebBucN90Oy6yoZD2w0JXPlnJ/QMgCgEhHwALOwTNGT6/2NLippcZJwQs9w1vbco3Rjxd58lkV71QkZAFA5CJgAecwc6vR2gPSyz2Cv2XcyDaO/J406h+e4mxpXBdCFgBEKgIWcA65qzx1qm+pW2r1BJ0H2zk67kk/+8yTz5Z+dyUhCwAiEQELOIuv9xvN3Gr0ak5gtmaorEcudeQ30pglFZWs33Ryqu25AQCBQcACzuLp1Z5Sa0q3Na/+CtJ/XlZxuvAXn1f0ZP13R0IWAEQSAhZwBvvLjV5d5+nRy2wlOKE5RfdYh4qQ9ctlFacLH+tAyAKASEHAAs7gxQJPfk/6cZvQXg/9vzs6Om4qKllxtvTzywhZABAJCFjAd/g9o2fXeLqjhaXUxNA3mP+mY8U+Wf+5pKKS9XB7QhYAhDsCFvAd0zcZbS4LnyBjWZZ+d6Wt4570yOKKnqwH24XH2AAAZ0bAAr4jd5Wn7qmWOtYPffXqW5Zl6Q9d7BP7ZPnsin2zAADhiYAFnGTFN0YLdxq93Tv8wotlWXqya8W1Cyt2fLd0f+vQ9ogBAM6MgAWcJHeVq4wk6aam4VO9OpllWcr9wanXLrwnk5AFAOGGgAX8y67DRlOKjP73Sls+OzwDllQRsp692pbfGN0735VjSUNbEbIAIJwQsIB/ySvw5FjS8KzwDyu2ZemF7o78nqt75ldUsu5oEf7jBoBYQcACJB1zjZ5f4+meTFspNcK3enUy27I0IbsiZA2d58pnSUOaE7IAIBwQsABJb20w2nlEGt0usgKKY1ua2NORa1zdObeiknVT08j6GQAgGnEkRswzxih3lae+l1hqUy8yqlcnc2xLr+Q4GtzM0m1zXP19kxfqIQFAzCNgIeYt3m209Bujh9tH7p+Dz7Y0+RpHAzMsDfnE1YwthCwACKXIXVGAAMld5SmzrnRd48irXp0szrY0pZej6xpbunm2q1lbCVkAECoELMS0LWVG7240eqidLduK7IAlSfGOpbd6O+pzsaUbZ7mas42QBQChQMBCTPvrGk+1fNI9UbSPVIJj6Z0+jnqmWRo401X+dkIWAFS36FlVgPN02G80vsDT/Vm2asdHfvXqZDV8lt671tHVjSzdMNPVwh2ELACoTgQsxKzX1xvtL5dGRdjWDJVV02dpel9HXRpYun6mq3/sImQBQHWJzpUF+B4VWzO4GtTEUvM60VW9Olmiz9Lf+znqeJGl62a4WrKbkAUA1SHoAeuhhx5Sdna2nnjiidPuO378uG677Tbl5OToT3/6kySpuLhYaWlpysnJ0d133x3s4SFGzd1utHq/InprhsqqFWfpg36OLkux1G+Gq6V7CFkAEGxBXV2WLl0qn8+nhQsXavny5dq1a9cp90+bNk3dunVTfn6+li1bpj179kiSBgwYoPz8fL366qvBHB5iWO4qT5elSDlp0Vu9OlnteEsfXeeoTbKlvjNcrfjGhHpIABDVghqwlixZol69ekmSevbsqWXLlp1y/8aNG3XppZdKktq0aaPPP/9ckjRz5kxlZ2fr9ddfD+bwEKOKSo0+2Gw0up0jKwq2ZqisOvGWPu7vqGUdS30+8uuLvYQsAAiW0wLW6NGj1bRpU1mWpZUrV55y37p169StWzdlZmaqc+fOWr169Tm/eUlJierUqSNJql27tkpKSk65PysrS/Pnz5cxRgsXLtSBAweUlpamwsJCzZo1S3l5edq7d29Vf0bgFM+s9pSSIN3VMnbC1bfqxlua2d9R09pSn4/8WrWPkAUAwXBawBoyZIgWLVqkJk2anPbgkSNHasSIEVq7dq3Gjh2rYcOGSZLWrFmjnJycUz7GjRun5ORklZaWSpIOHjyo5OTkU77foEGDtGfPHl177bVq0KCBGjZsqISEBCUmJqpmzZrKzs5WUVHRWQdfVlam0tLSEx/l5eVV+b9ADCg9ZvRyoaeRbWzV9MVewJKkegmWZvf36ZJEqfdHfn29n5AFAIFmGWPOeHRt2rSp3n//fXXo0EGStHv3brVs2VL79u2Tz+eTMUZpaWlatGiRWrZsecZv/vnnn2vKlCn6y1/+oltvvVXPPvusUlNTT3ucMUb33HOPXnjhBXmep6SkJBljdN1112nSpElKS0s75fHLly9Xp06dTvs+jz76qMaOHVupH3z//v2qV69epR6L6JG3Pk6/+jJBK/of0sU1AxcsInE+7S23dNOCmtp7zNL0HofVqjZBK1xE4nxC+GI+VY+UlJRTPvdV9gu3bNmitLQ0+XwVX2JZljIyMrR58+azBqzOnTtr0qRJys7OVr9+/ZSamqqdO3fqtdde05gxY7R582bdfffdsm1bjzzyiBITEzVr1iw9/vjjiouL05AhQ04LVyebP3/+iQAoSQkJCUpISKjsj3Tafwaim+sZvbTRryHNLbW/JPAHm0ibTymS5g0yuuYDvwYvStL8G3xqWTc2q3rhKNLmE8Ib86n6VTpgXajnnnvulM8bNWqkMWPGSJIyMjKUn59/yv19+/ZV3759K/W9k5KSTvR4Ad/noy1GRaXS5Jzo35qhshrWtDRngE85H/h1zYd+zb/BF9X7ggFAdan0StO4cWPt2LFDfr9fUsVpvc2bNysjIyNogwMCKXeVpy4NLHVNJWCdrFGipbkDfKrpSNd86FfxQU4VAkBVVXqladiwoTp27KjJkydLkt59912lp6ef9fQgEE5W7TOas93ExMaiF+LiWhUhy2dJvT70a3MZIQsAquK01WbkyJFKT0/X1q1b1a9fv1MCVF5envLy8pSZmalx48Zp4sSJ1TpY4EI9vdpVWqI0pBmnv84mPcnSvBt8MqYiZG07RMgCgAt1Wg9WXl7eWR+clZWlxYsXB3VAQKDtPWr02jqj/+pgK94hYJ1Lxr9CVs8P/LrmA7/mD/QpLZH/MwA4X5wvQdSbUODJSBrZhuleGU1rV5wuPOJWVLJ2HaaSBQDnixUHUe24Z/TcGk93tbDUoCaVmMpqUaciZJUeqwhZu48QsgDgfBCwENXe22i09ZD0cHsn1EOJOK3qVoSsfeUVl9X55ighCwAqi4CFqJa72lNOmqXLL6J6dSGykiv2ydp1RLr2I7/2EbIAoFIIWIhaS/d4+scutmaoqrb1LM253qeth6S+M1yVlBOyAOD7sPIgauWu8tSstjQwg+pVVbVPsfTJ9T5tPGjUb4arA8cIWQBwLgQsRKUdh42mbjAa1daWYxOwAuHyiyzNvt6ntQeMrpvhqpSQBQBnRcBCVHrha0/xtnRfFlM8kDrWtzT7ekdflxhd/7GrsuOELAA4E1YfRJ1y1+iFrz0Ny7SVnED1KtCubGBrZn9HX+4zGvCxq0OELAA4DQELUefNIqPdR6SH2jG9g+WqhrY+7u9o+V6jQbNcHfYTsgDgZKxAiCrGGOWuctW/saWsZKpXwdQt1dZH/Rx9ttvoplmujhKyAOAEAhaiyqKdRiv2Sg9TvaoW2Wm2PuznaNFOo5tnuyp3CVkAIBGwEGVyV3tqnSz1Tad6VV1yLrb1936O8ncYDfnE1TFCFgAQsBA9Nh00eq/YaHQ7W5ZFwKpOvS+xNb2vo9nbjG6b4+q4R8gCENsIWIgaz63xVCdOursV0zoU+qbbmtbH0YwtRnfOJWQBiG2sRIgKh44bTSjwNLy1rVpxVK9C5foMW+/0cfS3TUY/mufKT8gCEKMIWIgKr63zVHpcGtWWKR1qA5vYmtrL0Tsbje7Jd+USsgDEIFYjRDxjjJ5e7emmJpaa1KZ6FQ5ubmZrSi9HUzcY3beAkAUg9vhCPQCgqmZvM/q6RHqhO68XwsmtzW25RvrhPFc+y9WEHo5s3nwAIEYQsBDxcld56nCRlN2IxTvc3NHClt+T7s535bM9Pd/dJmQBiAkELES0tSVGH20xmtjTYWuGMDW0VUUl6975rmo4Um43J9RDAoCgI2Ahoj292lODGtIdzQlX4eyeTFuH/UY//dTTVQ0t3dWS07kAohtHOUSsknKjSWs9/biNrRo+Ala4+3EbW0NbWhq5yNW6AzS9A4huBCxErJcLPR3zpJ+wNUNEsCxLf73aUVpN6fY5fq5bCCCqsTIhIrme0TOrPd3W3FJaItWrSFE73tJbfXxaUyKNWeKFejgAEDQELESkv282Ki6THm7PFI40HS6y9ORVtp5Z7en9YkIWgOjE6oSIlLvK0w8aWurcgCkciX7a1tbgppbune9q00FOFQKIPqxOiDhf7jXK32H0CNWriGVZll7q4ahuvLgwNICoxAqFiJO7ylV6LenmZvReRbLkBEtv9nL0+R6jXy7lVCGA6ELAQkTZc8To9SKjB9vairMJWJGua6qtJzrb+sMXnj7eQsgCEBp7jhi9t9GTMYGrphOwEFHGF3iyJT3QmqkbLX5+ma3+jS3dne9q+yFOFQKoPoUlRj9e6Cpjil/D5rsqORa4780qhYhx3DP66xpPQ1tZuqgG1atoYVuWXunpKM6Whs5z5dKPBSCIjDGav8PToJl+tX7br/c3efqvDraK7vCpXkLg1hYCFiLGOxuMth+WRrfjWnbRpkFNS2/0cjR/p9HvV3KqEEDg+T2jN4s8dXnfVc4HropKjV7MdlR8h0//3dFR/QC/cOdahIgYuas99b7YUvsUqlfRqGearV93NPqf5Z56plnqmcbrPwBVV3rM6KVCT0+t8rS5TOp9saWPrnPUL92SbQVvPSFgISIs2e1pyW6jv/WlehXN/quDrfztRnfNdbVysKUGNQnTAC7MljKjp1d7Gv+1p8N+6Y4Wln5+maMOF1XPcYWAhYiQu8pTizrSgAwW3Gjm2JYmX+OowzS/7pnv6oN+TlBfYQKIPsu/MXryS1dvbTCqFVdxofmH2tlKT6reYwkBC2Fv2yGjtzcY/bmrzWIbAy6uZenVHEf9P3b15JeexlxO1RLAuXnGaMYWoye/9DRvh1GTJOnPXW3dl2mrdnxo1g0CFsLe82s81fRJ92bSkxMrrmtsa+zlRo9/7im7kaWuqfzuAZzuqN9o8nqjv3zl6usSqUsDS1N7ORrczJIvxHslErAQ1o76jfIKPN2baatOiF6FIDT+90pbC3YY3THX1YrBVkDfPg0gsn1ztGLbnufWeNpzRBrUxNL4bFtXp1qywuRMBwELYe2NIqO9R6WH2lHBiDVxtqUpvSr6sYYvcPVOHydsDpwAQmNtidH/rfI0aa0nS9KwTFuPtLeVmRx+xwYCFsKWMUa5q1wNyLDUsm74/fEg+JrUtjSxp6ObZ7v66xpPD7IHGhBzjDFauNPoya88/X2TUYOa0uMdbP2krR3wvasCiYCFsDV/h9GX+6Qnu1K9imU3NbX1UDuj//jM09WN7Gp7izWA0PJ7Ru9urAhWn+8xapssTch29MOWlmr4wv84QMBC2Mpd5aldvYpN4RDb/nSVrU93ebrtE7+W3ewL2buCAATfwWNGLxZ6yl3laVM1bgwaaAQshKWNpUbTNxm90J2+G0gJjqU3e/nU8T2/fvKpq9dymBdAtNn6r41B807aGPT9atwYNNAIWAhLz67xVC9BGtoqMv+wEHit6loa393RXfNc9b7Y6N4s5gYQDVZ8Y/TkV66mFoV2Y9BAI2Ah7JQdr7hu1I/b2EqMgPPsqD53trQ1Z7unBz91dVVDS23rMT+ASOQZo4+3VPRXzd1esTHon66ydX9W6DYGDTQCFsLOK2s9lR2XHmxLcztO93Q3R4t3+XX7HL+W3OQjhAMR5Lsbg3YOo41BA40VDGHFM0a5qzwNbmqpcYSXhxEciT5Lb/X2qahUemSxG+rhAKiEb44a/e9yV03e9GvEQleZdS0tuMHRkhsd3dbCjrpwJVHBQpj5eIvRulJpYk+yP86uXYqlZ7o5Gr7QVa+LPd3RgvkChKNvNwZ9Za0nKbw3Bg00AhbCSu4qT53qW+qWGv1/fKia+7IszdluacRCV1fWZzNaIFwYY7ToXxuD/u1fG4P+IgI2Bg00XvYhbHy932jWNqOH29u8BR/fy7IsvdDdUWpN6Y65rspdE+ohATHN7xlNLfJ01XRXPT5wtfaA0YRsR5vu8OmXHZ2YClcSAQth5OnVnlJrSrc1j60/Qly4OvGWpvb26at9RmP/6YV6OEBMOnjM6KmvXLWc6tcdc13VjpM+7Odo1RCf7m9tR8Su68HAKUKEhf3lRq+u8zT2clsJTmz+MeLCdKxv6c9X2Rq92FOviy0NasLrRqA6fLsx6PgCT4eOV2wM+t6ljq6ozzFcImAhTLxY4MnvSSNbszji/I1qZ2vudqNh812tHGwpg3egAkGzcq/Rk1+6evNfG4OOaG1rdBRsDBporGYIOb9n9OwaT3e2sJSayB8ozp9lWXq5p6PacdKdc10d9+jHAgLJM0YfbfbU+0O/rpjm18KdRn+6ytaWO33641UO4eoMCFgIuembjDaXSQ+3d0I9FESwegmW3uzlaMluo18vox8LCISjfqOXCjxd+o5fA2a6OnhcmtrL0frbfXrkUidqdl0PBk4RIuRyV3nKbmRx3h5V9oNUW7+/0uixzz3lpFnqm85rSOBCfHPU6Pk1np5d42nPEWlQE0svdLfVvZHFu7wriYCFkFrxjdHCnUbv9KF6hcAYc7mteTuMfpRf0Y+VxmlnoNLWHTD6v688TVrryahiY9CfxcjGoIHGyzuEVO4qVxlJ0o1N+ONFYNiWpVdzHDmWNHSeK5d+LOCcKjYG9XTTLL+y3vLr3WJPv+hga8tdPj3f3SFcXSACFkJm12GjKUVGo9pG53WoEDoNa1p6/RpH87Yb/b+V9GMBZ+L3jN4q8tR1uqvsv7MxaKBxihAhk1fgyWdLw9maAUFwzcW2ftXR6NfLPfVIs9QjjXkGSBUbg75U6Cl3lafiMqnXxZY+7OfousaWbPqrAoaAhZA45lY0UN7dyla9BP6gERy/vMJW/g6jO+e6+uIWi1fkiGlby4yeWe0p76SNQaexMWjQELAQEm9tMNp5RBrdjqoCgsexK04Vdpjm17B8V3/r5/AKHTHnqxJbL37h15tFRok+aWQbNgatDqxuqHbGGOWu8tT3Ektt6vEHjuC6pFZF0/uHWyreHQXEinLXaPQ/XOXMqXViY9Ctd7ExaHWhgoVqt3i30dJvjD7sx9YMqB79G9sac5nRY/+s2HOtS0NeWyK6rS0xumOuX2tKpP93+VH955VJvJmomnGUQbXLXeUps650XWP+2FF9ft/ZVqcGlm6f46qknK0bEL0mr/PU8T2/Dvmlz270aUTL44SrECBgoVptKTN6d6PRQ+1semFQreLsikvp7D8mPbDQlTGELESXsuNGw/L9+lG+q1uaWVp2s08dLuI4GyoELFSrv67xVMsn3dOKqYfq17S2pZd7OHpno9ELX9OPhejxxV6jK9/z652NRq/0dPRKjk9JcYSrUGKVQ7U57DcaX+Dp/iybC4QiZAY3s/VgW1s/+8zTF3upYiGyGWP01zWurpruVw1HWnazT3dnsrSHA34LqDavrzfaXy6NYmsGhNifr7LVJlm6bY5fZccJWYhMJeVGt85x9eCnnh7IsvXZjT5lcVmbsMFKh2pRsTWDq0FNLDWvwwEAoVXDZ2lqL5+2HZJ+usgN9XCA8/bZLk8dpvk1Z5vRtD6OnrnaUQ0fx9ZwQsBCtZi73Wj1funh9kw5hIfMZEt52Y5eW2/0ylr6sRAZPGP0xy8qrh14SS1LKwf7dHMzjqvhiH2wUC1yV3m6LEXKSeMVFsLHD1vamrPN008/ddWlARvfIrztPmJ0d76rWVuNHutg63862Ypj+4WwRexF0K0/YPTBZqPR7RxZbM2AMPNMN0dNkqTb5/p1xE8/FsLTnG2eLn/XrxV7jWb2d/REZ4dwFeYIWAi6Z1Z7SkmQ7mrJwQDhp1acpbd6+7TugPSzxZwqRHjxe0b//bmraz9y1T7F0heDfbo2naU7EvBbQlCVHjOauNbTyDa2atKAiTDVPsXS090c5RV4equIkIXwsKXMKOcDV+O+8PT7zrZm9nfUKJHjaKSgBwtBNXGtpyN+6adtyfIIb8OzLM3ZZumBha6ubMC7XRFaf9vkadh8V0k+af4Njq5uxDE00vAbQ9C4ntEzqz3d2tzSJbVYrBDeLMvS+GxHDWpKt89xdcylHwvVr9w1evgfrm6c5apno4p3CRKuIlNQf2v79u1Tp06dlJSUdNbHPPTQQ8rOztYTTzxxztsQeT7aYlRUytYMiBx14iuuV/jFPqPH/smpQlSvdQeMfjDdrxe+9vRMN1vTrnWUUoMXp5EqqCtf7dq1NXv2bHXt2vWM9y9dulQ+n08LFy7U8uXLtWvXrjPehsiUu8rTVQ0tXdWQgIXIcWUDW3+6ytb/rfL0902ELFSP19d76vieX2XHpc9u9GkU77qOeEFd+eLi4pSSknLW+5csWaJevXpJknr27Klly5ad8TZEnlX7jOZsN3qYy+IgAo1uZ2tQE0vD5rvaUsapQgTPoeNG9833a+g8Vzc3sbTsZp+uqE+wiganrX6jR49W06ZNZVmWVq5cecp969atU7du3ZSZmanOnTtr9erVVXrykpIS1alTR1JFtaukpOSMtyHyPL3a1cWJ0pDmHCgQeSzL0sQejmr5pDvnuvJ7hCwE3pd7ja5836+pG4wm9XT06jU+1Y7nmBktTgtYQ4YM0aJFi9SkSZPTHjxy5EiNGDFCa9eu1dixYzVs2DBJ0po1a5STk3PKx7hx4773yZOTk1VaWipJOnjwoJKTk89429mUlZWptLT0xEd5eXmlfmgE196jRq+tM/ppW3YZRuRKqWFpSi9Hn+02+s0yThUicIwxemGNqy7T/Yq3peU3+3RPJtX+aHPaNg09evQ44wN3796tpUuXatasWZKkW265RaNGjdLxQ3r4AAAgAElEQVT69evVtm1b5efnn/eTd+nSRVOmTNHAgQO1YMEC3XbbbWrQoMFpt33rueee03PPPacjR45IqjiFeLJHH31UY8eOrdRz79+//7zHi8rJLYyXMfG6tVGp9u2LjVf+zKfo1CZe+kXbeP1+Zbw6JpUpJ7V6LgzNfIpeB45Jjyyvob9ti9N9zY/pfy8rVw1P2rcveM/JfKoe322JqvQ+WFu2bFFaWpp8voovsSxLGRkZ2rx5s1q2bHnWr+vTp49WrFihPn366KmnnlL9+vX12muvacyYMercubMmTZqk7Oxs9evXT6mpqUpNTT3ttm89+OCDevDBB7V8+XJ16tRJ8+fPV4cOHU7cn5CQoISEhAv+z0DVHfeMJm7064ctLWWm1Qv1cKoV8yk6/c8PjJaUuPrpskStHOyrto0emU/RZ8luT3fMc7X/mPROH0e3NKslqVa1PDfzqfoFfaPRTz755LTbxowZc+Lfzz333Gn3n+m2M0lKSjrRr4Xw8N5Go62HpIfbO6EeChAQtmXptRxHHab59aN8VzP7O7J5dxfOg2eMnvzS0+Ofe7qygaV5NzhqWps5FO0qfdK3cePG2rFjh/x+v6SKc8ibN29WRkZG0AaHyJO72lNOmqXLLuLggeiRmmhp8jWO5mwzGreSfixU3u4jRgM+dvXoPz39/DJbCwYSrmJFpQNWw4YN1bFjR02ePFmS9O677yo9Pf2cpwcRW5bu8fSPXYaNRRGVel9i67+usPXLZZ4W7iBk4fvN3eapwzS/ln1j9PF1jsZ1cXjjTww5bSUcOXKk0tPTtXXrVvXr1++UAJWXl6e8vDxlZmZq3LhxmjhxYrUOFuEtd5WnZrWlgRkcQBCdft3R1tWplu6a52rv0dh4AwfOn98z+tVSV30+ctU22dIXt/jUrzEvPGPNaT1YeXl5Z31wVlaWFi9eHNQBITLtOGw0dYPRuM62HF6hIUr5bEtvXFPRjzVsvqu/9WW3bZxqa5nRXfNc/WOX0e+utDX2co6JsYpIjYB44WtP8bZ0XxZTCtEtPcnSKzmOPths9NQqThXi3/6+ydPl0/wqPmiUf4Ojx69wCFcxjNUQVVbuGr3wtadhmbaSEziYIPoNyLD180ttjf2np8/3ELJiXblr9LPFrgbNctW9kaWVg33q3ojlNdYxA1BlbxYZ7T4iPcR1BxFDnuhs64qLLN0+x9WBY/Rjxar1B4y6/c2vv67x9PQPbL1/raOUGrzQBAELVWSMUe4qV/0bW8pK5qCC2BHvWHqzl6N95dKIha6MIWTFminrPXV8z6/SY9LiG316qD09efg3AhaqZNFOoxV7pYepXiEGNatj6cVsR29tMBpfwKnCWHHouNH98/26a56rQU0sLb/Zp471CVY4VdB3ckd0y13tqXWy1Dedgwti05Dmtn6y3ejhxZ5+0NBmk90o99U+o9vn+LWpTJrY09E9rSyqVjgjyg64YJsOGr1XbDS6nc0BBjHtL11tZdWVbp/rV9lxThVGI2OM8r521eV9v3y2tPQmn4ZlcuzD2RGwcMGeW+OpTpx0dyumEWJbDZ+lqb192lImjfrUDfVwEGAHjhndPsfVjxd5ujfT1pIbfWpTj2CFc2NlxAU5dNxoQoGn4a1t1YrjQAO0Trb016sdvbLO6NW19GNFi3/u9nTFNL9mbTN6u7ejv3Z3VNPHMQ/fj4CFC/LaOk+lx6VRbZlCwLfuzrR1TytLP/3UVWEJpwojmWeM/vylq6v/5qphTUsrbvZpSHOOd6g8ZgvOm2eMcld7uqmJpSZcFR44xbNXO2qcJN02x68jfkJWJNpzxOiGma7GLPH0H5faWjjQUbM6HOtwfghYOG+ztxoVlEiPtGf6AN+VFGdpai+fCg9IP/+MU4WRZt72isvdLN1jNOM6R3+4ylEcl7vBBWCFxHnLXe3piouk7o046ABnctlFlnJ/YOv5rz29s4GQFQn8ntGvl7nq/aGrNskVl7u5rjFLJC4c+2DhvBSWGM3YYjSpJzsWA+cyorWtOduM7l/gqmN9S805xRS2tpYZ/XCeq0W7jH7bydYvOthcpBlVRjzHeXlmtaeGNaU7WnDwAc7FsixN6OHoohrSHXNdHXPpxwpHH2zy1GGaXxsOGuUPcPTfHR3CFQKCgIVKKyk3mrTW04/b2EpwOAAB36duvKWpvR2t3Gv0+OecKgwnx1yj/1jsauAsV91SK04JZqexJCJwmE2otJcLPR3zpJ+0YdoAldW5ga0/dLH15FeePtxMyAoHRaVG3f7m6tk1np76ga3pfR1dVIMXjQgsVkpUiusZPbPa0+3NLTVK5EAEnI9H2tu6IcPSPfmutpZxqjCU3iyq2Dj0wDGjxYN8erg9/aQIDgIWKuXvm42Ky6SH2ZoBOG+WZWlST0c1fdJd81z5PUJWdTvsN3pggV93znU1MMPSspt96tSAYIXgYbVEpeSu8tQt1dKVDZgywIW4qIalN65x9Okuo98u51RhdVq1z6jze369vt7opR6OJl/jqE484QrBxWqJ7/XlXqP8HUYPt2O6AFWRnWbrt51s/W6FpznbCFnBZkzFNVM7v++XbUnLbvbpviybU4KoFqyY+F65q1yl15JubsZBCaiqxy631fsSS0Pnudp1mFOFwXLgmNGdc12NWOjqnla2/nmTT23qcQxD9SFg4Zz2HDF6vcjowbY2l4sAAsCxLb2W48iT9KN8V54hZAXa53s8dZzm14wtRm/1dvRCtqOaPo5fqF4ELJzT+AJPtqQHWjNVgEBplGjp9WscfbLN6A9fcKowUDxj9OSXrrpNd1W/RsXeVrc259iF0GDm4ayOe0Z/XeNpaCuLPWKAAOtzScUlWX651NOnOwlZVfXNUaOBM1395xJPj1xqa+FAR824PBFCiICFs3png9H2w9Lodk6ohwJEpf/pZKtrQ0t3znW17yinCi/U/B2eLn/Xr3/uMfroOkd/uspRPFebQIgRsHBWz67x1PtiS+1TOFABweCzLU3p5eiQX7p3gStDP9Z5cT2j3yxz1etDV1l1LX0x2Kf+jVnWEB58oR4AwtPqfUb/2GX0dm+qV0AwNU6q2IR00CxXT6/29HD72P2bcz2jI6502H/yh/nO5/++/f1NRgt3Gv2mo63HO9hcpBlhhYCFM5pQ6KlhTWlQEw5YQLANbGLrkfZGY5Z4ujrVUvMwy1jGGB09LfhIR9wzhZ/zvO2k71vuVm48lqREn9Q4SZo3wFEPLtKMMETAwmmO+o1eXefpgSybPgagmvyhi61FO43umOvqkxwppZJfd9wLQMg5+eMMoemIX6rsycsaTkX4+fajpiMl+qwTn9ev8e199imPq/iwznzbd75ngiM2C0XYI2DhNO8WG+0vl4azNQNQbeIdS2/2dtRxml+3LkpUmxT/KdWdbwPRke+EH38lk49jSbXOEWjqxEuNEr8n+DjnDkM1HHGaDvgXAhZOM6HA0zVpllrV5UAJVKcWdSy90cvRrz93VVz27+CSkvCdqs9JQafm2ao+37mNjYKB6kXAwikKS4zm7zB645owawIBYsSADFs/SDqslJTKniQEEI44B4RTvFjo6aIE6eamvNoFAOBCEbBwQrlrNGmtp7tb2arBdbsAALhgBCycML3Y6JujXHcQAICqYiXFCRMKPXVPtdSmHtUrAACqgoAFSVJRqdEn2wzVKwAAAoDVFJKklwo9JcdLtzanegUAQFURsKDjntHLhZ6GtrRVk+Z2AACqjIAFfbDJaNcRaUQbpgMAAIHAigqNL/B0VUNLl6ZQvQIAIBAIWDFu00GjmVuNHshiKgAAECisqjHu5bWekuKk21tQvQIAIFAIWDHM7xm9VOjprha2kuIIWAAABAoBK4Z9vMVo2yGa2wEACDRW1hg2vsBTx/pSx/pUrwAACCQCVozadsjowy1GI9i5HQCAgGN1jVEvF3qq4Uh3tmAKAAAQaKyuMcgzFc3tdzS3VCee04MAAAQaASsGzd5qtKmM5nYAAIKFFTYGjS/wdGmK1KUB1SsAAIKBgBVjdh42+tumiuZ2yyJgAQAQDASsGDNprSefLf2wJb96AACChVU2hnjG6MVCT7c1t1QvgeoVAADBQsCKIfO2GxWVSg+w9xUAAEHFShtDJhR4apMsXZ1K9QoAgGAiYMWIPUeM3is2eoDmdgAAgo6AFSNeXedJkn7Uil85AADBxmobA4wxmlDg6ZZmlurXoHoFAECwEbBiwMKdRoUHaG4HAKC6sOLGgAkFnlrWkXLSqF4BAFAdCFhRbt9Ro7c30twOAEB1ImBFucnrPbmedA/N7QAAVBtW3Sj2bXP7TU0tpSZSvQIAoLoQsKLYZ7uNVu2nuR0AgOrGyhvFJhR4apok9bmE6hUAANWJgBWlDhwzerPIaHhrWzbN7QAAVCsCVpR6Y72nY550bya/YgAAqhurbxQyxmh8gacbMixdXIvqFQAA1Y2AFYWWfWO0ci/N7QAAhAorcBSaUGCUXku6Lp3qFQAAoUDAijJlx43eKPJ0f5YtxyZgAQAQCgSsKPNmkdFhv3RfFr9aAABChVU4yowv8HRduqWMJKpXAACECgErinyx1+jzPUYjaG4HACCkWImjyIQCT2mJ0oAMqlcAAIQSAStKHPYbTV7v6d5MWz6a2wEACKmgBqx9+/apU6dOSkpKOutjHnroIWVnZ+uJJ56QJBUXFystLU05OTm6++67gzm8qPL2BqMDx6ThNLcDABByQV2Na9eurdmzZ6tr165nvH/p0qXy+XxauHChli9frl27dkmSBgwYoPz8fL366qvBHF5UGV/g6dpLLDWrQ/UKAIBQC2rAiouLU0pKylnvX7JkiXr16iVJ6tmzp5YtWyZJmjlzprKzs/X6668Hc3hRY/U+o3/sorkdAIBwcdqKPHr0aDVt2lSWZWnlypWn3Ldu3Tp169ZNmZmZ6ty5s1avXl2lJy8pKVGdOnUkVVS7SkpKlJaWpsLCQs2aNUt5eXnau3dvlZ4jFkwo9NSghjSoCdUrAADCwWkBa8iQIVq0aJGaNGly2oNHjhypESNGaO3atRo7dqyGDRsmSVqzZo1ycnJO+Rg3btz3PnlycrJKS0slSQcPHlRycrISEhKUmJiomjVrKjs7W0VFRWf9+rKyMpWWlp74KC8vr+zPHTWO+o1eW1fR3B7vELAAAAgHvu/e0KNHjzM+cPfu3Vq6dKlmzZolSbrllls0atQorV+/Xm3btlV+fv55P3mXLl00ZcoUDRw4UAsWLNBtt92msrIyJSUlyRijpUuXatSoUWf9+p49e57y+aOPPqqxY8dW6rn3799/3uMNR+9s9mlfeU0NaVSqfftMqIcTs6JlPiE8MJ8QSMyn6vHdlqjTAtbZbNmyRWlpafL5Kr7EsixlZGRo8+bNatmy5Vm/rk+fPlqxYoX69Omjp556SvXr19drr72mMWPGqHPnzpo0aZKys7PVr18/paamatasWXr88ccVFxenIUOGKC0t7azfe/78+erQocOJzxMSEpSQkFDZH+mc/WGR4o1/+JWTJnVuUi/UQ4l50TCfED6YTwgk5lP1q3TAulCffPLJabeNGTPmxL+fe+65U+7r27ev+vbtW6nvnZSUdKKHKxatLTGav8PojWucUA8FAACcpNJvO2vcuLF27Nghv98vSTLGaPPmzcrIyAja4HBuEwo9pSRINzel9woAgHBS6YDVsGFDdezYUZMnT5Ykvfvuu0pPTz/n6UEET7lrNGmtp3ta2arhI2ABABBOTgtYI0eOVHp6urZu3ap+/fqdEqDy8vKUl5enzMxMjRs3ThMnTqzWweLfphcbfXNUeoC9rwAACDun9WDl5eWd9cFZWVlavHhxUAeEyplQ6Kl7qqU29aheAQAQbih/RKANpUafbDNUrwAACFOs0BHoxUJPdeOlIc2pXgEAEI4IWBHmuGc0sdDTj1raSqS5HQCAsETAijAfbDLaeYTmdgAAwhmrdISZUOjpqoaWLruI6hUAAOGKgBVBNh00+niL0QNZ/NoAAAhnrNQR5OW1npLipNtbUL0CACCcEbAihN8zernQ010tbCXFEbAAAAhnBKwI8fEWo62HaG4HACASsFpHiAmFnjrWlzo1oHoFAEC4I2BFgG2HjD7YTHM7AACRghU7Akws9FTDke5qya8LAIBIwIod5jxj9GKhpzuaW6oTz+lBAAAiAQErzM3earSpTBrRhl8VAACRglU7zE0o9HRpitSF5nYAACIGASuM7TpsNL24orndsghYAABECgJWGJu01pPPloa24tcEAEAkYeUOU54xmlDo6dZmluolUL0CACCSELDCVP52o6JSmtsBAIhErN5hanyBpzbJ0tWpVK8AAIg0BKww9M1Ro/eKjR5oTXM7AACRiIAVhl5d60mSfkRzOwAAEYkVPMwYYzS+wNPgZpbq16B6BQBAJCJghZlFO40KD0gjWvOrAQAgUrGKh5nxBZ5a1pFy0qheAQAQqQhYYWTfUaO3N9LcDgBApCNghZHJ6z25nnQPze0AAEQ0VvIwYYzRhAJPNzW1lJpI9QoAgEhGwAoTS3YbrdovPUBzOwAAEY/VPEyML/DUNEnqcwnVKwAAIh0BKwwcOGY0dYPR8Na2bJrbAQCIeASsMPDGek/lrnRvJr8OAACiASt6iH27c/sNGZYurkX1CgCAaEDACrFl3xit3EtzOwAA0YRVPcQmFBil15KuS6d6BQBAtCBghVDZcaM3ijzdn2XLsQlYAABECwJWCL1ZZHTouHRfFr8GAACiCSt7CE0o8NS/saWMJKpXAABEEwJWiHyx1+ifewzN7QAARCFW9xCZUOApLVEakEH1CgCAaEPACoHDfqPJ6z3dm2krjuZ2AACiDgErBN7eYHTgmDSc5nYAAKISK3wITCjwdO0llprVoXoFAEA0ImBVs9X7jD7dRXM7AADRjFW+mr1Y6KlBDenGJlSvAACIVgSsanTUb/TqOk/DMm3FOwQsAACiFQGrGk0rNtpXTnM7AADRjpW+Go0v8JSTZikzmeoVAADRjIBVTdaWGM3fYTSC5nYAAKIeq301ebHQU0qCdHNTqlcAAEQ7AlY1OOYaTVrr6e5Wtmr4CFgAAEQ7AlY1mL7JaM9RsfcVAAAxghW/Gowv8HR1qqW29aheAQAQCwhYQbah1OiTbTS3AwAQS1j1g+zFQk9146UhzaleAQAQKwhYQXTcM5pY6OlHLW0l0twOAEDMIGAF0YebjXYeobkdAIBYw8ofROMLPHVpYOmyi6heAQAQSwhYQbK5zOjjLTS3AwAQi1j9g+SlQk+14qTbW1C9AgAg1hCwgsDvGb1c6OmHLWwlxRGwAACINQSsIPh4i9HWQzS3AwAQq0gAQTCh0FPH+lKnBlSvAACIRQSsANt2yOjDzUYPZPFfCwBArCIFBNjEQk8JjnRXS/5rAQCIVaSAAPKM0YuFnu5obqlOPKcHAQCIVQSsAJq91WhTGc3tAADEOpJAAE0o9HRpinRVQ6pXAADEMgJWgOw6bDS9uKK53bIIWAAAxDICVoBMWuvJZ0tDW/FfCgBArCMNBMC3ze23NrNUL4HqFQAAsY6AFQD5243Wl9LcDgAAKpAIAmBCoafWyVL3RlSvAAAAAavKvjlqNG0jze0AAODfCFhV9OpaT5J0dyb/lQAAoAKpoAqMMRpf4GlwM0v1a1C9AgAAFQhYVbBop1HhAWkEze0AAOAkJIMqmFDgqWUdKSeN6hUAAPg3AtYF2l9u9PZGo+E0twMAgO8IasDat2+fOnXqpKSkpPO6/6GHHlJ2draeeOKJYA6vSiav8+T3pGE0twMAgO8IajqoXbu2Zs+era5du1b6/qVLl8rn82nhwoVavny5du3aFcwhXpBvm9tvbGopNZHqFQAAOFVQA1ZcXJxSUlLO6/4lS5aoV69ekqSePXtq2bJlwRziBVmy22jVfprbAQDAmZ2WEEaPHq2mTZvKsiytXLnylPvWrVunbt26KTMzU507d9bq1asDPqCSkhLVqVNHUkWFq6SkJODPUVXjCzw1TZL6XEL1CgAAnO60gDVkyBAtWrRITZo0Oe3BI0eO1IgRI7R27VqNHTtWw4YNkyStWbNGOTk5p3yMGzfuggaUnJys0tJSSdLBgweVnJx81seWlZWptLT0xEd5efkFPef5KD1mNHWD0fDWtmya2wEAwBn4vntDjx49zvjA3bt3a+nSpZo1a5Yk6ZZbbtGoUaO0fv16tW3bVvn5+QEZUJcuXTRlyhQNHDhQCxYs0G233XbWx/bs2fOUzx999FGNHTu2Us+zf//+CxrfxA1xKncTdFPDUu3bZy7oeyD6XOh8As6E+YRAYj5Vj++2PJ0WsM5my5YtSktLk89X8SWWZSkjI0ObN29Wy5Ytz/p1ffr00YoVK9SnTx899dRTql+/vl577TWNGTPmjPd37txZkyZNUnZ2tvr166fU1NSzfu/58+erQ4cOJz5PSEhQQkJCZX+kc/aHnc3r+cc1IMNSu0vqnffXIrpdyHwCzob5hEBiPlW/SgesC/XJJ5+cdtu34eps9z/33HOV+t5JSUkn+rWqw7I9Riv2Sv97Jc3tAADg7CqdFBo3bqwdO3bI7/dLqtiqYPPmzcrIyAja4MLN+AJP6bWk69LpvQIAAGdX6YDVsGFDdezYUZMnT5Ykvfvuu0pPTz/n6cFoUnbc6I0iT/dn2XJsAhYAADi70wLWyJEjlZ6erq1bt6pfv36nBKi8vDzl5eUpMzNT48aN08SJE6t1sKH0ZpHRoePSfVmcHgQAAOd2Wg9WXl7eWR+clZWlxYsXB3VA4WpCgaf+jS1lJFG9AgAA50Y5phK+2Gv0zz1GD7BzOwAAqAQSQyVMKPDUqKY0IIPqFQAA+H4ErO9x2G80eb2n+7JsxdHcDgAAKoGA9T3e3mB04Jh0P83tAACgkkgN32NCgadrL7HUvA7VKwAAUDkErHNYvc/o0100twMAgPNDcjiHFws9Nagh3diE6hUAAKg8AtZZHPUbvbrO07BMW/EOAQsAAFQeAessphUb7SuXhtPcDgAAzhPp4SwmFHjKSbOUmUz1CgAAnB8C1hmsLTHK30FzOwAAuDAkiDN4sdBTSoI0uCnVKwAAcP4IWN9xzDWatNbT3a1s1fARsAAAwPkjYH3H9E1Ge46K04MAAOCCkSK+Y0KBp6tTLbWtR/UKAABcGALWSTaUGs3eRnM7AACoGpLESV4q9FQ3Xrq1OdUrAABw4QhY/3LcM3q50NPQlrYSaW4HAABVQMD6lw83G+08QnM7AACoOtLEv4wv8NSlgaXLL6J6BQAAqoaAJWlzmdHHW2huBwAAgUGikPRyoadacdIdLaheAQCAqov5gOV6Ri8Verqrha2kOAIWAACoupgPWB9vNdp6SBrB6UEAABAgMZ8qxhd4uuIiqVMDqlcAACAwYjpgbTtk9OFmQ/UKAAAEVEwni4mFnhIc6a6WMf3fAAAAAixmk4VnKi6Nc3tzS3XiOT0IAAACJ2YD1vzdjorLaG4HAACBF7Pp4pWNcWpfT7qqIdUrAAAQWDEZsHYdNpqx3acRrW1ZFgELAAAEVkwGrElrPfksaWirmPzxAQBAkMVkwvjnHqNbM46rXgLVKwAAEHi+UA8gFF6/xtGBknJJtUI9FAAAEIVisoJVw2cpLiZ/cgAAUB2IGQAAAAFGwAIAAAgwAhYAAECAxWTAKi8v1x/+8AeVl5eHeiiIAswnBBLzCYHEfAodyxhjQj2I87V8+XJ16tRJy5YtU8eOHc/760tLS1W3bl0dOHBAderUCcIIEUuYTwgk5hMCifkUOjFZwQIAAAgmAhYAAECAReRGo0eOHJEkff311xf09WVlZZKklStXKikpKWDjQmxiPiGQmE8IJOZT9WrdurUSExMlRWgP1uuvv66hQ4eGehgAAAAnnNwbHpEB65tvvtHMmTPVtGlT1axZM9TDAQAAiPwKFgAAQDijyR0AACDACFgAAAABFlMB6+jRo7rpppuUmZmpyy+/XNdee63Wr18f6mEhCkycOFGWZen9998P9VAQwcrLyzVq1Ci1atVKl156KW/mQZV99NFH6tixozp06KD27dvrlVdeCfWQYkZEbtNQFSNGjFD//v1lWZaeffZZDR8+XPn5+aEeFiJYcXGxJkyYoK5du4Z6KIhwjz32mCzL0tq1a2VZlnbu3BnqISGCGWM0dOhQ5efn67LLLlNxcbFat26twYMHq3bt2qEeXtSLqQpWjRo1dP3118uyLElS165dVVxcHNpBIaJ5nqfhw4frmWeeUUJCQqiHgwh26NAhvfTSS/r9739/4hjVqFGjEI8Kkc6yLJWUlEiquGzORRddxLGqmsRUwPqu3Nxc3XjjjaEeBiLYX/7yF1199dXq1KlTqIeCCFdUVKSUlBQ98cQTuvLKK5Wdna05c+aEeliIYJZlaerUqRo8eLCaNGmi7t2765VXXlF8fHyohxYTYu4U4beeeOIJrV+/ngMYLtiqVav07rvvasGCBaEeCqKA3+/Xpk2b1LZtW40bN04rVqzQtddeq9WrVys1NTXUw0ME8vv9+t3vfqdp06apR48e+vzzzzVo0CB99dVXql+/fqiHF/VisoL15z//WdOmTdOMGTNObAgGnK+FCxequLhYrVq1UtOmTfXZZ59pxIgRev7550M9NESgjIwM2batH/7wh5KkK664Qs2aNdNXX30V4pEhUq1cuVLbt29Xjx49JEmdO3dWenq6VqxYEeKRxYaYC1h/+ctfNGXKFM2ePVvJycmhHg4i2E9+8hPt2LFDxcXFKi4uVteuXTV+/Hj95Cc/CfXQEIHq16+v3r17a+bMmZKkjRs3auPGjWrTpk2IR4ZI1bhxY+3YsePEdXvXr1+voqIiZWVl/f/27pgGghCKougLDRIoJpkCJ5gYA0igwAgF4Agf6JiMCTZkwz0KbvmTV/zNZWc4aohaltAAAABaSURBVCKccyqlJO+9QgiSJGutxhibywBAqrUqxqics4wxaq3puq7dWfhTzjn13vU8j4wxet9XpRTd97077Qi8ygEAAFjsuIkQAADg1ziwAAAAFuPAAgAAWOwD7gMQqDSUtPkAAAAASUVORK5CYII="
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# there is no exponential gradient vanishing\n",
    "# please compare it with the result in notebook `variational_qantum_circuit.ipynb`\n",
    "pyplot()\n",
    "ax = plot(2:maxN, grad_vars, yscale=:log10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Two qubits simulating a four-qubit systems\n",
    "\n",
    "Hadamard Gate with dashed box is applied only when measuring on the X basis.\n",
    "<div style=\"display:inline-block\"><img style=\"float: left\", src=\"images/hgate.png\" width=\"50\"/></div>\n",
    "\n",
    "The following gates Measure and Reset a qubit to 0.\n",
    "\n",
    "<div style=\"display:inline-block\"><img style=\"float:left\", src=\"images/mreset.png\" width=\"100\"></div>\n",
    "\n",
    "The goal of this section is to build the MPS-inspired sampler as our ansatz\n",
    "\n",
    "<div style=\"display:inline-block\"><img style=\"float:left\", src=\"images/twoqubit.png\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference: Variational Quantum Eigensolver with Fewer Qubits\n",
    "\n",
    "*Jin-Guo Liu, Yi-Hong Zhang, Yuan Wan, Lei Wang*\n",
    "\n",
    "https://arxiv.org/abs/1902.02663"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MeasureAndReset (generic function with 1 method)"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Hint\n",
    "MeasureAndReset(nbit::Int, i::Int) = Measure{nbit, 1, AbstractBlock}(Z, (i,), 0, false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ansatz_circuit_2qubit"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    ansatz_circuit_2qubit(nbit::Int, d::Int)\n",
    "\n",
    "Construct the above ansatz, `nrepeat` is the number of measure operations,\n",
    "`d` is the depth of each block between two measurements.\n",
    "\"\"\"\n",
    "function ansatz_circuit_2qubit(nbit::Int, d::Int)\n",
    "    Nv = 1\n",
    "    nbit_used = 1+Nv\n",
    "    nrepeat = nbit - Nv\n",
    "    circuit = chain(nbit_used)\n",
    "\n",
    "    for i=1:nrepeat\n",
    "        unit = chain(nbit_used)\n",
    "        for j=1:d\n",
    "            push!(unit, put(nbit_used, 1=>Rxz()))\n",
    "            push!(unit, put(nbit_used, 2=>Rxz()))\n",
    "            push!(unit, control(nbit_used, 1, 2=>shift(0.0)))\n",
    "        end\n",
    "        push!(unit, MeasureAndReset(nbit_used, 1))\n",
    "        if i==nrepeat\n",
    "            for k=2:nbit_used\n",
    "                push!(unit, MeasureAndReset(nbit_used, k))\n",
    "            end\n",
    "        end\n",
    "        push!(circuit, unit)\n",
    "    end\n",
    "    circuit\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice_size = 4\n",
    "vqe = VQE(TFI(lattice_size; h=0.5, periodic=false),\n",
    "    ansatz_circuit_2qubit(lattice_size, 2);\n",
    "    nshots=4000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nparameters(vqe.circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gensample"
      ]
     },
     "execution_count": 201,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    gensample(circuit, operator; nbatch=1024) -> Vector of Measure\n",
    "\n",
    "Generate samples from MPS-inspired circuit. Here, `nbatch` means nshot.\n",
    "`operator` is the operator to measure.\n",
    "This function returns a vector of `Measure` gates, results are stored in `m.results`.\n",
    "\"\"\"\n",
    "function gensample(circuit, operator; nbatch=1024)\n",
    "    mblocks = collect_blocks(Measure, circuit)\n",
    "    for m in mblocks\n",
    "        m.operator = operator\n",
    "    end\n",
    "    reg = zero_state(nqubits(circuit); nbatch=nbatch)\n",
    "    reg |> circuit\n",
    "    mblocks\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 0, E/N = 0.02584375\n",
      "Iter 10, E/N = -0.24709375\n",
      "Iter 20, E/N = -0.25425\n",
      "Iter 30, E/N = -0.27565625\n",
      "Iter 40, E/N = -0.27749999999999997\n",
      "Iter 50, E/N = -0.2753125\n",
      "Iter 60, E/N = -0.27965625\n",
      "Iter 70, E/N = -0.27903124999999995\n",
      "Iter 80, E/N = -0.27603125\n",
      "Iter 90, E/N = -0.27775\n",
      "Iter 100, E/N = -0.27890624999999997\n",
      "100.322192 seconds (1.15 G allocations: 130.021 GiB, 34.81% gc time)\n"
     ]
    }
   ],
   "source": [
    "dispatch!(vqe.circuit, :random)\n",
    "@time loss_history, circuit = train!(vqe; maxiter=100, α=0.3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CuYao"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gensample"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    gensample(circuit, operator; nbatch=1024) -> Vector of Measure\n",
    "\n",
    "Generate samples from MPS-inspired circuit. Here, `nbatch` means nshot.\n",
    "`operator` is the operator to measure.\n",
    "This function returns a vector of `Measure` gates, results are stored in `m.results`.\n",
    "\"\"\"\n",
    "function gensample(circuit, operator; nbatch=1024)\n",
    "    mblocks = collect_blocks(Measure, circuit)\n",
    "    for m in mblocks\n",
    "        m.operator = operator\n",
    "    end\n",
    "    reg = zero_state(nqubits(circuit); nbatch=nbatch) |> cu\n",
    "    reg |> circuit\n",
    "    mblocks\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 0, E/N = -0.0021562499999999984\n",
      "Iter 10, E/N = -0.24275000000000002\n",
      "Iter 20, E/N = -0.2655625\n",
      "Iter 30, E/N = -0.276625\n",
      "Iter 40, E/N = -0.277625\n",
      "Iter 50, E/N = -0.27696875\n",
      "Iter 60, E/N = -0.27659374999999997\n",
      "Iter 70, E/N = -0.28025\n",
      "Iter 80, E/N = -0.27690625\n",
      "Iter 90, E/N = -0.27653125\n",
      "Iter 100, E/N = -0.27893749999999995\n",
      " 81.006458 seconds (783.97 M allocations: 60.921 GiB, 22.79% gc time)\n"
     ]
    }
   ],
   "source": [
    "dispatch!(vqe.circuit, :random)\n",
    "@time loss_history, circuit = train!(vqe; maxiter=100, α=0.3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip5400\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip5400)\" points=\"\n",
       "0,1600 2400,1600 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip5401\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip5400)\" points=\"\n",
       "252.496,1503.47 2321.26,1503.47 2321.26,47.2441 252.496,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip5402\">\n",
       "    <rect x=\"252\" y=\"47\" width=\"2070\" height=\"1457\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  311.046,1503.47 311.046,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  798.962,1503.47 798.962,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1286.88,1503.47 1286.88,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1774.79,1503.47 1774.79,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2262.71,1503.47 2262.71,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.496,1472.02 2321.26,1472.02 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.496,1093.41 2321.26,1093.41 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.496,714.794 2321.26,714.794 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.496,336.182 2321.26,336.182 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,1503.47 2321.26,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,1503.47 252.496,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  311.046,1503.47 311.046,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  798.962,1503.47 798.962,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1286.88,1503.47 1286.88,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1774.79,1503.47 1774.79,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2262.71,1503.47 2262.71,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,1472.02 283.528,1472.02 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,1093.41 283.528,1093.41 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,714.794 283.528,714.794 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.496,336.182 283.528,336.182 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 311.046, 1557.47)\" x=\"311.046\" y=\"1557.47\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 798.962, 1557.47)\" x=\"798.962\" y=\"1557.47\">25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1286.88, 1557.47)\" x=\"1286.88\" y=\"1557.47\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1774.79, 1557.47)\" x=\"1774.79\" y=\"1557.47\">75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2262.71, 1557.47)\" x=\"2262.71\" y=\"1557.47\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 228.496, 1489.52)\" x=\"228.496\" y=\"1489.52\">-0.3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 228.496, 1110.91)\" x=\"228.496\" y=\"1110.91\">-0.2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 228.496, 732.294)\" x=\"228.496\" y=\"732.294\">-0.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 228.496, 353.682)\" x=\"228.496\" y=\"353.682\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 130.096, 775.359)\" x=\"130.096\" y=\"775.359\">Energy</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#009af9; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  311.046,88.4582 330.563,463.28 350.079,724.869 369.596,916.671 389.113,998.013 408.629,1071.5 428.146,1131.58 447.663,1169.02 467.179,1206.91 486.696,1244.35 \n",
       "  506.213,1240.65 525.729,1254.98 545.246,1267.46 564.762,1299.81 584.279,1257.75 603.796,1257.29 623.312,1320.15 642.829,1349.26 662.346,1319.68 681.862,1327.54 \n",
       "  701.379,1338.17 720.896,1336.32 740.412,1336.32 759.929,1356.2 779.445,1361.28 798.962,1354.81 818.479,1368.21 837.995,1385.78 857.512,1372.83 877.029,1352.5 \n",
       "  896.545,1347.88 916.062,1387.62 935.579,1333.55 955.095,1343.72 974.612,1400.57 994.128,1374.68 1013.65,1367.75 1033.16,1371.45 1052.68,1341.41 1072.2,1360.36 \n",
       "  1091.71,1340.02 1111.23,1347.42 1130.74,1408.42 1150.26,1394.09 1169.78,1374.68 1189.29,1373.3 1208.81,1389.47 1228.33,1395.94 1247.84,1400.57 1267.36,1378.84 \n",
       "  1286.88,1388.55 1306.39,1377.46 1325.91,1362.67 1345.43,1410.73 1364.94,1375.61 1384.46,1380.69 1403.98,1391.32 1423.49,1355.27 1443.01,1383.93 1462.53,1372.83 \n",
       "  1482.04,1386.24 1501.56,1364.52 1521.08,1407.5 1540.59,1399.18 1560.11,1389.01 1579.63,1374.68 1599.14,1361.28 1618.66,1361.28 1638.18,1412.12 1657.69,1414.43 \n",
       "  1677.21,1406.57 1696.73,1384.85 1716.24,1360.36 1735.76,1377.46 1755.28,1381.15 1774.79,1390.4 1794.31,1389.47 1813.83,1387.62 1833.34,1385.31 1852.86,1382.54 \n",
       "  1872.38,1389.01 1891.89,1404.72 1911.41,1407.04 1930.93,1414.43 1950.44,1372.83 1969.96,1380.69 1989.48,1424.14 2008.99,1409.81 2028.51,1381.62 2048.03,1396.41 \n",
       "  2067.54,1390.86 2087.06,1372.37 2106.58,1364.52 2126.09,1420.9 2145.61,1370.52 2165.13,1403.8 2184.64,1374.22 2204.16,1379.77 2223.68,1367.29 2243.19,1401.03 \n",
       "  2262.71,1380.69 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5402)\" style=\"stroke:#e26f46; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  311.046,1462.26 330.563,1462.26 350.079,1462.26 369.596,1462.26 389.113,1462.26 408.629,1462.26 428.146,1462.26 447.663,1462.26 467.179,1462.26 486.696,1462.26 \n",
       "  506.213,1462.26 525.729,1462.26 545.246,1462.26 564.762,1462.26 584.279,1462.26 603.796,1462.26 623.312,1462.26 642.829,1462.26 662.346,1462.26 681.862,1462.26 \n",
       "  701.379,1462.26 720.896,1462.26 740.412,1462.26 759.929,1462.26 779.445,1462.26 798.962,1462.26 818.479,1462.26 837.995,1462.26 857.512,1462.26 877.029,1462.26 \n",
       "  896.545,1462.26 916.062,1462.26 935.579,1462.26 955.095,1462.26 974.612,1462.26 994.128,1462.26 1013.65,1462.26 1033.16,1462.26 1052.68,1462.26 1072.2,1462.26 \n",
       "  1091.71,1462.26 1111.23,1462.26 1130.74,1462.26 1150.26,1462.26 1169.78,1462.26 1189.29,1462.26 1208.81,1462.26 1228.33,1462.26 1247.84,1462.26 1267.36,1462.26 \n",
       "  1286.88,1462.26 1306.39,1462.26 1325.91,1462.26 1345.43,1462.26 1364.94,1462.26 1384.46,1462.26 1403.98,1462.26 1423.49,1462.26 1443.01,1462.26 1462.53,1462.26 \n",
       "  1482.04,1462.26 1501.56,1462.26 1521.08,1462.26 1540.59,1462.26 1560.11,1462.26 1579.63,1462.26 1599.14,1462.26 1618.66,1462.26 1638.18,1462.26 1657.69,1462.26 \n",
       "  1677.21,1462.26 1696.73,1462.26 1716.24,1462.26 1735.76,1462.26 1755.28,1462.26 1774.79,1462.26 1794.31,1462.26 1813.83,1462.26 1833.34,1462.26 1852.86,1462.26 \n",
       "  1872.38,1462.26 1891.89,1462.26 1911.41,1462.26 1930.93,1462.26 1950.44,1462.26 1969.96,1462.26 1989.48,1462.26 2008.99,1462.26 2028.51,1462.26 2048.03,1462.26 \n",
       "  2067.54,1462.26 2087.06,1462.26 2106.58,1462.26 2126.09,1462.26 2145.61,1462.26 2165.13,1462.26 2184.64,1462.26 2204.16,1462.26 2223.68,1462.26 2243.19,1462.26 \n",
       "  2262.71,1462.26 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip5400)\" points=\"\n",
       "1700.29,312.204 2249.26,312.204 2249.26,130.764 1700.29,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1700.29,312.204 2249.26,312.204 2249.26,130.764 1700.29,130.764 1700.29,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#009af9; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  1724.29,191.244 1868.29,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1892.29, 208.744)\" x=\"1892.29\" y=\"208.744\">QMPS-Reuse</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip5400)\" style=\"stroke:#e26f46; stroke-width:12; stroke-opacity:1; fill:none\" points=\"\n",
       "  1724.29,251.724 1868.29,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip5400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1892.29, 269.224)\" x=\"1892.29\" y=\"269.224\">Exact</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = length(loss_history)\n",
    "gr()\n",
    "plot(0:M-1, [loss_history, fill(EG, M)], labels=[\"QMPS-Reuse\", \"Exact\"], lw=3, ylabel=\"Energy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}