SBRG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initializations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%reload_ext snakeviz\n",
    "%matplotlib inline\n",
    "%run SBRG2.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Core\n",
    "##Pauli Algebra\n",
    "Pauli matrix class. A Pauli matrix is denoted as\n",
    "$$\\sigma=\\mathrm{i}^{x\\cdot z}\\prod_{i\\in Xs}X_i\\prod_{i\\in Zs}Z_i=\\mathrm{i}^{x\\cdot z}\\prod_{i}X_i^{x_i}Z_i^{z_i},$$\n",
    "where $x_i=\\delta(i\\in Xs)$, $z_i=\\delta(i\\in Zs)$ are vectors of 0,1, and $X_i$, $Z_i$ are the Pauli matrix on site $i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pauli matrices can be created by `mkMat`\n",
    "- from Xs, Zs sets: `mkMat(Xs, Zs)`,\n",
    "- or from dict of indices: `mkMat({i:mu, ...})`,\n",
    "- or from list of indices: `mkMat([m0, mu1, ...])`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<Xs:[0, 1] Zs:[1, 2]>, <Xs:[0, 1] Zs:[1, 2]>, <Xs:[0, 1] Zs:[1, 2]>]"
      ]
     },
     "execution_count": 252,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[mkMat({0,1},{1,2}), mkMat({0:1,1:2,2:3}), mkMat([1,2,3])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Mat` is hashable and can be compared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6795291966596493067\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(hash(mkMat()))\n",
    "mkMat() == mkMat([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The commutativity of two Pauli operators can be checked by `is_commute`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 249,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_commute(mkMat([0,1,0]),mkMat([1,2,3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge two Pauli operators by indices. Coefficients are not calculated here. The coefficient can be restored from the commutation relations later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Xs:[0] Zs:[1, 2]>"
      ]
     },
     "execution_count": 250,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdot(mkMat([0,1,0]),mkMat([1,2,3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Pauli Monomial\n",
    "A term of Pauli operator `Term.mat` with a constant multiple `Term.val`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0 <Xs:[] Zs:[]>, 1.0 <Xs:[0, 1] Zs:[1]>, -3.2 <Xs:[0, 1] Zs:[1]>]"
      ]
     },
     "execution_count": 255,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[Term(), Term(mkMat([1,2])),Term(mkMat([1,2]),-3.2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider $\\sigma_A=\\mathrm{i}^{x_A\\cdot z_A} X_A Z_A$, $\\sigma_B=\\mathrm{i}^{x_B\\cdot z_B} X_B Z_B$, then the product reads:\n",
    "$$\\sigma_A\\sigma_B=\\mathrm{i}^{x_A\\cdot z_A+x_B\\cdot z_B+2 z_A\\cdot x_B}\\quad (X_A X_B) (Z_A Z_B)=\\mathrm{i}^{x_A\\cdot z_A+x_B\\cdot z_B-x_C\\cdot z_C+2 z_A\\cdot x_B}\\quad \\mathrm{i}^{x_C\\cdot z_C}X_C Z_C=\\mathrm{i}^{n^{C}_{AB}}\\; \\sigma_C,$$\n",
    "where $X_C=X_AX_B$, $Z_C=Z_AZ_B$ and $z_C = (z_A+z_B)\\mod 2$. The power of $\\mathrm{i}$ is\n",
    "$$n^{C}_{AB}=x_A\\cdot z_A+x_B\\cdot z_B-x_C\\cdot z_C+2 z_A\\cdot x_B$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 <Xs:[0, 1] Zs:[0]>"
      ]
     },
     "execution_count": 256,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dot(Term(mkMat([1,2])),Term(mkMat([3,3])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0 <Xs:[0] Zs:[1]>"
      ]
     },
     "execution_count": 257,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idot(Term(mkMat([2,0])),Term(mkMat([3,3])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Hamiltonian\n",
    "Hamiltonian `Ham` is a sum of Pauli monomial terms. It can be created from a list of terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0 <Xs:[0] Zs:[1]>, 0.8 <Xs:[] Zs:[1, 2]>, 0.5 <Xs:[1, 2] Zs:[]>, 0.2 <Xs:[] Zs:[0]>]"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = Ham(\n",
    "    [Term(mkMat([0,1,1]),0.5),\n",
    "     Term(mkMat([3,0,0]),0.2),\n",
    "     Term(mkMat([0,3,3]),0.8),\n",
    "     Term(mkMat([1,3,0]),1.0)])\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0 <Xs:[] Zs:[0]>, 0.8 <Xs:[] Zs:[1, 2]>, -0.5 <Xs:[0, 1, 2] Zs:[0, 1]>, -0.2 <Xs:[0] Zs:[1]>]"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.forward([Term(mkMat([0],[0,1]),-1.)])\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0 <Xs:[] Zs:[0]>, 0.8 <Xs:[] Zs:[1, 2]>, -0.5 <Xs:[0, 1, 2] Zs:[0, 1]>, -0.1 <Xs:[0] Zs:[1]>]"
      ]
     },
     "execution_count": 260,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.extend([Term(mkMat([1,3,0]),0.1)])\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0 <Xs:[] Zs:[0]>, -0.1 <Xs:[0] Zs:[1]>, -0.5 <Xs:[0, 1, 2] Zs:[0, 1]>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.remove(H.terms[1])\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[True, True, True, True]\n",
      "[[True, True, True], [True, True, True], [True, True]]\n",
      "[0, 1, 2, 3]\n"
     ]
    }
   ],
   "source": [
    "print([H.mats[term.mat] == term for term in H])\n",
    "print([[i in (term.mat.Xs | term.mat.Zs) for term in terms] for (i, terms) in H.imap.items()])\n",
    "print([term.pos for term in H])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SBRG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profiling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file '/var/folders/tl/lwpcq5qj049ftcj7pvhkzv_h0000gn/T/tmpnzp9nz'. \n"
     ]
    }
   ],
   "source": [
    "system = SBRG(TFIsing(256,J=4.,K=1.,h=1.,alpha=0.5))\n",
    "% snakeviz system.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benchmark with SBRG1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<__main__.SBRG at 0x106025410>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.seed(2)\n",
    "model = TFIsing(8,J=2.,K=1.,h=1.,alpha=1.)\n",
    "system = SBRG(model)\n",
    "system.tol = 0.\n",
    "system.max_rate = 1000.\n",
    "system.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Energy Coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "49\n"
     ]
    }
   ],
   "source": [
    "from IPython.display import clear_output\n",
    "def span(mat,L): # get the span of Zs of a Mat\n",
    "    mus = list(mat.Zs)\n",
    "    mus.sort()\n",
    "    diff = [snd-fst for fst, snd in zip(mus, mus[1:] + [mus[0]+L])]\n",
    "    i = diff.index(max(diff))\n",
    "    mu0 = mus[(i+1)%len(mus)]\n",
    "    mu1 = mus[i]\n",
    "    return (mu1-mu0)%L\n",
    "fname = 'ECdat_J2K1h1_a10'\n",
    "realization = 50\n",
    "L = 256\n",
    "try:\n",
    "    ECdat = load(fname)\n",
    "except:\n",
    "    ECdat = []\n",
    "for i in range(realization):\n",
    "    clear_output(wait=True)\n",
    "    print(i)\n",
    "    system = SBRG(TFIsing(L,J=2.,K=1.,h=1.,alpha=1.)).run()\n",
    "    ECdat.extend([[span(term.mat,L),len(term.mat.Zs),term.val] for term in system.Heff])\n",
    "dump(fname, ECdat)\n",
    "export(fname, ECdat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Statistics of Energy Coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "random.seed(2)\n",
    "system = SBRG(TFIsing(14,J=2.,K=1.,h=1.,alpha=1.)).run()\n",
    "from itertools import chain, combinations\n",
    "def subsets(iterable):\n",
    "    s = list(iterable)\n",
    "    return map(set, chain.from_iterable(combinations(s, r) for r in range(len(s)+1)))\n",
    "emap = [[s,0] for s in subsets(range(0,system.size))];\n",
    "for term in system.Heff:\n",
    "    Zs = term.mat.Zs\n",
    "    v = term.val\n",
    "    for em in emap:\n",
    "        if len(Zs & em[0])%2 == 1:\n",
    "            em[1] += v\n",
    "e0 = sum(term.val for term in system.Heff)\n",
    "es = [e0-2*em[1] for em in emap]\n",
    "export('Es_J2K1h1_a10_14',es)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Wave Function Overlap\n",
    "Benchmark with ED fidelity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def tolst(term):\n",
    "    return [list(term.mat.Xs),list(term.mat.Zs),term.val]\n",
    "random.seed(3)\n",
    "system = SBRG(TFIsing(8,J=1.,K=1.,h=1.,alpha=0.2))\n",
    "RGdat = []\n",
    "while system.phybits:\n",
    "    (H0, Rs, offdiag) = system.nextstep()\n",
    "    RGdat.append([tolst(H0),[tolst(R) for R in Rs],[tolst(S) for S in offdiag]])\n",
    "Hbdy = [tolst(term) for term in system.Hbdy]\n",
    "Heff = [tolst(term) for term in system.Heff]\n",
    "export('WFO_J1K1h1_a02_8',{'bits':system.size,'RGdat': RGdat,'Hbdy':Hbdy,'Heff':Heff})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###RG Flow\n",
    "Collect histogram of Hamiltonian coefficients (in logrithmic scale)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "realization = 512\n",
    "bins = [x/10 for x in range(-1,501)]\n",
    "hdic = {}\n",
    "for i in range(realization):\n",
    "    system = SBRG(TFIsing(512,J=1.,K=0.,h=1.,alpha_J=1.,alpha_K=1.,alpha_h=1.))\n",
    "    n = len(system.phybits)\n",
    "    while n>=4:\n",
    "        hs = [-math.log(abs(term.val)) for term in system.H]\n",
    "        hist, bin_edges = np.histogram(hs,bins)\n",
    "        try:\n",
    "            hdic[n] += hist\n",
    "        except:\n",
    "            hdic[n] = hist\n",
    "        system.flow(n/2)\n",
    "        n = len(system.phybits)\n",
    "hdic = {n: hist.tolist() for n, hist in hdic.items()}\n",
    "export('RG_J1K0h1_a101010_512', {'bins':bins,'hdic':hdic})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However the Hamiltonian coefficient distribution does not directly reflect how good is the perturbation. The small parameter of perturbation is off-diagonal to diagonal ratio (ODR). The following collects this statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1999\n"
     ]
    }
   ],
   "source": [
    "from IPython.display import clear_output\n",
    "import math\n",
    "import numpy as np\n",
    "def logODR(sys):\n",
    "    if not sys.H: # if H is empty\n",
    "        return[]\n",
    "    H0 = sys.H.terms[0]\n",
    "    E0 = abs(H0.val)\n",
    "    if not E0: # if E0 = 0, return\n",
    "        return []\n",
    "    Rs, pbit = sys.findRs(H0.mat)\n",
    "    sys.H.forward(Rs)\n",
    "    lr = [-math.log(abs(term.val/E0)) for term in sys.H.imap[pbit] if pbit in term.mat.Xs]\n",
    "    sys.H.backward(Rs)\n",
    "    return lr\n",
    "realization = 2000\n",
    "bins = [x/10 for x in range(-1,301)]\n",
    "fname = 'ODR_J1K1h1_a20_512'\n",
    "try:\n",
    "    hdic = load(fname)\n",
    "except:\n",
    "    hdic = {}\n",
    "for i in range(realization):\n",
    "    clear_output(wait=True)\n",
    "    print(i)\n",
    "    system = SBRG(TFIsing(512,J=1.,K=1.,h=1.,alpha=2.0))\n",
    "    n = len(system.phybits)\n",
    "    while n>=4:\n",
    "        hs = logODR(system)\n",
    "        hist, bin_edges = np.histogram(hs,bins)\n",
    "        try:\n",
    "            hdic[n] += hist\n",
    "        except:\n",
    "            hdic[n] = hist\n",
    "        system.flow(n/2)\n",
    "        n = len(system.phybits)\n",
    "dump(fname, hdic)\n",
    "hdic = {n: hist.tolist() for n, hist in hdic.items()}\n",
    "export(fname, {'bins':bins,'hdic':hdic})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dump(fname, hdic)\n",
    "hdic = {n: hist.tolist() for n, hist in hdic.items()}\n",
    "export(fname, {'bins':bins,'hdic':hdic})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After discussion with Kevin, I realize that it is better to use the $g$ statistics.\n",
    "$$g=\\frac{||\\Sigma||}{||H_0||}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import clear_output\n",
    "import math\n",
    "import numpy as np\n",
    "system = SBRG(TFIsing(512,J=1.,K=0.,h=2.,alpha=1.))\n",
    "hs = []\n",
    "while system.phybits:\n",
    "    system.nextstep()\n",
    "    hs.append(system.ODR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Entanglement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Half-system-size entanglement entropy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = 1.0\n",
    "N = 256\n",
    "filename = 'PD_a10_256'\n",
    "def EE(J0, K0, h0, sampling = 1):\n",
    "    try:\n",
    "        entropy_rec = load(filename)\n",
    "    except:\n",
    "        entropy_rec = {}\n",
    "    samples = []\n",
    "    L = int(N/2)\n",
    "    for s in range(sampling):\n",
    "        system = SBRG(TFIsing(N,J=J0,K=K0,h=h0,alpha=a))\n",
    "        system.max_rate = 1.8\n",
    "        system.run()\n",
    "        S = bipartite_entropy(system)\n",
    "        print('%3d: %f'%(s+1,S))\n",
    "        samples.append(S)\n",
    "    pt = (int(J0),int(K0),int(h0))\n",
    "    try:\n",
    "        entropy_rec[pt].extend(samples)\n",
    "    except:\n",
    "        entropy_rec[pt] = samples\n",
    "    dump(filename,entropy_rec)\n",
    "    export(filename,entropy_rec)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "EE(2,1,1,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Anderson Correlator\n",
    "First use the random cliford circuit to forward transform the operators to the holographic space. Then two operators are correlated if their product does not involves any flipping of the emergent qubit (such that can have expectation value $\\pm1$ in the product state of emergent qubits).\n",
    "\n",
    "Snapshot of the boundary operators mapped into the bulk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1,) ----------\n",
      "<Xs:[0, 2, 7] Zs:[0, 1, 7]>\n",
      "<Xs:[1, 2, 7] Zs:[0, 1, 7]>\n",
      "<Xs:[1, 3] Zs:[2, 4, 5, 6]>\n",
      "<Xs:[3] Zs:[]>\n",
      "<Xs:[4, 7] Zs:[0, 1, 2, 4, 5, 6, 7]>\n",
      "<Xs:[4, 7] Zs:[0, 1, 2, 4, 6, 7]>\n",
      "<Xs:[6, 7] Zs:[0, 1, 2, 4, 5, 6, 7]>\n",
      "<Xs:[0, 6, 7] Zs:[0, 1, 2, 4, 5, 6]>\n",
      "(1, 1) ----------\n",
      "<Xs:[0, 1] Zs:[]>\n",
      "<Xs:[2, 3, 7] Zs:[0, 1, 2, 4, 5, 6, 7]>\n",
      "<Xs:[1] Zs:[2, 4, 5, 6]>\n",
      "<Xs:[3, 4, 7] Zs:[0, 1, 2, 4, 5, 6, 7]>\n",
      "<Xs:[] Zs:[5]>\n",
      "<Xs:[4, 6] Zs:[5]>\n",
      "<Xs:[0] Zs:[7]>\n",
      "<Xs:[2, 6] Zs:[2, 4, 5, 6, 7]>\n"
     ]
    }
   ],
   "source": [
    "system = SBRG(XYZ(8, Jx=1,Jy=1,Jz=0,alpha=0.5))\n",
    "system.run()\n",
    "for mus in [(1,), (1,1)]:\n",
    "    print(mus, '----------')\n",
    "    ops = Ham([Term(mkMat({(i+j)%system.size:mus[i] for i in range(len(mus))}))\n",
    "               for j in range(system.size)])\n",
    "    ops.forward(system.RCC)\n",
    "    for term in ops:\n",
    "        print(term.mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collect Anderson-Edward correlator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999\n"
     ]
    }
   ],
   "source": [
    "from IPython.display import clear_output\n",
    "def AEC_update(system, mus, correlator):\n",
    "    L = system.size\n",
    "    terms = [Term(mkMat({(i+j)%L:mus[i] for i in range(len(mus))})) for j in range(L)]\n",
    "    new_corr = system.correlate(terms) # terms will be destroyed\n",
    "    for d, c in new_corr.items():\n",
    "        correlator[d] = correlator.get(d,0) + c\n",
    "realization = 1000\n",
    "fname = 'AEC_X1Y1Z0_a05_256'\n",
    "try:\n",
    "    AEC = load(fname)\n",
    "except:\n",
    "    AEC = {mus:{} for mus in [(1,),(2,),(3,),(1,1),(2,2),(3,3)]}\n",
    "for i in range(realization):\n",
    "    clear_output(wait=True)\n",
    "    print(i)\n",
    "    system = SBRG(XYZ(256, Jx=1,Jy=1,Jz=0,alpha=0.5))\n",
    "    system.run()\n",
    "    for mus, correlator in AEC.items():\n",
    "        AEC_update(system, mus, correlator)\n",
    "    if i%10 == 0:\n",
    "        dump(fname, AEC)\n",
    "        export(fname, AEC)\n",
    "dump(fname, AEC)\n",
    "export(fname, AEC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Find the case of long-range ZZ correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11\n",
      "found!\n"
     ]
    }
   ],
   "source": [
    "for i in range(1000):\n",
    "    clear_output(wait=True)\n",
    "    print(i)\n",
    "    system = SBRG(XYZ(64, Jx=1,Jy=1,Jz=0,alpha=0.5))\n",
    "    system.run()\n",
    "    opsX = Ham([Term(mkMat({i:1})) for i in range(system.size)])\n",
    "    opsY = Ham([Term(mkMat({i:2})) for i in range(system.size)])\n",
    "    opsZ = Ham([Term(mkMat({i:3})) for i in range(system.size)])\n",
    "    opsX.forward(system.RCC)\n",
    "    opsY.forward(system.RCC)\n",
    "    opsZ.forward(system.RCC)\n",
    "    corZ = {}\n",
    "    L = system.size\n",
    "    for (i,j) in combinations(range(len(opsZ)),2):\n",
    "        if len(opsZ.terms[i].mat.Xs ^ opsZ.terms[j].mat.Xs) == 0:\n",
    "            d = int(abs((j - i + L/2)%L - L/2))\n",
    "            corZ[d] = corZ.get(d,0) + 1\n",
    "    if any(d>3 for d in corZ):\n",
    "        print('found!')\n",
    "        export_Ham('taus',system.taus)\n",
    "        export('ops',[[i,list(opsX.terms[i].mat.Xs),\n",
    "                   list(opsY.terms[i].mat.Xs),\n",
    "                   list(opsZ.terms[i].mat.Xs)] for i in range(system.size)])\n",
    "        break"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}